求最大李雅普诺夫指数的matlab程序

程序一functiondx=Lorenz(t,x);dx(1,1)=10*(x(2)-x(1));dx(2,1)=x(1)*(30-x(3))-x(2);dx(3,1)=x(1)*x(2)-8/3*x(3);dx(4,1)=0;dx(5,1)=0;dx(6,1)=0;functionlambda_1=lyapunov_wolf1(data,N,m,tau,P)%该函数用来计算时间序列的最大Lyapunov指数--Wolf方法%m:嵌入维数%tau:时间延迟%data:时间序列%N:时间序列长度%P:时间序列的平均周期,选择演化相点距当前点的位置差，即若当前相点为I，则演化相点只能在|I－J|P的相点中搜寻%lambda_1:返回最大lyapunov指数值%**************************************************************************%ode计算整数阶系统的时间序列%******************************************************************delt_t1=0.001;t1=0:delt_t1:60;[tt1,y1]=ode45(@lorenz,t1,[-1,0,1]);xx1=y1(:,1)';x1=spline(tt1,xx1,t1);data=x1(20000:10:60000);%采样N=length(data);m=3;tau=11;%*****************************************************%FFT计算平均周期%**********************************************************x=data;xPower=abs(fft(x)).^2;NN=length(xPower);xPower(1)=[];%去除直流分量NN=floor(NN/2);xPower=xPower(1:NN);freq=(1:NN)/NN*0.5;[mP,index]=max(xPower);P=index;%*************************************************************min_point=1;%&&要求最少搜索到的点数MAX_CISHU=5;%&&最大增加搜索范围次数%FLYINGHAWK%求最大、最小和平均相点距离max_d=0;%最大相点距离min_d=1.0e+100;%最小相点距离avg_dd=0;Y=reconstitution(data,N,m,tau);%相空间重构M=N-(m-1)*tau;%重构相空间中相点的个数fori=1:(M-1)forj=i+1:Md=0;fork=1:md=d+(Y(k,i)-Y(k,j))*(Y(k,i)-Y(k,j));endd=sqrt(d);ifmax_ddmax_d=d;endifmin_ddmin_d=d;endavg_dd=avg_dd+d;endendavg_d=2*avg_dd/(M*(M-1));%平均相点距离dlt_eps=(avg_d-min_d)*0.02;%若在min_eps～max_eps中找不到演化相点时，对max_eps的放宽幅度min_eps=min_d+dlt_eps/2;%演化相点与当前相点距离的最小限max_eps=min_d+2*dlt_eps;%&&演化相点与当前相点距离的最大限%从P+1～M-1个相点中找与第一个相点最近的相点位置(Loc_DK)及其最短距离DKDK=1.0e+100;%第i个相点到其最近距离点的距离Loc_DK=2;%第i个相点对应的最近距离点的下标fori=(P+1):(M-1)%限制短暂分离，从点P+1开始搜索d=0;fork=1:md=d+(Y(k,i)-Y(k,1))*(Y(k,i)-Y(k,1));endd=sqrt(d);if(dDK)&(dmin_eps)DK=d;Loc_DK=i;endend%以下计算各相点对应的李氏数保存到lmd()数组中%i为相点序号，从1到(M-1)，也是i-1点的演化点；Loc_DK为相点i-1对应最短距离的相点位置，DK为其对应的最短距离%Loc_DK+1为Loc_DK的演化点，DK1为i点到Loc_DK+1点的距离，称为演化距离%前i个log2（DK1/DK）的累计和用于求i点的lambda值sum_lmd=0;%存放前i个log2（DK1/DK）的累计和fori=2:(M-1)%计算演化距离DK1=0;fork=1:mDK1=DK1+(Y(k,i)-Y(k,Loc_DK+1))*(Y(k,i)-Y(k,Loc_DK+1));endDK1=sqrt(DK1);old_Loc_DK=Loc_DK;%保存原最近位置相点old_DK=DK;%计算前i个log2（DK1/DK）的累计和以及保存i点的李氏指数if(DK1~=0)&(DK~=0)sum_lmd=sum_lmd+log(DK1/DK)/log(2);endlmd(i-1)=sum_lmd/(i-1);%以下寻找i点的最短距离：要求距离在指定距离范围内尽量短，与DK1的角度最小point_num=0;%&&在指定距离范围内找到的候选相点的个数cos_sita=0;%&&夹角余弦的比较初值——要求一定是锐角zjfwcs=0;%&&增加范围次数while(point_num==0)%*搜索相点forj=1:(M-1)ifabs(j-i)=(P-1)%&&候选点距当前点太近，跳过！continue;end%*计算候选点与当前点的距离dnew=0;fork=1:mdnew=dnew+(Y(k,i)-Y(k,j))*(Y(k,i)-Y(k,j));enddnew=sqrt(dnew);if(dnewmin_eps)|(dnewmax_eps)%&&不在距离范围，跳过！continue;end%*计算夹角余弦及比较DOT=0;fork=1:mDOT=DOT+(Y(k,i)-Y(k,j))*(Y(k,i)-Y(k,old_Loc_DK+1));endCTH=DOT/(dnew*DK1);ifacos(CTH)(3.14151926/4)%&&不是小于45度的角，跳过！continue;endifCTHcos_sita%&&新夹角小于过去已找到的相点的夹角，保留cos_sita=CTH;Loc_DK=j;DK=dnew;endpoint_num=point_num+1;endifpoint_num=min_pointmax_eps=max_eps+dlt_eps;zjfwcs=zjfwcs+1;ifzjfwcsMAX_CISHU%&&超过最大放宽次数，改找最近的点DK=1.0e+100;forii=1:(M-1)ifabs(i-ii)=(P-1)%&&候选点距当前点太近，跳过！continue;endd=0;fork=1:md=d+(Y(k,i)-Y(k,ii))*(Y(k,i)-Y(k,ii));endd=sqrt(d);if(dDK)&(dmin_eps)DK=d;Loc_DK=ii;endendbreak;endpoint_num=0;%&&扩大距离范围后重新搜索cos_sita=0;endendend%取平均得到最大李雅普诺夫指数lambda_1=sum(lmd)/length(lmd)functionX=reconstitution(data,N,m,tau)%该函数用来重构相空间%m为嵌入空间维数%tau为时间延迟%data为输入时间序列%N为时间序列长度%X为输出,是m*n维矩阵M=N-(m-1)*tau;%相空间中点的个数forj=1:M%相空间重构fori=1:mX(i,j)=data((i-1)*tau+j);endend以上是计算最大李氏指数的程序，可以运行。问题是里面的时间延迟tau和嵌入维数m该如何取值，程序里我是随便取的3和11*******************************************************************************程序二现我做了一个基于RHR算法的李氏指数计算方法，收敛速度很快，精确度也还可以，现传上：Lya1=[];Lya2=[];Lya3=[];V=eye(3);S=V;b1=0;a=0.4;c=0.2;gama=3.5;b=4.0;h=0.01;x(1)=0.1;y(1)=0;z(1)=0;n=0;whilez=200n=n+1;k1=h*y(n);m1=h*(-sin(x(n))-a*y(n)+b*cos(gama*z(n)).*sin(x(n))+c);k2=h*(y(n)+m1/2);m2=h*(-sin(x(n)+k1/2)-a*(y(n)+m1/2)+b*cos(gama*(z(n)+h/2)).*sin(x(n)+k1/2)+c);k3=h*(y(n)+m2/2);m3=h*(-sin(x(n)+k2/2)-a*(y(n)+m2/2)+b*cos(gama*(z(n)+h/2)).*sin(x(n)+k2/2)+c);k4=h*(y(n)+m3);m4=h*(-sin(x(n)+k3)-a*(y(n)+m3)+b*cos(gama*(z(n)+h)).*sin(x(n)+k3)+c);x(n+1)=x(n)+(k1+2*k2+2*k3+k4)/6;y(n+1)=y(n)+(m1+2*m2+2*m3+m4)/6;z(n+1)=n*h;J=[010;b*cos(gama*z(n+1))*cos(x(n+1))-cos(x(n+1))-a-b*gama*sin(gama*z(n+1))*sin(x(n+1));000];J=eye(3)+h*J;B=J*V*S;[V,S,U]=svd(B);a_max=max(diag(S));S=(1/a_max)*S;b1=b1+log(a_max);Lyapunov1=(log(diag(S))+b1)/(n*h);Lya1=[Lya1,Lyapunov1(1,];Lya2=[Lya2,Lyapunov1(2,];Lya3=[Lya3,Lyapunov1(3,];endLyapunov1n=1:20001;plot(n,Lya1,'k',n,Lya2,'k',n,Lya3,'k')%gridonaxis([0,30001,-0.8,0.5])title('LyapunovexponentsofWarship')xlabel('n'),ylabel('LCE')