计算Lorenz系统的Lyapunov指数

Lorenz2.m

function dX = Lorenz2(t,X) % % Lorenz吸引子，用来计算Lyapunov指数 % dx/dt = a(y-x), % dy/dt = x(r-z)-y, % dz/dt = xy-bz, global a; % 变量不放入参数表中 global b; global r; 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*(r-z)-y; dX(3) = x*y-b*z; % Lorenz吸引子的Jacobi矩阵 Jaco = [-a a 0; r-z -1 -x; y x -b]; dX(4:12) = Jaco*Y;

llp.m

function llp global a; global b; global r; a=16;b=4;r=45.92; y=[1;1;1;1;0;0;0;1;0;0;0;1]; lp=0; for k=1:300 [T,Y] = ode45('Lorenz2', 1, y); y = Y(size(Y,1),:); y0 = [y(4) y(7) y(10); y(5) y(8) y(11); y(6) y(9) y(12)]; y0=GS(y0); mod(1)=norm(y0(:,1)); mod(2)=norm(y0(:,2)); mod(3)=norm(y0(:,3)); lp = lp+log(abs(mod)); y0(:,1)=y0(:,1)/mod(1); y0(:,2)=y0(:,2)/mod(2); y0(:,3)=y0(:,3)/mod(3); y(4:12) = y0'; end lp=lp/300

>> llp
 
lp =
 
    1.5096   -0.5995  -21.8105
 

相连的程序有GS.m