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

李震波， 唐駕時

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

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

(1)

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

(2)

(3)

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

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

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

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

(4)

(5)

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

若存在有理分式函數(shù)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]。

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

(7)

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

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

考慮如下形式的自治振動系統(tǒng)

(8)

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

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

(9)

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

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

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

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

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

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

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

(10)

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

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

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

(11)

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

Z=α0-b0

(12)

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

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

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

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

α0=b0

(13)

其中b0可由式(14)確定

(14)

步驟3令如下函數(shù)為目標(biāo)函數(shù)

(15)

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

(16)

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

(17)

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

3 算 例

3.1 Helmholtz-Duffing振子的周期解

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

(18)

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

? ?

其中初值a0=a，c1=0。當(dāng)系統(tǒng)參數(shù)c1=-2，c2=-5和c3=10時，令初始角頻率ω0=0，增量Δω=0.000 1；令逼近式(7)中的L=M=3，式(12)中的精度控制指數(shù)N=6，即可利用第一類改進的廣義Padé逼近方法來逼近冪級數(shù)解式(9)。當(dāng)初值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)與數(shù)值解的絕對誤差曲線如圖3所示，本文所采用的數(shù)值解均為四階Runge-Kutta法所得。從圖1～圖3不難看出，所得之解有著較高的精度。

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

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

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

當(dāng)初值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+

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

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)與數(shù)值解的絕對誤差曲線如圖6所示。從圖4～圖6不難看出，所得之解有著很高的精度。

當(dāng)初值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)與數(shù)值解的絕對誤差曲線如圖9所示。

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

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

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

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

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

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

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

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

表2 系統(tǒng)式(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)的冪級數(shù)為式(9)所示，將式(9)代入式(22)即可求得式(9)的各項系數(shù)為

? ?

其中初值a0=a，a1=0。當(dāng)λ=-1，μ=5和v=2，初值a=0.6時，仍可用第一類改進的廣義Padé逼近方法來求解，令初始角頻率ω0=0，增量Δω=0.000 01；令逼近式(7)中的L=M=3，式(12)中的精度控制指數(shù)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)與數(shù)值解的絕對誤差曲線如圖12所示。從上述圖像可看出，解式(23)有著較高的精度，也進一步驗證了該方法的有效性。

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

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

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

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

根據(jù)式(17)，將系統(tǒng)式(22)的解描述為如下分段函數(shù)

(25)

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

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

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

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

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

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

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

3.3 勢能函數(shù)為無理函數(shù)的振子

本小節(jié)考慮如下形式無理振子

(26)

? ?

其中初值a0=a，a1=0。當(dāng)系統(tǒng)參數(shù)c1=-2，c2=1和c3=4時，系統(tǒng)的大致相圖如圖16所示。令初始角頻率ω0=0，增量Δω=0.000 01。令式中的L=M=3，式(12)中的精度控制指數(shù)N=6，即可利用第一類改進的廣義Padé逼近方法來逼近冪級數(shù)解式(9)。當(dāng)初值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)與數(shù)值解的絕對誤差曲線如圖19所示。

圖16 系統(tǒng)式(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不難看出，所得之解有著很高的精度。下面討論系統(tǒng)在不同初值時，利用第一類改進的廣義Padé逼近方法所求得的閉軌周期與數(shù)值方法得到的周期之間的關(guān)系。如表5所示，其中T1為本文求得的周期，T2為數(shù)值方法得到的周期，Te為二者的誤差。從上表中不難看出，該方法對于不同初值的周期軌均有著較好的精度。從本小節(jié)的內(nèi)容可以看出，本節(jié)提出的改進的廣義Padé逼近方法在無理振子的周期解求解中依然有著較高的精度，也進一步說明了該方法有著較廣的適用范圍。

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

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

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

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

4 結(jié) 論

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