劉耀峰，楊 喜，舒 婷

(吉首大學(xué)物理與機電工程學(xué)院，湖南 吉首 416000)

基于QR分解的Duffing系統(tǒng)Lyapunov指數(shù)求解方法*

劉耀峰，楊 喜，舒 婷

(吉首大學(xué)物理與機電工程學(xué)院，湖南 吉首 416000)

Lyapunov指數(shù)衡量了非線性軌跡的穩(wěn)定性和非線性系統(tǒng)的動力學(xué)特性，常作為判斷Duffing系統(tǒng)混沌態(tài)和大尺度周期態(tài)的依據(jù).根據(jù)Lyapunov指數(shù)的特征,Duffing系統(tǒng)的策動項采用正弦函數(shù)替代余弦函數(shù)方式，提出了一種基于QR分解的Duffing系統(tǒng)Lyapunov指數(shù)的求解方法.Matlab仿真結(jié)果表明了該算法的正確性、可靠性和有效性.

Duffing系統(tǒng)；Lyapunov指數(shù)；QR分解

近年來，隨著弱信號檢測技術(shù)和混沌理論的快速發(fā)展，一些研究者已經(jīng)成功地將混沌理論用于信號檢測領(lǐng)域，并取得了一些成果.Lyapunov指數(shù)作為混沌系統(tǒng)的一個重要參量，它既是混沌系統(tǒng)的一個重要判據(jù)，也是信號檢測模型中的重要參考量.目前，大多數(shù)關(guān)于Lyapunov指數(shù)的求解采用QR分解方法.Benettin等[1]最早利用GS標(biāo)準(zhǔn)正交化過程對系統(tǒng)進(jìn)行QR 分解，求解了系統(tǒng)的Lyapunov指數(shù).Lai等[2]提出了用Jacobian法求解系統(tǒng)的Lyapunov指數(shù)，這種方法適應(yīng)于有噪聲的環(huán)境，與p-范數(shù)法相比，它的計算量比較小.然而重復(fù)正交化會帶來龐大地計算量，影響計算效率，Rangaraian等提出了QR法的改進(jìn)算法RHR算法，有效地避免了重復(fù)正交化.Udwadia等[3-4]針對RHR算法，提出了RHR改進(jìn)算法，提高了計算效率.張賓[5]對低維系統(tǒng)Lyapunov指數(shù)的求解做了深入研究，并把Lyapunov指數(shù)引入到微弱信號檢測領(lǐng)域.宋春云等[6]把最大Lyapunov指數(shù)作為混沌判據(jù)引入到Duffing系統(tǒng)信號檢測模型中，金天等[7]分析了Lyapunov指數(shù)存在統(tǒng)計特性，確定了系統(tǒng)檢測概率與誤警概率等計算方法，李琳等[8]利用Lyapunov指數(shù)來確定系統(tǒng)的檢測閾值.

筆者根據(jù)QR分解的基本思想，提出了基于QR分解的Duffing系統(tǒng)Lyapunov指數(shù)方法，詳細(xì)分析了該方法設(shè)計思路，并根據(jù)該算法編寫程序并仿真，仿真結(jié)果驗證了該算法的正確性、可靠性和有效性，為下一步研究弱信號檢測提供了平臺.

1 Lyapunov指數(shù)

混沌系統(tǒng)的動力學(xué)特征可以通過系統(tǒng)的Lyapunov指數(shù)[5]和相軌跡來描述.當(dāng)相空間中相鄰軌跡f(t,x0)和f(t,x0+Δx)隨著時間推移時，其軌跡將按f(t,x0)=f(t,x0+Δx)eλt規(guī)律相互吸引或離開，把這種相互吸引或離開的平均變化率稱為Lyapunov指數(shù).該指數(shù)衡量了非線性軌跡的穩(wěn)定性，并從統(tǒng)計特性上反映非線性系統(tǒng)的動力學(xué)特性.同時，Lyapunov指數(shù)也量度了混沌系統(tǒng)對初始值敏感這一特性，即一個微小擾動，都將使混沌系統(tǒng)的Lyapunov指數(shù)發(fā)生改變，因此可以利用Lyapunov指數(shù)作為判斷系統(tǒng)混沌態(tài)和大尺度周期態(tài)的依據(jù).

對于N維相空間中的連續(xù)動力學(xué)系統(tǒng)，考慮相圖中第i個維度方向上相鄰的2個點.在時間0處和時間t處，設(shè)該2個點的間距分別為Pi(0)和Pi(t)，那么該系統(tǒng)在相圖的第i個維度方向上的無量綱值Lyapunov指數(shù)為

其中：λi為Lyapunov指數(shù)；t為演變時間.

Lyapunov指數(shù)與相圖中隨時間演化得到的軌線之間收縮和擴張是息息相關(guān)地，在Lyapunov指數(shù)為負(fù)值方向上的軌道是收縮的，運動相對來說是穩(wěn)定地且有規(guī)則的，對初始條件不敏感；相反，在Lyapunov指數(shù)為正值方向上的軌道是分離的，對初始條件敏感.N維系統(tǒng)擁有N個Lyapunov指數(shù)，系統(tǒng)的所有Lyapunov指數(shù)組成一個集合，稱為Lyapunov指數(shù)譜.對于Duffing系統(tǒng)是否存在動力學(xué)混沌，可以先求解Duffing系統(tǒng)的Lyapunov指數(shù)譜，再借助于指數(shù)譜中最大Lyapunov指數(shù)是否大于零來直接判斷[9]，當(dāng)系統(tǒng)最大Lyapunov指數(shù)大于零時，系統(tǒng)隨時間演化得到的相軌跡上的相鄰點有排斥離開趨勢，最終處于不穩(wěn)定狀態(tài)，即混沌狀態(tài)；當(dāng)系統(tǒng)最大Lyapunov指數(shù)小于零時，系統(tǒng)隨時間演化得到的相軌跡上的相鄰點有吸引收縮趨勢，最終系統(tǒng)會處于穩(wěn)定狀態(tài)，即周期狀態(tài)；當(dāng)系統(tǒng)的最大Lyapunov指數(shù)等于零時，系統(tǒng)將處于混沌和周期狀態(tài)之間的臨界狀態(tài)[10].

2 基于QR分解的Duffing系統(tǒng)Lyapunov指數(shù)求解方法

QR分解法即Jacobian方法，其基本思想是將已知動力系統(tǒng)的基本解矩陣，分解成正交矩陣Q和對角線元素都為正的上三角矩陣R的乘積，其中矩陣R是正交矩陣Q和系統(tǒng)Jacobian矩陣的函數(shù)，最終根據(jù)上三角矩陣R的對角線元素可以求得系統(tǒng)的Lyapunov指數(shù).

考慮混沌Duffing振子的系統(tǒng)方程為

(1)

其中：k為阻尼比；-ax(t)+bx3(t)為非線性恢復(fù)力；γsin(ωt)為周期策動力；γ為周期策動力振幅.

(2)

其中令y1(t),y2(t),y3(t)的初始值分別為y1(0)=y10,y2(0)=y20,y3(0)=0.則(2)式的線性變分方程為

(3)

其中：I3是3階的單位陣；J(t)是三維自治系統(tǒng)的Jacobian矩陣；Y(t)是(3)式的基本解矩陣，且為

(4)

對矩陣Y(t)進(jìn)行QR分解，記作：Y(t)=Q(t)R(t).經(jīng)過推導(dǎo)可以得到文獻(xiàn)[6，10，12]所示的滿足(3)式的Lyapunov指數(shù)表達(dá)式為

(5)

其中Rii(t)是上三角矩陣R(t)的主對角線元素.利用文獻(xiàn)[12]中的連續(xù)系統(tǒng)中QR分解算法可得

(6)

對ln(Rii(t))求導(dǎo)可得

(7)

(7)式兩邊同時做積分運算后，結(jié)合(6)式可得

因此(5)式可進(jìn)一步化為

(8)

其中(QT(t)J(t)Q(t))ii就是矩陣QT(t)J(t)Q(t)的主對角線元素.

從上述推導(dǎo)過程可知，利用QR分解算法求解系統(tǒng)的Lyapunov指數(shù)，則按照下列步驟求解：

(1) 將(1)式轉(zhuǎn)化成三維自治系統(tǒng)(2)，求出Duffing系統(tǒng)的Jacobian矩陣；

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

(3) 對(8)式進(jìn)行數(shù)值積分，從而求得Lyapunov指數(shù).

3 仿真實驗與討論

混沌Duffing振子方程(1)在k=0.5，a=b=1，ω=1時，系統(tǒng)狀態(tài)將隨著策動力振幅γ的變化而變化，同時Lyapunov指數(shù)也會發(fā)生相應(yīng)變化.從文獻(xiàn)[6-8，13]可知，Duffing系統(tǒng)存在一個閾值γd，該閾值是系統(tǒng)從混沌狀態(tài)變化到周期狀態(tài)的臨界值.文中混沌Duffing振子的系統(tǒng)閾值為γd=0.826 0，當(dāng)系統(tǒng)的策動力幅值γ<γd時，可以通過QR分解算法求得系統(tǒng)最大Lyapunov指數(shù)大于零，此時可以判斷系統(tǒng)處于混沌狀態(tài)，從相應(yīng)的相軌跡圖中可得到驗證；當(dāng)γ>γd時，同樣可以用QR分解算法求得系統(tǒng)最大Lyapunov指數(shù)小于零，此時可以判斷系統(tǒng)處于周期狀態(tài)，同樣也可以從相應(yīng)的相軌跡圖中得到驗證.

設(shè)仿真的時間步長為0.01，每次迭代步數(shù)為10，總的循環(huán)次數(shù)為100 000次.利用庫塔算法得到Duffing系統(tǒng)的Lyapunov指數(shù)譜.為了減少誤差，提高系統(tǒng)的檢測概率，應(yīng)刪除不穩(wěn)定迭代，因為仿真得到的指數(shù)譜中均存在一定的過渡帶[14-15]，在這些過渡帶中的Lyapunov指數(shù)的數(shù)值，呈現(xiàn)正負(fù)交差狀態(tài).為了準(zhǔn)確的判定系統(tǒng)的狀態(tài)，要選取穩(wěn)定的Lyapunov指數(shù)值.文中選取i=10 000點處的穩(wěn)定的Lyapunov指數(shù)值.仿真結(jié)果如圖1～6所示.當(dāng)γ=0.825 9時，相軌跡仿真結(jié)果為圖1的混沌狀態(tài).從圖2的仿真結(jié)果可知，穩(wěn)定時的最大Lyapunov指數(shù)值約為0.078 5；當(dāng)γ=0.826 0時，相軌跡仿真結(jié)果為圖3的間歇混沌狀態(tài).從圖4的仿真結(jié)果可知，穩(wěn)定時的最大Lyapunov指數(shù)值約為0.003 1；當(dāng)γ=0.826 1時，相軌跡仿真結(jié)果為圖5的周期狀態(tài).從圖6的仿真結(jié)果可知，穩(wěn)定時最大Lyapunov指數(shù)值約為-0.012 2，與文獻(xiàn)[9]中的結(jié)論一致.

圖2 γ=0.825 9時的混沌狀態(tài)的相軌跡曲線

圖3 γ=0.826 0時的Lyapunov指數(shù)演化曲線

圖4 γ=0.8260時的臨界狀態(tài)的相軌跡曲線

圖5 γ=0.826 1時的Lyapunov指數(shù)演化曲線

圖6 γ=0.826 1時的周期狀態(tài)的相軌跡曲線

4 結(jié)語

本文分析了用QR分解求解Duffing系統(tǒng)Lyapunov指數(shù)的方法，仿真結(jié)果表明，當(dāng)Duffing系統(tǒng)處于臨界狀態(tài)時，系統(tǒng)再加入一個微弱的信號后，系統(tǒng)狀態(tài)將由混沌態(tài)躍遷到周期狀態(tài)，且圖形變化明顯.

[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特性指數(shù)的算法研究及其在弱信號混沌檢測中的應(yīng)用[D].吉林:吉林大學(xué),2002.

[6] 宋春云.最大Lyapunov特性指數(shù)在微弱信號檢測中的應(yīng)用[J].聲學(xué)技術(shù)，2007,26(1):126-129.

[7] 金 天,張 驊.基于統(tǒng)計方法的混沌Duffing振子弱信號檢測與估計[J].中國科學(xué):信息科學(xué)，2011,41(10):1 184-1 199.

[8] 李 琳,劉春剛,石 碩,等.基于混沌振子和Lyapunov指數(shù)的微弱信號檢測方法[J].黑龍江大學(xué)學(xué)報:自然科學(xué)版,2012,29(4):556-560.

[9] 李 月,楊寶俊.混沌振子系統(tǒng)(L-Y)與檢測[M].北京:科學(xué)出版社,2005:83-101.

[10] 王曉亮.基于多小波變換與Duffing振子的微弱信號檢測與估計[D].哈爾濱:哈爾濱工程大學(xué),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混沌系統(tǒng)的頻譜感知算法[J].北京理工大學(xué)學(xué)報,2011,31(3):329-332.

[14] 楊紅英,葉 昊,王桂增,等.Duffing振子的Lyapunov指數(shù)與Floquet指數(shù)研究[J].儀器儀表學(xué)報,2008,29(5):927-931.

[15] 劉海波,吳德偉,金 偉,等.Duffing振子微弱信號檢測方法研究[J].物理學(xué)報,2013,62(5):42-47.

(責(zé)任編輯 陳炳權(quán))

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-)，男，湖南婁底人，吉首大學(xué)物理與機電工程學(xué)院無線電物理碩士研究生，主要從事信號檢測研究

楊 喜(1978-)，男，湖南湘陰人，吉首大學(xué)信息科學(xué)與工程學(xué)院副教授，博士，碩士生導(dǎo)師，主要從事認(rèn)知無線電研究.

TP391.41

A

10.3969/j.issn.1007-2985.2014.01.014