mylyapunov.m

Lyapunov 数值计算

% ??????Lyapunov?

%               ??lorenz??

%               dx/dt = SIGMA*(y - x)

%               dy/dt = BETA*x - y -x*z

%               dz/dt= x*y - GAMA*z

%   In this demo, SIGMA = 16,BETA = 45.92, GAMA = 4

%      ???: 

%      LE1 = 1.5, LE2 = 0.00, LE3 = -22.5

%      Author's E_mail: [e-mail]
clc;
sigma=16;beta=45.92;gama=4;
h=0.005;k=6000;
x=zeros(1,k);
y=zeros(1,k);
z=zeros(1,k);
x(1)=20;y(1)=10;z(1)=50;V=eye(3);S=V;I=V;b=0;



for i=1:k

    k11=sigma*(-x(i)+y(i));

    k21=beta*x(i)-x(i)*z(i)-y(i);

    k31=x(i)*y(i)-gama*z(i);

    

    k12=sigma*(-(x(i)+0.5*h*k11)+(y(i)+0.5*h*k21));

    k22=beta*(x(i)+0.5*h*k11)-(x(i)+0.5*h*k11)*(z(i)+0.5*h*k31)-(y(i)+0.5*h*k21);

    k32=(x(i)+0.5*h*k11)*(y(i)+0.5*h*k21)-gama*(z(i)+0.5*h*k31);

    

    k13=sigma*(-(x(i)+0.5*h*k12)+(y(i)+0.5*h*k22));

    k23=beta*(x(i)+0.5*h*k12)-(x(i)+0.5*h*k12)*(z(i)+0.5*h*k32)-(y(i)+0.5*h*k22);

    k33=(x(i)+0.5*h*k12)*(y(i)+0.5*h*k22)-gama*(z(i)+0.5*h*k32);

    

    k14=sigma*(-(x(i)+h*k13)+(y(i)+h*k23));

    k24=beta*(x(i)+h*k13)-(x(i)+h*k13)*(z(i)+h*k33)-(y(i)+h*k23);

    k34=(x(i)+h*k13)*(y(i)+h*k23)-gama*(z(i)+h*k33);

    

    J1=[-sigma sigma 0;beta-z(i) -1 -x(i);y(i) x(i) -gama];

    J2=[-sigma sigma 0;

        beta-(z(i)+0.5*h*k31) -1 -(x(i)+0.5*h*k11);

        y(i)+0.5*h*k21 x(i)+0.5*h*k11 -gama];

    J3=[-sigma sigma 0;

        beta-(z(i)+0.5*h*k32) -1 -(x(i)+0.5*h*k12);

        y(i)+0.5*h*k22 x(i)+0.5*h*k12 -gama];

    J4=[-sigma sigma 0;

        beta-(z(i)+h*k33) -1 -(x(i)+h*k13);

        y(i)+h*k23 x(i)+h*k13 -gama];

    

    J=I+h*(J1+2*J2*(I+0.5*h*J1)+2*J3*(I+0.5*h*J2*(I+0.5*h*J1))+J4*(I+h*J3*(I+0.5*h*J2*(I+0.5*h*J1))))/6;

        

    B=J*V*S;

    [V,S,U]=svd(B);

    am=max(diag(S));

    S=S/am;

    b=b+log(am);
x(i+1)=x(i)+h*(k11+2*k12+2*k13+k14)/6;
y(i+1)=y(i)+h*(k21+2*k22+2*k23+k24)/6;
z(i+1)=z(i)+h*(k31+2*k32+2*k33+k34)/6;

    

end

Lyapunov1=(log(diag(S))+b)/(k*h)