基于Tikhonov正則化改進的IHB法求解Mathieu-Duffing系統多重解*

王德亮， 劉濟科， 劉廣

1.中山大學航空航天學院， 廣東 深圳 518107

2.深圳市智能微小衛星星座技術與應用重點實驗室， 廣東 深圳 518107

現實中的各種振動系統都含有非線性因素（陳予恕，1992；陳樹輝，2007；Amabili，2008）。然而，早期的數學工具和振動理論通常會忽略掉系統的非線性因素而當作線性系統來求解（Devore，1998）。某些時候線性近似可能會超過容許誤差，甚至導致定性上的錯誤（Dowell et al.，1988）。隨著求解技術尤其是計算技術的發展，非線性振動理論開始逐步完善，其研究方法也日益完善。非線性振動研究可以分為定性和定量兩種方法，本文主要關注定量方法。定量方法主要有數值解法和解析方法，此外還有二者相結合的半數值半解析方法。數值方法，如Runge-Kutta法（徐加初等，2008）、Newmark-β法（劉廣等，2017）和精細積分法（富明慧等，2008）等。半解析方法有時域最小殘值法TMRM（Liu et al.，2021a；Liu et al.，2021b；Liu et al.，2022；Zheng et al.，2022a），同倫分析HAM（Liao et al.，2004）和增量諧波平衡法IHB（張丹偉等，2018；趙倩等，2018）。

IHB 法最早由Lau 和Cheung et al.于20 世紀80年代初提出（Ferri，1986），該方法在應用時不受小參數的限制，是一種適用于弱非線性系統和強非線性系統的半數值半解析方法（Cheung et al.，1990）。IHB 法具有概念直觀、易于應用的特點，廣泛地應用在非線性振動的各種研究中。近30 年來，涌現了大量的基于IHB 法或改進IHB 法的研究。即便如此，原始的IHB 法仍有以下一些局限：一方面，IHB 法在本質上是諧波平衡法和Newton-Raphson 迭代法的結合，其收斂性強烈地依賴于迭代初值的選取（黃建亮等，2021）；另一方面，對于參數激勵的非線性振動系統，可能同時具有多個周期解或極限環解，其中還包含穩定解和不穩定解。在IHB 法中，由于線性增量方程的雅可比矩陣可能是病態的，如果參數的初始猜測值遠離精確解，計算很可能不收斂（Liu et al.，2018；錢征文等，2011）。而系統同時存在多重解時，IHB法可能難以獲得全部的半解析解（Zheng et al.， 2022b）。

在后續的基于IHB 法的振動分析中，近些年發展了一些選取迭代初值的思路，比如：結合數值方法獲得數值響應，通過快速傅里葉變換（FFT）獲得系統的近似頻率和振幅，從而確定迭代初值（Ju et al.， 2020）。但這種思路對不穩定解無能為力。另外，可通過求解系統的近似線性系統的級數解作為迭代初值（Chen et al.， 2009）。但某些系統具有豐富的非線性特性，低階近似解無法準確揭示這些特性，從而造成方法失效。

本文擬通過在IHB 法的線性增量方程中引入Tikhonov 正則化，來提高IHB 法的收斂性能。并通過改進的IHB 法，求解了具有多重解的非線性參數激勵Mathieu-Duffing 系統。數值結果表明，改進的TIHB 不僅可以快速獲得高精度的半解析解，而且能通過選取不同的迭代初值獲得Mathieu-Duffing系統不同的穩定或不穩定周期解。

1 基于Tikhonov正則化的TIHB法

1.1 IHB法

考慮如下形式的非線性Mathieu-Duffing 系統（Shen et al.，2008）：

通過引入時間變換τ=ωt，方程（1）變為

其中上標“'”表示對無量綱時間τ的導數。對于系統的穩態周期響應，其半解析解為

接下來，通過IHB 法求解方程（3）所示的半解析解。IHB法的第一步是增量過程。假設x0為方程（1）的精確解，通過引入增量Δx和Δω，可以得到其鄰近點x=x0+ Δx和ω=ω0+ Δω， 代入式（2）得

引入方程

方程（4）可簡化為

其中R(x0)表示原始方程的殘差，ο(Δx，Δω)表示關于增量Δx和Δω的高階小量。僅保留關于Δx和Δω的一階小量，可得增量方程

截取Fourier 級數的前N項作為系統的近似解，得

將方程（8）代入方程（6）中，并引入Galerkin 平均過程

其中

積分式（9）并化簡為代數方程組，即可得到關于Δa、Δω的矩陣形式的增量方程

其中Ja和Jω為方程關于諧波系數和頻率的雅可比矩陣，Rˉ為殘差（5）積分之后的結果。方程（1）所示的Mathieu-Duffing 系統為已知周期的參數激勵系統，其穩態周期解的頻率為ω，因此無需求解。則式（9）可化為

其中

通過選取合適的迭代初值，并通過方程（12）求解諧波系數的迭代更新量Δa，然后根據以下方程進行更新

直到

則終止迭代。其中c為預先給定的收斂誤差。

1.2 引入Tikhonov正則化的改進IHB法

在某些的參數條件下，方程（11）中的Jacobi矩陣Ja和Jω可能有較大的條件數，即諧波系數和頻率的Jacobi矩陣可能是病態的。此時，如果諧波系數的初值選擇不當，上述原始的IHB 法的迭代可能不收斂。如何確定IHB 法的迭代初值一直是一個難題，常用的方法可以首先獲取系統的穩定數值解，然后通過對數值解進行快速傅里葉變換（FFT）來獲得系統的迭代初值。這個過程無疑是低效且繁瑣的，此外對于系統的不穩定解，這種方式還會失效。Tikhonov 正則化經常被用于反問題中來解決可能出現的病態問題，可嘗試在方程（12）中引入Tikhonov正則化，則方程（12）變為

進一步得

方程（16）中，‖ ·‖2表示2-范數，I為單位矩陣，λ為正則化參數。從方程（17）中可以看出，不同的正則化參數λ會獲得不同的迭代更新量Δaλ，一個合適的正則化參數可以平衡計算量和精度之間的矛盾。

常用的確定正則化參數λ的方法有L 曲線法（Ju et al.，2020）。L 曲線法的思想是，選取不同的正則化參數來計算獲得方程（16）中與，并且分別以和作為橫縱坐標繪制L 曲線，一般認為L 曲線最大曲率點對應的正則化參數為最優參數λOPT.令則曲率κ的表達式為

其中“'”表示對λ的導數。上述通過引入Tikhonov正則化改進的IHB法被稱為TIHB法。

2 數值算例

2.1 TIHB法求解Mathieu-Duffing系統

對于方程（1）所示的非線性參數激勵Mathieu-Duffing系統，其近似周期解可以通過方程（8）所示的截斷Fourier級數表示，即

當穩態周期解的頻率為參激頻率的一半，即響應周期為激勵周期的2倍，此時的穩態周期解稱之為2-周期解。Shen et al.（2008）的研究表明，當選取參數ζ= 0.125，β= 1，ω= 2， 3.6 ＜α＜5.264 566時，系統還同時存在1-周期解。上述研究表明，隨著參數的變化，該非線性參數激勵Mathieu-Duffing 系統同時存在多個周期解，且周期解中還有穩定解與不穩定解。

現有的方法在求解這些周期解時存在各種困難。對于數值方法，其無法求得不穩定解；而原始的IHB 法由于其收斂性強烈依賴于諧波初值的選擇，因此在追蹤多個不同解時，需要較為精準地給出不同的迭代初值，這無疑是一項困難的工作。接下來將嘗試通過本文提出的TIHB 法來解決上述多解問題。進一步選取α= 4，N= 40.

本文所提的TIHB 法和4 階Runge-Kutta（R-K）法獲得的2-周期對應的相圖，如圖1所示。從圖中可以看到，TIHB 法所得半解析解和R-K 法獲得的結果符合得非常好，且即便考慮諧波階數N= 40時，僅僅需要40 步迭代就可以收斂。圖2 給出了TIHB 法獲得的半解析解對應的殘差，也即方程（5）對應的系統殘差。從圖中可以看到，一個周期內的系統殘差在1011量級，這表明TIHB 法對應的半解析解具有非常高的精度。

圖1 α = 4時，2-周期響應的相平面圖Fig.1 Phase diagram of the 2-period response with α = 4

圖2 α = 4時，2-周期響應對應的殘差圖Fig.2 The residual diagram corresponding to 2-periodic response with α = 4

2.2 TIHB法的收斂性能

為了進一步驗證本文所提的TIHB 法相較于原始IHB 法的收斂性能，對方程（19）所示的半解析解的一階諧波系數賦初值，剩余的諧波系數初值全部設定為0。例如，令一階諧波系數a1和b1在參數區域[-2.5，2.5]內遍歷取值，然后分別通過原始IHB 法和本文所提的TIHB 法來計算結果。最終的收斂情況如圖3所示，其中空心圓表示該迭代初值收斂到系統的“0”解，實心點則表示收斂周期解，空白點則表示該初值不收斂。上述收斂的迭代初值都在40步之內使系統的殘差降到10-11以內；空白處的點則是經過500次迭代仍舊無法將殘差降到10-10以下的。對比圖3（a）和圖3（b）可以發現，本文提出的TIHB 法在此區域內每一個點都是可以收斂的，但是原始IHB 法只是在部分點收斂。上述結果表明，Tikhonov 正則化的引入，確實能夠顯著地增大傳統IHB法的收斂域。

圖3 原始IHB法和TIHB法的收斂性Fig.3 The convergence of the original IHB method and the TIHB method

結合前文的計算以及對應的收斂性分析，可以發現：相比于傳統的IHB 法，TIHB 法在不損失精度的前提下，其收斂性得到了極大的改善，或者說其收斂域得到了顯著的增大。事實上，在當前參數下，系統除了2-周期周期解還同時存在2個1-周期周期解。通過引入約束a2k-1=b2k-1= 0到方程（12）中，就可以進一步求得系統的1-周期解。圖4 為系統的1-周期解sol1 和sol2，其中周期解sol1為穩定周期解，周期解sol2為不穩定周期解。

圖4 系統的1-周期解對應的相圖Fig.4 Phase diagram of the 1-period response

同樣地，令一階諧波系數a1和b1在[-2.5，2.5]內遍歷取值，得到原始IHB 法和TIHB 法關于1-周期解的收斂結果，如圖5所示。此時系統的收斂情況分為3 種，分別是收斂到穩定周期解sol1、不穩定周期解sol2 和“0 解”。即便此時的情況比前文中2-周期周期解更復雜，仍舊可以看到TIHB 法各解的收斂域得到了顯著的擴大。從圖5（a）中還可以發現，原始IHB 法收斂域中，每個迭代初值似乎是隨機收斂到某個解的。

圖5 原始IHB法和TIHB法的收斂性Fig.5 The convergence of the original IHB method and the TIHB method

保持Mathieu-Duffing 系統中其他參數不變，取α= 5.1。在當前參數設定下，此時系統存在5個周期解，其中兩個是不穩定周期解，3 個是穩定周期解。通過TIHB 法并結合初值遍歷過程，可以依次獲得這5 個周期解。TIHB 法獲得的半解析解對應的相圖，如圖6 所示。其中，周期解sol2 和sol5為不穩定解，用虛線表示；其余解為穩定周期解，用實線表示。

圖6 α = 5.1時，5個2-周期解的相圖Fig.6 The phase diagrams of five 2-period response with α = 5.1

進一步，通過原始IHB 法和TIHB 法對α= 5.1的系統進行求解，結果如表1所示。表中的周期解sol1-sol5 和圖6 中的結果相對應。表1 顯示，原始IHB 法的收斂效果遠不如TIHB 法。原始IHB 法在所給的步長和區域內，無法迭代獲得解sol1 和解sol5，只得到了3個穩定周期解和一個“0”解；且各解的收斂域非常小，在2 601 個初值點域上，僅有36 個非零收斂點和321 個“0”解收斂點。相比之下，TIHB 法可以很好地計算出5 個不同的周期解，在同樣的遍歷區域內，未能收斂的初值僅有80個。

表1 α = 5.1時，IHB和TIHB法的收斂情況Table 1 The convergence of the IHB and the TIHB method with α = 5.1

3 結 論

通過在原始的增量諧波平衡法的迭代過程中引入Tikhonov 正則化，提出了一種改進的IHB 法，并借助該方法成功求解了具有多重解的Mathieu-Duffing 系統。和原始的IHB 法相比，可得出以下結論：（1）TIHB 法可以快速獲得具有高精度的半解析解，其收斂性能得到了顯著的改善；（2）在多組不同的參數設定下，結合初值遍歷的方式，本文所提的TIHB 法能夠準確找出Mathieu-Duffing 方程中同時存在的多個周期解，而原始的IHB 法僅能獲得部分周期解。