基于QR分解的Duffing系统Lyapunov指数求解方法*

刘耀峰，杨 喜，舒 婷

(吉首大学物理与机电工程学院，湖南 吉首 416000)

基于QR分解的Duffing系统Lyapunov指数求解方法*

刘耀峰，杨 喜，舒 婷

(吉首大学物理与机电工程学院，湖南 吉首 416000)

Lyapunov指数衡量了非线性轨迹的稳定性和非线性系统的动力学特性，常作为判断Duffing系统混沌态和大尺度周期态的依据.根据Lyapunov指数的特征,Duffing系统的策动项采用正弦函数替代余弦函数方式，提出了一种基于QR分解的Duffing系统Lyapunov指数的求解方法.Matlab仿真结果表明了该算法的正确性、可靠性和有效性.

Duffing系统；Lyapunov指数；QR分解

近年来，随着弱信号检测技术和混沌理论的快速发展，一些研究者已经成功地将混沌理论用于信号检测领域，并取得了一些成果.Lyapunov指数作为混沌系统的一个重要参量，它既是混沌系统的一个重要判据，也是信号检测模型中的重要参考量.目前，大多数关于Lyapunov指数的求解采用QR分解方法.Benettin等[1]最早利用GS标准正交化过程对系统进行QR 分解，求解了系统的Lyapunov指数.Lai等[2]提出了用Jacobian法求解系统的Lyapunov指数，这种方法适应于有噪声的环境，与p-范数法相比，它的计算量比较小.然而重复正交化会带来庞大地计算量，影响计算效率，Rangaraian等提出了QR法的改进算法RHR算法，有效地避免了重复正交化.Udwadia等[3-4]针对RHR算法，提出了RHR改进算法，提高了计算效率.张宾[5]对低维系统Lyapunov指数的求解做了深入研究，并把Lyapunov指数引入到微弱信号检测领域.宋春云等[6]把最大Lyapunov指数作为混沌判据引入到Duffing系统信号检测模型中，金天等[7]分析了Lyapunov指数存在统计特性，确定了系统检测概率与误警概率等计算方法，李琳等[8]利用Lyapunov指数来确定系统的检测阈值.

笔者根据QR分解的基本思想，提出了基于QR分解的Duffing系统Lyapunov指数方法，详细分析了该方法设计思路，并根据该算法编写程序并仿真，仿真结果验证了该算法的正确性、可靠性和有效性，为下一步研究弱信号检测提供了平台.

1 Lyapunov指数

混沌系统的动力学特征可以通过系统的Lyapunov指数[5]和相轨迹来描述.当相空间中相邻轨迹f(t,x0)和f(t,x0+Δx)随着时间推移时，其轨迹将按f(t,x0)=f(t,x0+Δx)eλt规律相互吸引或离开，把这种相互吸引或离开的平均变化率称为Lyapunov指数.该指数衡量了非线性轨迹的稳定性，并从统计特性上反映非线性系统的动力学特性.同时，Lyapunov指数也量度了混沌系统对初始值敏感这一特性，即一个微小扰动，都将使混沌系统的Lyapunov指数发生改变，因此可以利用Lyapunov指数作为判断系统混沌态和大尺度周期态的依据.

对于N维相空间中的连续动力学系统，考虑相图中第i个维度方向上相邻的2个点.在时间0处和时间t处，设该2个点的间距分别为Pi(0)和Pi(t)，那么该系统在相图的第i个维度方向上的无量纲值Lyapunov指数为

其中：λi为Lyapunov指数；t为演变时间.

Lyapunov指数与相图中随时间演化得到的轨线之间收缩和扩张是息息相关地，在Lyapunov指数为负值方向上的轨道是收缩的，运动相对来说是稳定地且有规则的，对初始条件不敏感；相反，在Lyapunov指数为正值方向上的轨道是分离的，对初始条件敏感.N维系统拥有N个Lyapunov指数，系统的所有Lyapunov指数组成一个集合，称为Lyapunov指数谱.对于Duffing系统是否存在动力学混沌，可以先求解Duffing系统的Lyapunov指数谱，再借助于指数谱中最大Lyapunov指数是否大于零来直接判断[9]，当系统最大Lyapunov指数大于零时，系统随时间演化得到的相轨迹上的相邻点有排斥离开趋势，最终处于不稳定状态，即混沌状态；当系统最大Lyapunov指数小于零时，系统随时间演化得到的相轨迹上的相邻点有吸引收缩趋势，最终系统会处于稳定状态，即周期状态；当系统的最大Lyapunov指数等于零时，系统将处于混沌和周期状态之间的临界状态[10].

2 基于QR分解的Duffing系统Lyapunov指数求解方法

QR分解法即Jacobian方法，其基本思想是将已知动力系统的基本解矩阵，分解成正交矩阵Q和对角线元素都为正的上三角矩阵R的乘积，其中矩阵R是正交矩阵Q和系统Jacobian矩阵的函数，最终根据上三角矩阵R的对角线元素可以求得系统的Lyapunov指数.

考虑混沌Duffing振子的系统方程为

(1)

其中：k为阻尼比；-ax(t)+bx3(t)为非线性恢复力；γsin(ωt)为周期策动力；γ为周期策动力振幅.

(2)

其中令y1(t),y2(t),y3(t)的初始值分别为y1(0)=y10,y2(0)=y20,y3(0)=0.则(2)式的线性变分方程为

(3)

其中：I3是3阶的单位阵；J(t)是三维自治系统的Jacobian矩阵；Y(t)是(3)式的基本解矩阵，且为

(4)

对矩阵Y(t)进行QR分解，记作：Y(t)=Q(t)R(t).经过推导可以得到文献[6，10，12]所示的满足(3)式的Lyapunov指数表达式为

(5)

其中Rii(t)是上三角矩阵R(t)的主对角线元素.利用文献[12]中的连续系统中QR分解算法可得

(6)

对ln(Rii(t))求导可得

(7)

(7)式两边同时做积分运算后，结合(6)式可得

因此(5)式可进一步化为

(8)

其中(QT(t)J(t)Q(t))ii就是矩阵QT(t)J(t)Q(t)的主对角线元素.

从上述推导过程可知，利用QR分解算法求解系统的Lyapunov指数，则按照下列步骤求解：

(1) 将(1)式转化成三维自治系统(2)，求出Duffing系统的Jacobian矩阵；

(2) 求解正交矩阵Q(t)；

(3) 对(8)式进行数值积分，从而求得Lyapunov指数.

3 仿真实验与讨论

混沌Duffing振子方程(1)在k=0.5，a=b=1，ω=1时，系统状态将随着策动力振幅γ的变化而变化，同时Lyapunov指数也会发生相应变化.从文献[6-8，13]可知，Duffing系统存在一个阈值γd，该阈值是系统从混沌状态变化到周期状态的临界值.文中混沌Duffing振子的系统阈值为γd=0.826 0，当系统的策动力幅值γ<γd时，可以通过QR分解算法求得系统最大Lyapunov指数大于零，此时可以判断系统处于混沌状态，从相应的相轨迹图中可得到验证；当γ>γd时，同样可以用QR分解算法求得系统最大Lyapunov指数小于零，此时可以判断系统处于周期状态，同样也可以从相应的相轨迹图中得到验证.

设仿真的时间步长为0.01，每次迭代步数为10，总的循环次数为100 000次.利用库塔算法得到Duffing系统的Lyapunov指数谱.为了减少误差，提高系统的检测概率，应删除不稳定迭代，因为仿真得到的指数谱中均存在一定的过渡带[14-15]，在这些过渡带中的Lyapunov指数的数值，呈现正负交差状态.为了准确的判定系统的状态，要选取稳定的Lyapunov指数值.文中选取i=10 000点处的稳定的Lyapunov指数值.仿真结果如图1～6所示.当γ=0.825 9时，相轨迹仿真结果为图1的混沌状态.从图2的仿真结果可知，稳定时的最大Lyapunov指数值约为0.078 5；当γ=0.826 0时，相轨迹仿真结果为图3的间歇混沌状态.从图4的仿真结果可知，稳定时的最大Lyapunov指数值约为0.003 1；当γ=0.826 1时，相轨迹仿真结果为图5的周期状态.从图6的仿真结果可知，稳定时最大Lyapunov指数值约为-0.012 2，与文献[9]中的结论一致.

图2 γ=0.825 9时的混沌状态的相轨迹曲线

图3 γ=0.826 0时的Lyapunov指数演化曲线

图4 γ=0.8260时的临界状态的相轨迹曲线

图5 γ=0.826 1时的Lyapunov指数演化曲线

图6 γ=0.826 1时的周期状态的相轨迹曲线

4 结语

本文分析了用QR分解求解Duffing系统Lyapunov指数的方法，仿真结果表明，当Duffing系统处于临界状态时，系统再加入一个微弱的信号后，系统状态将由混沌态跃迁到周期状态，且图形变化明显.

[1] BENETTIN G,GALGANI L,GIORGILLI A,et al.Lyapunov Characteristic Exponents for Smooth Dynamical Systems and for Hamiltonian Systems;A Method For Computing All of Them.Part 1:Theory[J].Meccanica,1980,15(1):9-20.

[2] LAI D,CHEN G.Statistical Analysis of Lyapunov Exponents from Time Series:A Jacobian Approach[J].Mathematical and Computer Modelling,1998,27(7):1-9.

[3] UDWADIA F E,VON BREMEN H F.An Efficient and Stable Approach for Computation of Lyapunov Characteristic Exponents of Continuous Dynamical Systems[J].Applied Mathematics and Computation,2001,121(2):219-259.

[4] UDWADIA F E,VON BREMEN H F.Computation of Lyapunov Characteristic Exponents for Continuous Dynamical Systems[J].Zeitschrift FüR Angewandte Mathematik Und Physik ZAMP,2002,53(1):123-146.

[5] 张 宾.Lyapunov特性指数的算法研究及其在弱信号混沌检测中的应用[D].吉林:吉林大学,2002.

[6] 宋春云.最大Lyapunov特性指数在微弱信号检测中的应用[J].声学技术，2007,26(1):126-129.

[7] 金 天,张 骅.基于统计方法的混沌Duffing振子弱信号检测与估计[J].中国科学:信息科学，2011,41(10):1 184-1 199.

[8] 李 琳,刘春刚,石 硕,等.基于混沌振子和Lyapunov指数的微弱信号检测方法[J].黑龙江大学学报:自然科学版,2012,29(4):556-560.

[9] 李 月,杨宝俊.混沌振子系统(L-Y)与检测[M].北京:科学出版社,2005:83-101.

[10] 王晓亮.基于多小波变换与Duffing振子的微弱信号检测与估计[D].哈尔滨:哈尔滨工程大学,2012.

[11] DIECI L,RUSSELL R D,VAN VLECK E S.On the Compuation of Lyapunov Exponents for Continuous Dynamical Systems[J].SIAM Journal on Numerical Analysis,1997,34(1):402-423.

[12] MCDONALD E J,HIGHAM D J.Error Analysis of QR Algorithms for Computing Lyapunov Exponents[J].Electronic Transactions on Numerical Analysis,2001(12):234-251.

[13] 杨 淼,安建平,陈 宁,等.基于Duffing混沌系统的频谱感知算法[J].北京理工大学学报,2011,31(3):329-332.

[14] 杨红英,叶 昊,王桂增,等.Duffing振子的Lyapunov指数与Floquet指数研究[J].仪器仪表学报,2008,29(5):927-931.

[15] 刘海波,吴德伟,金 伟,等.Duffing振子微弱信号检测方法研究[J].物理学报,2013,62(5):42-47.

(责任编辑 陈炳权)

QRDecompositionBasedMethodfortheComputationofLyapunovExponentinDuffingSystems

LIU Yaofeng,YANG Xi,SHU Ting

(College of Physics and Mechanical & Electrical Eengineering,Jishou University,Jishou 416000,Hunan China)

In order to compute the Lyapunov exponent of Duffing system,this paper adopts a kind of classical method based on QR decomposition to calculate the Lyapunov exponent of Duffing system.According to the basic characteristics of the Lyapunov exponent,we presents a kind of Lyapunov exponent algorithm based on QR decomposition for the known Duffing system.By using software simulation,we know that the Duffing system Lyapunov exponent algorithm based on QR decomposition is correctness.

Duffing system；Lyapunov exponent；QR decomposition

1007-2985(2014)01-0058-05

2013-10-15

刘耀峰(1985-)，男，湖南娄底人，吉首大学物理与机电工程学院无线电物理硕士研究生，主要从事信号检测研究

杨 喜(1978-)，男，湖南湘阴人，吉首大学信息科学与工程学院副教授，博士，硕士生导师，主要从事认知无线电研究.

TP391.41

A

10.3969/j.issn.1007-2985.2014.01.014