Lyapunov指数计算

Lyapunov指数的计算与仿真摘要：Lyapunov指数在研究动力系统的分岔、混沌运动特性中起着重要的作用，是衡量系统动力学特性的一个重要定量指标，它表征了系统在相空间中相邻轨道间收敛或发散的平均指数率。对于系统是否存在动力学混沌，可以从最大Lyapunov指数是否大于零非常直观的判断出来：一个正的Lyapunov指数意味着在相空间中，无论初始两条轨线的间距多么小，其差别都会随时间的演化而成为指数率的增加以致达到无法预测，这就是混沌现象。本文介绍了Lyapunov指数的定义、性质及计算原理和数值计算的实现方法，应用Matlab软件、编写了计算Lyapunov指数程序，计算了Logistic映射系统、Henon系统及时间序列的Lyapunov指数。实例的计算机仿真表明Lyapunov指数是研究分岔、混沌运动，解决工程实际问题的有效方法之一。关键词：非线性动力系统；混沌；Lyapunov指数；时间序列1.Lyapunov指数的定义及性质1.1Lyapunov指数的定义对于一维离散映射系统)(1nnxFx=+，在每次迭代过程中，初始两点是相互分离还是相互靠拢由下式决定：1dxdF，映射使两点分开1dxdF，映射使两点靠拢但是在不断地迭代过程中，dxdF的值要不断变化。为了表示从整体看两相邻初始状态分离的情况，必须对时间（或迭代次数）取平均。为此，设平均每次迭代所引起的指数分离中的指数为λ，则初始两点0x与ε+0x相距ε，经n次迭代后，点变为)(0xFn与)(0ε+xFn，相距λεne，即：)(0xneλε=)(0ε+xFn-)(0xFn（1）取极限0→ε，∞→n，由上式得到：0)(ln1lim)()(ln1limlim)(0000xxnnnnndxxdFnxFxFnx=∞→→∞→=−+=εελε（2）上式取极限后便与初始点0x无关，并根据复合函数求导链式法则有：∑−=∞→==10'0)(ln1)(limniinxFnxλλ（3）由式（1）可知，当0λ时，相邻两点终归要靠拢合并成一点，这对应于定点（不动点）或周期运动；反之，当0λ时，相邻两点要分开，这对应于混沌运动。对于连续微分系统，设其运动由下述自治方程描述：nixxfxnii,,2,1),,(1==（4）又设它有已知解),,(020,100nxxxx，令nixxxiii,,2,10=+=δ（5）为0x邻域内的另一解，将式（5）代入式（4）得：∑==njjijixxLdtxd10)(δδ（6）矩阵0)()(0xjiijxfxL∂∂=（7）为线性化演化算子或李雅普诺夫矩阵。我们知道，只要)(0xLij有一个特征值l的实部是正的，ixδ就会指数发散：ltiiexx0δδ=（8）即空间中两相邻轨道将分开。但是这种线性化处理方法仅在解得邻域内有效，其结果不能说明沿轨道长时间整体结果如何。为此必须仿照一维离散映射系统中多次迭代求平均的方法求长时间（t）过程中单位时间内xδ的平均变化。对于一维（n=1）情形，设相邻两点间距离为ε，推广（或替代）式（8）应有：tetλεε)0()(=（9）式（8）与式（9）的区别在于，式（8）是线性情形下导出的，它只在局域成立，因为特征值il是x的函数，从而式（8）也只有在t很小时才能成立；式（9）虽然也利用了局域式（8），但它是极长时间的平均式，λ值是取平均的结果。因此，式（9）是对整个吸引子整体而言的，λ是一个整体量。此时的Lyapunov指数的定义可写成：)0()(ln1limlim0)0(εελεttt→∞→=（10）由此二式可知：0λ，轨道收缩；0λ，轨道分离。在n维相空间中，相邻两点间距离（位移）εδ=x是n维矢量。随着时间演化，各方向伸长或压缩不一样。试想在0=t时以0x为中心，以小量)0(ε为半径作一n维球面。经时间t后，此球面将演化为n维椭球面，如图1所示（图是n=3的情形）。此椭球第i个坐标轴方向的半长轴为)()(ttiiεδ=。因此这时系统应该有n个Lyapunov指数，其中第i个Lyapunov指数应为：图1三维球面随时间的演化)0()(ln1limlim0)0(εελεttiti→∞→=（11）通常将此n个Lyapunov指数按下述方式排列：nλλλ≥≥≥21（12）1.2Lyapunov指数的性质当相空间是一维（单变量）时，吸引子只可能是稳定定点，其周围的点都要趋于此定点，此时必有0λ。反之，如果不存在定点或定点不稳定，即不存在吸引子，则0λ。对于二维情形，吸引子或者是稳定定点或者是极限环。对于稳定定点任意方向的ixδ都要收缩，故这时两个Lyapunov指数都是负的：01λ，02λ。至于极限环，取ixδ沿环线的法线方向，它一定要收缩，故此法线方向的Lyapunov指数总是小于零；当取ixδ沿环线的切线方向时，它既不会增大也不会缩小，因此沿切线方向应有0=λ。对于三维情形（n=3），类似的分析结果是：),,(),,(321−−−=λλλ：稳定定点),,0(),,(321−−=λλλ：极限环（周期运动）),0,0(),,(321−=λλλ：二维环面（准周期振荡）)0,,(),,(321++=λλλ：不稳极限环)0,0,(),,(321+=λλλ：不稳二维环面),0,(),,(321−+=λλλ：奇怪吸引子（混沌运动）我们可以看出，Lyapunov指数可以表征系统运动特征：所有Lyapunov指数取值的集合决定系统在相空间轨线的性质；沿某一方向Lyapunov指数iλ取值的正负和大小表示长时间系统在相空间中相邻轨线沿该方向平均发散（0iλ）或收敛（0iλ）的快慢程度；最大Lyapunov指数1λ决定相邻轨线是（01≤λ）否（01λ）能靠拢形成稳定轨道或稳定定点；最小Lyapunov指数nλ则表示相空间中所有轨线能（0nλ）否（0nλ）收缩形成稳定吸引子。也可按如下方式理解iλ的意义：如图1所示，相空间中一个半径为ε的超球体随着时间的演化要变为超椭球体。各iλ分别决定超椭球体的各半长轴的大小。故1λ决定椭球体的t时刻后最大半长轴的长度（te1λε）；（21λλ+）则决定最大和第二大两长轴所构成的椭圆的面积（te)(221λλπε+）。依此类推，前面k个长轴所构成的超球体的体积由kλλλ+++21决定。而∑=nii1λ则决定整个n维超球体的体积收缩的快慢。由以上分析可以看出：（1）任何吸引子，不论其奇怪与否，至少有一个Lyapunov指数小于零，否则轨线就不可能收缩为吸引子。（2）稳定定态和周期运动（以及准周期运动）不可能有正的Lyapunov指数。稳定定态的Lyapunov指数都是负的；周期运动的01=λ，其余的Lyapunov指数也都是负的。（3）混沌运动至少有一个正的Lyapunov指数。反之，如果计算得知系统至少有一个正的Lyapunov指数，则可肯定系统做混沌运动。（4）当4≥n时，还可能有不止一个Lyapunov指数大于零，这时的运动称为超混沌。2.Lyapunov指数的计算仿真2.1一维离散映射系统的Lyapunov指数对于一维离散映射系统，我们采用式（3）计算系统的Lyapunov指数。当然，实际计算中n不可能取无穷大，通常选取n为200。为了保证计算结果准确，我们选取稳态过程中的某点为初值。举例Logistic映射系统：nnnxxux)1(1−=+（13）其中]1,0[∈nx，]4,1[∈u根据（3）式编写matlab程序计算Logistic映射系统的Lyapunov指数谱，即Lyapunov指数随参数u的变化情况。图2Logistic映射分岔图及Lyapunov指数谱从其系统的分岔图及对应的Lyapunov指数谱我们可以看出，随着参数u的增加，系统由倍周期分岔进入混沌，而且在倍周期分岔点处，系统的Lyapunov指数为零。当系统进入混沌后，其Lyapunov指数为变为正。从混沌图中我们可以看出在混沌区有一个周期窗口出现，在对应的Lyapunov指数谱中可以看到此点处Lyapunov指数为负。2.2多维离散映射系统的Lyapunov指数在离散动力系统中，我们经常以1−=ttTTT的形式来处理映射xTxtt=，其中x属于紧密相连的流行M，映射是从有限的n维空间映射到n维空间，参数t代表一个非负整数或一个实数。切空间序列txdT可以通过考虑1−=txxTtxdTdTdTT迭代得到。我们可以考虑这些运算在标准基中的矩阵表示，事实上，任何标准正交基都是有效地。然后，Lyapunov指数可以通过对这些切空间矩阵表示（决定于流行M上某合适的一点x）的乘积进行QR分解得到。我们记与it=时对应的切空间的矩阵表示为iJ。为了得到所有Lyapunov指数的较好估计，需要对矩阵乘积11JJJmm−进行QR分解。其过程可按如下方式进行，其中IQ=0，我们有RQRRRQRRRRQJJJqrRRQJJJJqrRQJJJJQJJJJqrJJJqrmmmiiiimmmmmmmmmm========−−−−−−−−][])][([])][([])][([)]([][1212211112234111231012111（14）其中][⋅qr代表QR分解。从1J开始，在第i步中，我们给1−iQ左乘矩阵iJ得1−=iiiQJB，然后对其进行QR分解iiiRQB=，mi,,2,1=。矩阵R是矩阵乘积12RRRm。进一步我们知道，矩阵R对角线上的每个元素就是所有矩阵iR对角线上对应元素的乘积。因此，所有Lyapunov指数的近似值可以通过下式得到：∑==miikkkRmx1),(ln1，nk,,2,1=（15）算法可表述为下：初始化初始化Q为nn×的单位矩阵初始化LCEvector为n维零向量fori=1tom_iterationsQJBi=计算B的QR分解QRB=LCEvector=LCEvector+log(diagR)EndLCEvector=LCEvector/m_iteration则系统的n个Lyapunov指数由向量LCEvector给出。举例Henon系统：=+−=++nnnnnbxyyaxxT1211:其中]1,0[,∈nnyx，其Jacobie矩阵为−=012baxJn。固定b为0.3，参数a从0到1.4变化，用matlab编写计算程序，得到如图3所示的系统分岔图及Lyapunov指数谱。图3Henon二维映射系统分岔图及Lyapunov指数谱同样，从系统的分岔图我们可以看出，随着参数a的增加，系统由倍周期分岔进入混沌，在倍周期分岔点上，最大Lyapunov指数为零。当系统进入混沌后，最大Lyapunov指数为正。在参数a稍大于1.2时出现了周期窗口，我们可以看出此时的最大Lyapunov指数为负。2.3时间序列的Lyapunov指数从Lyapunov指数的性质中我们可以看出，最重要的是最大Lyapunov指数1λ，由其可直接判断系统运动是不是混沌。因此，对于测得的时间序列我们通常只需计算其最大Lyapunov指数即可。G.Benettin等人提出了一个计算动力系统的1λ的方法。选取两个很靠近的初始状态0x和0y，随着时间的演化，它们的状态分别是),,(21xxx和),,(21yyy，则第n时刻（以时间间隔τ为单位）两轨道上的状态为)(τnx和)(τny，它们之间的距离是：)()(ττnxnydn−=（16）于是在n+1时刻两轨道之间的距离变为：])1[(])1[(1ττ+−+=+nxnydn（17）为避免计算时出现溢出，选取一个新的计算起点])1[('τ+ny代替])1[(τ+ny，])1[('τ+ny是在])1[(τ+ny和])1[(τ+nx的连线上。但距离])1[(τ+nx保持最初的0d。以])1[(τ+nx和])1[('τ+ny为起点再求τ)2(+n时刻两状态和其间距2+nd。这样，每次都是从距离为0d的两状态出发，得到一系列的距离1d、2d、…。当τ很小时，由式（11）可求