lorenz_lyapunov.m

小数据量法求Lyapunov指数_龙格库塔

clear
clc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
h=0.01;
k=8000;
x(1)=0;
y(1)=1;
z(1)=0;
for i=1:k
    kx1=fx(x(i),y(i));
    ky1=fy(x(i),y(i),z(i));
    kz1=fz(x(i),y(i),z(i));
    
    kx2=fx(x(i)+(1/2)*h*kx1,y(i)+(1/2)*h*ky1);
    ky2=fy(x(i)+(1/2)*h*kx1,y(i)+(1/2)*h*ky1,z(i)+(1/2)*h*kz1);
    kz2=fz(x(i)+(1/2)*h*kx1,y(i)+(1/2)*h*ky1,z(i)+(1/2)*h*kz1);
   
    kx3=fx(x(i)+(1/2)*h*kx2,y(i)+(1/2)*h*ky2);
    ky3=fy(x(i)+(1/2)*h*kx2,y(i)+(1/2)*h*ky2,z(i)+(1/2)*h*kz2);
    kz3=fz(x(i)+(1/2)*h*kx2,y(i)+(1/2)*h*ky2,z(i)+(1/2)*h*kz2);
    
    kx4=fx(x(i)+h*kx3,y(i)+h*ky3);
    ky4=fy(x(i)+h*kx3,y(i)+h*ky3,z(i)+h*kz3);
    kz4=fz(x(i)+h*kx3,y(i)+h*ky3,z(i)+h*kz3);

    x(i+1)=x(i)+(h/6)*(kx1+2*kx2+2*kx3+kx4);
    y(i+1)=y(i)+(h/6)*(ky1+2*ky2+2*ky3+ky4);
    z(i+1)=z(i)+(h/6)*(kz1+2*kz2+2*kz3+kz4);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Data = x(1000:1:6000);
m=8;
tau=10;
N=5001;
P=ave_T(Data,N);
[Y_num,Y]=small_data_sets_Lyapunov1(Data,N,m,tau,P,h);
XX=1:Y_num;
plot(XX,Y(XX));hold on;
linear=80;%线性长度的变化引起最大Lapunov指数的变化，线性长度和最大Lapunov指数也受步长h的影响
XX1=1:linear;
pp=polyfit(XX1,Y(XX1),1);
lambda_1=pp(1)
yp=polyval(pp,XX1);
plot(XX1,yp,'r');