matlab 计算雅可比矩阵的lorenz李雅普图

最新推荐文章于 2023-12-07 12:22:08 发布

studyer_domi

最新推荐文章于 2023-12-07 12:22:08 发布

阅读量384

点赞数

分类专栏： Matlab系列案例文章标签： matlab 雅可比矩阵 lorenz 李雅普图

本文链接：https://blog.csdn.net/qingfengxd1/article/details/128266833

版权

Matlab系列案例专栏收录该内容

485 篇文章 244 订阅

订阅专栏

1、内容简介

略
629-可以交流、咨询、答疑

matlab 计算雅可比矩阵的lorenz李雅普图

2、内容说明

略

3、仿真分析

clear
clc
close all
x0_span = -5:0.01:5;
Z = [];
for i = 1:length(x0_span)
    i
    ly(i,:) = ly_spectra(x0_span(i));
end

figure
plot(ly(:,[2 3 4 5]),'DisplayName','ly')
figure
plot(sum(ly,2))

function ly = ly_spectra(x1)
% clear 
% clc
T = 13.5;
dt = 0.01;
a = 3.39;
b = 1.95;
c = 9.36;
r = 0.42;
g = 1.42;
A = 0.0676;
B = 0.3682;
N=T/dt; 
x2=0;x3=0;x4=0; x5=0;
[t,x] = ode45('func',[0:1:T],[x1;x2;x3;x4;x5]);
n=length(t);
x1=x(n,1); x2=x(n,2); x3=x(n,3); x4=x(n,4); x5=x(n,5);
[t,x] = ode45('func',[0:dt:T],[x1;x2;x3;x4;x5]);

JN = eye(5);
w = eye(5);
J = eye(5);
for k=1:N
    x1 = x(k,1);
    x2 = x(k,2);
    x3 = x(k,3);
    x4 = x(k,4);
    x5 = x(k,5);
    W1 = A+B*abs(x4);
    W2 = A+B*abs(x5);
    J = (eye(5)+[-c*r-c/W2 c c 0 0; -a a*g a*g 0 0; -b b*g b*g-b*W1 0 0;0 0 1 0 0;1/W2 0 0 0 0]*dt);
    [~,w] = qr(J*w);

    w(:,1) = w(:,1)/norm(w(:,1));
    w(:,2) = w(:,2)/norm(w(:,2));
    w(:,3) = w(:,3)/norm(w(:,3));
    w(:,4) = w(:,4)/norm(w(:,4));
    w(:,5) = w(:,5)/norm(w(:,5));
end