張錦濤，王 昕，呂小紅，金 花

（蘭州交通大學 機電工程學院，蘭州 730070）

Duffing 系統(tǒng)是為描述機械系統(tǒng)中彈簧硬化效應而建立的一類典型非線性系統(tǒng)，很多工程中的振動問題可以使用該系統(tǒng)來研究。Duffing 系統(tǒng)雖然是比較簡單的非線性系統(tǒng)，但是卻包含豐富的非線性特性，如周期運動的跳變、周期解的共存、混沌激變等。目前，對于Duffing 系統(tǒng)的研究已有很多成果。文獻[1]中提出了一種簡單的周期參數(shù)切換方法，該方法可以找到廣義Duffing系統(tǒng)中任何可以數(shù)值逼近的穩(wěn)定極限環(huán)。文獻[2]中研究了一類具有不連續(xù)非線性阻尼的Duffing方程周期解問題，證明其周期解只能是諧波解和次諧波解，然后利用龐加萊-波爾定理研究了該方程的解存在唯一性。文獻[3]中研究了雙勢阱Duffing 方程混沌與周期吸引子的對稱破缺激變。文獻[4]中運用諧波平衡法計算系統(tǒng)的頻率-振幅響應曲線，再通過預估-校正算法找到了系統(tǒng)隨參數(shù)變化的跳躍區(qū)間，研究了其跳躍機理。文獻[5]中采用數(shù)值方法計算雙參數(shù)平面上的最大Lyapunov 指數(shù)，得到分岔曲線，分析了系統(tǒng)在雙參數(shù)平面上的分岔及其全局特性。文獻[6]中基于非線性振動系統(tǒng)的主要特征提出一種新的識別Duffing 非線性系統(tǒng)剛度的方法。文獻[7]中應用諧波平衡法求解Duffing系統(tǒng)的振動方程，研究了常數(shù)激勵與簡諧激勵聯(lián)合作用下系統(tǒng)的骨架曲線及幅頻響應特性。文獻[8]中通過多尺度法得到系統(tǒng)的近似解析解，研究在多頻激勵下同時發(fā)生主共振和1/3次亞諧共振的動力學行為與穩(wěn)定性。文獻[9]中運用多尺度法求解一類含五次非線性恢復力的Duffing系統(tǒng)的幅頻響應方程，分析系統(tǒng)的共振與分岔特性。以上文獻多運用不同解析法對Duffing 方程的周期解進行計算，進而研究了系統(tǒng)中存在的多種非線性特性，但這些研究對于系統(tǒng)中共存周期解的產(chǎn)生原因及其分岔演化過程卻鮮有提及。因此，系統(tǒng)中共存周期解的存在區(qū)域及其分岔演化規(guī)律沒有完全被揭示。

本文將延續(xù)算法和打靶法相結(jié)合，延拓追蹤系統(tǒng)中共存的周期解并應用Floquet 理論判斷其穩(wěn)定性及分岔。最后應用簡單胞映射法計算系統(tǒng)共存周期解的吸引域，并結(jié)合相圖和Poincaré 映射圖研究周期解的吸引域隨系統(tǒng)參數(shù)變化的侵蝕演變過程。

1 Duffing方程的延拓打靶法

式中，無量綱參數(shù)α為阻尼系數(shù)，β為剛度系數(shù)，ξ為非線性項系數(shù)，f為外激勵幅值，ω為激勵頻率。為方便數(shù)值計算，令=v，則方程可降階為：

為了研究系統(tǒng)參數(shù)變化時周期解的穩(wěn)定性及其分岔演化規(guī)律，本文將延續(xù)算法和一般打靶法結(jié)合求解系統(tǒng)的周期解，并對其進行追蹤分析[10]。

對于Duffing系統(tǒng)(2)，用X表示狀態(tài)向量[x,v]T,用μ表示參數(shù)向量[α,β,ξ,ω]T,可簡記為：

其中，求解系統(tǒng)的周期解需要構(gòu)造系統(tǒng)的Poincaré映射，以此將周期解求解問題轉(zhuǎn)化為滿足下列條件的映射不動點曲線問題，即：

式中：P(μ,X)是X的Poincaré映射。

求解式(4)這樣的非線性方程組，一般使用Newton-Raphson 方法進行迭代求解。延拓打靶法的基本思想是將參數(shù)空間[μ0,μ] 進行離散，當μ=μ0時求解式(4)可以得到X(μ0)，記為X0。將X0作為μ=μ1時的近似值X01，由式(4)進一步迭代可以得到X(μ1)，記為X1。以此類推得到一個近似的不動點X0k。此時可以用Newton-Raphson 方法進一步迭代求解，求得此時的不動點為Xk。由式(4)可得：

牛頓迭代格式如下：

針對部分老化管線結(jié)垢、腐蝕嚴重，注水生產(chǎn)時壓力損失大、刺漏狀況頻發(fā)等問題，2017年完成17條腐蝕老化注水管線、1條注水干線的改造，將管線材質(zhì)由復合管改為無縫鋼管，改造總長度16.3 km，涉及17口注水井，對應注配間柱塞泵降壓 0.5 MPa，節(jié)電 4.2×104kWh。

式中：i為迭代次數(shù)，實際計算中P(Xk)和DP(Xk)可由式(3)和式(7)以[X0k,I]為初始值積分一個Poincaré映射周期T得到。

為了提高算法的效率，對延續(xù)過程中的近似值先進行預估，再迭代校正。將式(4)轉(zhuǎn)化為一般常微分方程的初值問題：

求解得到Xk時，運用Euler 法可以對μ=μk+1時的不動點進行預估得到X0k+1。

式中：Δμ為系統(tǒng)參數(shù)的變化量，?G(μk,Xk)/?X可由式(3)、式(5)、式(7)聯(lián)合求解得到。?G(μk,Xk)/?μ可由式(3)和式(10)以[Xk,0]為初始值積分一個Poincaré 映射周期T得到。本文計算中選取激勵頻率ω為延拓參數(shù)，其余參數(shù)均為常數(shù)。

上述計算過程中，在求得不動點的同時得到了Poincaré 映射的Jacobi 矩陣DP(Xk)的特征值，根據(jù)Floquet 理論可判斷周期解的穩(wěn)定性。當所有Floquet 乘子都在復平面上的單位圓內(nèi)時周期解穩(wěn)定，存在一個Floquet乘子在單位圓外時周期解不穩(wěn)定。隨著參數(shù)的變化，當存在一個特征乘子從單位圓穿出時，周期解發(fā)生分岔而失穩(wěn)，可從乘子穿出單位圓的方式識別其分岔類型。

2 基于延拓打靶法的分岔分析

Duffing 系統(tǒng)存在豐富的動力學行為。為研究周期解的分岔演化過程，取參數(shù)：α=0.83，β=1.0，ξ=0.5，f=2.0，數(shù)值計算系統(tǒng)的多初值分岔圖如圖1所示。

從圖1中可以看出，在ω∈[0.8，0.965]時系統(tǒng)中只有穩(wěn)定的周期一解P1a（灰線）。隨著ω的增加，當ω=0.965 238 6 時，系統(tǒng)出現(xiàn)了新的穩(wěn)定周期一解P1b（黑線），此時P1b的Floquet 乘子為(0.004 5，0.9997)，由Floquet 理論可知，P1b發(fā)生了鞍結(jié)分岔。在ω∈[0.965，1.032]內(nèi)，穩(wěn)定的周期一解P1a和P1b共存。隨后，周期一解P1b經(jīng)過叉式分岔產(chǎn)生兩個穩(wěn)定的反對稱周期一解P1c和P1d。在ω∈[1.032，1.041]內(nèi)，反對稱的周期一解P1c和P1d以及周期一解P1a共存。當ω=1.040 942 6 時，周期一解P1a的Floquet乘子為(0.999 9，0.006 7)，由Floquet 理論可知，周期一解P1a發(fā)生鞍結(jié)分岔。當ω＞1.041 時，周期一解P1a消失，系統(tǒng)中只有反對稱的周期一解P1c和P1d共存。P1c和P1d隨著ω的增加經(jīng)倍周期分岔序列進入混沌，此時，兩個反對稱的混沌吸引子共存。當ω=1.085 時，兩個反對稱混沌吸引子融合激變，生成一個全新的混沌吸引子。

圖1 多初值分岔圖

由圖1可見穩(wěn)定周期三解P3a的窗口。當ω減小時，P3a經(jīng)叉式分岔產(chǎn)生兩個反對稱的周期三解P3b（灰線）和P3c（黑線），然后經(jīng)倍周期分岔序列嵌入混沌。當ω增加時，P3a直接轉(zhuǎn)遷為混沌。繼續(xù)增加ω，當ω=1.833時，混沌吸引子再次激變?yōu)閮刹糠郑串a(chǎn)生兩個反對稱的混沌吸引子。隨后混沌吸引子經(jīng)倍周期倒分岔退出混沌，最終進入穩(wěn)定的反對稱周期一解的窗口。

為了深入研究系統(tǒng)周期解的穩(wěn)定性與分岔，應用延拓打靶法追蹤周期一解在ω∈[0.9，1.05]區(qū)間內(nèi)的分岔演化。延拓計算之前，首先在系統(tǒng)的狀態(tài)空間選擇一個考察區(qū)域Ω=(x，v)=(-3≤x≤3，-3≤v≤3)，然后將狀態(tài)平面Ω劃分為100×100=10 000 個小網(wǎng)格。以每個小網(wǎng)格的中心點為初始值應用打靶法計算系統(tǒng)的周期解，則共有10 000 個初始值需要考察。當ω=1 時，求得系統(tǒng)有3 個共存的周期解（2 個穩(wěn)定的周期一解(0.031 239 28，1.746 523 36)和(-1.028 217 48，0.657 523 80)，1個不穩(wěn)定的周期一解(-1.066 474 56，1.417 145 48)）。分別以每個周期解為初值進行延拓計算，結(jié)果如圖2所示。圖中實線表示穩(wěn)定的周期解，“+”標記線表示不穩(wěn)定的周期解。

圖2 周期一解的延拓追蹤

由圖2可知，在ω∈[0.965，1.032]內(nèi)，系統(tǒng)存在兩個穩(wěn)定的周期一解P1a和P1b，以及一個不穩(wěn)定的周期一解UP1a。在ω∈[1.032，1.041]內(nèi)，3 個穩(wěn)定的周期一解P1a、P1c和P1d，以及兩個不穩(wěn)定的周期一解UP1a和UP1b共存。在ω∈[1.041，1.05]內(nèi)，兩個穩(wěn)定的周期一解P1c和P1d，以及一個不穩(wěn)定的周期一解UP1b共存。當ω=1.041時，穩(wěn)定的周期一解P1a與不穩(wěn)定的周期一解UP1a相撞后消失，由前面的分析可知，此時系統(tǒng)發(fā)生鞍結(jié)分岔。當ω=0.965時，穩(wěn)定的周期一解P1b與不穩(wěn)定的周期一解UP1a經(jīng)鞍結(jié)分岔消失。由此可見，鞍結(jié)分岔是Duffing系統(tǒng)中共存吸引子消失的重要原因。當ω=1.032 時，穩(wěn)定的周期一解P1b失穩(wěn)產(chǎn)生不穩(wěn)定的周期一解UP1b，同時通過叉式分岔產(chǎn)生兩個穩(wěn)定的反對稱周期一解P1c和P1d。

以ω=1.41時的8個共存周期三解為初值應用延拓打靶法延拓追蹤每個周期三解在ω∈[1.41，1.70]內(nèi)的分岔演化如圖3所示。圖中穩(wěn)定的周期解用實線表示，不穩(wěn)定的周期解用不同符號的虛點區(qū)分。由圖3可知，在ω∈[1.444，1.632]內(nèi)，一個穩(wěn)定的周期三解P3i與5個不穩(wěn)定的周期三解P3d至P3h共存。圖3(b)為圖3(a)黑框中不穩(wěn)定周期三解的局部放大。當ω減小至1.444 時，周期三解P3i失穩(wěn)轉(zhuǎn)遷為不穩(wěn)定的周期三解P3c，同時通過叉式分岔產(chǎn)生兩個穩(wěn)定的反對稱周期三解P3a和P3b，在此分岔過程中，周期三解P3d至P3h的模式及穩(wěn)定性保持不變。因此，在ω∈[1.41,1.444]內(nèi)，系統(tǒng)存在兩個穩(wěn)定的反對稱周期三解P3a和P3b，以及6個不穩(wěn)定的周期三解P3c至P3h。圖4(a)至圖4(h)依次給出了ω=1.41 時周期三解P3a至P3h的相圖和Poincaré 映射圖。當ω增加至1.631 978 7時，穩(wěn)定的周期三解P3i與不穩(wěn)定的周期三解P3d碰撞之后消失，此時P3i的Floquet 乘子為(0.999 6，0.000 1)，結(jié)合圖1并由Floquet 理論可知，周期三解P3i通過鞍結(jié)分岔直接進入混沌。因此，在ω∈[1.632，1.700]內(nèi)，4個不穩(wěn)定的周期三解P3e至P3h與混沌吸引子共存。在周期三解的窗口中，由于穩(wěn)定的周期三解和多個不穩(wěn)定的周期三解共存，導致系統(tǒng)的穩(wěn)定性降低。

圖3 周期三解的延拓追蹤和局部放大

圖4 ω=1.41時共存周期三解的相圖和Poincaré映射圖

圖4 ω=1.41時共存周期三解的相圖和Poincaré映射圖

由以上分析可知，Duffing系統(tǒng)在分岔圖的周期解窗口中，穩(wěn)定和不穩(wěn)定周期解的共存現(xiàn)象普遍存在，而不穩(wěn)定周期解經(jīng)過較長時間后會趨于穩(wěn)定的周期解。當穩(wěn)定周期解共存時，系統(tǒng)的響應與初值條件密切相關。采用周期解的吸引域能夠更為清晰地描述系統(tǒng)周期的初值響應和全局特性。因此，下文將根據(jù)簡單胞映射法計算穩(wěn)定周期解的吸引域，進一步揭示系統(tǒng)隨參數(shù)變化時周期解吸引域的演化規(guī)律以及系統(tǒng)的全局分岔特性。

3 共存吸引子的吸引域演化

當多個吸引子共存時，系統(tǒng)的動力學行為受初始條件的影響較大。選用第2 節(jié)所用參數(shù)，對系統(tǒng)的全局動力學特性進行研究。在系統(tǒng)的相空間中選取一個考察區(qū)域H={-4.0≤x≤4.0，-4.0≤v≤4.0}，將其劃分為800×800 個胞。以每個胞的中心點為初值，應用簡單胞映射法求解圖2中穩(wěn)定周期一解的吸引域，揭示激勵頻率ω增大時系統(tǒng)周期解的吸引域演化與分岔過程，如圖5所示。

由第2 節(jié)中的分析可知，在ω∈[0.9，0.965]時系統(tǒng)只存在一個穩(wěn)定的周期一解P1a。圖5(a)為ω=0.950 時的吸引域，其中白色“×”號表示周期一吸引子P1a，黑色區(qū)域為其吸引域，此時系統(tǒng)在考察區(qū)域H 內(nèi)只存在一個周期解。對應的相圖和Poincaré 映射圖如圖6(a)所示。

隨著ω增大，當ω=0.965 時系統(tǒng)經(jīng)鞍結(jié)分岔產(chǎn)生新的穩(wěn)定周期一解P1b，并與P1a共存。圖5(b)給出了ω=0.968時共存吸引子的吸引域，其中白色“×”號表示周期一吸引子P1a（94.09%），白色“+”號表示周期一吸引子P1b，灰色區(qū)域為P1b的吸引域（5.91%）。對比圖5(a)，可見在H 內(nèi)P1a原本完整的黑色吸引域中出現(xiàn)了P1b的灰色吸引域，即P1a在部分初值范圍內(nèi)發(fā)生周期跳躍轉(zhuǎn)遷為P1b，P1a的吸引域被P1b的吸引域侵蝕。當ω=1.0時共存吸引子P1a和P1b的吸引域如圖5(c)所示，可見隨著ω增大，系統(tǒng)中灰色吸引域的面積不斷增大，P1a的吸引域被進一步侵蝕。此時P1b的穩(wěn)定性增強，但是其吸引域面積（15.79%）遠小于P1a的吸引域面積（84.21%），P1a的穩(wěn)定性仍然最強。P1a和P1b對應的相圖和Poincaré映射圖如圖6(b)所示，兩運動軌跡的振動幅值不同，P1b的振幅相對較小。進一步增加ω至1.025時，P1a（85.71 %）和P1b（14.29 %）的吸引域如圖5(d)所示。由圖可見，P1a在H內(nèi)的右側(cè)被P1b部分侵蝕，而P1b集中分布的灰色吸引域被P1a侵蝕，即原來的P1b在部分初值范圍內(nèi)轉(zhuǎn)遷為P1a，兩個吸引域嵌套在一起，其初始條件的微小擾動有可能使系統(tǒng)出現(xiàn)不同的周期響應，P1b的穩(wěn)定性明顯降低。

繼續(xù)增大ω，當ω=1.032時，P1b經(jīng)叉式分岔產(chǎn)生兩個反對稱的周期一吸引子P1c和P1d。當ω＞1.032時，系統(tǒng)中3個吸引子P1a、P1c和P1d共存。圖5(e)給出了ω=1.035時不同吸引子的吸引域，其中“×”號表示P1a，黑色為其吸引域（65.22%），“◇”號表示P1c，灰色為其吸引域（19.92%），“△”號表示P1d，白色為其吸引域（14.86%）。圖中兩個反對稱吸引子P1c和P1d的吸引域互相嵌套并纏繞在一起，兩個吸引子的穩(wěn)定性較弱，而P1a的吸引域面積較大，穩(wěn)定性較強。3個共存吸引子的相圖和Poincaré映射圖如圖6(c)所示。隨著ω增大，P1a的吸引域（56.76%）被纏繞在一起的P1c（21.27%）和P1d（21.97%）的吸引域分別從下往上和從右往左逐漸侵蝕，P1a的穩(wěn)定性逐漸降低，如圖5(f)所示。

當ω增加至1.041時，P1a發(fā)生鞍結(jié)分岔而消失，系統(tǒng)中只剩下P1c和P1d共存。ω=1.042 時，兩個共存吸引子的吸引域如圖5(g)所示。由圖可見，系統(tǒng)在H 內(nèi)只存在P1c（48.99%）和P1d（51.01%）的吸引域，原本P1a的吸引域被完全侵蝕。圖5(h)給出了ω=1.05 時兩個反對稱吸引子P1c（50.81 %）和P1d（49.19%）的吸引域。對比圖5(g)可見，隨著ω增大，P1c的吸引域和P1d的吸引域面積逐漸接近，并且纏繞在一起，具有相似的分形結(jié)構(gòu)。此時系統(tǒng)對初值極度敏感，初值的微小變化會使系統(tǒng)運動到不同的周期軌道上。P1c和P1d對應的相圖和Poincaré映射圖如圖6(d)所示。

圖5 隨激勵頻率ω增大吸引域的演化與岔過程

圖6 相圖和Poincaré映射圖

由以上分析可知，系統(tǒng)中出現(xiàn)共存吸引子時，隨著ω增大，新吸引子的穩(wěn)定性會逐漸增強，如圖5(b)和圖5(c)所示，P1b的吸引域從中心向外擴展，面積逐漸增大；但是在分岔點附近吸引子的穩(wěn)定性會減弱，如圖5(d)和圖5(f)所示，P1b和P1a在發(fā)生分岔前原本完整的吸引域被與其共存的其余吸引子的吸引域所侵蝕。鞍結(jié)分岔會引起系統(tǒng)新的周期解的出現(xiàn)或消失，從而導致考察區(qū)域內(nèi)周期解類型的突變及吸引域拓撲結(jié)構(gòu)的改變，其穩(wěn)定性也隨之發(fā)生變化，如圖5(a)和圖5(b)以及圖5(f)和圖5(g)所示。當系統(tǒng)中由于叉式分岔出現(xiàn)、反對稱吸引子共存時，反對稱吸引子的吸引域拓撲結(jié)構(gòu)相似，互相纏繞，導致系統(tǒng)的穩(wěn)定性降低。吸引域的演化和侵蝕對系統(tǒng)周期解的穩(wěn)定性有重要的影響。

4 結(jié)語

本文應用打靶法和Floquet 理論計算得到Duffing系統(tǒng)穩(wěn)定和不穩(wěn)定的周期解，并對其進行延拓追蹤，分析了共存周期解的分岔演化過程及其全局動力學特性。得到了以下結(jié)論：

(1)Duffing 系統(tǒng)中多個穩(wěn)定及不穩(wěn)定周期解的共存現(xiàn)象較為普遍，并具有一定反對稱性。當反對稱吸引子通過倍周期分岔進入或退出混沌時，會發(fā)生混沌吸引子的內(nèi)部激變。

(2)穩(wěn)定周期解與不穩(wěn)定周期解碰撞而發(fā)生鞍結(jié)分岔是共存吸引子出現(xiàn)和消失的重要原因，并且鞍結(jié)分岔會導致吸引子的吸引域拓撲結(jié)構(gòu)發(fā)生改變。當系統(tǒng)由于叉式分岔產(chǎn)生、兩個反對稱吸引子共存時，其吸引域拓撲結(jié)構(gòu)相似，彼此嵌套纏繞在一起，會導致系統(tǒng)的穩(wěn)定性降低。Duffing系統(tǒng)的周期吸引子主要通過倍周期分岔或鞍結(jié)分岔進入混沌。

(3)當系統(tǒng)中出現(xiàn)共存吸引子時，穩(wěn)定周期解受初始條件的影響較大。隨著系統(tǒng)參數(shù)的變化，新出現(xiàn)吸引子的穩(wěn)定性會逐漸增強，但在分岔點附近，即將發(fā)生分岔的吸引子穩(wěn)定性會降低。