非線性保守系統(tǒng)周期運(yùn)動(dòng)的Hermite插值解法*

朱金文 楊德慶

（上海交通大學(xué)船舶海洋與建筑工程學(xué)院，上海 200240）

引言

很早就有學(xué)者［1］認(rèn)識(shí)到動(dòng)力系統(tǒng)微分方程的解可展開為時(shí)間的冪級(jí)數(shù)形式.該級(jí)數(shù)收斂半徑通常很小，因此不能構(gòu)成系統(tǒng)有效的解.通過解析延拓可以擴(kuò)大解的收斂區(qū)域，然而這既不實(shí)用，也無法給出運(yùn)動(dòng)系統(tǒng)的一般性質(zhì).

研究人員一直在尋求整個(gè)時(shí)間段內(nèi)收斂的解的表達(dá)式.1884年，Poincaré［1］認(rèn)識(shí)到轉(zhuǎn)換時(shí)間變量的必要性，提出了以下時(shí)間轉(zhuǎn)換公式：

其中τ為新的時(shí)間變量，t0和h為常數(shù).由于無窮的時(shí)間被轉(zhuǎn)換為有限的時(shí)間，系統(tǒng)運(yùn)動(dòng)冪級(jí)數(shù)解的收斂要求降低到｜τ｜＝1的圓內(nèi).然而，一段時(shí)間后τ將趨于常數(shù)，運(yùn)動(dòng)隨之停止.1996年，Qaisi［2］意識(shí)到將獨(dú)立時(shí)間變量t轉(zhuǎn)換為諧振時(shí)間變量τ更為合理：

無窮的時(shí)間區(qū)域0≤t≤∞變?yōu)橛邢薜臅r(shí)間區(qū)域-1≤τ≤+1，且τ以頻率ω做簡(jiǎn)諧振動(dòng).基于該變換，Qaisi提出了研究周期運(yùn)動(dòng)的冪級(jí)數(shù)法［2］.該方法將運(yùn)動(dòng)系統(tǒng)對(duì)諧振時(shí)間τ進(jìn)行零點(diǎn)冪級(jí)數(shù)展開，由于τ以頻率ω在［-1，1］之間做簡(jiǎn)諧振動(dòng)，故只要級(jí)數(shù)在｜τ｜＝1的圓內(nèi)收斂，則原系統(tǒng)可解.然而，由于冪級(jí)數(shù)的固有特性，往往無法保證［-1，1］內(nèi)的收斂性；即使收斂，其收斂速度也比較慢.因此，該方法在應(yīng)用中受到了限制.

本文在Qaisi諧振時(shí)間變換的基礎(chǔ)上提出周期運(yùn)動(dòng)的Hermite插值解法，與冪級(jí)數(shù)法不同，采用兩點(diǎn)Hermite插值函數(shù)代替一點(diǎn)冪級(jí)數(shù)展開，解決了收斂性問題并提高收斂速度.

1 Hermite插值解法

為了闡述，考慮三次Duffing振子的自由振動(dòng)

其初始條件為 x（0）＝A與 ˙x（0）＝0，點(diǎn)表示對(duì)時(shí)間t求導(dǎo)，x表示系統(tǒng)位移，ε為系統(tǒng)非線性參數(shù).吳曉和黃翀［3］在研究功能梯度材料橢圓板的非線性熱振動(dòng)及屈曲中用攝動(dòng)法對(duì)（3）進(jìn)行了求解.

為了用Hermite插值法求解方程（3），首先將獨(dú)立時(shí)間變量t轉(zhuǎn)換為諧振時(shí)間變量τ

對(duì)應(yīng)的微分方程（3）變?yōu)?/p>

其中撇表示對(duì)時(shí)間τ求導(dǎo).

由于系統(tǒng)的初始狀態(tài)點(diǎn)由τ＝0唯一確定，振子在“振動(dòng)時(shí)間τ”的半個(gè)振動(dòng)周期內(nèi)必定已經(jīng)完成了一次振動(dòng)，因此振動(dòng)時(shí)間τ的頻率ω必定是振子振動(dòng)頻率Ω（A）的一半：

另外，由方程（4）容易得出x（τ）是偶函數(shù).

令（4）中，由初始條件可以得到

另一方面，易知該振子做等幅值振動(dòng)，故x（1）＝-A；令（4）中 τ＝1，可以得到

事實(shí)上，對(duì)方程（4）求導(dǎo)，可依次得到0或1點(diǎn)處的高階導(dǎo)數(shù).

由上述條件可以得到τ∈［0，1］中的任意Hermite插值函數(shù).可以斷言，只要方程（4）的解在τ∈［0，1］是光滑的，就可以得到這個(gè)區(qū)間上足夠精確的逼近解.又因x（τ）為偶函數(shù)，故可方便的得到x（τ）在 τ∈［-1，1］上的近似解.只要將 τ替換為sinωt就得到了原方程的解.

例如，為了得到 x（τ）在 τ∈［-1，1］上的近似解，可以構(gòu)造插值函數(shù)如下

其中ci（i＝1，2）為待定常數(shù).進(jìn)一步，為了確定常數(shù) ci（i＝1，2），令 X（1）＝x（1）和 X′（1）＝x′（1），得到兩個(gè)線性方程，解得 ci（i＝1，2）如下

顯然式（8）使用了0點(diǎn)直到2階的導(dǎo)數(shù)和1點(diǎn)直到1階的導(dǎo)數(shù)，故（8）可稱為［2，1］階 Hermite插值函數(shù).

將（8）中的τ替換為sinωt并進(jìn)行三角函數(shù)化簡(jiǎn)，可得方程（3）的近似解為

可見方程（3）的［2，1］階 Hermite插值解，即是在諧波解的基礎(chǔ)上加上了一個(gè)修正項(xiàng).

另一方面，振子（3）的頻率可以通過能量積分得到.對(duì)應(yīng)（3）的 Hamilton量［4］為

其中為動(dòng)能，為勢(shì)能.由初始條件

易知

振子的周期為

令 x＝A cosθ可將（13）化為 T＝4K（m）／f（A），其中

則頻率為

事實(shí)上，對(duì)于非線性振子，并不總能得到其頻率解析表達(dá)式，然而經(jīng)常能夠得到其周期的積分表達(dá)式，如式（13）所示，由數(shù)值積分即可求得其自然頻率.

2 一類非線性振子的近似通解

注意到由［2，1］階 Hermite插值得到的方程（3）的近似解由一次諧波項(xiàng)及另外一個(gè)修正項(xiàng)組成，如果將方程（3）改寫為

其中那么（10）可改寫為

可以推斷，對(duì)于一般的非線性振子

其中 f（x）為偶函數(shù)，初始條件為 x（0）＝x0，˙x（0）＝˙x0，其近似通解為

這個(gè)結(jié)論很容易證明：首先不考慮相位θ，則（19）轉(zhuǎn)化為（17）；為了得到（17），只需用 Hermite插值法求解方程（18）的［2，1］階解，其過程與求解（3）的過程完全類似.

另一方面，考慮修正項(xiàng)中所起的作用.在一個(gè)周期內(nèi)cosΩt-cos3Ωt隨時(shí)間的變化如圖1所示

由圖1容易看出當(dāng)修正項(xiàng)系數(shù)為正時(shí)，由于修正項(xiàng)的“削波”作用，振子的振動(dòng)波形將趨于平坦；反之，若修正項(xiàng)系數(shù)為負(fù)，振子的振動(dòng)波形將趨于尖銳.圖2表示了修正項(xiàng)對(duì)振子影響.

圖1 cosΩt-cos3Ωt在一個(gè)周期內(nèi)隨時(shí)間的變化Fig.1 The variation of cosΩt-cos3Ωt with time in a period

為了方便，可將f（A）稱為“修正頻率”.當(dāng)修正頻率小于自然頻率時(shí)，振子隨時(shí)間變化的波形趨于平坦，表現(xiàn)出希望在峰值處滯留的特性；當(dāng)修正頻率大于自然頻率時(shí)，波形趨于尖銳，表現(xiàn)出疏離峰值的特性；而當(dāng)修正頻率恰好等于自然頻率時(shí)，振子的振動(dòng)為簡(jiǎn)諧振動(dòng)，振子為線性振子.

圖 2 cosΩt+1／8（cosΩt-cos3Ωt）與 cosΩt-1／8（cosΩt-cos3Ωt）在一個(gè)周期內(nèi)隨時(shí)間的變化Fig.2 The variation of cosΩt+1／8（cosΩt-cos3Ωt）and cosΩt-1／8（cosΩt-cos3Ωt）with time in a period

例如，單擺在做大振幅擺動(dòng)時(shí)總是趨于在最大振幅處滯留；而振子彈簧為硬彈簧時(shí)，總是趨于疏離最大振幅的位置.

當(dāng)然，由Hermite插值解法可以求解更高階的解，只需要利用0和1處的高階導(dǎo)數(shù)信息.

3 高階Hermite插值

由于方程（3）有精確解［5］

其中 f（A）與 m如（14）所示.因此可以由（20）來檢驗(yàn)Hermite插值解的精度.

取非線性參數(shù)ε＝10，幅值A(chǔ)＝1，求解Duffing方程（3）的［6，4］階 Hermite插值解.當(dāng) A＝1時(shí)，可以得到 Ω（1）＝2.86664，因此 ω＝1.43332.則可以求得τ＝0處x（τ）的冪級(jí)數(shù)為

圖3 非線性參數(shù)取 ε＝10時(shí)，振幅 A＝0.1，1，10的［6，4］階Hermite插值解與精確解的對(duì)比Fig.3 The comparison between［6，4］th Hermite interpolation solution and exact solution when nonlinear parameterε＝10

為了得到（3）的［6，4］階 Hermite插值解，構(gòu)造插值函數(shù)

其中 ci（i＝1，2，3，4）為待定常數(shù).

令 X（1）＝x（1），X′（1）＝x′（1），X＂（1）＝x＂（1）和 X（3）（1）＝x（3）（1），得到 4個(gè)線性方程組.求解該方程組，得到 ci（i＝1，2，3，4）

將（22）中替換為sinωt并進(jìn)行三角函數(shù)的化簡(jiǎn)，得到（3）的近似解

為了檢驗(yàn)解（24）的精度，將精確解（20）進(jìn)行Fourier展開

可以看到對(duì)應(yīng)的諧波系數(shù)彼此十分接近.

為了進(jìn)一步考察Hermite插值解的收斂性，可以進(jìn)一步求解［8，4］階、［10，5］階和［12，6］階解，取其前4個(gè)諧波的諧波系數(shù)，列表如下

表 1 ［8，4］階、［10，5］階和［12，6］階Hermite插值解的前4個(gè)諧波系數(shù)Table 1 The first four harmonic coefficients of［8，4］th，［10，5］th，and［12，6］th solutions

可以看出Hermite插值方法具有很好的收斂性.

圖3分別給出了A＝0.1，1，10的Hermite插值解與精確解的對(duì)比.

由圖3可知，Hermite插值方法在求解非線性保守系統(tǒng)周期運(yùn)動(dòng)時(shí)簡(jiǎn)單且高效.

4 討論

第3節(jié)給出了Duffing方程的高階Hermite插值解，而事實(shí)上對(duì)很多非線性振子，［2，1］階Hermite插值解的精度已經(jīng)足夠.例如：廣義Duffing簡(jiǎn)諧振子［6］

若不考慮相位θ，由式（17）可直接寫出（26）的解為

其頻率也容易由能量積分得到

同樣的，對(duì)于Duffing簡(jiǎn)諧振子［7-10］

若不考慮相位θ，由式（17）可直接寫出（29）的解為

容易得到其頻率為

可以驗(yàn)證Duffing簡(jiǎn)諧振子在幅值分別取0.01，0.1，1，10時(shí)解（30）的精確性，如圖 4所示，其中由數(shù)值積分所得頻率依次為 0.00847179，0.0843886，0.63678，0.990916.

5 小結(jié)

本文提出了非線性保守系統(tǒng)周期運(yùn)動(dòng)的Hermite插值解法，由［2，1］階 Hermite插值給出了一類非線性振子的近似通解

分析了保守非線性振子的一些振動(dòng)特性，介紹了振子的高階Hermite插值.

圖4 Duffing簡(jiǎn)諧振子在幅值分別為0.01，0.1，1，10時(shí)解（30）與數(shù)值解的對(duì)比Fig.4 The comparison between the［2，1］th Hermite interpolation solution and numerical solution for the Duffing harmonic oscillator when the amplitudes are 0.01，0.1，1，10

研究表明，Hermite插值解法具有很好的收斂性，在求解高階解時(shí)，相當(dāng)于Fourier展開，一般而言［6，4］階近似解可滿足精度要求.對(duì)很多非線性振子［2，1］階近似解就已經(jīng)具有很高的精度.故Hermite插值法具有很好的收斂性和較高的收斂速度.

需要說明的是，對(duì)于非等振幅的振子，同樣可以使用該方法進(jìn)行求解.只需先由能量積分公式確定出兩個(gè)幅值，再通過數(shù)值積分得到振子頻率.因此，本文提出的Hermite插值解法對(duì)于求解非線性保守系統(tǒng)周期運(yùn)動(dòng)是一種簡(jiǎn)單有效的數(shù)值解析方法.

1 Whittaker E T.A treatise on the analytical dynamics of particles and rigid bodies.Cambridge：Cambridge University Press，England，1965

2 QaisiM I.A power series approach for the study of periodic motion.Journal of Sound and Vibration，1996，196（4）：247～252

3 吳曉，黃翀.功能梯度材料橢圓板的非線性熱振動(dòng)及屈曲.動(dòng)力學(xué)與控制學(xué)報(bào)，2013，11（2）：165～171（Wu X，Huang C.Nonlinear thermal vibration and buckling of functionally graded elliptical plate.Journal of Dynamics and Control，2013，11（2）：165～171（in Chinese））

4 Akbarzade M，Kargar A.Application of the Hamiltonian approach to nonlinear vibrating equations.Mathematical and Computer Modelling，2011，54（9-10）：2504～2514

5 Mickens R E.Oscillations in planar dynamic systems.Singapore：World Scientific，1996

6 Wang H L，Chung K Y.Analytical solutions of a generalized Duffing-harmonic oscillator by a nonlinear time transformationmethod.Physics Letters A，2012，376（12-13）：1118～1124

7 Mickens R E.Mathematical and numerical study of the Duffing-Harmonic oscillator.Journal of Sound and Vibration，2001，244（3）：563～567

8 Lim CW，Wu B S.A new analytical approach to the Duffing-harmonic oscillator.Physics Letters A，2003，311（4-5）：365～373

9 Lim CW，Wu B S，Sun W P.Higher accuracy analytical approximations to the Duffing-harmonic oscillator.Journal of Sound and Vibration，2006，296（4-5）：1039～1045

10 Fesanghary M，Pirbodaghi T，Asghari M，et al.A new analytical approximation to the Duffing-harmonic oscillator.Chaos，Solitons＆Fractals，2009，42（1）：571～576