基于HB-AFT算法的陣發(fā)性混沌振動研究新方法*

李磊 張智勇,? 芮筱亭 陳予恕

(1.南京理工大學 理學院, 南京 210094) (2.南京理工大學 發(fā)射動力學研究所, 南京 210094)(3.哈爾濱工業(yè)大學 航天學院, 哈爾濱 150001)

引言

隨著科學技術(shù)的發(fā)展,混沌理論有了長足的發(fā)展. 陣發(fā)性混沌理論從20 世紀80 年代起步,在等離子體[1]、非線性電路[2]、流體力學[3]、細胞力學[4]和機器人動力學[5]等分析中取得了大量理論和實驗成果,由此受到了廣泛的關(guān)注.陣發(fā)性混沌的主要特征是規(guī)則和不規(guī)則運動偽隨機交替出現(xiàn).常見的陣發(fā)性混沌經(jīng)過折疊分岔、亞臨界Hopf分岔和亞臨界倍周期分岔伴隨系統(tǒng)全局結(jié)構(gòu)的重整化產(chǎn)生,Pomeau和Manneville[6]把這三類陣發(fā)失穩(wěn)嵌入混沌的方式分別定義為陣發(fā)Ⅰ、Ⅱ和Ⅲ型混沌,顯然這三種分岔并不是陣發(fā)混沌產(chǎn)生的充分條件[7].近20年來,陣發(fā)性混沌理論得到了進一步豐富[8],陣發(fā)V型、X型混沌、開關(guān)(on-off)陣發(fā)、激變陣發(fā)等多種新的陣發(fā)混沌響應(yīng)模式被理論和實驗發(fā)現(xiàn).陣發(fā)性混沌演化涉及系統(tǒng)全局動力學特性的轉(zhuǎn)變,因此理論上研究困難,常采用數(shù)值方法進行計算分析.然而,陣發(fā)性混沌往往發(fā)生在失穩(wěn)不動點附近,由于臨界慢化現(xiàn)象[9]的發(fā)生,對數(shù)值積分模擬要求高精度和較長的計算時間.綜上所述,就陣發(fā)性混沌研究而言,無論在理論發(fā)展、工程應(yīng)用還是研究手段領(lǐng)域,都存在大量開放性問題,方興未艾.

有許多經(jīng)典的數(shù)值方法、定量方法和解析方法可對非線性動力系統(tǒng)的特性進行分析.一般把所求問題看成常微分方程的初值問題,采用數(shù)值積分方法來求解系統(tǒng)的漸近穩(wěn)態(tài)響應(yīng)[10],并采用Poincare截面、頻譜圖、Lyapunov指數(shù)和數(shù)值分岔圖等分析系統(tǒng)的混沌特征.但是,由于陣發(fā)性混沌響應(yīng)的陡變特性,在數(shù)值積分求解過程中需要頻繁改變步長,且臨界慢化現(xiàn)象使得接近失穩(wěn)位置的穩(wěn)態(tài)響應(yīng)計算非常耗時[11].小參數(shù)攝動法、平均法、多尺度法等各種經(jīng)典的定量方法受到小參數(shù)條件、可積條件、級數(shù)展開條件等的限制[12],而諧波平衡法對于復(fù)雜非線性系統(tǒng)諧波平衡本身就是極為困難的.另外,上述定量解析方法在高次諧波解求解時,由計算工作量原因,通常僅求解一次或二次諧波解,可是相關(guān)研究表明高次、小諧波分量的忽略可能會對系統(tǒng)的穩(wěn)定性帶來極大的誤差甚至錯誤[13].

HB-AFT方法是Yamauchi[14]提出來的一種半解析半數(shù)值的隱式諧波平衡法.該方法在求解系統(tǒng)周期響應(yīng)過程中把響應(yīng)和非線性函數(shù)同時設(shè)為諧波解形式,根據(jù)系統(tǒng)的離散時頻特性建立諧波系數(shù)之間的關(guān)系,對于復(fù)雜非線性項無需級數(shù)展開、積分處理等近似過程,因此該方法具有一定的普遍適用性[15]. 隨后,Kim等[16]引入廣義的DFT變換技術(shù),使HB-AFT方法可自動求解系統(tǒng)的準周期解. 最近,Didier等[17]提出的隨機HB-AFT方法能很好的用于包含隨機強非線性的轉(zhuǎn)子系統(tǒng)的非線性響應(yīng)分析. Zhang等[18]將同倫延拓技術(shù)嵌入HB-AFT方法,結(jié)合Hsu求解Floquet單值矩陣的離散方法,給出了一整套追蹤系統(tǒng)周期響應(yīng)和穩(wěn)定性分析的全頻域快速方法.

Duffing振子系統(tǒng)已經(jīng)成為分析非線性動力學行為的一種經(jīng)典模型,而且在許多情況下可以使用不同形式的Duffing方程對一些工程非線性問題進行定性分析[19,20].Duffing系統(tǒng)是研究陣發(fā)性混沌振動的一種主要模型,陣發(fā)性混沌現(xiàn)象在各類型Duffing系統(tǒng)中被大量研究[21-23].簡諧激勵Duffing系統(tǒng)(1)式的陣發(fā)性混沌現(xiàn)象被廣泛研究[11,24],不過缺少對其動態(tài)演化的相關(guān)分析.

(1)

鑒于以上背景,本研究從簡諧激勵Duffing系統(tǒng)陣發(fā)性混沌現(xiàn)象出發(fā),基于非線性動力學理論,采用HB-AFT方法結(jié)合Floquet理論,對非線性動力系統(tǒng)的典型陣發(fā)性混沌行為的演化展開研究,擬提供一種陣發(fā)性混沌研究的新方法.

1 理論與方法

1.1 HB-AFT方法求解周期響應(yīng)

不失一般性,求解響應(yīng)周期為2π的解,對于非線性系統(tǒng)

(2)

首先,進行響應(yīng)和非線性函數(shù)諧波平衡化過程：

(3)

把等式(3)代入系統(tǒng)(2),各階諧波平衡可得：

g(P,Q)=0

(4)

其中,P、Q分別表示響應(yīng)和非線性力諧波系數(shù).

把Q記為已知,可采用Newton-Raphson迭代求解不動點P,

J(i)(P(i+1)-P(i))+g(i)=0

(5)

其中,迭代Jacobian矩陣,

J=dg(P,Q)/dP

=?g(P,Q)/?P+?g(P,Q)/?Q·dQ/dP

(6)

迭代求解(5)式,只有(6)式中dQ/dP是未知的,下面通過AFT變換[22]給出該關(guān)系.

系統(tǒng)響應(yīng)x(t)的時域離散信息可由反有限傅里葉變換(Inverse Discrete Fourier transform, IDFT)給出：

(7)

這里Pk=ak+ibk,n=0,…,N-1,x(n)為x(t)在第n個時間點的值,N為時間離散點數(shù).

(8)

則Q可由有限傅里葉變換(Discrete Fourier transform, DFT)給出：

(9)

其中Qk=ck+idk,當n=0,φ為1, 否則φ等于2.

從(7)～(9)式可見,Qk是Pk的函數(shù),因此可以求得dQ/dP的顯式表達式.

1.2 嵌入弧長延拓的HB-AFT方法

以λ為控制參數(shù),由(5)式迭代求解會在轉(zhuǎn)向點失效[25],對此問題可用弧長延拓法把(4)式轉(zhuǎn)化為：

(10)

在初始條件(P(s0),λ(s0))=(P,λ)0下的求解問題.式(10)中λ為系統(tǒng)參數(shù)變量,s為曲線(P(s),λ(s))的弧長變量,(P,λ)0為(10)式一個已知解.

常用牛頓迭代校正法[25]求解(10)式,具體求解過程為：

首先,由已知解(P,λ)0預(yù)估一個解

(11)

然后,對(P,λ)1進行Newton-Raphson迭代,

(12)

1.3 穩(wěn)定性分析

對由HB-AFT方法求得的系統(tǒng)周期解,基于Floquet理論[7],由Hsu階躍函數(shù)法求得系統(tǒng)Floquet單值矩陣[17].

(13)

對于(13)式系統(tǒng),用ΔU擾動HB-AFT過程求得的T為周期的解U*(τ)如下：

d(U*(t)+ΔU)/dt=F(t,U*(t)+ΔU)

(14)

可得(14)式的線性化表達式

dΔU/dt=?F(t,U*(t))/?U(t)*·ΔU

=A(t,U*(t))ΔU

(15)

周期解U*(τ)的局部線性穩(wěn)定性可通過Floquet理論,由周期變系數(shù)微分方程(15)式ΔU的穩(wěn)定性來判斷.由Hsu階躍函數(shù)法,給出單值矩陣的計算公式：

(16)

就周期解分岔特性而言,周期為T周期解穩(wěn)定的條件為所有的Floquet乘子(即單值矩陣B的特征值λm)都位于復(fù)平面上的單位圓內(nèi),根據(jù)λm通過單位圓的方式可知周期解的失穩(wěn)的三種基本形式[7]：

(1)λm在+1處越出單位圓,系統(tǒng)周期響應(yīng)可能發(fā)生跨臨界分岔、對稱破缺分岔、折疊分岔；

(2)λm在-1處越出單位圓,則周期響應(yīng)發(fā)生倍周期分岔；

(3)若兩個共軛乘子離開單位圓,則周期響應(yīng)發(fā)生二次Hopf分岔.

一般來說,若系統(tǒng)發(fā)生折疊分岔、亞臨界Hopf分岔或亞臨界倍周期分岔后,系統(tǒng)的原穩(wěn)定周期吸引子并未發(fā)生轉(zhuǎn)遷[24],伴隨系統(tǒng)結(jié)構(gòu)的重整化將分別發(fā)生陣發(fā)Ⅰ、Ⅱ和Ⅲ混沌現(xiàn)象[7].

2 算例

對于(1)式Duffing系統(tǒng),選取文獻[11,24]參數(shù)c=0.25,k=-1,α=1,ω=1進行分析.當激勵幅值A(chǔ)=0時,對無阻尼未擾系統(tǒng)(1)式進行平衡點穩(wěn)定性分析,可知系統(tǒng)具有兩個中心(±1,0)和一個鞍點(0,0).

圖1 全局周期1解分支Fig.1 The trajectory of global period-1 branch

圖2 數(shù)值分岔圖,其中黑點、紅點分別為隨A向前、向后掃頻的數(shù)值積分結(jié)果Fig.2 Numerical bifurcation diagram, the black and red dots denotes A sweeping up and down

以激勵幅值A(chǔ)為控制參數(shù)(即式(10)中λ=A),取諧波次數(shù)K=10,式(12)迭代誤差取10-15,采用嵌入弧長延拓的HB-AFT方法追蹤系統(tǒng)的周期1解分支(相對激勵周期而言).如圖2所示,系統(tǒng)受擾后從平衡點出現(xiàn)三條周期解分支,與未擾系統(tǒng)類似,兩中心產(chǎn)生的周期1分支是穩(wěn)定的,鞍點受擾產(chǎn)生的周期1分支是不穩(wěn)定的,這與文獻[26]的解析分析結(jié)果一致.兩條穩(wěn)定的周期解在A1、A2點通過折疊分岔失穩(wěn),如表1、表2所示,周期1解分支的Floquet乘子在轉(zhuǎn)向點穿過+1軸離開單位圓.隨后,不穩(wěn)定的周期1解分支在A3、A4通過折疊分岔產(chǎn)生半穩(wěn)定[27]的周期解分支A3-A5、A4-A5.另一方面,從鞍點產(chǎn)生的周期解分支在A6經(jīng)折疊分岔產(chǎn)生較大振幅的穩(wěn)定周期軌線A6-A7.上述結(jié)果與數(shù)值分岔圖是吻合的.隨控制參數(shù)A的改變,如圖3所示系統(tǒng)的全局周期解共存特性及其穩(wěn)定性變化顯著.

表1 周期1分支在轉(zhuǎn)向點A1附近的Floquet乘子λmTable 1 Floquet multipliers λmof period-1 branch around the turning point A1

表2 周期1分支在轉(zhuǎn)向點A2附近的Floquet乘子λmTable 2 Floquet multipliers λmof period-1 branch around the turning point A2

圖3 系統(tǒng)在A激勵幅值取(a) 0.1,(b) 0.2,(c) 0.4和(d) 0.8時的穩(wěn)定(實線)、不穩(wěn)定(虛線)周期1解共存特性Fig.3 Coexistence of stable (solid) and unstable (dotted) period-1 solutions, excitation amplitude A taken as (a) 0.1, (b) 0.2, (c) 0.4 and (d) 0.8

如表1、表2所示,隨著控制參數(shù)A加大,兩條穩(wěn)定的周期解分支在AT=0.26497948248附近同時發(fā)生折疊分岔,這是由于系統(tǒng)的對稱性導(dǎo)致的.如圖4所示,系統(tǒng)折疊分岔失穩(wěn)后系統(tǒng)的穩(wěn)態(tài)響應(yīng)并未發(fā)生危險分岔跳躍到其他吸引子上[24],且在臨近不動點A1、A2位置經(jīng)歷逆切分岔(inverse tangent bifurcation)[8]緩慢通過,由此可以判斷系統(tǒng)此時觸發(fā)陣發(fā)Ⅰ型混沌振動模式,可以認為AT是此系統(tǒng)陣發(fā)Ⅰ型混沌的閾值.

圖4 激勵幅值A(chǔ)在0.16到0.32之間的全局周期1解分支和數(shù)值分岔圖Fig.4 The trajectory of global period-1 branch and numerical bifurcation diagram for excitation amplitude A taken 0.16 to 0.32

如圖5所示,隨著控制參數(shù)越過陣發(fā)I型混沌閾值A(chǔ)T,系統(tǒng)的時域響應(yīng)出現(xiàn)了規(guī)則與不規(guī)則偽隨機交替.當A值繼續(xù)增加,不規(guī)則運動所占區(qū)間增大,最后系統(tǒng)完全進入不規(guī)則運動,由圖6 Poincare映射可見,系統(tǒng)響應(yīng)在此過程中逐漸遠離分岔不動點A1、A2(圖6局部放大圖a3和a4),這也是陣發(fā)Ⅰ型混沌與陣發(fā)Ⅱ型、Ⅲ型混沌不同特征之一[8].另外,系統(tǒng)存在穩(wěn)定的周期1解與混沌解共存的狀態(tài),見圖6(c).

圖5 激勵幅值A(chǔ)取(a) 0.26498,(b) 0.26499,(c) 0.26501和(d) 0.27000時的時間歷程Fig.5 Time series for excitation amplitude A taking (a) 0.26498,(b) 0.26499,(c) 0.26501 and (d) 0.27000

圖6 激勵幅值A(chǔ)取(a) 0.26498,(b) 0.27000和(c) 0.45000時的Poincare映射,其中藍線為相軌線Fig.6 Poincare maps for excitation amplitude A taking (a) 0.26498,(b) 0.27000 and (c) 0.45000, where the blue line is phase trajectory

3 結(jié)論

以Duffing系統(tǒng)為研究對象,對經(jīng)典的單頻激勵Duffing系統(tǒng)的分岔特性進行了深入研究.發(fā)現(xiàn)該系統(tǒng)包含豐富的周期解分支共存和失穩(wěn)特性,且在一定條件下存在陣發(fā)I型混沌行為,具有一定的基礎(chǔ)理論價值.陣發(fā)性混沌演化由于涉及系統(tǒng)全局動力學特性的轉(zhuǎn)變,理論上分析困難.本研究針對該問題,從非線性動力學角度出發(fā),采用嵌入弧長延拓的HB-AFT方法結(jié)合Floquet理論,提供了一套陣發(fā)性混沌演化的研究方法.本文方法本質(zhì)上屬于頻域方法,可以快速而精確地自動追蹤系統(tǒng)周期解分支并判斷其穩(wěn)定性.通過系統(tǒng)全局周期解特性的分析,能夠?qū)ο到y(tǒng)經(jīng)典陣發(fā)類型的演化進行全局分析.