lorenz_correlation.m

来自「lorenz_lyapunov 用于混沌动力学的计算」· M 代码 · 共 46 行

clear;close all;
%Lorenz;m=5;tau=11;
space=0.01;tau=11;
X0=[3.0;4.0;5.0];[t,x]=ode45('DYdt',[0:space:30],X0);
x=x(1001:end,1);N=length(x);max_m=8;
for m=1:max_m 
    k=N-(m-1)*tau;
    for i=1:m
        for j=1:k
            y(i,j)=x(j+(i-1)*tau);
        end
    end
    [L,W]=size(y);
    for i=1:W
        for j=i:W
            d(i,j)=norm(y(:,i)-y(:,j),2);
            d(j,i)=d(i,j);
        end
    end
    rmin=min(min(d));rmax=max(max(d));
    R=[];CR=[];
    for r=rmin:(rmax-rmin)/20:rmax
        H=sum(sum((r-d>0)))./(W^2);
        Cr=H/(W.^2);
        R=[R r];
        CR=[CR Cr];
    end
    eval(['R' num2str(m) '=R;']);
    eval(['CR' num2str(m) '=CR;']);
end
figure
hold on
D=[];
for m=1:max_m
    x1=eval(['R' num2str(m)]);
    y1=eval(['CR' num2str(m)]);
    x2=log(x1(2:end));y2=log(y1(2:end));
    plot(x2,y2,'.');plot(x2,y2);
    x3=x2(1:5);y3=y2(1:5);
    beta = polyfit(x3,y3,1);
    D=[D beta(1)];
end
xlabel('log(R)');ylabel('log(CR)');
figure
plot(1:max_m,D);