lorenzlp.m

Lorenz系统的李雅谱诺夫指数

% Lyapunov exponents of Lorenz
function lorenzlp
Z1=[];
Z2=[];
Z3=[];
global a;
global b;
global r;
a=10;b=8/3;
for r=linspace(0,500,80);
    y=[1;1;1;1;0;0;0;1;0;0;0;1];
    lp=0;
    for k=1:100
        [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);
     y(4:12) = y0';
    end
    lp=lp/100;    
    Z1=[Z1 lp(1)];     
    Z2=[Z2 lp(2)];    
    Z3=[Z3 lp(3)];
end
r=linspace(0,500,80);
plot(r,Z1,'-',r,Z2,'-',r,Z3,'-');
title('Lyapunov exponents of Lorenz')
xlabel('parameter r');
ylabel('lyapunov exponents');
grid on