《求最大李雅普诺夫指数的matlab程序.docx》由会员分享，可在线阅读，更多相关《求最大李雅普诺夫指数的matlab程序.docx（7页珍藏版）》请在三一文库上搜索。
1、程序一function dx=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;function lambda_1=lyapunov_wolf1(data,N,m,tau,P)%该函数用来计算时间序列的最大 Lyapunov指数-Wolf方法% m:嵌入维数% tau:时间延迟% data :时间序列% N:时间序列长度% P:时间序列的平均周期，选择演化相点距当前点的位置差，即若当前相点为I，则演化相点只能在
2、|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).A2;NN=length(xPower);xPower(1)=;%
3、去除直流分量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;for i = 1 : (M-1)for j = i+1 :
4、 M%最大相点距离%最小相点距离%相空间重构%重构相空间中相点的个数d = 0;for k = 1 : md = d + (Y(k,i)-Y(k,j)*(Y(k,i)-Y(k,j);endd = sqrt(d);if max_d < d max_d = d;endif min_d > dmin_d = d;endavg_dd = avg_dd + d;endendavg_d = 2*avg_dd/(M*(M-1);%平均相点距离dlt_eps = (avg_d - min_d) * 0.02 ;对max_eps的放宽幅度min_eps = min_d + dlt_eps / 2 ;
5、max_eps = min_d + 2 * dlt_eps ;%若在min_epsmax_eps中找不到演化相点时,%演化相点与当前相点距离的最小限%&&演化相点与当前相点距离的最大限% 从P+1M-1个相点中找与第一个相点最近的相点位置(Loc_DK)及其最短距离 DKDK = 1.0e+100;%第i个相点到其最近距离点的距离Loc_DK = 2;%第i个相点对应的最近距离点的下标for i = (P+1):(M-1)%限制短暂分离，从点 P+1开始搜索d = 0;for k = 1 : md = d + (Y(k,i)-Y(k,1)*(Y(k,i)-Y(k,1);endd
6、 = sqrt(d);if (d < DK) & (d > min_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 )的累计和for i = 2 : (
7、M-1)% 计算演化距离DK1 = 0;for k = 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点的最短距离：要求距离在指定
8、距离范围内尽量短，与 DK1的角度最小point_num = 0 ; % &&在指定距离范围内找到的候选相点的个数cos_sita = 0 ; %&&夹角余弦的比较初值一一要求一定是锐角zjfwcs=0 ;%&&增加范围次数while (point_num = 0)% *搜索相点for j = 1 : (M-1)if abs(j-i) <=(P-1)%&&候选点距当前点太近，跳过！continue;end%*计算候选点与当前点的距离dnew = 0;for k = 1 : mdnew = dnew + (Y(k,i)-Y(k
9、,j)*(Y(k,i)-Y(k,j);enddnew = sqrt(dnew);if (dnew < min_eps)|( dnew > max_eps )%&& 不在距离范围，跳过！continue;end%*计算夹角余弦及比较DOT = 0;for k = 1 : mDOT = DOT+(Y(k,i)-Y(k,j)*(Y(k,i)-Y(k,old_Loc_DK+1);endCTH = DOT/(dnew*DK1);if acos(CTH) > (3.14151926/4)%&& 不是小于 45 度的角，跳过!continue;endif C
10、TH > cos_sita %&&新夹角小于过去已找到的相点的夹角，保留 cos_sita = CTH;Loc_DK = j;DK = dnew;endpoint_num = point_num +1;endif point_num <= min_pointmax_eps = max_eps + dlt_eps;zjfwcs =zjfwcs +1;if zjfwcs > MAX_CISHU %&&超过最大放宽次数，改找最近的点 DK = 1.0e+100;for ii = 1 : (M-1)if abs(i-ii) <= (P-1)%&a
11、mp;&候选点距当前点太近，跳过！continue;endd = 0;for k = 1 : md = d + (Y(k,i)-Y(k,ii)*(Y(k,i)-Y(k,ii);endd = sqrt(d);if (d < DK) & (d > min_eps)DK = d;Loc_DK = ii;endendbreak;endpoint_num = 0;%&&扩大距离范围后重新搜索cos_sita = 0;endendend%取平均得到最大李雅普诺夫指数lambda_1=sum(lmd)/length(lmd)function X=reconstit
12、ution(data,N,m,tau)%该函数用来重构相空间% m为嵌入空间维数% tau为时间延迟% data为输入时间序列% N为时间序列长度% X为输出，是m*n维矩阵M=N-(m-1)*tau;% 相空间中点的个数for j=1:M%相空间重构for i=1:mX(i,j)=data(i-1)*tau+j);endend以上是计算最大李氏指数的程序，可以运行。问题是里面的时间延迟tau和嵌入维数 m该如何取值，程序里我是随便取的3和11*程序二现我做了一个基于RHF法的李氏指数计算方法，收敛速度很快，精确度也还可 以，现传上：Lya1=;Lya2=;Lya3=;V=eye(3);S=V
13、;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;while z<=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*c
14、os(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 = 0 1 0;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);0 0 0;J=eye(3)+
15、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);,ILya1=Lya1,Lyapunov1(1, a"Lya2=Lya2,Lyapunov1(2,少;Lya3=Lya3,Lyapunov1(3, W； endLyapunov1n=1:20001;plot(n,Lya1,'k',n,Lya2,k,n,Lya3,k)%grid onaxis(0,30001,-0.8,0.5)title('Lyapunov exponents of Warship')xlabel('n'),ylabel('LCE')