劉耀峰，楊 喜，舒 婷

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

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

劉耀峰，楊 喜，舒 婷

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

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

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

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

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

1 Lyapunov指數(shù)

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

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

其中：λi為L(zhǎng)yapunov指數(shù)；t為演變時(shí)間.

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

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

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

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

(1)

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

(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)

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

(5)

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

(6)

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

(7)

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

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

(8)

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

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

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

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

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

3 仿真實(shí)驗(yàn)與討論

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

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

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

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

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

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

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

4 結(jié)語(yǔ)

本文分析了用QR分解求解Duffing系統(tǒng)Lyapunov指數(shù)的方法，仿真結(jié)果表明，當(dāng)Duffing系統(tǒng)處于臨界狀態(tài)時(shí)，系統(tǒng)再加入一個(gè)微弱的信號(hào)后，系統(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ù)的算法研究及其在弱信號(hào)混沌檢測(cè)中的應(yīng)用[D].吉林:吉林大學(xué),2002.

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

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

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

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

[10] 王曉亮.基于多小波變換與Duffing振子的微弱信號(hào)檢測(cè)與估計(jì)[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é)報(bào),2011,31(3):329-332.

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

[15] 劉海波,吳德偉,金 偉,等.Duffing振子微弱信號(hào)檢測(cè)方法研究[J].物理學(xué)報(bào),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é)物理與機(jī)電工程學(xué)院無(wú)線電物理碩士研究生，主要從事信號(hào)檢測(cè)研究

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

TP391.41

A

10.3969/j.issn.1007-2985.2014.01.014