余弦廣義Padé逼近法及其在強非線性振子周期解求解中的應用

李震波， 唐駕時

(1.南華大學 數理學院，湖南 衡陽 421001；2.湖南大學 機械與運載工程學院，長沙 410082)

近年來，非線性振動系統的定量研究一直是學者們研究的熱門課題。隨著計算機技術的發展和對自然界認識的不斷加深，對振動系統的定量研究逐步從弱非線性上升到了強非線性，并相繼提出了很多有效的定量研究方法，如橢圓函數攝動法[1]、橢圓函數Lindsted-Poincare法[2-3]、非線性時間變換法[4-7]、同倫分析法[8-10]、變分迭代法[11]以及攝動-增量法[12]等。此外，利用Padé逼近來求解強非線性振子同、異宿軌也取得了一些成果。Emaci等[13]基于如下的經典Padé逼近式，構造了非線性薛定諤方程的同宿解。

(1)

由于經典Padé逼近式形式固定，在求解某些非線性振子時，精度不理想。因此，學者們在經典Padé逼近式的基礎上進行了改進，構造了類Padé逼近式，有效地提高了求解精度。如Mikhlin在文獻[14]中構造了如式(2)所示的一類同時含有多項式函數和指數函數的類Padé逼近式，并基于該類Padé逼近式求解了非線性薛定諤方程、洛倫茲方程以及Duffing方程的同宿解。

(2)

(3)

張琪昌等在文獻[15-16]中提出了一種改進的類Padé逼近方法，用指數函數代替了式(2)中的多項式函數，構造了如式(3)所示的類Padé逼近式，并求解了具有立方非線性的Duffing系統同、異宿解。上述各種類Padé逼近雖然對逼近式的形式做出了各種改進，但卻沒有對Padé逼近的定義做出推廣，而且已有的各類Padé逼近式本身均為非周期函數且通常在無窮遠處存在極限。因此，從現有的文獻來看，沒有直接利用Padé逼近來求解微分方程周期解的相關報道。

李震波等在文獻[17]中提出了一種廣義Padé逼近方法，該方法將經典Padé逼近式的分子和分母由多項式函數推廣至任意連續函數構成函數級數，該推廣不僅從理論上統一了現有的各種Padé和類Padé逼近，也為構造更多廣義形式的Padé逼近式來逼近微分方程的周期解提供了思路和方法。本文正是基于廣義Padé逼近方法，構造了一類新形式的廣義Padé逼近式。該逼近式的分子和分母均為余弦函數構成的函數級數，其本身具有周期性，利用該式來直接逼近一個表達式未知但冪級數形式已知的周期函數，如強非線性振子的周期解，取得了較理想的逼近結果，不僅得到了被逼近函數的高精度解析表達式，同時也較準確的預測了周期，從而擴大了Padé逼近方法在振動領域中的適用范圍?？偟膩砜?，該方法具有以下特點：① 繼承了Padé逼近方法求解過程簡單直接的優點，便于進行計算機編程計算；②所得之解的精度不受非線性項系數大小和初始振幅大小的影響，適用于強非線性；③ 該方法并不局限于某一個特定的系統，而是有著較廣的適用范圍。因此，對于廣義Padé逼近方法的研究和推廣具有一定的實際意義和理論價值。

1 一類新的廣義Padé逼近式

早在1821年，Cauch就提出了Padé逼近的理論。由于在計算數學、量子力學、量子場論、原子和分子物理及控制理論等領域中有重要應用，Padé逼近自20世紀70年代以來在理論、方法及應用上得到了很大的發展，出版了許多專著[18-20]。Padé逼近方法是從冪級數出發獲得有理函數逼近式的一種十分簡潔而且非常有效的方法。其基本思想就是對于一個給定的形式冪級數，構造一個有理函數，稱其為Padé逼近式，使其Taylor展開有盡可能多的項與原來的冪級數相吻合。在經典Padé逼近理論的基礎上，李震波等提出了一種廣義Padé逼近方法，首先定義如下記號

(4)

(5)

定義1(廣義Padé逼近) 設f(z)為由下述形式冪級數所定義的函數

若存在有理分式函數PL(z)/QM(z)∈G(L,M)滿足

f(z)-PL(z)/QM(z)=O(zL+M+1)

(6)

則稱PL(z)/QM(z)為f(z)在G(L,M)中的廣義Padé逼近式，記為GPA[L/M]f(z)，或簡記為GPA[L/M]。

與經典Padé逼近定義相比，定義1將PL(z)和QM(z)由多項式函數推廣至任意函數組成的函數級數，使得在構造廣義Padé逼近式時具有更廣泛的選擇范圍，也使得經典Padé逼近定義更完備，擴大了Padé逼近的適用范圍，從理論上統一了現有的Padé和各種類Padé逼近方法。當選取PL(z)和QM(z)為多項式函數時，廣義Padé逼近退化為經典Padé逼近。眾所周知，經典Padé逼近式由于本身為非周期函數且通常在無窮遠處存在極限，因此，從現有的文獻來看，沒有直接利用Padé逼近來求解微分方程周期解的相關報道，然而定義1的提出，為構造周期性的廣義Padé逼近式提供了新的思路和方法，基于該定義，我們通過引入余弦函數，構造了如下形式的廣義Padé逼近式

(7)

式中：α2i和β2i為待定的廣義Padé系數；ω為廣義Padé角頻率，周期為π/ω。利用上式來逼近微分方程的周期解，取得了較理想的逼近結果，詳見下文的求解和算例。

2 強非線性自治振子的周期解

考慮如下形式的自治振動系統

(8)

當g(x)為非線性函數時，上述方程需要利用攝動法求解，但攝動法的不足之處在于只能求解弱非線性的情形。當g(x)為高階多項式函數或有理函數時，求解上述方程的周期解仍然是非常值得研究的課題。為此，本文基于廣義Padé逼近方法，研究了當g(x)為任意非線性函數時，式(8)的周期解。

用如下形式的冪級數來表示式(8)的解

(9)

將式(9)代入式(8)，即可求得上述冪級數解的前N項系數關于初值的表達式。令周期軌道與x軸的交點為初值點，即a0=a，a1=0。則冪級數解的系數ai可表述為如下形式的代數多項式

a2n=Pn(a0,cm),a2n+1=0,n≥0

接下來要利用上述冪級數的解來確定廣義Padé逼近式(7)的系數α2i和β2i。由于廣義Padé角頻率ω未知，若直接利用Padé逼近方法，容易造成求解困難。因此，本文在構造了一類新的廣義Padé逼近式的同時，也對經典Padé逼近方法的步驟進行了相應的改進，使之更有效的逼近一類強非線性振子的周期解，并較準確的預測了周期。

以求解初值為a0=a，a1=0的周期軌為例，該方法的具體步驟如下：

步驟1令廣義Padé角頻率ω為主動變量，而不再是未知數，其增量為Δω。

步驟2根據經典Padé逼近的方法，將式(7)右端在原點處泰勒展開，與式(9)中的同階項進行比較，從中得到如下形式的L+M+1個方程

Hk(α2i,β2i;ω0+Δωn,a0)=0

(10)

式中：α2i和β2i為待定的廣義Padé系數；a0為所求周期軌的初值；n為迭代次數；ω0為ω的迭代初值。

步驟3根據具體情況設置ω0的值(也可直接從零開始迭代)，并根據所求結果的精度要求設置增量Δω的值，通常情況下可令Δω=0.000 1。設置好ω0和Δω的值后，代入式(10)，并迭代求解α2i和β2i，每迭代一次可得一組解。

步驟4確定目標函數，即何時終止迭代。由于求解過程中選取的初值點為周期軌道與x軸交點，則當軌道到達與x軸的另一個交點時，正好經過半個周期，即π/2ω。令該周期軌道與x軸的另一個交點為B(b0,0)，則b0可由式(11)確定

(11)

將t=π/2ω代入式(7)可得其值為α0，因此令如下函數為目標函數

Z=α0-b0

(12)

將第三步所求得的α0代入式(12)，若使得Z<1×10-N(N∈Z+為精度控制指數)，則當次迭代所求得的α2i和β2i為該周期軌道的廣義Padé系數，代入式(7)即得其廣義Padé解，周期T=π/(ω0+Δωn)(其中n為迭代次數)；否則，繼續迭代。我們稱上述方法為第一類改進的廣義Padé逼近方法。

通過對Padé逼近方法的求解過程進行上述改進，可較理想的求得g(x)為多項式函數時系統式(8)的周期解。然而當g(x)為有理函數時的，對于部分周期較大的周期軌，上述方法所得之結果不能令人滿意，因此我們對上述求解過程進行了再一次改進，通過引入新的目標函數并對求解結果進行分段描述等手法，得到了當系統式(8)為對稱有理振子時的高精度周期解。其具體步驟為：

步驟1仍然令廣義Padé角頻率ω為主動變量且以ω0為初值，Δω為增量；并利用相同的方法得到式(10)。

步驟2將t=π/2ω代入式(7)可得其值為α0，因此令

α0=b0

(13)

其中b0可由式(14)確定

(14)

步驟3令如下函數為目標函數

(15)

h(0,c0)-h(a0,0)=0

(16)

步驟4將相關數據代入式(10)并進行迭代求解。將每次迭代所求得的α2i和β2i代入目標函數式(15)進行判別，若使得Z<1×10-N(N∈Z+為精度控制指數)，則當次迭代所求得的α2i和β2i為該周期軌道的廣義Padé系數，并將其表述出如下分段函數即得其廣義Padé解，周期T=π/(ω0+Δωn)(其中n為迭代次數)；否則，繼續迭代。

(17)

通過對廣義Padé逼近方法的上述二次改進，可求得系統式(8)為對稱有理振子時的高精度周期解，擴大了該方法的有效范圍，算例證明如下。

3 算 例

3.1 Helmholtz-Duffing振子的周期解

考慮如下形式的Helmholtz-Duffing振子

(18)

當c2c3=0時，該振子的周期解可用橢圓函數進行精確描述；本文則僅討論c2c3≠0時的情形。首先，令式(18)的冪級數為式(9)所示，將式(9)代入式(18)即可求得式(9)的各項系數為

? ?

其中初值a0=a，c1=0。當系統參數c1=-2，c2=-5和c3=10時，令初始角頻率ω0=0，增量Δω=0.000 1；令逼近式(7)中的L=M=3，式(12)中的精度控制指數N=6，即可利用第一類改進的廣義Padé逼近方法來逼近冪級數解式(9)。當初值a=1時，求得廣義Padé角頻率ω=1.133 8，周期T=2.770 8，逼近式為

GPA[3/3]=(0.361 4+1.160 9cos(1.133 8t)2+

0.508 3cos(1.133 8t)4-0.007 5cos(1.133 8t)6)/(1+

2.239 1cos(1.133 8t)2-1.203 0cos(1.133 8t)4-

0.012 9cos(1.133 8t)6)

(19)

其相圖如圖1所示，其時間歷程如圖2所示，解式(19)與數值解的絕對誤差曲線如圖3所示，本文所采用的數值解均為四階Runge-Kutta法所得。從圖1～圖3不難看出，所得之解有著較高的精度。

—表示數值解； …表示本文所得之解圖1 解式(19)的相圖Fig.1 The phase portrait of solution (19)

—表示數值解； …表示本文所得之解圖2 解式(19)的時間歷程曲線 Fig.2 Time response curve of solution (19)

圖3 解式(19)與數值解的絕對誤差曲線Fig.3 The absolute error curve between solution (19) and numerical solution

當初值a=-0.36時，求得廣義Padé角頻率ω=0.701 5，周期T=4.478 4，逼近式為

GPA[3/3]=(-0.112 3-1.182 7cos(0.701 5t)2-

1.003 6cos(0.701 5t)4-0.035 5cos(0.701 5t)6)/(1+

朝向主要考慮東南西北是個方向，東邊選百葉窗、南窗配深色窗簾、西窗選有褶簾、北邊選藝術窗簾。當建筑本身處于不規則形狀，使得窗戶所在位置不是以上四個朝向時，則需根據實際情況考慮。

9.201 2cos(0.701 5t)2-3.888 2cos(0.701 5t)4+

0.170 6cos(0.701 5t)6)

(20)

其相圖如圖4所示，其時間歷程如圖5所示，解式(20)與數值解的絕對誤差曲線如圖6所示。從圖4～圖6不難看出，所得之解有著很高的精度。

當初值a=1.15時，所求周期軌為大周期軌，令增量Δω=0.000 01，求得廣義Padé角頻率ω=0.845 25，周期T=3.716 8，逼近式為

GPA[3/3]=(-0.659 3+1.838 7cos(0.845 25t)2-

0.499 4cos(0.845 25t)4+0.185 3cos(0.845 25t)6)/(1+

1.200 6cos(0.845 25t)2-1.474 4cos(0.845 25t)4+

(21)

其相圖如圖7所示，其時間歷程如圖8所示，解式(21)與數值解的絕對誤差曲線如圖9所示。

—表示數值解； …表示本文所得之解圖4 解式(20)的相圖Fig.4 The phase portrait of solution (20)

—表示數值解； …表示本文所得之解圖5 解式(20)的時間歷程曲線Fig.5 Time response curve of solution (20)

圖6 解式(20)與數值解的絕對誤差曲線Fig.6 The absolute error curve between solution (20) and numerical solution

—表示數值解； …表示本文所得之解圖7 解式(21)的相圖Fig.7 The phase portrait of solution (21)

—表示數值解； …表示本文所得之解圖8 解式(21)的時間歷程曲線Fig.8 Time response curve of solution (21)

圖9 解式(21)與數值解的絕對誤差曲線Fig.9 The absolute error curve between solution (21) and numerical solution

接下來討論系統式(18)在不同參數時，該方法的可靠性。表1描述了在c1，c2和初值不變，最高階非線性項系數c3由小增大時，利用第一類改進的廣義Padé逼近方法所求得的閉軌周期與數值方法得到的周期之間的關系。表2則描述了在參數不變，初始振幅由小增大時，利用第一類改進的廣義Padé逼近方法所求得的閉軌周期與數值方法得到的周期之間的關系。其中T1為本文求得的周期，T2為數值方法得到的周期，Te為二者的誤差。從上述二表中不難看出，該方法的誤差不會隨著非線性項系數的增大而增大，精度亦不隨著初始振幅的增大而減小，因此該方法在強非線性和大振幅下均有較高的可靠性。

表1 系統式(18)在不同參數下的周期比較Tab.1 The comparisons of the periodic with different parameters of system (18)

表2 系統式(18)在不同初始振幅下的周期比較Tab.2 The comparison of periodic with different initial amplitude of system (18)

3.2 廣義Duffing-Harmonic振子的周期解

廣義Duffing-Harmonic振子可由式(22)描述

(22)

考慮μv≠0的情形。Duffing-Harmonic振子常被用來檢驗一些近似求解算法的精度和有效性。本文亦選取該振子來驗證改進的廣義Padé逼近方法之有效性和精度。首先，令式(22)的冪級數為式(9)所示，將式(9)代入式(22)即可求得式(9)的各項系數為

? ?

其中初值a0=a，a1=0。當λ=-1，μ=5和v=2，初值a=0.6時，仍可用第一類改進的廣義Padé逼近方法來求解，令初始角頻率ω0=0，增量Δω=0.000 01；令逼近式(7)中的L=M=3，式(12)中的精度控制指數N=6，求得廣義Padé角頻率ω=0.547 62，周期T=5.736 8，逼近式為

GPA[3/3]=(0.247 37+0.161 3cos(0.547 62t)2+

0.009 9cos(0.547 62t)4+0.013 0cos(0.547 62t)6)/(1-

0.380 7cos(0.547 62t)2+0.110 5cos(0.547 62t)4-

0.010 4cos(0.547 62t)6)

(23)

其相圖如圖10所示，時間歷程如圖11所示，解式(23)與數值解的絕對誤差曲線如圖12所示。從上述圖像可看出，解式(23)有著較高的精度，也進一步驗證了該方法的有效性。

—表示數值解； …表示本文所得之解圖10 解(23)的相圖Fig.10 The phase portrait of solution (23)

—表示數值解； …表示本文所得之解圖11 解(23)的時間歷程曲線Fig.11 Time response curve of solution (23)

圖12 解式(23)與數值解的絕對誤差曲線Fig.12 The absolute error curve between solution (23) and numerical solution

當初值a=0.9時，式(22)的周期軌為大周期軌，即在相圖上包圍了所有奇點的軌道。此時，利用第一類改進的廣義Padé逼近方法所求得的結果精度不太理想，由于式(22)為對稱振子，因此可采用第二類改進的廣義Padé逼近方法來求解，令初始角頻率ω0=0，增量Δω=0.000 001；令逼近式(7)中的L=M=4，式(12)中的精度控制指數N=6，可求得廣義Padé角頻率ω=0.441 812，周期T=7.110 7，逼近式為

35.690 9cos(0.441 812t)4+41.584 9cos(0.441 812t)6-

22.422 4cos(0.441 812t)8)/(1-12.535 7cos(0.441 812t)2+

15.439 6cos(0.441 812t)4-12.019 9cos(0.441 812t)6+

2.142 1cos(0.441 812t)8)

(24)

根據式(17)，將系統式(22)的解描述為如下分段函數

(25)

式中：ω為所求得的廣義Padé角頻率。其相圖如圖13所示，時間歷程如圖14所示，解式(25)與數值解的絕對誤差曲線如圖15所示。由上述圖像可知，第二類改進的廣義Padé逼近方法仍然有著較高的精度和可靠性。為進一步說明，我們繪制了相應的表格。表3描述了在λ，v和初值不變，最高階非線性項系數μ由小增大時，利用第二類改進的廣義Padé逼近方法所求得的閉軌周期與數值方法得到的周期之間的關系。表4則描述了在參數不變，初始振幅由小增大時，利用第二類改進的廣義Padé逼近方法所求得的閉軌周期與數值方法得到的周期之間的關系。其中T1為本文求得的周期，T2為數值方法得到的周期，Te為二者的誤差。

表3 系統式(22)在不同參數下的周期比較Tab.3 The comparisons of the periodic with different parameters of system (22)

—表示數值解； …表示本文所得之解圖13 解式(25)的相圖Fig.13 The phase portrait of solution (25)

—表示數值解； …表示本文所得之解圖14 解式(25)的時間歷程曲線Fig.14 Time response curve of solution (25)

圖15 解式與(25)數值解的絕對誤差曲線Fig.15 The absolute error curve between solution (25) and numerical solution

表4 系統式(22)在不同初始振幅下的周期比較Tab.4 The comparison of periodic with different initial amplitudeof system (22)

從表3和表4中可見，在強非線性和大振幅下，本文所求得的廣義Duffing-Harmonic振子的周期均保持著較高的精度，進一步表明了該方法的有效性和可靠性。

3.3 勢能函數為無理函數的振子

本小節考慮如下形式無理振子

(26)

? ?

其中初值a0=a，a1=0。當系統參數c1=-2，c2=1和c3=4時，系統的大致相圖如圖16所示。令初始角頻率ω0=0，增量Δω=0.000 01。令式中的L=M=3，式(12)中的精度控制指數N=6，即可利用第一類改進的廣義Padé逼近方法來逼近冪級數解式(9)。當初值a=0.5時，求得廣義Padé角頻率ω=0.419 08，周期T=7.496 4，逼近式為

GPA[3/3]=(0.960 3-1.058 0cos(0.419 05t)2+

0.304 5cos(0.419 05t)4-0.018 0cos(0.419 05t)6)/

(1-0.531 4cos(0.419 05t)2-0.086 1cos(0.419 05t)4-

0.004 6cos(0.419 05t)6)

(27)

其相圖如圖17所示，其時間歷程如圖18所示，解式(27)與數值解的絕對誤差曲線如圖19所示。

圖16 系統式(26)在c1=-2，c2=1和c3=4時的大致相圖Fig.16 The rough phase portrait of system (26) when c1=-2,c2=1,c3=4

從圖17～圖19不難看出，所得之解有著很高的精度。下面討論系統在不同初值時，利用第一類改進的廣義Padé逼近方法所求得的閉軌周期與數值方法得到的周期之間的關系。如表5所示，其中T1為本文求得的周期，T2為數值方法得到的周期，Te為二者的誤差。從上表中不難看出，該方法對于不同初值的周期軌均有著較好的精度。從本小節的內容可以看出，本節提出的改進的廣義Padé逼近方法在無理振子的周期解求解中依然有著較高的精度，也進一步說明了該方法有著較廣的適用范圍。

—表示數值解； …表示本文所得之解圖17 解式(27)的相圖Fig.17 The phase portrait of solution (27)

—表示數值解； …表示本文所得之解圖18 解式(27)的時間歷程曲線Fig.18 Time response curve of solution (27)

圖19 解式(27)與數值解的絕對誤差曲線.Fig.19 The absolute error curve between solution (27) and numerical solution

表5 系統式(26)在不同初始振幅下的周期比較Tab.5 The comparison of periodic with different initial amplitude of system (26)

4 結 論

本文基于廣義Padé逼近方法，構造了一類由余弦函數級數組成的新的廣義Padé逼近式；同時，對經典Padé逼近的求解過程進行了合理的改進，并求得了一類勢能函數為高階多項式、有理函數和無理函數振子的高精度解析周期解及其周期，為非線性振動的定量研究提供了新的思路和參考方法，也擴大了Padé逼近方法的適用范圍，因此對該方法的研究具有一定的理論價值和實際意義?？偟膩碚f該方法具有以下特點：①繼承了Padé逼近方法求解過程簡單直接的優點，便于進行計算機編程計算；②在強非線性和大振幅下，所得結果均保持較高精度，保證了該方法的可靠性；③該方法并不局限于某一個固定的系統，而是有著較廣的適用范圍。將該方法進行進一步的改進，亦可用于求解阻尼振動和受迫振動的情形，這也將是我們下一步的工作。