李雅普诺夫指数
一、李雅普诺夫指数的提出与历史
1961年冬季的一天，为了考察一条更长的序列，洛伦兹走了一条捷径。他在进行天气模式计算时没有从头开始运行，而是从中途开始。作为计算的初值，他直接输入了上次运算的输出结果，然后他穿过大厅下楼，清净的去喝一杯咖啡。一个小时之后他回来时，看到了出乎意料的事。从几乎相同出发点开始，洛伦兹看到他的计算机产生的天气模式差别愈来愈大，终至毫无相似之处。就是这件事播下了一门新科学的种子。
稳定体系的相轨线相应于趋向某个平衡点，如果出现越来越远离平衡点，则系统是不稳定的。系统只要有一个正值就会出现混沌运动。判断一个非线性体统是否存在混沌运动时，需要检查它的李雅普诺夫指数λ是否为正值。
二、李雅普诺夫指数的定义
Lyapunov 指数是描述时序数据所生成的相空间中两个极其相近的初值所产生的轨道，随时间推移按指数方式分散或收敛的平均变化率。任何一个系统，只要有一个Lyapunov 大于零，就认为该系统为混沌系统。
李雅普诺夫指数是指在相空间中相互靠近的两条轨线随着时间的推移，按指数分离或聚合的平均变化速率。李雅普诺夫指数的定义: 首先考虑一维映射假设初始位置附近有一点，则经过一次迭代后，这两点之间的距离为：
（1）
并利用微分中值定理有：
（2）
n次迭代后，并利用微分中值定理，这两点之间的距离为：
（3）
由（3）式可得：（4）又由复合函数的微分规则有：
其中
那么式（4）就变为：
（5）
则称（6）为Lyapunov指数。
一维映射就对应一个李雅普诺夫指数，而且当时，该系统具有混沌特性。当时，
对应着分岔点或系统的周期解，既系统出现周期现象。时，系统有稳定的不动点，
即此时对应的是一个点。而对于多维系统则有多个李雅普诺夫指数。
Lyapunov 特性指数沿某一方向取值的正负和大小表示长时间系统在吸引子中相邻轨线
沿该方向平均发散或收敛i的快慢程度，仅从数学角度考虑，Lyapunov特性指数无量纲。n 维系统具有 n 个 Lyapunov 特性指数，形成指数谱。其中数值最大的被称为最大 Lyapunov 特性指数。最大 Lyapunov 指数定义为
其中，表示时刻最邻近零点间的距离；M为计算总步数。最大 Lyapunov指数不仅是区别混沌吸引子的重要指标，也是混沌系统对于初始值敏感性的定量描述。其中一维系统只有一个指数，二维系统有两个指数来表征。在实际计算中,要计算所有的Lyapunov指数,计算量较大,尤其当系统维数L较大时更为突出.所以注意力集中在计算系统的最大Lyapunov指数λm上.
三、李雅普诺夫指数的物理意义
系统的Lyapunov指数谱可有效地表征变量随时间演化时，系统对初值的敏感性。指数小于零说明体系的相体积在该方向上是收缩的，此方向的运动是稳定的；
而正的指数值则表明了体系的相体积在该方向上不断膨胀和折叠，以致吸引子中本来邻近的轨线变得越来越不相关，从而使初态对任何不确定性的系统的长期行为成为不可预测，即所谓的初值敏感性。
进一步意义：设某一系统的指数谱为
（从大到小排列），若该系统具有混沌吸引子，则必须同时满足以下条件
（1）至少存在一个正李雅普诺夫指数
（2）至少存在某一指数为0
（3）指数谱之和为负。
混沌运动的基本特点是运动状态对初始条件的高度敏感性。两个极为靠近的初值所产生的轨道,随时间推移按指数形式分离,Lyapunov指数是定量描述这一现象的量。
对所讨论的Duffing振子,若它的Lyapunov指数均小于零, :若存在一个Lyapunov特性指数大于零,就说明系统是处于混沌状态。这种判别方法计算简单,物理意义明确,误差小。四、计算此指数的几种方法
用Logistic映射产生的模拟时间序列数据，采用两种从实验数据时间序列恢复动力学的方法，计算混沌吸引子的Lyapunov指数。
一种方法是S.J.Chang和J.Wright提出的混合嫡法〔8一〕，这种方法特别适合一维的实验系统。另一种方法是A.wolf提出的重构吸引子法〔7〕，这种方法可以推广到相空间维数及
动力学规律都不知道的更普遍的实验系统，在原则上可以计算系统的全部正Lyapunov指数谱。
具体的方法还有
1.从动力学规律计算Lyapunov指数
2. Chang一Wright混合嫡法
Chang一wright混合嫡法仅适用于可化为一维凸映射的情况。一般来说，总点数N及盒子数凡越大，所得结果越精确。
3. Wolf重构法
Wolf重构法以Takens的延迟坐标重构相空间技术为基础，对于一个由观测得到的实验数据时间序列x(t)以延迟坐标重构m维相空间中的一条轨道，
计算最大Lyapunov指数的Wolf程序，一般适合于嵌入维m>1的重构吸引子的时间序列。
常微分大作业--李雅普诺夫稳定性 11091059洪一洲 从19世纪末以来，李雅普诺夫稳定性理论一直指导着关于稳定性的研究和应用。不少学者遵循李雅普诺夫所开辟的研究路线对第二方法作了一些新的发展。一方面，李雅普诺夫第二方法被推广到研究一般系统的稳定性。例如，1957年，В．И．祖博夫将李雅普诺夫方法用于研究度量空间中不变集合的稳定性。随后，J.P.拉萨尔等又对各种形式抽象系统的李雅普诺夫稳定性进行了研究。在这些研究中，系统的描述不限于微分方程或差分方程，运动平衡状态已采用不变集合表示，李雅普诺夫函数是在更一般意义下定义的。1967年，D.布肖对表征在集合与映射水平上的系统建立了李雅普诺夫第二方法。这时，李雅普诺夫函数已不在实数域上取值，而是在有序定义的半格上取值。另一方面，李雅普诺夫第二方法被用于研究大系统或多级系统的稳定性。此时，李雅普诺夫函数被推广为向量形式，称为向量李雅普诺夫函数。用这种方法可建立大系统稳定性的充分条件。 1.李雅普诺夫稳定性概念 忽略输入后，非线性时变系统的状态方程如下 ),(t x f x = （1） 式中，x 为n 维状态向量；t 为时间变量；),(t x f 为n 维函数，其展开式为 12(,,,,)i i n x f x x x t = n i ,,1 = 假定方程的解为 ),;(00t x t x ，x 0和t 0 分别为初始状态向量和初始时刻，0000),;(x t x t x =。 平衡状态 如果对于所有t ，满足 0),(==t x f x e e （2） 的状态x e 称为平衡状态（又称为平衡点）。平衡状态的各分量不再随时间变化。若已知状态 方程，令0=x 所求得的解x ，便是平衡状态。 对于线性定常系统Ax x = ，其平衡状态满足0=e Ax ，如果A 非奇异，系统只有惟一的零解，即存在一个位于状态空间原点的平衡状态。至于非线性系统，0),(=t x f e 的解可能有多个，由系统状态方程决定。 控制系统李雅普诺夫意义下的稳定性是关于平衡状态的稳定性，反映了系统在平衡状态附近的动态行为。鉴于实际线性系统只有一个平衡状态，平衡状态的稳定性能够表征整个系统的稳定性。对于具有多个平衡状态的非线性系统来说，由于各平衡状态的稳定性一般并不相同，故需逐个加以考虑，还需结合具体初始条件下的系统运动轨迹来考虑。
Lorenz系统 文档分两个文件方程m文件和计算L指数m文件分开写，复制粘贴即可运行matlab2012a，改写方程文件和参数即可算自己的系统，其中最大L指数用的是经典的柏内庭（G.Benettin）计算方法，准确快速无误！附计算结果图！！ 方程m文件： function dX = Loren(t,X) global a; %变量不放入参数表中 global b; global c; x=X(1); y=X(2); z=X(3); % Y的三个列向量为相互正交的单位向量 %输出向量的初始化 dX = zeros(6,1); % Lorenz吸引子 dX(1)=a*(y-x); dX(2)=x*(b-z)-y; dX(3)=x*y-c*z; end 计算最大L指数文件 Z=[]; global a; global b;
global c; a=10; c=8/3; d0=1e-7; for b=linspace(0,500,500) lsum=0; x=1;y=1;z=1; x1=1;y1=1;z1=1+d0; for i=1:100 [T1,Y1]=ode45('Loren',1,[x;y;z;16;b;4]); [T2,Y2]=ode45('Loren',1,[x1;y1;z1;16;b;4]); n1=length(Y1);n2=length(Y2); x=Y1(n1,1);y=Y1(n1,2);z=Y1(n1,3); x1=Y2(n2,1);y1=Y2(n2,2);z1=Y2(n2,3); d1=sqrt((x-x1)^2+(y-y1)^2+(z-z1)^2); x1=x+(d0/d1)*(x1-x); y1=y+(d0/d1)*(y1-y); z1=z+(d0/d1)*(z1-z); if i>50 lsum=lsum+log(d1/d0); end
1. 关于连续系统Lyapunov指数的计算方法 连续系统LE的计算方法主要有定义方法、Jacobian方法、QR分解方法、奇异值分解方法，或者通过求解系统的微分方程，得到微分方程解的时间序列，然后利用时间序列（即离散系统）的LE求解方法来计算得到。 关于连续系统LE的计算，主要以定义方法、Jacobian方法做主要介绍内容。 （1）定义法
关于定义法求解的程序，和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;
多变量时间序列最大李雅普诺夫指数的计算 作者：卢山， 王海燕， Lu Shan， Wang Hai-Yan 作者单位：东南大学经济管理学院,南京,210096 刊名： 物理学报 英文刊名：ACTA PHYSICA SINICA 年，卷(期)：2006,55(2) 被引用次数：5次 参考文献(14条) 1.Liu W D;Ren K F;Meunier S查看详情 2003 2.徐莉梅;胡岗;史朋亮查看详情 2000 3.游荣义;陈忠;徐慎初基于小波变换的混沌信号相空间重构研究[期刊论文]-物理学报 2004(9) 4.肖方红;阎桂荣;韩宇航混沌时序相空间重构参数确定的信息论方法[期刊论文]-物理学报 2005(2) 5.Cao L Y;Mees A;Judd K查看详情 1998 6.Boccaletti S;Valladares D L;Louis M查看详情 2002 7.Zhang H;Ma X K;Yang Y查看详情[期刊论文]-Chin Phys 2005 8.Rosenblum M G;Pikovsky A S;Kurths J查看详情 1996 9.王海燕;盛昭瀚;张进多变量时间序列复杂系统的相空间重构[期刊论文]-东南大学学报(自然科学版) 2003(1) 10.杨绍清;贾传荧两种实用的相空间重构方法[期刊论文]-物理学报 2002(11) 11.Rosenstei MT;Collins J J;De L C J查看详情 1993 12.Zou Y L;ZhuJ;Chen G R查看详情[期刊论文]-Chin Phys 2005 13.谢勇;徐健学;杨红军皮层脑电时间序列的相空间重构及非线性特征量的提取[期刊论文]-物理学报 2002(2) 14.Abarbanel H Analysis of Observed Chaotic Data 1996 引证文献(5条) 1.聂春燕.王祝文.李泽.崔炳民储集层测井信号的非线性混沌特性[期刊论文]-吉林大学学报（地球科学版）2011(1) 2.刘立霞.苗海峰多变量时间序列最大Lyapunov指数的噪声估计[期刊论文]-计算机工程与应用 2010(22) 3.徐威.郭静波混沌直扩信号检测的最大Lyapunov指数方法[期刊论文]-应用科学学报 2009(2) 4.刘志平.何秀凤.何习平基于多变量最大Lyapunov指数高边坡稳定分区研究[期刊论文]-岩石力学与工程学报2008(z2) 5.赵敏.FAN Yin-hai.孙辉电力推进船舶电力负荷的多变量混沌局部预测[期刊论文]-系统仿真学报 2008(11)本文链接：https://www.doczj.com/doc/3d3463901.html,/Periodical_wlxb200602018.aspx
第六章 李雅普诺夫稳定性分析 在反馈控制系统的分析设计中，系统的稳定性是首先需要考虑的问题之一。因为它关系到系统是否能正常工作。 经典控制理论中已经建立了劳斯判据、Huiwitz 稳定判据、Nquist 判据、对数判据、根轨迹判据等来判断线性定常系统的稳定性，但不适用于非线性和时变系统。分析非线性系统稳定性及自振的描述函数法，则要求系统的线性部分具有良好的滤除谐波的性能；而相平面法则只适合于一阶、二阶非线性系统。 1892年俄国学者李雅普诺夫（Lyapunov ）提出的稳定性理论是确定系统稳定性的更一般的理论，它采用状态向量来描述，不仅适用于单变量、线性、定常系统，还适用于多变量、非线性、时变系统。 §6-1 外部稳定性和内部稳定性 系统的数学模型有输入输出描述（即外部描述）和状态空间描述（即内部描述），相应的稳定性便分为外部稳定性和内部稳定性。 一、外部稳定性 1、定义（外部稳定性）： 若系统对所有有界输入引起的零状态响应的输出是有界的，则称该系统是外部稳定的。 (外部稳定性也称为BIBO （Bounded Input Bounded Output ）稳定性) 说明： （1）所谓有界是指如果一个函数)(t h ，在时间区间],0[∞中，它的幅值不会增至无穷，即存在一个实 常数k ，使得对于所有的[]∞∈0 t ，恒有∞<≤k t h )(成立。 （2）所谓零状态响应，是指零初始状态时非零输入引起的响应。 2、系统外部稳定性判据 线性定常连续系统 ∑),,(C B A 的传递函数矩阵为 Cx y Bu Ax x =+= BU A sI X BU X A sI CX Y BU AX sX 1)()(--==-=+= B A sI C s G 1 )()(--= 当且仅当)(s G 极点都在s 的左半平面内时，系统才是外部稳定（或BIBO 稳定）的。
求解系统的Lyapunov指数谱程序 Lyapunov 指数是描述时序数据所生成的相空间中两个极其相近的初值所产生的轨道，随时间推移按指数方式分散或收敛的平均变化率。任何一个系统，只要有一个Lyapunov 大于零，就认为该系统为混沌系统。 李雅普诺夫指数是指在相空间中相互靠近的两条轨线随着时间的推移，按指数分离或聚合的平均变化速率。 一 chen系统的Lyapunov指数谱 function dX = Chen2(t,X) % Chen吸引子，用来计算Lyapunov指数 % dx/dt=a*(y-x) % dy/dt=(c-a)*x+c*y-x*z % dz/dt=x*y-b*z global a; % 变量不放入参数表中 global b; global c; 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); % Lorenz吸引子 dX(1) = a*(y-x); dX(2) = (c-a)*x+c*y-x*z; dX(3) = x*y-b*z; % Lorenz吸引子的Jacobi矩阵 Jaco = [-a a 0; c-a-z c -x; y x -b]; dX(4:12) = Jaco*Y; Z1=[];
Z2=[]; Z3=[]; global a; global b; global c; b=3;c=28; for a=linspace(32,40,100); y=[1;1;1;1;0;0;0;1;0;0;0;1]; lp=0; for k=1:200 [T,Y] = ode45('Chen2', 1, y); y = Y(size(Y,1),:); y0 = [y(4) y(7) y(10); y(5) y(8) y(11); y(6) y(9) y(12)]; y0=GS(y0); mod(1)=norm(y0(:,1)); mod(2)=norm(y0(:,2)); mod(3)=norm(y0(:,3)); lp = lp+log(abs(mod)); y0(:,1)=y0(:,1)/mod(1); y0(:,2)=y0(:,2)/mod(2); y0(:,3)=y0(:,3)/mod(3); y(4:12) = y0'; end lp=lp/200; Z1=[Z1 lp(1)]; Z2=[Z2 lp(2)]; Z3=[Z3 lp(3)]; end a=linspace(32,40,100); plot(a,Z1,'-',a,Z2,'-',a,Z3,'-'); title('Lyapunov exponents of Chen') xlabel('b=3,c=28,parameter a'),ylabel('lyapunov exponents') grid on
【总结】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 0x-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:iteratetimes tspan = tstart:tstep:(tstart + tstep*steps); [T,Y] = ode45('Rossler_ly', tspan, y); % 取积分得到的最后一个时刻的值 y = Y(size(Y,1),:); % 重新定义起始时刻 tstart = tstart + tstep*steps;
李雅普诺夫稳定性方法 李雅普诺夫第一方法又称间接法，它是通过系统状态方程的解来判断系统的稳定性。如果其解随时间而收敛，则系统稳定；如果其解随时间而发散，则系统不稳定。 李雅普诺夫第二方法又称直接法，它不通过系统状态方程的解来判断系统的稳定性，而是借助李雅普诺夫函数对稳定性作出判断，是从广义能量的观点进行稳定性分析的。例如有阻尼的振动系统能量连续减小（总能量对时间的导数是负定的），系统会逐渐停止在平衡状态，系统是稳定的。由于李雅普诺夫第一方法求解通常很烦琐，因此李雅普诺夫第二方法获得更广泛的应用。李雅普诺夫第二方法的难点在于寻找李雅普诺夫函数。迄今为止，尚没有通用于一切系统的构造李雅普诺夫函数的方法。 对于系统[]t ,f x x = ，平衡状态为,0e =x 满足()0f e =x 。如果存在一个标量函数()x V ,它满足()x V 对所有x 都具有连续的 一阶偏导数；同时满足()x V 是正定的；则 （1）若()x V 沿状态轨迹方向计算的时间导数()dt /)(dV V x x = 为半负定，则平衡状态稳定； （2） 若()x V 为负定，或虽然()x V 为半负定，但对任意初始状 态不恒为零，则平衡状态渐近稳定。进而当∞→∞→)(V x x 时，，则系统大范围渐近稳定； （3） 若()x V 为正定，则平衡状态不稳定。 判断二次型 x x x P )(V τ=的正定性可由赛尔维斯特
（Sylvester ）准则来确定，即正定（记作V(x)>0）的充要条件为P 的所有主子行列式为正。如果P 的所有主子行列式为非负，为正半定（记作V(x)≥0）；如果-V(x)为正定,则V(x)为负定（记作V(x)<0）；如果-V(x)为正半定,则V(x)为负半定（记作V(x)≤0）。 例： []正定。 则)(V 0 1121412110 ,0411 10,010x x x 1121412110x x x )(V 321321x x >---->>----=??? ????????????? 例： )x x (x x x ) x x (x x x 2 2212122221121+--=+-= (0,0）是唯一的平衡状态。设正定的标量函数为 ∞→∞→<+-=+--++-=+=??+??=+=)V(,且当0 )x 2(x )]x (x x x [2x )]x (x x [x 2x x 2x x 2x dt dx x V dt dx x V )(V x x )V(2222122212122221121221122112 2 21x x x x 故系统在坐标原点处为大范围渐近稳定。
李雅普诺夫指数 一、李雅普诺夫指数的提出与历史 1961年冬季的一天，为了考察一条更长的序列，洛伦兹走了一条捷径。他在进行天气模式计算时没有从头开始运行，而是从中途开始。作为计算的初值，他直接输入了上次运算的输出结果，然后他穿过大厅下楼，清净的去喝一杯咖啡。一个小时之后他回来时，看到了出乎意料的事。从几乎相同出发点开始，洛伦兹看到他的计算机产生的天气模式差别愈来愈大，终至毫无相似之处。就是这件事播下了一门新科学的种子。 稳定体系的相轨线相应于趋向某个平衡点，如果出现越来越远离平衡点，则系统是不稳定的。系统只要有一个正值就会出现混沌运动。判断一个非线性体统是否存在混沌运动时，需要检查它的李雅普诺夫指数λ是否为正值。 二、李雅普诺夫指数的定义 Lyapunov 指数是描述时序数据所生成的相空间中两个极其相近的初值所产生的轨道，随时间推移按指数方式分散或收敛的平均变化率。任何一个系统，只要有一个Lyapunov 大于零，就认为该系统为混沌系统。 李雅普诺夫指数是指在相空间中相互靠近的两条轨线随着时间的推移，按指数分离或聚合的平均变化速率。李雅普诺夫指数的定义: 首先考虑一维映射假设初始位置附近有一点，则经过一次迭代后，这两点之间的距离为： （1） 并利用微分中值定理有： （2） n次迭代后，并利用微分中值定理，这两点之间的距离为： （3） 由（3）式可得：（4）又由复合函数的微分规则有： 其中 那么式（4）就变为： （5） 则称（6）为Lyapunov指数。 一维映射就对应一个李雅普诺夫指数，而且当时，该系统具有混沌特性。当时， 对应着分岔点或系统的周期解，既系统出现周期现象。时，系统有稳定的不动点，
物理学报 ACTA PHYSICA SINICA 2000 Vol.49 No.4 P.636-640 一种最大李雅普诺夫指数估计的稳健算法 杨绍清章新华赵长安 最大李雅普诺夫指数是诊断和描述动态系统混沌的重要参数.在深入研究相空间重构技术和轨道跟踪法的基础上，提出了一种从标量混沌时间序列中估计最大李雅普诺夫指数的新算法.该算法能够克服现有算法的不足，主要有以下三个优点：1）很高的精度；2）几乎不受噪声的影响；3）所需的计算时间和存贮空间小, 能进行在线计算. PACC： 0545 A ROBUST METHOD FOR ESTIMATING THE LARGEST LYAPUNOV EXPONENT YANG 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) ABSTRACT The 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.
Lorenz 系统 文档分两个文件方程m文件和计算L指数m文件分开写，复制粘贴即可运行matlab2012a，改写方程文件和参数即可算自己的系统，其中最大L指数用的是经典的柏内庭（G.Benettin）计算方法，准确快速无误！附计算结果图！！ 方程m文件： function dX = Loren(t,X) global a; % 变量不放入参数表中 global b; global c; x=X(1); y=X(2); z=X(3); % Y的三个列向量为相互正交的单位向量 % 输出向量的初始化 dX = zeros(6,1); % Lorenz吸引子 dX(1)=a*(y-x); dX(2)=x*(b-z)-y; dX(3)=x*y-c*z; end 计算最大L指数文件 Z=[]; global a; global b; global c; a=10; c=8/3; d0=1e-7; for b=linspace(0,500,500) lsum=0; x=1;y=1;z=1; x1=1;y1=1;z1=1+d0; for i=1:100 [T1,Y1]=ode45('Loren',1,[x;y;z;16;b;4]); [T2,Y2]=ode45('Loren',1,[x1;y1;z1;16;b;4]); n1=length(Y1);n2=length(Y2); x=Y1(n1,1);y=Y1(n1,2);z=Y1(n1,3); x1=Y2(n2,1);y1=Y2(n2,2);z1=Y2(n2,3); d1=sqrt((x-x1)^2+(y-y1)^2+(z-z1)^2); x1=x+(d0/d1)*(x1-x); y1=y+(d0/d1)*(y1-y); z1=z+(d0/d1)*(z1-z); if i>50
基于正定二次型的 李雅普诺夫稳定性分析 张俊超 （控制科学与工程、控制理论与控制工程、2010010215） 摘要：李雅普诺夫稳定性理论以状态向量描述为基础，不仅适用于单变量、 线性、定常系统，而且适用于多变量、非线性、时变系统。但要应用李氏判据判断系统稳定性，就要涉及到系统矩阵A特征值的求解以及根据系统状态方程构造正定二次型的李雅普诺夫函数来判断系统稳定性。 1.问题的提出 我们在处理实际工程问题时，经常需要判断系统稳定性，一般稳定性判据都有一定局限性，李雅普诺夫稳定性理论是确定系统稳定性的一般的理论，不仅适用于单变量、线性、定常系统，而且适用于多变量、非线性、时变系统，它以状态向量描述为基础，结合正定二次型的相关知识对系统稳定性进行判断。 2.问题的求解 李雅普诺夫稳定性理论分析系统稳定性的两种方法： （1）利用线性系统微分方程的解来判断系统的稳定性 ——李雅普诺夫第一法（间接法） 李雅普诺夫第一法的主要内容 1）用一次近似式表示状态方程，即：X=AX+B(x) 如果A的全部特征值都具有负实部，则系统在平衡点xe=0处是稳定的， 且系统的稳定性与高阶项B(x)无关。 2）如果X=AX+B(x)的A的特征值中至少有一个具有正实部，则无论B(x)如何，系统在平衡点xe=0处为不稳定的。 3）如果X=AX+B(x)的A的含有等于零的特征值，则系统的稳定性由B(x)决定。李雅普诺夫第一法是根据系统矩阵A的特征值来判断系统的稳定性的。 （2）构造李雅普诺夫函数，利用构造的李氏函数判断系统稳定性 ——李雅普诺夫第二法（直接法） 观察振动现象，若系统能量(含动能和位能)随时间推移而衰减，则系统迟早会达到平衡状态。基本思想：若系统内部能量随时间↑而↓,最终到达静止状态，系统稳定。虚构一个能量函数（李雅普诺夫函数) V(x，t)=f（x 1，x 2 ， (x) n ，t） V(x)=f（x 1，x 2 ， (x) n ） V(x,t)或V(x)是一个标量函数。能量总大于零，故为正定函数。能量随随时间增加而衰减，即：V(x，t)或V(x)的导数小于零。
目录 摘要 (1) 关键词 (1) Abstract (1) Keywords (1) 前言 (1) 1．预备知识 (1) 1.1李雅普诺夫第一法 (5) 1.2李雅普诺夫第二法 (1) 1.3线性系统的特征 (2) 2．李雅普诺夫意义下的稳定性 (2) 2.1稳定与一致稳定 (2) 2.2 渐进稳定和一致渐近稳定 (3) 2.3 不稳定 (3) 3．李雅普诺夫稳定性定理 (3) 4．线性系统的李雅普诺夫稳定性分析 (4) 小结 (7)
参考文献 (7)
李雅普诺夫方法在线性系统中的应用 摘要：在判定线性系统稳定性时，李雅普诺夫方法的优点在于无须求解系统方程的解，就能对系统的稳定性进行分析．文章介绍了李雅普诺夫稳定性分析在线性系统中的应用． 关键词：正定矩阵；标量函数；渐近稳定 Application of Lyapunov’s method in linear system Abstract:In determining the stability of linear systems,the advantages of the Lyapunov’s method is without solving the system equation,which can analyze the stability of the systems.we introduce the application in linear system analysis in Lyapunov stability in the paper. Keywords:positive definite matrix;Scalar function; asymptotic stability 前言 自动控制系统最重要的特性之一是稳定性．系统的稳定性，表示在遭受外界扰动偏离原来的平衡状态，而在扰动消失后，系统自身仍有能力恢复到原来平衡状态的一种“顽性”[]1．本文中，我们把研究对象集中到线性系统上，来讨论线性系统的稳定性问题．对于这个问题的讨论，都是建立在李雅普诺意义的稳定性的基本概念之上的．１．预备知识 1.1李雅普诺夫第一法 李雅普诺夫第一法又称间接法，它的基本思路是通过系统状态方程的解来判别系统的稳定性．对于线性定常系统，只需解出特征方程的解即可作出稳定性判断． 1.2李雅普诺夫第二法 李雅普诺夫第二法又称直接法，是通过构造一个类似于“能量”的李雅普诺夫函数，并分析它和其一次导数的定号性，直接对系统平衡状态的稳定性作出判断．
李雅普诺夫指数 ? 1.李雅普诺夫指数的定义 ? 2. 李雅普诺夫指数的划分意义 ? 3. 李雅普诺夫指数用在混沌中,如何应用 一李雅普诺夫指数的定义 李雅普诺夫指数是指在相空间中相互靠近的两条轨线随着时间的推移,按指数分离或聚合的平均变化速率。李雅普诺夫指数的定义: 首先考虑一维映射假设初始位置附近有一点,则经过一次迭代后,这两点之间的距离为： （1） 并利用微分中值定理有： （2） n次迭代后,并利用微分中值定理,这两点之间的距离为： （3） 由（3）式可得：（4）又由复合函数的微分规则有： 其中 那么式（4）就变为： （5） 则称（6）为Lyapunov指数。 一维映射就对应一个李雅普诺夫指数,而且当时,该系统具有混沌特性。当时,对应着分岔点或系统的周期解,既系统出现周期现象。时,系统有稳定的不动点,即此时对应的是一个点。而对于多维系统则有多个李雅普诺夫指数。 Lyapunov 特性指数沿某一方向取值的正负和大小表示长时间系统在吸引子中相邻轨线沿该方向平均发散或收敛i的快慢程度,仅从数学角度考虑,Lyapunov特性指数无量纲。n 维系统具有n 个Lyapunov 特性指数,形成指数谱。 其中数值最大的被称为最大Lyapunov 特性指数。最大Lyapunov 指数定义为 其中,表示时刻最邻近零点间的距离；M为计算总步数。最大Lyapunov指数不仅是区别混沌吸引子的重要指标,也是混沌系统对于初始值敏感性的定量描述。其中一维系统只有一个指数,二维系统有两个指数来表征。在实际计算中,要计算所有的Lyapunov指数,计算量较大,尤其当系统维数L较大时更为突出.所以注意力集中在计算系统的最大Lyapunov指数λm上. 二李雅普诺夫指数的物理意义 系统的Lyapunov指数谱可有效地表征变量随时间演化时,系统对初值的敏感性。指数小于零说明体系的相体积在该方向上是收缩的,此方向的运动是稳定的； 而正的指数值则表明了体系的相体积在该方向上不断膨胀和折叠,以致吸引子中本来邻近的轨线变得越来越不相关,从而使初态对任何不确定性的系统的长期行为成为不可预测,即所谓的初值敏感性。
求取lyapunov指数的小数据量方法，采用混合编程lylorenz.f90的程序如下： program lylorenz parameter(n=3,m=12,st=100) integer::i,j,k real y(m),z(n),s(n),yc(m),h,y1(m),a,b,r,f(m),k1,k2,k3 y(1)=10. y(2)=1. y(3)=0. a=10. b=8./3. r=28. t=0. h=0.01 !!!!!initial conditions do i=n+1,m y(i)=0. end do do i=1,n y((n+1)*i)=1. s(i)=0 end do open(1,file='lorenz1.dat') open(2,file='ly lorenz.dat') do 100 k=1,st !!!!!!!!st iterations call rgkt(m,h,t,y,f,yc,y1) !!!!normarize vector 123 z=0. do i=1,n do j=1,n z(i)=z(i)+y(n*j+i)**2 enddo if(z(i)>0.)z(i)=sqrt(z(i)) do j=1,n y(n*j+i)=y(n*j+i)/z(i) enddo end do !!!!generate gsr coefficient k1=0.
k2=0. k3=0. do i=1,n k1=k1+y(3*i+1)*y(3*i+2) k2=k2+y(3*i+3)*y(3*i+2) k3=k3+y(3*i+1)*y(3*i+3) end do !!!!conduct new vector2 and 3 do i=1,n y(3*i+2)=y(3*i+2)-k1*y(3*i+1) y(3*i+3)=y(3*i+3)-k2*y(3*i+2)-k3*y(3*i+1) end do !!!generate new vector2 and 3,normarize them do i=2,n z(i)=0. do j=2,n z(i)=z(i)+y(n*j+i)**2 enddo if(z(i)>0.)z(i)=sqrt(z(i)) do j=2,n y(n*j+i)=y(n*j+i)/z(i) end do end do !!!!!!!update lyapunov exponent do i=1,n if(z(i)>0)s(i)=s(i)+log(z(i)) enddo 100 continue do i=1,n s(i)=s(i)/(1.*st*h*1000.) write(2,*)s(i) enddo end !!!!!!!!!!!!!!!!!!!!!!!!! subroutine rgkt(m,h,t,y,f,yc,y1) real y(m),f(m),y1(m),yc(m),a,b,r integer::i,j do j=1,1000 call df(m,t,y,f)
【总结】 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];
计算最大lyapunov指数的编程思想 通过计算最大的李雅谱诺夫指数来检验系统是否混沌，一个正的最大李雅谱诺夫指数标致了该系统是否混沌，当你已经获得已知产生混沌的方程，这是相对容易做的。当你仅仅只有试验数据记录时，这种计算的困难程度是几乎不可能的，这里我们不考虑这种情形。计算的一般思路是跟踪两个较近的轨道，计算它们分离的平均对数率，无论它们分离得多远，其中之一的轨道将最终回到另一条轨道的附近。在每次迭代过程中，一程序将实现这个工作，完整的程序过程如下：1：在吸引域内，开始于任何初值条件。 最好开始于吸引子上的一个已知点，这种情形，第2步可以省略掉。 2：迭代直至轨道落到吸引子上。 这需要自己的判断或通过研究得到的有关该系统的一些知识，对绝大多数系统，只要迭代几百次，我们就认为足够了。通常情况下，除非你选的初值点刚好靠近分岔点，要不在几百次后它不会远离吸引子的。 3：选择（几乎任意的）靠近的点（设两点的距离为d0） 适合的d0的选择是至少要比计算机中的浮点数的1000倍大，如：在（8字节）的双精度的（对这些计算的最低推荐），变量有52－比特的尾数，故，精度是：2-52=2.22x10-16，因此，我们只要有d0=10-12 就足够了。 4：两轨道向前迭代一次，并计算分离量d1 分离量的计算：两轨道的每个分量上的差的平方和再开平方，如：对二维的系统有变量x和y，那么分离量是: d =[(xa-xb)2+(ya-yb)2]1/2, 这里下标（a和b）分别标出了两轨道。 5：以任何方便的底数估计log|d1/d0| 通常是采用自然对数（e为底数），但对映射，在每次少许的迭代中，李雅谱诺夫指数被引用，在这种情形下，我们需要以2做为底数（注：log2x=1.4427logex）。如果d1近似等于0，那么在估计对数时，将会得到错误，在这种情形下，尝试用更大的d0,如果这还不够，那么你不得不忽略这里发生的值。但是这样做的话，你的李雅谱诺夫指数计算可能有错误。 6：调整一条轨道，使得在方向d1上的分离量为d0 这也许是最难的一步，且易出错，如（在二维的情形），假如轨道b是要调整的，在一次调整后它的值是(xb1,yb1),那么它将重新初始化为：xb0=xa1+d0(xb1-xa1)/d1 和yb0=ya1+d0(yb1-ya1)/d1. 7：重复4－6步，计算第5步的平均值
§6.4 李雅普诺夫第二方法上一节我们介绍了稳定性概念,但是据此来判明系统解的稳定性,其应用范围是极其有限的. 李雅普诺夫创立了处理稳定性问题的两种方法:第一方法要利用微分方程的级数解,在他之后没有得到大的发展;第二方法是在不求方程解的情况下,借助一个所谓的李雅普诺夫函数)(x V 和通过微分方程所计算出来的导数 dt x dV )(的符号性质, 就能直接推断出解的稳定性,因此又称为直接法.本节主要介绍李雅普诺夫第二方法. 为了便于理解,我们只考虑自治系统 ) (x F d t d x =n R x ∈ (6.11) 假设T n x F x F x F ))(,),(()(1 =在{}K x R x G n ≤∈=上连续,满足局部利普希茨条件,且 O O F =)(. 为介绍李雅普诺夫基本定理,先引入李雅普诺夫函数概念. 定义6.3 若函数 R G x V →:)( 满足0)(=O V ,)(x V 和 i x V ??) ,,2,1(n i =都连续,且若存在K H ≤<0,使在 {}H x x D ≤=上)0(0)(≤≥x V ,则称)(x V 是常正(负)的;若在D 上除O x ≠外总有 )0(0)(<>x V ,则称)(x V 是正(负)定的;既不是常正又不是常负的函数称为变号函数. 通常我们称函数)(x V 为李雅普诺夫函数.易知: 函数2 22 1x x V +=在),(21x x 平面上为正定的; 函数 )(2 22 1x x V +-=在),(21x x 平面上为负定的; 函数222 1x x V -=在),(21x x 平面上为变号函数;
函数 21x V =在),(21x x 平面上为常正函数. 李雅普诺夫函数有明显的几何意义. 首先看正定函数),(21x x V V =. 在三维空间),,(21V x x 中, ),(21x x V V =是一个位于坐标面21Ox x 即0=V 上方的曲面.它与坐标面21Ox x 只在一个点,即原点)0,0,0(O 接触(图6-1(a)).如果用水平面 C V =(正常数)与),(21x x V V =相交,并将截口垂直投影到21Ox x 平面上,就得到一组一 个套一个的闭曲线族C x x V =),(21 (图6-1(b)),由于),(21x x V V =连续可微,且 0)0,0(=V ,故在021==x x 的充分小的邻域中, ) ,(21x x V 可以任意小.即在这些邻域中 存在C 值可任意小的闭曲线C V =. 对于负定函数),(21x x V V =可作类似的几何解释,只是曲面),(21x x V V =将在坐标面21Ox x 的下方. 对于变号函数),(21x x V V =,自然应对应于这样的曲面,在原点O 的任意邻域,它既有在21Ox x 平面上方的点,又有在其下方的点. 定理6.1 对系统(6.11),若在区域D 上存在李雅普诺夫函数)(x V 满足 (1) 正定; (2) )(1 ) 11.5(x F x V dt dV i n i i ∑=??= 常负, (a) (b)