《求最大李雅普诺夫指数的matlab程序》由会员分享，可在线阅读，更多相关《求最大李雅普诺夫指数的matlab程序（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,则演化相点只能在IIJ
2、卜P的相点中搜寻% lambda_l:返回最大lyapunov指数值*1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1*1* *1* *1* *1* *1* *1* *1* *1* *1* *1*1*1* *1* *1*J/-. Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx
3、 Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx TxTx Tx Tx Tx% ode 计算整数阶系统的时间序列zy *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1*1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1*1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1*1* *1* *1* *1* *1* *1* *1* *1* *1* *1*
4、*1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1*1* *1* *1* *1* *1* *1*J/-. Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx delt_tl = 0.00l;tl = 0:delt_tl:60;
5、ttl，yl=ode45(lorenz，tl，-l，0，l);xxl = yl(:，l);xl = spline(ttl， xxl， tl);data= x1(20000:10:60000);% 采样N=length(data);m=3;tau=11;zyJ/-. Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx% FFT 计算平均周期
6、zy *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1*J/-. Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx T
7、x Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Txx=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 ; %
8、&最大增加搜索范围次数 %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 : Md = 0;for k = 1 : m%最大相点距离%最小相点距离%相空间重构%重构相空间中相点的个数d = d + (Y(k,i)-Y(k,j)*(Y(k,i)-Y(k,j); end d = sqrt(d); if max_d d min_d = d;endavg_dd = avg_dd
9、 + d; endendavg_d = 2*avg_dd/(M*(M-1);%平均相点距离%若在min_epsmax_eps中找不到演化相点时,%演化相点与当前相点距离的最小限%&演化相点与当前相点距离的最大限dlt_eps = (avg_d - min_d) * 0.02 ; 对 max_eps 的放宽幅度min_eps = min_d + dlt_eps / 2 ;max_eps = min_d + 2 * dlt_eps ;% 从P+1M-1个相点中找与第一个相点最近的相点位置(Loc_DK)及其最短距离DK DK = 1.0e+100;%第i个相点到其最近距离点的距离Loc_DK =
10、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 = sqrt(d);if (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点的距离，称
11、为演化距离 % 前i个log2(DK1/DK)的累计和用于求i点的lambda值summd = 0 ;% 存放前i个log2 (DK1/DK)的累计和for i = 2 : (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)su
12、m_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)% * 搜索相点 for j = 1 : (M-1) if abs(j-i) =(P-1)%&候选点距当前点太近，跳过！continue;end%*计算候选点与当前
13、点的距离dnew = 0; for k = 1 : m dnew = dnew + (Y(k,i)-Y(k,j)*(Y(k,i)-Y(k,j);end dnew = sqrt(dnew);if (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;en
14、dif CTH cos_sita %&新夹角小于过去已找到的相点的夹角，保留 cos_sita = CTH;Loc_DK = j;DK = dnew;endpoint_num = point_num +1;endif point_num MAX_CISHU %&超过最大放宽次数，改找最近的点 DK = 1.0e+100;for ii = 1 : (M-1)if abs(i-ii) = (P-1)%&候选点距当前点太近，跳过！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 min_eps)DK = d;L
