李亚普诺夫指数计算方法
克隆法计算李雅普诺夫指数
克隆法计算李雅普诺夫指数李雅普诺夫指数（Lyapunov exponent）是用来衡量动态系统稳定性的一个重要指标。
在混沌理论中，李雅普诺夫指数可以用来预测一个系统对初始条件的敏感性。
克隆法是一种计算李雅普诺夫指数的数值方法。
其基本思想是：对于一个给定的动态系统，我们首先生成两个几乎完全相同的初始条件，然后让它们分别演化。
随着时间的推移，这两个初始条件会逐渐分离，我们可以通过测量它们之间的距离来计算李雅普诺夫指数。
具体步骤如下：
生成两个几乎完全相同的初始条件。
让这两个初始条件分别演化。
计算它们之间的距离。
重复上述步骤多次，并取平均值。
将平均值与时间作图，并求出斜率，即为李雅普诺夫指数。
需要注意的是，克隆法是一种数值方法，其结果会受到初始条件的选择、时间步长的选择等因素的影响。
因此，在使用克隆法计算李雅普诺夫指数时，需要选择合适的参数，并进行多次模拟以获得更准确的结果。
最大李指数计算方法
物理学报ACTA PHYSICA SINICA2000 Vol.49 No.4 P.636-640一种最大李雅普诺夫指数估计的稳健算法杨绍清章新华赵长安最大李雅普诺夫指数是诊断和描述动态系统混沌的重要参数.在深入研究相空间重构技术和轨道跟踪法的基础上，提出了一种从标量混沌时间序列中估计最大李雅普诺夫指数的新算法.该算法能够克服现有算法的不足，主要有以下三个优点：1）很高的精度；2）几乎不受噪声的影响；3）所需的计算时间和存贮空间小, 能进行在线计算.PACC： 0545A ROBUST METHOD FOR ESTIMATING THE LARGESTLYAPUNOV EXPONENTYANG SHAO-QING(Harbin Institute of Technology,Harbin 150001,China)ZHANG XIN-HUA(Dalian Naval Academy,Dalian 116018,China)ZHAO CHANG-AN(Harbin Institute of Technology,Harbin 150001,China)ABSTRACTThe largest Lyapunov exponent is an important parameter of detecting and characterizing chaos produced from a dynamical system. In this paper, based on the technology of phase space reconstruction and the methods of trajectory tracing, a new algorithm is proposed for estimating the largest Lyapunov exponent from a scalar chaotic time series. This method, which can overcome the deficiencies of the existing methods, has three main advantages: (1) It has highly accurate results; (2) It is little affected by noise; (3) It only needs a little time of computation and small space of memory and can calculate the largest Lyapunov exponent on line.1 引言近几年来，混沌信号的诊断及其特性的描述已经广泛地应用于时间序列的分析中［1］.在诊断和描述混沌信号时，最大李雅普诺夫指数(λ1)不仅是一个很重要的不变量，而且是判断混沌存在的一个重要依据.因此，利用λ1去诊断和描述混沌仍然是一种主要方法.其他有些方法(如K-熵法等［2］)与这种方法没有本质的区别.目前估计λ1的方法主要有两种［3］：一种是分析法(analytic approach)，一种是轨道跟踪法(trajectory tracing method).前者是用一个函数（如局部多项式或神经网络等）来建立系统模型，然后估计系统的雅可比矩阵，进而求取λ1；后者是直接从λ1的定义出发跟踪系统的两条轨道，获取λ1.由于轨道跟踪法不像分析法那样易受系统拓扑结构的影响，因而受到许多学者的高度重视.自从1985年Wolf［4］提出轨道跟踪法以来，这种方法得到了较大的改进.最有代表性的是Rosenstein等［5］的工作,他们给出了一种比较优化的算法: 首先绘制<ln dk(j)> 相对j。
李雅普诺夫方法
满足渐近稳定的球域 S ( )只是状态空间中的有限部分，这时称平 衡状态 x e 为局部渐近稳定，并且称S ( ) 为渐近稳定吸引区，表示只 有从该区域出发的受扰运动才能被“吸引”至平衡状态x e 。
线性系统若是渐近稳定（且A非奇异），必为全局渐近稳定。非 线性系统一般只能是小范围渐近稳定。
一般情况下，李雅普诺夫函数与状态和时间有关，表示为V ( x , t )， 如果不显含时间 t ，则表示为V ( x ) 。
（一）预备知识
1．标量函数的定号性 设 V ( x ) 为关于n维向量 x 的标量函数，并且在 x 0 处，有V(x) 0，
则对于任意的非零向量 x 0 ，有：
① 若V(x) 0 ，V ( x ) 为正定； ② 若V(x) 0 ，V ( x ) 为正半定；
条件（2）表明该“能量”函数随着系统的运动不断衰减， 条件（3）表示了满足渐近稳定的条件可扩展至整个状态空间。
2．渐近稳定判定定理2 ：
系统及平衡状态同上，如果 V ( x , t ) 满足条件：
（1） V ( x , t ) 为正定；
（2） V&( x,t) 为负半定，但它在非零解运动轨线上不恒为零，即 对于 x 0 有V&(x,t) 0 ； 则系统的平衡状态x e 0 是渐近稳定的。同样，如果还满足
3) A的特征值的实部有一部分为0，其它均具负实部，非线性系统 在x e 0 的稳定性不能得出明确结论，而取决于 f ( x ) 的高阶导数 项。一般可通过其它方法（如找合适的Lyapunov函数)确定其稳 定性。
李雅普诺夫第一法需要求出系统的全部特征值，这对于高阶系统 存在一定的困难，经典控制理论中针对线性定常系统提出了一些有 效的工程方法，可视为该法在线性定常系统中的工程应用。
Lyapunov指数的计算方法
总结Lyapunov 指数的计算方法近期为了把计算LE的一些问题弄清楚,看了有7～9本书下面以吕金虎混沌时间序列分析及其应用、马军海复杂非线性系统的重构技术为主线,把目前已有的LE计算方法做一个汇总1. 关于连续系统Lyapunov指数的计算方法连续系统LE的计算方法主要有定义方法、Jacobian方法、QR分解方法、奇异值分解方法,或者通过求解系统的微分方程,得到微分方程解的时间序列,然后利用时间序列即离散系统的LE求解方法来计算得到;关于连续系统LE的计算,主要以定义方法、Jacobian方法做主要介绍内容;1定义法关于定义法求解的程序,和matlab板块的“连续系统LE求解程序”差不多;以Rossler 系统为例Rossler系统微分方程定义程序function dX = Rossler_lyt,X% Rossler吸引子,用来计算Lyapunov指数% a=,b=,c=% dx/dt = -y-z,% dy/dt = x+ay,% dz/dt = b+zx-c,a = ;b = ;c = ;x=X1; y=X2; z=X3;% Y的三个列向量为相互正交的单位向量Y = X4, X7, X10;X5, X8, X11;X6, X9, X12;% 输出向量的初始化,必不可少dX = zeros12,1;% Rossler吸引子dX1 = -y-z;dX2 = x+ay;dX3 = b+zx-c;% Rossler吸引子的Jacobi矩阵Jaco = 0 -1 -1;1 a 0;z 0 x-c;dX4:12 = JacoY;求解LE代码：% 计算Rossler吸引子的Lyapunov指数clear;yinit = 1,1,1;orthmatrix = 1 0 0;0 1 0;0 0 1;a = ;b = ;c = ;y = zeros12,1;% 初始化输入y1:3 = yinit;y4:12 = orthmatrix;tstart = 0; % 时间初始值tstep = 1e-3; % 时间步长wholetimes = 1e5; % 总的循环次数steps = 10; % 每次演化的步数iteratetimes = wholetimes/steps; % 演化的次数mod = zeros3,1;lp = zeros3,1;% 初始化三个Lyapunov指数Lyapunov1 = zerositeratetimes,1;Lyapunov2 = zerositeratetimes,1;Lyapunov3 = zerositeratetimes,1;for i=1:iteratetimestspan = tstart:tstep:tstart + tstepsteps;T,Y = ode45'Rossler_ly', tspan, y;% 取积分得到的最后一个时刻的值y = YsizeY,1,:;% 重新定义起始时刻tstart = tstart + tstepsteps;y0 = y4 y7 y10;y5 y8 y11;y6 y9 y12;%正交化y0 = ThreeGSy0;% 取三个向量的模mod1 = sqrty0:,1'y0:,1;mod2 = sqrty0:,2'y0:,2;mod3 = sqrty0:,3'y0:,3;y0:,1 = y0:,1/mod1;y0:,2 = y0:,2/mod2;y0:,3 = y0:,3/mod3;lp = lp+logabsmod;%三个Lyapunov指数Lyapunov1i = lp1/tstart;Lyapunov2i = lp2/tstart;Lyapunov3i = lp3/tstart;y4:12 = y0';end% 作Lyapunov指数谱图i = 1:iteratetimes;ploti,Lyapunov1,i,Lyapunov2,i,Lyapunov3程序中用到的ThreeGS程序如下：%G-S正交化function A = ThreeGSV % V 为33向量v1 = V:,1;v2 = V:,2;v3 = V:,3;a1 = zeros3,1;a2 = zeros3,1;a3 = zeros3,1;a1 = v1;a2 = v2-a1'v2/a1'a1a1;a3 = v3-a1'v3/a1'a1a1-a2'v3/a2'a2a2;A = a1,a2,a3;计算得到的Rossler系统的LE为————Wolf文章中计算得到的Rossler系统的LE为———— 0需要注意的是——定义法求解的精度有限,对有些系统的计算往往出现计果和理论值有偏差的现象;正交化程序可以根据上面的扩展到NN向量,这里就不加以说明了,对matlab用户来说应该还是比较简单的2Jacobian方法通过资料检索,发现论坛中用的较多的LET工具箱的算法原理就是Jacobian方法;基本原理就是首先求解出连续系统微分方程的近似解,然后对系统的Jacobian矩阵进行QR分解,计算Jacobian矩阵特征值的乘积,最后计算出LE和分数维;经过计算也证明了这种方法精度较高,对目前常见的混沌系统,如Lorenz、Henon、Duffing等的Lyapunov指数的计算精度都很高,而且程序编写有一定的规范,个人很推荐使用;虽然我自己要做的系统并不适用LET工具箱可以在网络上找到,这里就不列出了关于LET工具箱如果有问题,欢迎加入本帖讨论Jacobian法求解Lyapunov指数.JPG对离散动力系统,或者说是非线性时间序列,往往不需要计算出所有的Lyapunov指数,通常只需计算出其最大的Lyapunov指数即可;“1983年,格里波基证明了只要最大Lyapunov指数大于零,就可以肯定混沌的存在”;目前常用的计算混沌序列最大Lyapunov指数的方法主要有一下几种：1由定义法延伸的Nicolis方法2Jacobian方法3Wolf方法4P－范数方法5小数据量方法其中以Wolf方法和小数据量方法应用最为广泛,也最为普遍;下面对Nicolis方法、Wolf方法以及小数据量方法作一一介绍;1Nicolis方法这种方法和连续系统的定义方法类似,而且目前应用很有限制,因此只对其理论进行介绍,编程应用方面就省略了Nicolis方法求最大2Wolf方法Wolf方法求最大Wolf方法的Matlab程序如下：function lambda_1=lyapunov_wolfdata,N,m,tau,P% 该函数用来计算时间序列的最大Lyapunov 指数--Wolf 方法% m: 嵌入维数一般选大于等于10% tau:时间延迟一般选与周期相当,如我选2000% data:时间序列可以选1000；% N:时间序列长度满足公式：M=N-m-1tau=24000-10-11000=5000% P:时间序列的平均周期,选择演化相点距当前点的位置差,即若当前相点为I,则演化相点只能在|I－J|>P的相点中搜寻P=周期=600% lambda_1: 返回最大lyapunov指数值min_point=1 ; %&&要求最少搜索到的点数MAX_CISHU=5 ; %&&最大增加搜索范围次数%FLYINGHAWK% 求最大、最小和平均相点距离max_d = 0; %最大相点距离min_d = +100; %最小相点距离avg_dd = 0;Y=reconstitutiondata,N,m,tau; %相空间重构可将此段程序加到整个程序中,在时间循环内,可以保存时间序列的地方;见完整程序;M=N-m-1tau; %重构相空间中相点的个数for i = 1 : M-1for j = i+1 : Md = 0;for k = 1 : md = d + Yk,i-Yk,jYk,i-Yk,j;endd = sqrtd;if max_d < dmax_d = d;endif min_d > dmin_d = d;endavg_dd = avg_dd + d;endendavg_d = 2avg_dd/MM-1; %平均相点距离dlt_eps = avg_d - min_d ; %若在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 = +100; %第i个相点到其最近距离点的距离Loc_DK = 2; %第i个相点对应的最近距离点的下标for i = P+1:M-1 %限制短暂分离,从点P+1开始搜索d = 0;for k = 1 : md = d + Yk,i-Yk,1Yk,i-Yk,1;endd = sqrtd;if d < DK & d > min_epsDK = 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个log2DK1/DK的累计和用于求i点的lambda值sum_lmd = 0 ; % 存放前i个log2DK1/DK的累计和for i = 2 : M-1 % 计算演化距离DK1 = 0;for k = 1 : mDK1 = DK1 + Yk,i-Yk,Loc_DK+1Yk,i-Yk,Loc_DK+1;endDK1 = sqrtDK1;old_Loc_DK = Loc_DK ; % 保存原最近位置相点old_DK=DK;% 计算前i个log2DK1/DK的累计和以及保存i点的李氏指数if DK1 ~= 0& DK ~= 0sum_lmd = sum_lmd + logDK1/DK /log2;endlmdi-1 = sum_lmd/i-1; 此处可以保存不同相点i对应的李氏指数,见完整程序;; % 以下寻找i点的最短距离:要求距离在指定距离范围内尽量短,与DK1的角度最小 point_num = 0; % &&在指定距离范围内找到的候选相点的个数cos_sita = 0 ; %&&夹角余弦的比较初值——要求一定是锐角zjfwcs=0 ; %&&增加范围次数while point_num == 0% 搜索相点for j = 1 : M-1if absj-i <=P-1 %&&候选点距当前点太近,跳过continue;end%计算候选点与当前点的距离dnew = 0;for k = 1 : mdnew = dnew + Yk,i-Yk,jYk,i-Yk,j;enddnew = sqrtdnew;if dnew < min_eps| dnew > max_eps %&&不在距离范围,跳过continue;end%计算夹角余弦及比较DOT = 0;for k = 1 : mDOT = DOT+Yk,i-Yk,jYk,i-Yk,old_Loc_DK+1;endCTH = DOT/dnewDK1;if acosCTH > 3./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 = +100;for ii = 1 : M-1if absi-ii <= P-1 %&&候选点距当前点太近,跳过continue;endd = 0;for k = 1 : md = d + Yk,i-Yk,iiYk,i-Yk,ii;endd = sqrtd;if d < DK & d > min_epsDK = d;Loc_DK = ii;endendbreak;endpoint_num = 0 ; %&&扩大距离范围后重新搜索cos_sita = 0;endendend%取平均得到最大李雅普诺夫指数此处只有一个值,若为正说明体系是混沌的,若为负则说明体系是周期性的确定性运动lambda_1=sumlmd/lengthlmd;程序中用到的reconstitution函数如下：此段程序可直接放在时间循环内部,即可以保存时间序列的地方;见完整程序范例;function X=reconstitutiondata,N,m,tau%该函数用来重构相空间% m 为嵌入空间维数% tau 为时间延迟% data 为输入时间序列% N 为时间序列长度% X 为输出,是mn维矩阵M=N-m-1tau; %相空间中点的个数for j=1:M %相空间重构for i=1:mXi,j=datai-1tau+j;endend这里声明一下,这些程序并非我自己编写的,均是转载,其使用我已经验证过,绝对可以运行3小数据量方法说小数据量方法是目前最实用、应用最广泛的方法应该不为过吧,呵呵小数据量方法求最大Lyapunov指数.JPG上面两种方法,即Wolf方法和小数据量方法,在计算LE之前,都要求对时间序列进行重构相空间,重构相空间的优良对于最大LE的计算精度影响非常大因此重构相空间的几个参数的确定就非常重要;1时间延迟主要推荐两种方法——自相关函数法、C－C方法自相关函数法——对一个混沌时间序列,可以先写出其自相关函数,然后作出自相关函数关于时间t的函数图像;根据数值试验结果,当自相关函数下降到初始值的1－1/e时,所得的时间t即为重构相空间的时间延迟;C－C方法——可以同时计算出时间延迟和时间窗口,个人推荐使用这种方法2平均周期平均周期的计算可以采用FFT方法;在matlab帮助中有一个太阳黑子的例子,现摘录如下：load %装载数据文件year = sunspot:,1; %读取年份信息wolfer = sunspot:,2; %读取黑子活动数据plotyear,wolfer %绘制原始数据图title'Sunspot Data'Y = fftwolfer; %快速FFT变换N = lengthY; %FFT变换后数据长度Y1 = ; %去掉Y的第一个数据,它是所有数据的和power = absY1:N/2.^2; %求功率谱nyquist = 1/2;freq = 1:N/2/N/2nyquist; %求频率plotfreq,power, grid on %绘制功率谱图xlabel'cycles/year'title'Periodogram'period = 1./freq; %年份周期plotperiod,power, axis0 40 0 2e7, grid on ％绘制年份－功率谱曲线ylabel'Power'xlabel'PeriodYears/Cycle'mp,index = maxpower; %求最高谱线所对应的年份下标periodindex %由下标求出平均周期3嵌入维数目前嵌入维数的主要计算方法是采用Grassberger和Procaccia提出的G－P算法计算出序列的关联维数d,然后利用嵌入维数m>=2d+1,选取合适的嵌入维数;G—P算法程序如下：function ln_r,ln_C=G_Pdata,N,tau,min_m,max_m,ss% the function is used to calculate correlation dimention with G-P algorithm% 计算关联维数的G－P算法% data:the time series 时间序列% N: the length of the time series 时间序列长度% tau: the time delay 时间延迟% min_m:the least embedded dimention m 最小的嵌入维数% max_m:the largest embedded dimention m 最大的嵌入维数% ss:the stepsize of r r的步长%skyhawkfor m=min_m:max_mY=reconstitutiondata,N,m,tau;%reconstitute state spaceM=N-m-1tau;%the number of points in state spacefor i=1:M-1for j=i+1:Mdi,j=maxabsY:,i-Y:,j;%calculate the distance of each twoend %points in state space 计算状态空间中每两点之间的距离endmax_d=maxmaxd;%the max distance of all points 得到所有点之间的最大距离d1,1=max_d;min_d=minmind;%the min distance of all points 得到所有点间的最短距离 delt=max_d-min_d/ss;%the stepsize of r 得到r的步长for k=1:ssr=min_d+kdelt;Ck=correlation_integralY,M,r;%calculate the correlation integralln_Cm,k=logCk;%lnCrln_rm,k=logr;%lnrfprintf'%d/%d/%d/%d\n',k,ss,m,max_m;endplotln_rm,:,ln_Cm,:;hold on;endfid=fopen'','w';fprintffid,'% %\n',ln_r;fclosefid;fid = fopen'','w';fprintffid,'% %\n',ln_C;fclosefid;程序中的correlation_integral函数如下：function C_I=correlation_integralX,M,r%the function is used to calculate correlation integral%C_I:the value of the correlation integral%X:the reconstituted state space,M is a mM matrix%m:the embedding demention%M:M is the number of embedded points in m-dimensional sapce%r:the radius of the Heaviside function,sigma/2<r<2sigma%calculate the sum of all the values of Heaviside%skyhawksum_H=0;for i=1:M% fprintf'%d/%d\n',i,M;for j=i+1:Md=normX:,i-X:,j,inf;%calculat the distances of each two points in matris M with sup-normsita=heavisider,d;%calculate the value of the heaviside functionsum_H=sum_H+sita;endendC_I=2sum_H/MM-1;%the value of correlation integral发布于2007-08-03 20:35:01感谢octopussheng的总结这是现有的方法的一个总结,不知道你对这些方法有些什么体会,或者说他们的局限,现在还有作这方面的改进的吗先总结总结,具体的比较会逐步贴出来的,而且这也需要大家的一起努力啊光是我在这里说,也只是我自己的一点体会,希望用过这些方法的一起来参与这个总结工作以上的各种方法在实际应用的时候要根据具体情况来选择;一般地,如果已知系统方程当然系统不能太过复杂时,则计算Lyapunov指数采用定义法、Jacobian方法要精确、简单些而如果系统方程比较复杂如超维系统、或者为一时间序列时,则推荐采样Wolf 方法、小数据量方法;Wolf方法的特点是时间序列无噪声,空间中小向量的演变高度非线性,而Jacobian 方法则是噪声大,空间中小向量的演变接近线性;小数据量方法的优点在于：1对小数据组的计算可靠；2计算量较小,比wolf方法快很多；3编程、操作较为容易;而关于时间延迟、嵌入维数、平均周期的确定,还是推荐C－C方法和G－P算法,结果更为可靠一些你这已经做得很好了;这个具体的应用我也是做得很少的,还有请注意一下matlab版中的一个帖子, 早就注意了,我这里面一些代码就是从那边套用的试过了吗都没有问题吧计算起来,原来论坛里面一些程序不能运行,套用,那你就要注意,在这里重写就要有我们的特色;不能够是简单的重请详细读帖子,上面有说明,当然不是简单的重复啦,我这里是理论结合实践,呵呵程序都测试能算的,不过对不同的系统,精度不敢保证的太高,呵呵,只要能运行就可以参照自己来编写精度可以根据自己需要来选择呵呵,确实,保证能够使用希望这样总结一下对大家有所帮助吧非常有价值;谢谢。
分数阶duffing李雅普诺夫指数
分数阶Duffing 方程是经典的非线性动力系统方程之一，它是对Duffing 方程的一种扩展，将二阶导数改为分数阶导数。
分数阶Duffing 方程在许多领域，如力学、控制系统和电子工程中都有重要的应用。
利用李雅普诺夫指数（Lyapunov exponent），可以评估一个动力系统的混沌性质。
李雅普诺夫指数是描述非线性系统稳定性和混沌行为的重要工具。
对于分数阶Duffing 方程，李雅普诺夫指数可以用来判断系统的混沌性质。
具体来说，当系统的李雅普诺夫指数是正值时，表示系统是混沌的；当指数为负值时，表示系统是稳定的；当指数等于零时，表示系统是临界的。
通过计算系统的李雅普诺夫指数，可以了解分数阶Duffing 方程的动力学行为，包括混沌运动、吸引子的结构等。
李雅普诺夫指数的计算通常需要借助数值方法和计算机仿真进行。
需要注意的是，详细的数学推导和计算方法超出了本回答的范围。
如果你对分数阶Duffing 方程和李雅普诺夫指数有更多研究兴趣，建议参考相关的数学和动力系统的文献。
李雅普诺夫第一方法和第二方法
李雅普诺夫第一方法和第二方法李雅普诺夫第一方法，也称为不动点迭代法或迭代法。
这种方法基于一个重要的定理，即如果函数g(x)在给定区间[a,b]上连续，并且对于这个区间上任意的x，都有g(x)在[a,b]内的闭区间上有g(x)∈[a,b]，那么在这个区间上存在唯一的不动点c，满足c=g(c)。
所谓不动点，即函数g(x)的值等于其自变量x的值。
李雅普诺夫第一方法的核心思想是通过迭代计算不动点c的近似值，即x_n+1=g(x_n)，不断逼近真实的解。
迭代的过程中，从一个初始值x_0开始，通过将x_0代入g(x)得到新的近似值x_1，再将x_1代入g(x)得到新的近似值x_2，以此类推，直到达到预定的精度要求或者迭代次数时停止迭代。
具体地，李雅普诺夫第一方法的步骤如下：1.选择一个合适的初始值x_0。
2.根据迭代公式x_n+1=g(x_n)计算x_1,x_2,...,x_n，直到满足停止条件。
3.如果满足停止条件，则该迭代过程收敛于不动点c，并且c是非线性方程的实根的一个近似解。
李雅普诺夫第二方法，也称为牛顿迭代法，是一种更加高效的求解非线性方程的数值计算方法。
牛顿迭代法的基本思想是利用函数的局部线性近似来逼近非线性方程的根。
根据泰勒级数展开，可以将非线性方程f(x)=0在一些近似解x_n的邻域中展开成一个一次项和高阶项的级数。
利用一次项的值来逼近非线性方程的解，可以得到迭代公式x_n+1=x_n-f(x_n)/f'(x_n)，其中f'(x_n)代表函数f(x)在x_n处的导数值。
具体地，李雅普诺夫第二方法的步骤如下：1.选择一个合适的初始值x_0。
2.根据迭代公式x_n+1=x_n-f(x_n)/f'(x_n)计算x_1,x_2,...,x_n，直到满足停止条件。
3.如果满足停止条件，则该迭代过程收敛于非线性方程的实根，并且x_n是方程的一个近似解。
与李雅普诺夫第一方法相比，李雅普诺夫第二方法的优点在于收敛速度更快。
最大李雅谱诺夫数数值计算
最大李雅谱诺夫数数值计算————————————————————————————————作者：————————————————————————————————日期：最大李雅谱诺夫指数的数值计算------------------------------------------------------------------------最大李雅谱诺夫指数的数值计算J.C.SprottDepartment of Physics, University of Wisconsin, Madison, WI 53706, USAOctober 15, 1997(Revised August 31, 2004)袁来翻译/上海大学/05/03/22原文出处:/chaos/lyapexp.htm通常通过计算最大的李雅谱诺夫指数来检验系统是否混沌，一个正的最大李雅谱诺夫指数标致了该系统是否混沌，当你已经获得已知产生混沌的方程，这是相对容易做的。
当你仅仅只有试验数据记录时，这种计算的困难程度是几乎不可能的，这里我们不考虑这种情形。
计算的一般思路是跟踪两个较近的轨道，计算它们分离的平均对数率，无论它们分离得多远，其中之一的轨道将最终回到另一条轨道的附近。
在每次迭代过程中，一程序将实现这个工作，完整的程序过程如下：1：在吸引域内，开始于任何初值条件。
最好开始于吸引子上的一个已知点，这种情形，第2步可以省略掉。
2：迭代直至轨道落到吸引子上。
这需要自己的判断或通过研究得到的有关该系统的一些知识，对绝大多数系统，只要迭代几百次，我们就认为足够了。
通常情况下，除非你选的初值点刚好靠近分岔点，要不在几百次后它不会远离吸引子的。
3：选择（几乎任意的）靠近的点（设两点的距离为d0）适合的d0的选择是至少要比计算机中的浮点数的1000倍大，如：在（8字节）的双精度的（对这些计算的最低推荐），变量有52－比特的尾数，故，精度是：2-52=2.22x10-16，因此，我们只要有d0=10-12 就足够了。
-Lyapunov指数的计算方法
【总结】Lyapunov指数的计算方法非线性理论近期为了把计算LE的一些问题弄清楚，看了有7～9本书！下面以吕金虎《混沌时间序列分析及其应用》、马军海《复杂非线性系统的重构技术》为主线，把目前已有的LE计算方法做一个汇总！1. 关于连续系统Lyapunov指数的计算方法连续系统LE的计算方法主要有定义方法、Jacobian方法、QR分解方法、奇异值分解方法，或者通过求解系统的微分方程，得到微分方程解的时间序列，然后利用时间序列（即离散系统）的LE求解方法来计算得到。
关于连续系统LE的计算，主要以定义方法、Jacobian方法做主要介绍内容。
（1）定义法定义法求解Lyapunov指数.JPG关于定义法求解的程序，和matlab板块的“连续系统LE求解程序”差不多。
以Rossler系统为例Rossler系统微分方程定义程序function dX = Rossler_ly(t,X)% Rossler吸引子，用来计算Lyapunov指数% a=0.15,b=0.20,c=10.0% dx/dt = -y-z,% dy/dt = x+ay,% dz/dt = b+z(x-c),a = 0.15;b = 0.20;c = 10.0;x=X(1); y=X(2); z=X(3);% Y的三个列向量为相互正交的单位向量Y = [X(4), X(7), X(10);X(5), X(8), X(11);X(6), X(9), X(12)];% 输出向量的初始化，必不可少dX = zeros(12,1);% Rossler吸引子dX(1) = -y-z;dX(2) = x+a*y;dX(3) = b+z*(x-c);% Rossler吸引子的Jacobi矩阵Jaco = [0 -1 -1;1 a 0;z 0 x-c];dX(4:12) = Jaco*Y;求解LE代码：% 计算Rossler吸引子的Lyapunov指数clear;yinit = [1,1,1];orthmatrix = [1 0 0;0 1 0;0 0 1];a = 0.15;b = 0.20;c = 10.0;y = zeros(12,1);% 初始化输入y(1:3) = yinit;y(4:12) = orthmatrix;tstart = 0; % 时间初始值tstep = 1e-3; % 时间步长wholetimes = 1e5; % 总的循环次数steps = 10; % 每次演化的步数iteratetimes = wholetimes/steps; % 演化的次数mod = zeros(3,1);lp = zeros(3,1);% 初始化三个Lyapunov指数Lyapunov1 = zeros(iteratetimes,1); Lyapunov2 = zeros(iteratetimes,1); Lyapunov3 = zeros(iteratetimes,1);for i=1:iteratetimestspan = tstart:tstep:(tstart + tstep*steps); [T,Y] = ode45('Rossler_ly', tspan, y);% 取积分得到的最后一个时刻的值y = Y(size(Y,1),:);% 重新定义起始时刻tstart = tstart + tstep*steps;y0 = [y(4) y(7) y(10);y(5) y(8) y(11);y(6) y(9) y(12)];%正交化y0 = ThreeGS(y0);% 取三个向量的模mod(1) = sqrt(y0(:,1)'*y0(:,1));mod(2) = sqrt(y0(:,2)'*y0(:,2));mod(3) = sqrt(y0(:,3)'*y0(:,3));y0(:,1) = y0(:,1)/mod(1);y0(:,2) = y0(:,2)/mod(2);y0(:,3) = y0(:,3)/mod(3);lp = lp+log(abs(mod));%三个Lyapunov指数Lyapunov1(i) = lp(1)/(tstart);Lyapunov2(i) = lp(2)/(tstart);Lyapunov3(i) = lp(3)/(tstart);y(4:12) = y0';end% 作Lyapunov指数谱图i = 1:iteratetimes;plot(i,Lyapunov1,i,Lyapunov2,i,Lyapunov3)程序中用到的ThreeGS程序如下：%G-S正交化function A = ThreeGS(V) % V 为3*3向量v1 = V(:,1);v2 = V(:,2);v3 = V(:,3);a1 = zeros(3,1);a2 = zeros(3,1);a3 = zeros(3,1);a1 = v1;a2 = v2-((a1'*v2)/(a1'*a1))*a1;a3 = v3-((a1'*v3)/(a1'*a1))*a1-((a2'*v3)/(a2'*a2))*a2;A = [a1,a2,a3];计算得到的Rossler系统的LE为———— 0.063231 0.092635 -9.8924Wolf文章中计算得到的Rossler系统的LE为————0.09 0 -9.77需要注意的是——定义法求解的精度有限，对有些系统的计算往往出现计果和理论值有偏差的现象。