姜劍，王兆清，莊美玲

（山東建筑大學力學研究所，山東 濟南250101）

0 引言

過去幾十年，多自由度系統(tǒng)的振動得到了廣泛的研究。Moochhala和 Raynor采用近似方法（Approximatemethod）分析了不規(guī)則物體的多自由度振動［1］，Huang研究了雙自由度非線性系統(tǒng)的諧波振蕩［2］，Gilchrist分析了雙自由度保守擬線性系統(tǒng)的自由振蕩［3］。雙自由度系統(tǒng)的振動在物理工程和許多實際工程中都有著非常重要的運用。例如雙彈簧支撐的彈性梁、銑床的振動都可以采用雙自由度系統(tǒng)振動模型分析研究［4］。雙自由度非線性振動系統(tǒng)可歸結(jié)為兩個非線性微分方程，很多情況下，求解非線性方程組的準確解析解是極其困難甚至不可能。因此，求解非線性方程大體上有兩類辦法，近似解析解法和數(shù)值方法。近似解析解法運用最廣泛的是攝動方法，該方法不適用于求解強非線性方程并且運用時還有其他很多弱點。許多文獻研究了攝動法的改進方法。比如同倫攝動法、max-min approach method［5－6］。Ladygina等 采用多 尺度法（multiscalemethod）分析了雙自由度保守系統(tǒng)的非線性振動和它的自振頻率［7］。Cveticanin同時采用雅可比橢圓函數(shù)和三角函數(shù)得到了一個雙自由度非線性系統(tǒng)振動的解析解［8－9］。數(shù)值方法中，Chen采用廣義伽遼金方法分析了雙自由度非線性振動問題［10］。Al-Karak運用拉普拉斯分解法求解了非線性微分方程問題［11］。

配點法求解線性微分方程初值問題，可以一次性得到計算節(jié)點處的函數(shù)值，是一種高效率的計算方法［12］。Rashidinia采用配點法求解了非線性問題［13］。依據(jù)未知函數(shù)近似方法的不同，配點法主要有基于特征多項式插值的擬譜方法、基于Lagrange插值的微分求積法、基于重心Lagrange插值的重心插值配點法和基于有理Haar函數(shù)插值的配點法等［14－17］。擬譜方法是一種高精度數(shù)值方法，其具有指數(shù)收斂精度［14］。基于重心Lagrange插值的配點法在求解微分方程邊值問題和初值問題時，具有很高的計算精度［16，18］。

直接采用配點法離散非線性微分方程，得到非線性代數(shù)方程，通常采用Newton法求解非線性代數(shù)方程。本文先假設未知函數(shù)的初始值，將非線性微分方程組線性化，運用重心有理插值配點法求解線性化的微分方程初值問題，再迭代逼近求解非線性微分方程的數(shù)值解。將通過算例詳細介紹本文求解多自由度非線性振動方法的高效和高精度性。

1 重心有理插值配點法

1.1 重心有理插值及其微分矩陣

對于給定的插值節(jié)點0＝t1＜t2＜… ＜tn＝T和相應的函數(shù)值 fj＝f（tj），j＝1，2，…，n，選擇任意的整數(shù)d滿足0≤d≤n，重心有理插值可由式（1）表示為

函數(shù)f（t）的m階導數(shù)可由式（3）表示為

函數(shù)f（t）在節(jié)點t1，t2，…，tn處的m階導數(shù)可由式（4）表示為

用矩陣的形式（5）表示為

1.2 多自由度系統(tǒng)線性振動控制方程的離散公式

雙自由度耦合系統(tǒng)線性振動控制方程一般由式（6）表示為

式中：x、y分別為線性系統(tǒng)振動的物理量；t為時間。

將時間區(qū)域［0，T］離散為0＝t1＜t2＜… ＜tn＝T，運用重心有理插值近似未知函數(shù)可由式（7）得

將式代入方程式，利用微分矩陣，并令方程在離散時間點 t1，t2，…，tn上成立，方程式用矩陣的形式（8）表示為

式中：D（2）為重心有理插值在時間節(jié)點 t1，t2，…，tn上的二階微分矩陣；I為 n階單位矩陣；k1、k2、k3、k4為常 數(shù)；x、y、f1、f2分 別 為 函 數(shù) x（t）、y（t）、f1（t）、f2（t）在時間節(jié)點0＝t1，t2，…，tn＝T上函數(shù)值構(gòu)成的列向量。

初始條件x（0）＝A、˙x（0）＝B、y（0）＝M，˙y（0）＝N，其中，A、B、M、N為常數(shù)，采用重心有理插值配點法將初始條件離散由式（9）表示為

定義兩個索引符號，其中 eni（i＝1，2，…，n）表示 n階單位矩陣 I的第 i行向量；dni（i＝1，2，…，n）表示n階微分矩陣D（n）第i行向量，初始條件可以寫成矩陣的形式由式（10）表示為

采用附加法施加初始條件式到式中，可求解代數(shù)方程式（11），得到線性振動控制方程數(shù)值解

2 多自由度非線性系統(tǒng)的線性化迭代法

由于非線性公式千差萬別，文章以具體的算例來說明多自由度非線性耦合系統(tǒng)的線性化迭代法的運用過程。采用控制方程在離散節(jié)點上的數(shù)值解同解析解差的最大值Ea‖uc－ue‖∞＝maxi｜uc－ue｜說明方法的計算精度，uc是數(shù)值解，ue是解析解。

算例1：如圖1所示，k1為彈簧的剛度，k2為非線性彈簧的剛度，m為物體質(zhì)量。該系統(tǒng)控制方程式（12）為［6］

給定初始條件由式（13）表示為

圖1 由非線性彈簧連接的兩個小車圖

方程中含有立方非線性項，假定未知函數(shù)的一組初始解，，代入方程組中的非線性項，得到方程的線性化方程式（14）表示為

在給定的時間節(jié)點上 t＝［t1，t2，…，tn］T，采用重心有理插值配點法離散方程，可以得到方程的迭代式（15）為

采用附加法施加初始條件式（13）到式（15），求解方程得到修正值xk，yk。給定控制精度ε，如果‖xk－xk－1‖∞＜ε且‖yk－yk－1‖∞＜ε，則xk和yk為非線性微分方程組的數(shù)值解，否則計算xk－1和yk－1直至數(shù)值解滿足控制精度或達到設置的最大迭代次數(shù)為止。

3 數(shù)值算例

算例1的耦合系統(tǒng)參數(shù)為m＝1 kg，k1＝k2＝1 N／m，時間區(qū)域?。?，3］，時間節(jié)點t＝［t1，t2，…，tn］T分別采用等距節(jié)點和第二類Chebyshe點［18］。算例1解析解為［9］

式中，cn為雅可比橢圓函數(shù)，60個節(jié)點重心有理插值迭代配點法計算解與解析解比較如圖2所示。時間節(jié)點數(shù)量與數(shù)值解誤差Ea＝maxi｜ue－ue｜和迭代次數(shù)的關(guān)系見表1?？梢钥闯觯捎幂^少數(shù)量節(jié)點，即可達到很高的計算精度。

表1 算例1重心有理插值迭代配點法計算精度

算例2 如圖3所示，k1為線性彈簧剛度，k2、k3為非線性彈簧剛度。該系統(tǒng)控制方程式（18）為［6］

給定初始條件由式（19）表示為

圖2 算例1重心有理插值迭代配點法與解析解計算結(jié)果圖

圖3 雙自由度耦合系統(tǒng)圖

耦合系統(tǒng)參數(shù)為m＝1 kg，k1＝k2＝k3＝1 N／m，時間區(qū)域取［0，3］，時間節(jié)點 t＝［t1，t2，…，tn］T采用等距節(jié)點和第二類 Chebyshe點［18］。算例解析解為［8］

其中，時間節(jié)點數(shù)量與數(shù)值解誤差 Ea＝maxi｜uc－ue｜和迭代次數(shù)的關(guān)系見表2。

表2 算例2重心有理插值迭代配點法計算精度

4 結(jié)論

通過本研究可知：

（1）數(shù)值算例表明，重心有理插值迭代配點法求解多自由度系統(tǒng)非線性振動問題時具有非常高的計算精度，可快速得到時間節(jié)點處的高精度數(shù)值解。當處理大量節(jié)點尤其是等距節(jié)點時，有理插值仍能保持應有的計算穩(wěn)定性。采用Lagrange多項式插值確定未知函數(shù)時，當節(jié)點數(shù)量很大時，計算矩陣接近奇異，因此，Lagrange多項式計算時不能采用數(shù)量較多的計算節(jié)點，限制了其計算精度的提高。

（2）在求解大跨度時間問題時，可以將求解區(qū)間分割為若干子區(qū)間，在每個子區(qū)間上采用重心有理迭代配點法計算，將前一區(qū)間末的函數(shù)值，作為下一區(qū)間計算的初始值，依次計算得到整個區(qū)間上節(jié)點的數(shù)值解。

［1］Moochhala Y.E.，Raynor S..Free vibration ofmulti-degree-offreedom nonlinear systems［J］.International Journal of Non-Linear Mechanics，1972，7（6）：651－661.

［2］Huang T.C..Harmonic oscillations of nonlinear two-degree-of-freedom systems［J］.Journal of Applied Mechanics，1995，22：107－110.

［3］Gilchrist A.O..The free oscillations of conservative quasilinear systemswith two-degrees-of-freedom［J］.International Journal of Mechanical Science，1961，3：286－311.

［4］Dimarogonas A.D，Haddad S..Vibration for Engineers［M］.New Jersey：Prentice-Hall，Englewood Cliffs，1992.

［5］Bayat M.，Shahidi M.，Barari A.，et al.The approximate analysis of nonlinear behavior of structure under harmonic loading［J］.International Journal of Physical Science，2010，5（7）：1074－1080.

［6］Bayat M.，Pakar I.，Shahidi M..Analysis of nonlinear vibration of coupled systems with cubic nonlinearity［J］.MECHANIKA，2011，17（6）：620－629.

［7］Ladygina Y.V.，Manevich A.I..Free oscillations of a nonlinear cubic system with two-degrees-of-freedom and close natural frequencies［J］.Journal of Applied Mathematics and Mechanics，1993，57：257－266.

［8］Cveticanin L..Vibrations of a coupled two-degree-of-freedom system［J］.Journal of Sound and Vibration，2001，247：279－292.

［9］Cveticanin L..Themotion of a two-mass system with non-linear connection［J］.Journal of Sound and Vibration，2002，252：361－369.

［10］Chen G..Applications of a generalized Galerkin'smethod to nonlinear oscillations of two-degree-of-freedom systems［J］.Journal of Sound and Vibration，1987，119：225－242.

［11］Jordan A.K..Modified Laplace Decomposition Method［J］.World Applied Sciences Journal，2012，18（11）：1481－1486.

［12］王兆清，唐炳濤，李樹忱，等.求解間斷邊值問題的重心插值單元配點法［J］.山東建筑大學學報，2009，24（2）：115－118.

［13］Rashididinia J.，GhasemiM.，Jalilian R..A collocationmethod for the solution of nonlinear one-dimensional parabolic equations［J］.Mathematical Sciences，2010，4（1）：87－104.

［14］Yan J.P.，Guo B.Y..A collocation method for lnitial value problems of second-order ODEs by using laguerre functions［J］.Numer Math Theory Method Application，2011，4（2）：283－295.

［15］Geeta A.，Brajesh K.S..Numerical solution of Burgers'equation with modified cubic B-spline differential quadrature method［J］.Applied Mathematics and Computation，2013，224：166－177.

［16］段英峰，王兆清，林本芳.重心有理插值配點法分析矩形板自由振動［J］.山東建筑大學學報，2009，24（6）：434－437.

［17］Ordokhani Y..A collocation method for solving nonlinear differential equations via hybrid of rationalized haar functions［J］.Journal of Science，Tarbiat Moallem University，2007，7（3）：223－232.

［18］李樹忱，王兆清.高精度無網(wǎng)格重心插值配點法—算法、程序及工程應用［M］.北京：科學出版社，2012.