非线性保守系统周期运动的Hermite插值解法*

朱金文 杨德庆

（上海交通大学船舶海洋与建筑工程学院，上海 200240）

引言

很早就有学者［1］认识到动力系统微分方程的解可展开为时间的幂级数形式.该级数收敛半径通常很小，因此不能构成系统有效的解.通过解析延拓可以扩大解的收敛区域，然而这既不实用，也无法给出运动系统的一般性质.

研究人员一直在寻求整个时间段内收敛的解的表达式.1884年，Poincaré［1］认识到转换时间变量的必要性，提出了以下时间转换公式：

其中τ为新的时间变量，t0和h为常数.由于无穷的时间被转换为有限的时间，系统运动幂级数解的收敛要求降低到｜τ｜＝1的圆内.然而，一段时间后τ将趋于常数，运动随之停止.1996年，Qaisi［2］意识到将独立时间变量t转换为谐振时间变量τ更为合理：

无穷的时间区域0≤t≤∞变为有限的时间区域-1≤τ≤+1，且τ以频率ω做简谐振动.基于该变换，Qaisi提出了研究周期运动的幂级数法［2］.该方法将运动系统对谐振时间τ进行零点幂级数展开，由于τ以频率ω在［-1，1］之间做简谐振动，故只要级数在｜τ｜＝1的圆内收敛，则原系统可解.然而，由于幂级数的固有特性，往往无法保证［-1，1］内的收敛性；即使收敛，其收敛速度也比较慢.因此，该方法在应用中受到了限制.

本文在Qaisi谐振时间变换的基础上提出周期运动的Hermite插值解法，与幂级数法不同，采用两点Hermite插值函数代替一点幂级数展开，解决了收敛性问题并提高收敛速度.

1 Hermite插值解法

为了阐述，考虑三次Duffing振子的自由振动

其初始条件为 x（0）＝A与 ˙x（0）＝0，点表示对时间t求导，x表示系统位移，ε为系统非线性参数.吴晓和黄翀［3］在研究功能梯度材料椭圆板的非线性热振动及屈曲中用摄动法对（3）进行了求解.

为了用Hermite插值法求解方程（3），首先将独立时间变量t转换为谐振时间变量τ

对应的微分方程（3）变为

其中撇表示对时间τ求导.

由于系统的初始状态点由τ＝0唯一确定，振子在“振动时间τ”的半个振动周期内必定已经完成了一次振动，因此振动时间τ的频率ω必定是振子振动频率Ω（A）的一半：

另外，由方程（4）容易得出x（τ）是偶函数.

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

另一方面，易知该振子做等幅值振动，故x（1）＝-A；令（4）中 τ＝1，可以得到

事实上，对方程（4）求导，可依次得到0或1点处的高阶导数.

由上述条件可以得到τ∈［0，1］中的任意Hermite插值函数.可以断言，只要方程（4）的解在τ∈［0，1］是光滑的，就可以得到这个区间上足够精确的逼近解.又因x（τ）为偶函数，故可方便的得到x（τ）在 τ∈［-1，1］上的近似解.只要将 τ替换为sinωt就得到了原方程的解.

例如，为了得到 x（τ）在 τ∈［-1，1］上的近似解，可以构造插值函数如下

其中ci（i＝1，2）为待定常数.进一步，为了确定常数 ci（i＝1，2），令 X（1）＝x（1）和 X′（1）＝x′（1），得到两个线性方程，解得 ci（i＝1，2）如下

显然式（8）使用了0点直到2阶的导数和1点直到1阶的导数，故（8）可称为［2，1］阶 Hermite插值函数.

将（8）中的τ替换为sinωt并进行三角函数化简，可得方程（3）的近似解为

可见方程（3）的［2，1］阶 Hermite插值解，即是在谐波解的基础上加上了一个修正项.

另一方面，振子（3）的频率可以通过能量积分得到.对应（3）的 Hamilton量［4］为

其中为动能，为势能.由初始条件

易知

振子的周期为

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

则频率为

事实上，对于非线性振子，并不总能得到其频率解析表达式，然而经常能够得到其周期的积分表达式，如式（13）所示，由数值积分即可求得其自然频率.

2 一类非线性振子的近似通解

注意到由［2，1］阶 Hermite插值得到的方程（3）的近似解由一次谐波项及另外一个修正项组成，如果将方程（3）改写为

其中那么（10）可改写为

可以推断，对于一般的非线性振子

其中 f（x）为偶函数，初始条件为 x（0）＝x0，˙x（0）＝˙x0，其近似通解为

这个结论很容易证明：首先不考虑相位θ，则（19）转化为（17）；为了得到（17），只需用 Hermite插值法求解方程（18）的［2，1］阶解，其过程与求解（3）的过程完全类似.

另一方面，考虑修正项中所起的作用.在一个周期内cosΩt-cos3Ωt随时间的变化如图1所示

由图1容易看出当修正项系数为正时，由于修正项的“削波”作用，振子的振动波形将趋于平坦；反之，若修正项系数为负，振子的振动波形将趋于尖锐.图2表示了修正项对振子影响.

图1 cosΩt-cos3Ωt在一个周期内随时间的变化Fig.1 The variation of cosΩt-cos3Ωt with time in a period

为了方便，可将f（A）称为“修正频率”.当修正频率小于自然频率时，振子随时间变化的波形趋于平坦，表现出希望在峰值处滞留的特性；当修正频率大于自然频率时，波形趋于尖锐，表现出疏离峰值的特性；而当修正频率恰好等于自然频率时，振子的振动为简谐振动，振子为线性振子.

图 2 cosΩt+1／8（cosΩt-cos3Ωt）与 cosΩt-1／8（cosΩt-cos3Ωt）在一个周期内随时间的变化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

例如，单摆在做大振幅摆动时总是趋于在最大振幅处滞留；而振子弹簧为硬弹簧时，总是趋于疏离最大振幅的位置.

当然，由Hermite插值解法可以求解更高阶的解，只需要利用0和1处的高阶导数信息.

3 高阶Hermite插值

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

其中 f（A）与 m如（14）所示.因此可以由（20）来检验Hermite插值解的精度.

取非线性参数ε＝10，幅值A＝1，求解Duffing方程（3）的［6，4］阶 Hermite插值解.当 A＝1时，可以得到 Ω（1）＝2.86664，因此 ω＝1.43332.则可以求得τ＝0处x（τ）的幂级数为

图3 非线性参数取 ε＝10时，振幅 A＝0.1，1，10的［6，4］阶Hermite插值解与精确解的对比Fig.3 The comparison between［6，4］th Hermite interpolation solution and exact solution when nonlinear parameterε＝10

为了得到（3）的［6，4］阶 Hermite插值解，构造插值函数

其中 ci（i＝1，2，3，4）为待定常数.

令 X（1）＝x（1），X′（1）＝x′（1），X＂（1）＝x＂（1）和 X（3）（1）＝x（3）（1），得到 4个线性方程组.求解该方程组，得到 ci（i＝1，2，3，4）

将（22）中替换为sinωt并进行三角函数的化简，得到（3）的近似解

为了检验解（24）的精度，将精确解（20）进行Fourier展开

可以看到对应的谐波系数彼此十分接近.

为了进一步考察Hermite插值解的收敛性，可以进一步求解［8，4］阶、［10，5］阶和［12，6］阶解，取其前4个谐波的谐波系数，列表如下

表 1 ［8，4］阶、［10，5］阶和［12，6］阶Hermite插值解的前4个谐波系数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插值解与精确解的对比.

由图3可知，Hermite插值方法在求解非线性保守系统周期运动时简单且高效.

4 讨论

第3节给出了Duffing方程的高阶Hermite插值解，而事实上对很多非线性振子，［2，1］阶Hermite插值解的精度已经足够.例如：广义Duffing简谐振子［6］

若不考虑相位θ，由式（17）可直接写出（26）的解为

其频率也容易由能量积分得到

同样的，对于Duffing简谐振子［7-10］

若不考虑相位θ，由式（17）可直接写出（29）的解为

容易得到其频率为

可以验证Duffing简谐振子在幅值分别取0.01，0.1，1，10时解（30）的精确性，如图 4所示，其中由数值积分所得频率依次为 0.00847179，0.0843886，0.63678，0.990916.

5 小结

本文提出了非线性保守系统周期运动的Hermite插值解法，由［2，1］阶 Hermite插值给出了一类非线性振子的近似通解

分析了保守非线性振子的一些振动特性，介绍了振子的高阶Hermite插值.

图4 Duffing简谐振子在幅值分别为0.01，0.1，1，10时解（30）与数值解的对比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插值解法具有很好的收敛性，在求解高阶解时，相当于Fourier展开，一般而言［6，4］阶近似解可满足精度要求.对很多非线性振子［2，1］阶近似解就已经具有很高的精度.故Hermite插值法具有很好的收敛性和较高的收敛速度.

需要说明的是，对于非等振幅的振子，同样可以使用该方法进行求解.只需先由能量积分公式确定出两个幅值，再通过数值积分得到振子频率.因此，本文提出的Hermite插值解法对于求解非线性保守系统周期运动是一种简单有效的数值解析方法.

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 吴晓，黄翀.功能梯度材料椭圆板的非线性热振动及屈曲.动力学与控制学报，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