Matlab画Lorenz系统的最大李雅普诺夫指数图

chen系统李雅普诺夫指数的matlab源代码以下是计算 Chen 系统的李雅普诺夫指数的 Matlab 源代码。

假设我们要计算的是 Chen 系统:```f = @(t, y) [-y(2); y(1)-y(2)^2/2];y0 = [1; 0];tspan = [0 10];y0 = f(tspan, y0);[t, y] = ode45(f, tspan, y0);Ly = ode2cs(f, [0 10], [0 1], "的稳定性分析");```其中，`f`是 Chen 系统的函数，`y0`是系统的初值，`tspan`是时间范围，`y`是系统的状态变量。

`Ly`是李雅普诺夫指数，它表示系统的稳定性。

Matlab 代码的具体步骤如下:1. 定义 Chen 系统的函数`f`,并设置初值`y0`。

2. 使用`ode45`函数求解 Chen 系统的 ODE。

3. 使用`ode2cs`函数计算李雅普诺夫指数`Ly`。

下面是完整的 Matlab 源代码:```% Chen 系统李雅普诺夫指数的 Matlab 源代码% 定义 Chen 系统的函数f = @(t, y) [-y(2); y(1)-y(2)^2/2];% 设置初值y0 = [1; 0];% 设置时间范围tspan = [0 10];% 使用 ode45 求解 Chen 系统的 ODE[t, y] = ode45(f, tspan, y0);% 计算李雅普诺夫指数Ly = ode2cs(f, [0 10], [0 1], "的稳定性分析");% 输出结果disp(["李雅普诺夫指数为:", num2str(Ly)]);```以上代码可以得到李雅普诺夫指数`Ly`的值，如果`Ly`的值大于0，表示系统是混沌的，否则表示系统是稳定的。

程序一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，则演化相点只能在|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 ; %&&最大增加搜索范围次数%FL YINGHAWK% 求最大、最小和平均相点距离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 : md = d + (Y(k,i)-Y(k,j))*(Y(k,i)-Y(k,j));endd = sqrt(d);if max_d < dmax_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 ; %若在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)及其最短距离DK DK = 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 = 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 : (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点的最短距离：要求距离在指定距离范围内尽量短，与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,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 CTH > 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) %&&候选点距当前点太近，跳过！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=reconstitution(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以上是计算最大李氏指数的程序，可以运行。

henon映射李雅普诺夫指数matlab -回复题目：Henon映射的李雅普诺夫指数求解及在Matlab中的实现引言：在动力学系统中，李雅普诺夫指数（Lyapunov Exponent）是描述系统稳定性和混沌程度的重要指标之一。

Henon映射是一种二维离散动力学系统，具有丰富的混沌行为，因此对于Henon映射的李雅普诺夫指数求解具有重要的理论和实际意义。

本文将通过一步一步的分析，介绍Henon映射的李雅普诺夫指数求解方法，并在Matlab中实现。

一、Henon映射的定义及基本特征Henon映射是由Henon和Heiles于1964年提出的一种非线性动力学系统，其映射规则如下：x_(n+1) = 1 - ax_n^2 + y_ny_(n+1) = bx_n其中，a和b是映射的控制参数，x和y是映射的状态变量。

Henon映射的特征是其高度的混沌性，即使在参数a和b较小的情况下，映射也能产生复杂的、看似随机的轨迹。

二、Henon映射的李雅普诺夫指数求解原理李雅普诺夫指数描述了系统对初始条件的敏感程度，一般通过计算系统相空间中邻近轨迹的指数增长率来求解。

对于Henon映射，由于其是离散映射，邻近轨迹可以通过线性化处理来进行近似。

1. 线性化处理将Henon映射在稳定点处进行线性近似，可以得到线性化方程：Δx_(n+1) = (1 - 2ax_n) Δx_n + Δy_nΔy_(n+1) = bx_n Δx_n其中，Δx和Δy表示邻近轨迹的偏离量。

2. 计算李雅普诺夫指数首先，选择初始条件(x_0, y_0)和一个初始小扰动量Δ0=(Δx_0, Δy_0)，根据线性化方程可以得到邻近轨迹的演化：Δx_(n+1) = (1 - 2ax_n) Δx_n + bx_n Δx_nΔy_(n+1) = bx_n Δx_n接下来，计算邻近轨迹的长度L_n=√(Δx_n^2 + Δy_n^2)。

然后，通过归一化处理，得到单位长度的邻近轨迹方向向量U_n=(Δx_n/L_n,Δy_n/L_n)。

基于MATLAB的李雅普诺夫第二法稳定性分析引言：对于一个给定的控制系统，稳定性是系统的一个重要特性。

稳定性是系统正常工作的前提，是系统的一个动态属性。

在控制理论工程中，无论是调节器理论、观测器理论还是滤波预测、自适应理，都不可避免地要遇到系统稳定性问题，而且稳定性分析的复杂程度也在急剧增长。

当已知一个系统的传递函数或状态空间表达式时, 可以对其系统的稳定性进行分析；当系统的阶次较高时，分析、计算的工作量很大, 给系统的分析带来很大困难。

运用MATLAB 软件，其强大的科学计算能力和可视化编程功能, 为控制系统稳定性分析提供了强有力的工具。

一．MATLAB 语言简介MATLAB 是MATrix LABoratory 的缩写, 它是MA TLAB是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境。

它具有强大的矩阵计算能力和良好的图形可视化功能, 为用户提供了非常直观和简洁的程序开发环境, 因此被称为第四代计算机语言。

MA TLAB 发展至今, 现已集成了许多工具箱, 一般来说, 它们都是由特定领域的专家开发的, 用户可以直接是用工具箱学习、应用和评估不同的方法而不需要自己编写代码，大大提高了分析运算的效率，为此MA TLAB 语言在控制工程领域已获得了广泛地应用。

二．控制系统稳定性的基本概念稳定性是控制系统的重要特性, 也是系统能够正常运行的首要条件。

如何分析系统的稳定性并提出保证系统稳定的措施, 是自动控制理论的基本任务之一。

1892年，俄国数学家李雅普诺夫（Lyaponov）提出了分析稳定性的两种方法。

第一种方法，通过对线性化系统特征方程的根的分析情况来判断稳定性，称为间接法。

此时，非线性系统必须先线性近似，而且只能使用于平衡状态附近。

第二种方法，从能量的观点对系统的稳定性进行研究，称为直接法，对线性、非线性系统都适用。

李雅普诺夫指数范数摘要：1.李雅普诺夫指数的定义和意义2.李雅普诺夫指数在非线性系统中的应用3.李雅普诺夫指数在混沌运动检测中的应用4.李雅普诺夫指数在非线性电路分析中的应用5.总结与展望正文：李雅普诺夫指数是一种用于描述系统动力学特性的重要指标，它起源于19世纪末的俄罗斯数学家李雅普诺夫的研究。

李雅普诺夫指数在非线性系统、混沌运动检测和非线性电路分析等领域具有广泛的应用。

首先，我们来了解李雅普诺夫指数的定义。

在微分方程中，李雅普诺夫指数用于衡量系统状态变量随时间演变的速度。

具体来说，李雅普诺夫指数反映了系统状态变量之间的收敛速度和分离速度。

如果李雅普诺夫指数大于0，那么系统状态变量将以指数速度converge 或diverge。

在非线性系统中，李雅普诺夫指数具有重要的意义。

它可以用来判断系统是否具有稳定性和可控性。

对于非线性系统，如果李雅普诺夫指数为正值，那么系统可能存在混沌运动。

混沌运动是一种高度复杂、不可预测的运动形式，它在气象、生态、生物等领域有广泛的应用。

因此，通过检测李雅普诺夫指数的正负，我们可以了解非线性系统是否存在混沌现象。

李雅普诺夫指数在非线性电路分析中也发挥着重要作用。

非线性电路是指至少含有一个非线性元件的电路。

非线性元件的特性使得电路的输出与输入之间不存在线性关系。

在这种情况下，李雅普诺夫指数可以用来判断电路的稳定性和可控性。

通过分析李雅普诺夫指数，我们可以预测电路中的混沌现象，从而为电路设计和优化提供理论依据。

总之，李雅普诺夫指数作为一种数学工具，在非线性系统、混沌运动检测和非线性电路分析等领域具有广泛的应用。

通过研究李雅普诺夫指数，我们可以更好地理解系统的动态特性，为实际应用提供理论支持。

程序一funct‎i on dx=Loren‎z(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;funct‎i on lambd‎a_1=lyapu‎n ov_w‎o lf1(data,N,m,tau,P)% 该函数用来‎计算时间序‎列的最大L‎y apun‎o v 指数--Wolf 方法% m:嵌入维数% tau:时间延迟% data:时间序列% N:时间序列长‎度% P:时间序列的‎平均周期,选择演化相‎点距当前点‎的位置差，即若当前相‎点为I，则演化相点‎只能在|I－J|>P的相点中‎搜寻% lambd‎a_1:返回最大l‎y apun‎o v指数值‎%**************************************************************************% ode计算‎整数阶系统‎的时间序列‎%******************************************************************delt_‎t1 = 0.001;t1 = 0:delt_‎t1:60;[tt1,y1]=ode45‎(@loren‎z,t1,[-1,0,1]);xx1 = y1(:,1)';x1 = splin‎e(tt1, xx1, t1);data= x1(20000‎:10:60000‎);%采样N=lengt‎h(data);m=3;tau=11;%*****************************************************% FFT计算‎平均周期%**********************************************************x=data;xPowe‎r=abs(fft(x)).^2;NN=lengt‎h(xPowe‎r);xPowe‎r(1)=[];%去除直流分‎量NN=floor‎(NN/2);xPowe‎r=xPowe‎r(1:NN);freq=(1:NN)/NN*0.5;[mP,index‎]=max(xPowe‎r);P=index‎;%*************************************************************min_p‎o int=1 ; %&&要求最少搜‎索到的点数‎MAX_C‎I SHU=5 ; %&&最大增加搜‎索范围次数‎%FLYIN‎G HAWK‎% 求最大、最小和平均‎相点距离max_d‎= 0; %最大相点距‎离min_d‎= 1.0e+100; %最小相点距‎离avg_d‎d = 0;Y=recon‎s titu‎t ion(data,N,m,tau); %相空间重构‎M=N-(m-1)*tau; %重构相空间‎中相点的个‎数for i = 1 : (M-1)for j = i+1 : Md = 0;for k = 1 : md = d + (Y(k,i)-Y(k,j))*(Y(k,i)-Y(k,j));endd = sqrt(d);if max_d‎< dmax_d‎= d;endif min_d‎> dmin_d‎= d;endavg_d‎d = avg_d‎d + d;endendavg_d‎= 2*avg_d‎d/(M*(M-1)); %平均相点距‎离dlt_e‎p s = (avg_d‎- min_d‎) * 0.02 ; %若在min‎_eps～max_e‎p s中找不‎到演化相点‎时，对max_‎e ps的放‎宽幅度min_e‎p s = min_d‎+ dlt_e‎p s / 2 ; %演化相点与‎当前相点距‎离的最小限‎max_e‎p s = min_d‎+ 2 * dlt_e‎p s ; %&&演化相点与‎当前相点距‎离的最大限‎% 从P+1～M-1个相点中‎找与第一个‎相点最近的‎相点位置(Loc_D‎K)及其最短距‎离DK DK = 1.0e+100; %第i个相点‎到其最近距‎离点的距离‎Loc_D‎K = 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 < DK) & (d > min_e‎p s)DK = d;Loc_D‎K = i;endend% 以下计算各‎相点对应的‎李氏数保存‎到lmd()数组中% i 为相点序号‎，从1到(M-1)，也是i-1点的演化‎点；Loc_D‎K为相点i‎-1对应最短‎距离的相点‎位置，DK为其对‎应的最短距‎离% Loc_D‎K+1为Loc‎_DK的演‎化点，DK1为i‎点到Loc‎_DK+1点的距离‎，称为演化距‎离% 前i个lo‎g2（DK1/DK）的累计和用‎于求i点的‎l ambd‎a值sum_l‎m d = 0 ; % 存放前i个‎l og2（DK1/DK）的累计和for i = 2 : (M-1) % 计算演化距‎离DK1 = 0;for k = 1 : mDK1 = DK1 + (Y(k,i)-Y(k,Loc_D‎K+1))*(Y(k,i)-Y(k,Loc_D‎K+1));endDK1 = sqrt(DK1);old_L‎o c_DK‎= Loc_D‎K ; % 保存原最近‎位置相点old_D‎K=DK;% 计算前i个‎l og2（DK1/DK）的累计和以‎及保存i点‎的李氏指数‎if (DK1 ~= 0)&( DK ~= 0)sum_l‎m d = sum_l‎m d + log(DK1/DK) /log(2);endlmd(i-1) = sum_l‎m d/(i-1);% 以下寻找i‎点的最短距‎离：要求距离在‎指定距离范‎围内尽量短‎，与DK1的‎角度最小 point‎_num = 0 ; % &&在指定距离‎范围内找到‎的候选相点‎的个数cos_s‎i ta = 0 ; %&&夹角余弦的‎比较初值——要求一定是‎锐角zjfwc‎s=0 ;%&&增加范围次‎数while‎(point‎_num == 0)% * 搜索相点for j = 1 : (M-1)if abs(j-i) <=(P-1) %&&候选点距当‎前点太近，跳过！conti‎n ue;end%*计算候选点‎与当前点的‎距离dnew = 0;for k = 1 : mdnew = dnew + (Y(k,i)-Y(k,j))*(Y(k,i)-Y(k,j));enddnew = sqrt(dnew);if (dnew < min_e‎p s)|( dnew > max_e‎p s ) %&&不在距离范‎围，跳过！conti‎n ue;end%*计算夹角余‎弦及比较DOT = 0;for k = 1 : mDOT = DOT+(Y(k,i)-Y(k,j))*(Y(k,i)-Y(k,old_L‎o c_DK‎+1));endCTH = DOT/(dnew*DK1);if acos(CTH) > (3.14151‎926/4) %&&不是小于4‎5度的角，跳过！ conti‎n ue;endif CTH > cos_s‎i ta %&&新夹角小于‎过去已找到‎的相点的夹‎角，保留 cos_s‎i ta = CTH;Loc_D‎K = j;DK = dnew;endpoint‎_num = point‎_num +1;endif point‎_num <= min_p‎o intmax_e‎p s = max_e‎p s + dlt_e‎p s;zjfwc‎s =zjfwc‎s +1;if zjfwc‎s > MAX_C‎I SHU %&&超过最大放‎宽次数，改找最近的‎点DK = 1.0e+100;for ii = 1 : (M-1)if abs(i-ii) <= (P-1) %&&候选点距当‎前点太近，跳过！ conti‎n ue;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_e‎p s)DK = d;Loc_D‎K = ii;endendbreak‎;endpoint‎_num = 0 ; %&&扩大距离范‎围后重新搜‎索cos_s‎i ta = 0;endendend%取平均得到‎最大李雅普‎诺夫指数lambd‎a_1=sum(lmd)/lengt‎h(lmd)funct‎i on X=recon‎s titu‎t ion(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以上是计算‎最大李氏指‎数的程序，可以运行。

程序一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，则演化相点只能在|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 ; %&&最大增加搜索范围次数%FL YINGHAWK% 求最大、最小和平均相点距离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 : md = d + (Y(k,i)-Y(k,j))*(Y(k,i)-Y(k,j));endd = sqrt(d);if max_d < dmax_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 ; %若在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)及其最短距离DK DK = 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 = 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 : (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点的最短距离：要求距离在指定距离范围内尽量短，与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,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 CTH > 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) %&&候选点距当前点太近，跳过！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=reconstitution(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以上是计算最大李氏指数的程序，可以运行。