参考原理：

Lyapunov 指数的计算与仿真：https://wenku.baidu.com/view/ae8e4f80680203d8ce2f2476.html

Lorenz 系统的动力学方程

function dX = Lorenz(t,X,params) a = params(1);
b = params(2);
c = params(3);x=X(1);
y=X(2);
z=X(3);dX = zeros(3,1);
dX(1)=a*(y-x);
dX(2)=x*(b-z)-y;
dX(3)=x*y-c*z;end

使用 G. Benettin 计算方法

Z=[];
a=10;
c=8/3;
d0=1e-7;
bs = linspace(0,300,301);
transient = 50;
for b=bsparams = [a,b,c];lsum=0; x=1;y=1;z=1;   % #初始基准点x1=1;y1=1;z1=1+d0;  % #初始偏离点for i=1:100  [T1,Y1]=ode45(@(t,X) Lorenz(t,X,params),[0 1],[x;y;z]); [T2,Y2]=ode45(@(t,X) Lorenz(t,X,params),[0 1],[x1;y1;z1]); n1=length(Y1);n2=length(Y2);     x=Y1(n1,1);y=Y1(n1,2);z=Y1(n1,3);   x1=Y2(n2,1);y1=Y2(n2,2);z1=Y2(n2,3);   d1=sqrt((x-x1)^2+(y-y1)^2+(z-z1)^2);     % #新的偏离点在上一次计算的两轨迹末端的连线上，且距离仍等于d0x1=x+(d0/d1)*(x1-x);   y1=y+(d0/d1)*(y1-y);    z1=z+(d0/d1)*(z1-z); % #舍弃暂态过程的数据，因为初始基准点不一定在吸引子上if i> transientlsum=lsum+log(d1/d0);  end  end  Z=[Z lsum/(i-transient)];
end plot(bs,Z,'-k');
title('Lorenz System''s LLE v.s. parameter b')
xlabel('parameter b'),ylabel('Largest Lyapunov Exponents');
grid on;

代码二

clear;
yinit = [0.1,0.1,0.1];
orthmatrix = [1 0 0;0 1 0;0 0 1];y = zeros(12,1);
% 初始化输入
y(1:3) = yinit;
y(4:12) = orthmatrix;
tstart = 0; % 时间初始值
tstep = 1e-3; % 时间步长
wholetimes = 1e5; % 总的循环次数
steps = 10; % 每次演化的步数
iteratetimes = wholetimes/steps; % 演化的次数
mod = zeros(3,1);
lp = zeros(3,1);
% 初始化三个Lyapunov指数
Lyapunov1 = zeros(iteratetimes,1);
Lyapunov2 = zeros(iteratetimes,1);
Lyapunov3 = zeros(iteratetimes,1);
for i=1:iteratetimestspan = tstart:tstep:(tstart + tstep*steps);  [T,Y] = ode45(@(t,y) Lorenz_ly(t,y), tspan, y);% 取积分得到的最后一个时刻的值y = Y(size(Y,1),:);% 重新定义起始时刻tstart = tstart + tstep*steps;y0 = [y(4) y(7) y(10);y(5) y(8) y(11);y(6) y(9) y(12)];%正交化y0 = ThreeGS(y0);% 取三个向量的模mod(1) = sqrt(y0(:,1)'*y0(:,1));mod(2) = sqrt(y0(:,2)'*y0(:,2));mod(3) = sqrt(y0(:,3)'*y0(:,3));y0(:,1) = y0(:,1)/mod(1);y0(:,2) = y0(:,2)/mod(2);y0(:,3) = y0(:,3)/mod(3);lp = lp+log(abs(mod));%三个Lyapunov指数Lyapunov1(i) = lp(1)/(tstart);Lyapunov2(i) = lp(2)/(tstart);Lyapunov3(i) = lp(3)/(tstart);y(4:12) = y0';
end
% 作Lyapunov指数谱图
figure,
i = 1:iteratetimes;
plot(i,Lyapunov1,i,Lyapunov2,i,Lyapunov3)%G-S正交化
function A = ThreeGS(V) % V 为3*3向量
v1 = V(:,1);
v2 = V(:,2);
v3 = V(:,3);
a1 = zeros(3,1);
a2 = zeros(3,1);
a3 = zeros(3,1);
a1 = v1;
a2 = v2-((a1'*v2)/(a1'*a1))*a1;
a3 = v3-((a1'*v3)/(a1'*a1))*a1-((a2'*v3)/(a2'*a2))*a2;
A = [a1,a2,a3];
endfunction dX = Rossler_ly(t,X)
% Rossler吸引子，用来计算Lyapunov指数
%        a=0.15,b=0.20,c=10.0
%        dx/dt = -y-z,
%        dy/dt = x+ay,
%        dz/dt = b+z(x-c),
a = 0.20;
b = 0.20;
c = 5.7;
x=X(1); y=X(2); z=X(3);
% Y的三个列向量为相互正交的单位向量
Y = [X(4), X(7), X(10);X(5), X(8), X(11);X(6), X(9), X(12)];
% 输出向量的初始化，必不可少
dX = zeros(12,1);
% Rossler吸引子
dX(1) = -y-z;
dX(2) = x+a*y;
dX(3) = b+z*(x-c);
% Rossler吸引子的Jacobi矩阵
Jaco = [0 -1 -1;1 a 0;z 0 x-c];
dX(4:12) = Jaco*Y;
endfunction dX = Lorenz_ly(t,X)
% Lorenz 吸引子，用来计算Lyapunov指数
a = 10;
b = 28;
c = 8/3;
x=X(1); y=X(2); z=X(3);
% Y的三个列向量为相互正交的单位向量
Y = [X(4), X(7), X(10);X(5), X(8), X(11);X(6), X(9), X(12)];
% 输出向量的初始化，必不可少
dX = zeros(12,1);
% Lorenz 吸引子
dX(1)=a*(y-x);
dX(2)=x*(b-z)-y;
dX(3)=x*y-c*z;
% Lorenz 吸引子的Jacobi矩阵
Jaco = [-a a 0;b-z -1 -x;y x -c];
dX(4:12) = Jaco*Y;
end

matlab 计算 Lorenz 系统最大李雅普诺夫指数相关推荐

matlab混沌指数的计算,matlab实现混沌系统最大李雅普诺夫指数
李雅普诺夫指数(Lyapunov)是一个较为典型的判断一个系统是否具有混沌特性以及混沌的程度分析方法. 李指数:在相空间中初始时无限接近的两个轨道,随着时间的不断推移按指数收敛或发散的平均变化率,它可 ...
Matlab求解混沌系统最大李雅普诺夫指数
李雅普诺夫指数(Lyapunov)是一个较为典型的判断一个系统是否具有混沌特性以及混沌的程度分析方法. 李雅普诺夫指数:在相空间中初始时无限接近的两个轨道,随着时间的不断推移按指数收敛或发散的平均变化 ...
MATLAB计算气象水文要素年内分配指数
MATLAB计算气象水文要素年内分配指数降水/径流等水文要素年内分配指标 1 不均匀系数nuniformity coefficient (Cn) 1.1 原理 1.2 MATLAB代码 1.3 案例 ...
lorenz画图matlab,Matlab画Lorenz系统的最大李雅普诺夫指数图
Lorenz 系统文档分两个文件方程m文件和计算L指数m文件分开写,复制粘贴即可运行matlab2012a,改写方程文件和参数即可算自己的系统,其中最大L指数用的是经典的柏内庭(G.Benettin ...
应用matlab计算线性定常系统的矩阵指数
syms s t A = [5 0 3;0 5 2;9 1 3]; Fs = inv(seye(3)-A); %求预解矩阵FS=[(sI-A)]^(-1),eye(3)为33矩阵 eAt = ilap ...
李雅普诺夫指数对应的特征方向
洛伦兹系统的李雅普诺夫指数的符号为 (+,0,−)(+,0,-)(+,0,−) Any continuous time-dependent dynamical system without a fix ...
Matlab编程之混沌系统李雅普诺夫指数分析
简介李雅普诺夫指数是衡量混沌系统的一个重要参数,下面截图是对其具体解释. 代码实现: clc; clear; global kk;e=0ina1=0;final2=10;for kk=ina1:1: ...
Lorenz系统Lyapunov指数图
Lorenz系统Lyapunov指数图 Lyapunov指数是定量判断Lorenz系统混沌状态的重要指标,最大Lyapunov指数大于零则判定该系统是混沌系统. 结果图:
matlab 计算雅可比矩阵的lorenz李雅普图
1.内容简介略 629-可以交流.咨询.答疑 matlab 计算雅可比矩阵的lorenz李雅普图 2.内容说明略 3.仿真分析 clear clc close all x0_span = -5:0 ...
matlab绘制庞加莱截面_matlab庞加莱截面法画Lorenz系统分岔图(附图).doc
利用庞加莱截面法画的Lorenz 系统的分岔图,复制改成其他系统即可运行,Matlab12a可以运行,附分岔图,见下页! function Lorenz_bifur_r Z=[]; for r=li ...

matlab 计算 Lorenz 系统最大李雅普诺夫指数

参考原理：

Lorenz 系统的动力学方程

使用 G. Benettin 计算方法

代码二

matlab 计算 Lorenz 系统最大李雅普诺夫指数相关推荐

最新文章

热门文章