Duffing系統共存周期解的穩定性與分岔演化

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

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

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

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

1 Duffing方程的延拓打靶法

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

為了研究系統參數變化時周期解的穩定性及其分岔演化規律，本文將延續算法和一般打靶法結合求解系統的周期解，并對其進行追蹤分析[10]。

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

其中，求解系統的周期解需要構造系統的Poincaré映射，以此將周期解求解問題轉化為滿足下列條件的映射不動點曲線問題，即：

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

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

牛頓迭代格式如下：

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

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

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

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

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

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

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

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

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

圖1 多初值分岔圖

由圖1可見穩定周期三解P3a的窗口。當ω減小時，P3a經叉式分岔產生兩個反對稱的周期三解P3b（灰線）和P3c（黑線），然后經倍周期分岔序列嵌入混沌。當ω增加時，P3a直接轉遷為混沌。繼續增加ω，當ω=1.833時，混沌吸引子再次激變為兩部分，即產生兩個反對稱的混沌吸引子。隨后混沌吸引子經倍周期倒分岔退出混沌，最終進入穩定的反對稱周期一解的窗口。

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

圖2 周期一解的延拓追蹤

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

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

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

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

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

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

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

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

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

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

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

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

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

圖6 相圖和Poincaré映射圖

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

4 結語

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

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

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

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