張宇功，常 勝，王小斌

(蘭州交通大學(xué)數(shù)理學(xué)院，甘肅 蘭州730070)

從理論上對(duì)于非線性動(dòng)力系統(tǒng)的研究比較困難，特別是高維非線性系統(tǒng)，因此我們希望通過一定的方法對(duì)其進(jìn)行簡化，并希望簡化后的系統(tǒng)能夠保持原系統(tǒng)的動(dòng)力學(xué)特性，通過國內(nèi)外專家研究已經(jīng)產(chǎn)生了很多簡化系統(tǒng)的途徑，例如中心流形［1］理論、范式方法、Lyapunov-Schmidt方法等等.

Duffing方程是非線性理論中常用的代表性微分方程，盡管是從簡單物理模型中得出來的非線性振動(dòng)模型，但是其模型具有代表性.工程中的許多非線性振動(dòng)問題的數(shù)學(xué)模型都可以轉(zhuǎn)化為該模型，因此，人們對(duì)Duffing方程進(jìn)行了廣泛而深入的研究，并且對(duì)Duffing方程的混沌運(yùn)動(dòng)有了越來越深入的了解和認(rèn)識(shí)［2］.

本文首先描述了中心流形，用中心流形對(duì)一類比Duffing系統(tǒng)更一般的含有二次和三次非線性項(xiàng)的受迫振動(dòng)系統(tǒng)［3］(Duffing方程可以視為這一非線性動(dòng)力方程的特殊情形)降維，并研究其穩(wěn)定性及分岔特性，它代表著許多非線性振動(dòng)的控制方程［4］.最后用 Matlab 進(jìn)行仿真.

1 中心流形

考慮自治系統(tǒng)

其中f:D→Rn是連續(xù)可微的，且D?Rn是包含原點(diǎn)x=0的定義.

定義1［5］設(shè)M是一個(gè)連通的度量空間，滿足

1)M有開覆蓋{Uα}，即Uα是開集，M=Uα;

2)對(duì)任意的α∈A，Uα同胚于Rn中的單位開球:B={x|x∈ Rn，x＜ 1};

3)若Uα∩Uβ≠φ，hα和hβ分別是從Uα和Uβ到Rn中單位開球B的同胚映射，則h=hα·是從Rn中的開集hβ到Rn中開集hα的可微映射，且對(duì)任意x∈hβ，雅克比矩陣的行列式非零，即 Dh(x)≠0.則稱M是一個(gè)n維可微流形.

k維流形有嚴(yán)格的定義，見文獻(xiàn)［6］.

定義2［7］如果 η(x(0))=0?η(x(t))≡0，?t∈［0，t1)?R其中［0，t1)是解x(t)有定義的時(shí)間區(qū)間，則稱流形{η(x)=0}是方程(1)的不變流形.

現(xiàn)在假設(shè)f是2次連續(xù)可微的，則方程(1)可表示為

其中

是2次可微的，且

我們只考慮無法線性化處理的情況，故假設(shè)A有k個(gè)實(shí)部為0的特征值，m=n-k個(gè)特征值實(shí)部為負(fù).我們總可以找到一個(gè)相似變換矩陣T，將A轉(zhuǎn)換為分塊對(duì)角矩陣

其中A1的所有特征值實(shí)部為0，A2的所有的特征值實(shí)部為負(fù).顯然，A1是k×k矩陣，而A2是m ×m矩陣，應(yīng)用變量代換

將方程(1)轉(zhuǎn)換為

其中g(shù)1和g2具有～f的性質(zhì).具體地說，它們是2次連續(xù)可微的，且

其中i=1，2.若z=h(y)是方程(2)和(3)的不變流形，且h是光滑的，則如果

就稱z為中心流形，關(guān)于中心流形的一些概念見文獻(xiàn)［8］.

定理1［7］如果g1和g2是二次連續(xù)可微的，且滿足(4)，A1的所有特征值均實(shí)部為零，A2的所有特征值均實(shí)部為負(fù)，則存在一個(gè)常數(shù)δ＞0和對(duì)于所有y＜δ有定義的連續(xù)可微函數(shù)h(y)，使得z=h(y)是方程(2)和方程(3)的中心流形.

如果系統(tǒng)(2)～(3)的初始狀態(tài)位于中心流形，即z(0)=h(y(0))，那么對(duì)于所有的t≥0，解(y(t)，z(t))將位于該流形內(nèi)，即z(t)≡h(y(t)).在這種情況下，中心流形內(nèi)系統(tǒng)的運(yùn)動(dòng)可由k階方程

描述，這個(gè)方程稱為降階系統(tǒng).

定理2［7］在定理1的條件下，如果降階系統(tǒng)(5)的原點(diǎn)y=0是漸近穩(wěn)定的(或非穩(wěn)定的)，則整個(gè)系統(tǒng)(2)～(3)的原點(diǎn)也是漸近穩(wěn)定的(或非穩(wěn)定的).

如果能夠求出臨界情形下非線性系統(tǒng)(1)的中心流形，那么就可以通過研究其約化系統(tǒng)的穩(wěn)定性來研究系統(tǒng)(1)的臨界穩(wěn)定性問題.根據(jù)中心流形對(duì)流的不變性，中心流形z=h(y)的求解可考慮系統(tǒng)(1)的臨界穩(wěn)定性問題，根據(jù)中心流形的不變性，在應(yīng)用定理2時(shí)，需要求出中心流形z=h(y).

函數(shù)h為偏微分方程

的一個(gè)解，其邊界條件為

在大多數(shù)情況下，該方程對(duì)h不能準(zhǔn)確求解(如果能獲得準(zhǔn)確的解，則表明求出整個(gè)系統(tǒng)的解)，但該解能夠以y的泰勒級(jí)數(shù)任意逼近.

定理3［7］如果可以找到一個(gè)連續(xù)可微的函數(shù)φ(y)，且φ(y)=0，［?φ/?y］(0)=0，使得對(duì)于p ＞1，有N(φ(y))=O(yp)，則對(duì)于足夠小的y，有h(y)-φ(y)=O(yp)且降階系統(tǒng)可表示為y˙=A1y+g1(y，φ(y))+O(yp+1).

對(duì)于標(biāo)量方程y˙=ayp+O(yp+1)，其中 p為正整數(shù).如果p是奇數(shù)且a＜0，則原點(diǎn)是漸近穩(wěn)定的;如果p是奇數(shù)且a＞0，或者p是偶數(shù)且a≠0，則原點(diǎn)是不穩(wěn)定的.

2 質(zhì)量 －彈簧系統(tǒng)模型

圖1所示的質(zhì)量 -彈簧系統(tǒng)中［7］，在水平面上滑動(dòng)并通過彈簧連接到豎直表面的物體m受到一個(gè)外力F的作用.

定義物體距離參考點(diǎn)的位移為x，根據(jù)牛頓運(yùn)動(dòng)定律，有

其中Ff是摩擦阻力，這里我們考慮摩擦阻力為黏滯摩擦力，即當(dāng)物體黏滯介質(zhì)，如空氣或潤滑劑中運(yùn)動(dòng)時(shí)，會(huì)有由于黏滯性引起的摩擦力.這個(gè)力通常按照速度的非線性函數(shù)建模，即Ff=h(v)，h(0)=0.當(dāng)速度較小時(shí)可假設(shè)為Ff=αv.Fsp是彈簧的回復(fù)力.設(shè)Fsp只是位移x的函數(shù)，即Fsp=g(x)，同時(shí)假設(shè)參考點(diǎn)位于g(x)=0處，外力由我們?cè)O(shè)定.對(duì)于不同的F，F(xiàn)f和g，會(huì)出現(xiàn)幾個(gè)不同的自治和非自治二階系統(tǒng)模型.

對(duì)于其他幾種情我們不予考慮，我們只考慮彈簧的恢復(fù)力為

和線性黏滯摩擦力

時(shí)系統(tǒng)的特性.

可得到無驅(qū)動(dòng)力的含二次和三次非線性項(xiàng)的系統(tǒng):

其中α為阻尼系數(shù)，β，γ為常數(shù).

顯然當(dāng)β=0，F(xiàn)=Acos(w t)時(shí)，系統(tǒng)就是典型的Duffing系統(tǒng).

為研究方便，接下來假設(shè)當(dāng)m=1時(shí)，系統(tǒng)隨α變化時(shí)的分岔行為.

設(shè)x1=x，x2=，則(1)式可化為

其雅克比矩陣為

其特征值為λ1=0，和λ2=-α，A有實(shí)部為零和實(shí)部為負(fù)的特征值，因此原點(diǎn)為非雙曲平衡點(diǎn).

令M是矩陣A對(duì)應(yīng)的基解矩陣，

即，

系統(tǒng)在原點(diǎn)的雅克比矩陣為

令T=M-1

則，可得

運(yùn)用變量代換如下

為了理論研究方便，不妨取β=γ=1，則將系統(tǒng)化為

該系統(tǒng)有中心流形，

嘗試令 y=h(z)=az2+O(3).

由中心流形偏微分方程

且由

比較兩邊x的同次冪系數(shù)得系統(tǒng)的中心流形近似表達(dá)式為

因此將其代入

可得

化簡可得

因此原系統(tǒng)的穩(wěn)定性與該約化系統(tǒng)的穩(wěn)定性等價(jià)，即原點(diǎn)是不穩(wěn)定的.

3 數(shù)值仿真驗(yàn)證

首先考慮有阻尼無驅(qū)動(dòng)情況即對(duì)于F=0時(shí)的情況，此時(shí)w取何值并不影響系統(tǒng)模擬.初值取(0.1，0.2)，α =0.2，β =1，γ =0.3 時(shí)從時(shí)間響應(yīng)圖(圖1)可以看出雖然初值相差0.1但在有限時(shí)間內(nèi)2條曲線很快會(huì)成為1條重合的曲線，這說明了系統(tǒng)穩(wěn)定.當(dāng) α =0.2，β =1，γ =0.3 時(shí)的相圖2 和 α =0.5，β=1，γ=0.3時(shí)的相圖3可以看出隨阻尼系數(shù)增大，系統(tǒng)會(huì)很快恢復(fù)到平衡態(tài).從Lyapunov指數(shù)(圖4)及分岔(圖5)可以看出，系統(tǒng)隨阻尼系數(shù)α增大從開始的混沌狀態(tài)很快進(jìn)入到穩(wěn)定狀態(tài)，從相圖及Poincare截面(圖6)也可得出系統(tǒng)是穩(wěn)定的.

從初值取(0.1，0.2)，α =0.2，β =1，γ =1 時(shí)的相圖7和α =0.2，β=-1，γ=1時(shí)的相圖8可以看出，β只改變平衡點(diǎn)的位置而不改變系統(tǒng)的穩(wěn)定性.從初值取(0.1，0.2)，α =0.2，β =1，γ =1 時(shí)的相圖7和α =0.2，β=1，γ =-1時(shí)的相圖9及時(shí)序圖10可以看出，γ的正負(fù)可以影響系統(tǒng)平衡點(diǎn)的穩(wěn)定性，這與我們前面的結(jié)論剛好一致.

接下來我們考慮有阻尼有驅(qū)動(dòng)的情況即F≠0，我們?nèi)=1.初值取(0.1，0.2)，α =0.2，β =1，γ =0.3，w=1，F(xiàn)=1 時(shí)的相圖12 和 α =0.2，β =1，γ =0.3，w=1，F(xiàn)=10時(shí)的相圖13可以看出當(dāng)F增大時(shí)，系統(tǒng)不容易恢復(fù)到平衡態(tài)，并且隨F增大系統(tǒng)會(huì)出現(xiàn)分岔及混沌現(xiàn)象.從圖14～16同樣可以得出在其它因素不變時(shí)，隨F增大系統(tǒng)會(huì)出現(xiàn)分岔及混沌現(xiàn)象.

4 結(jié)語

本文首先描述了中心流形，其次用中心流形對(duì)比Duffing系統(tǒng)更一般的含有二次和三次非線性項(xiàng)的受迫振動(dòng)系統(tǒng)降維，并研究其穩(wěn)定性及分岔特性.通過數(shù)值模擬，我們發(fā)現(xiàn)各參數(shù)對(duì)系統(tǒng)影響程度有所不同，從線性化可以看出由于α＞0，若γ＞0則(，0)為一漸近穩(wěn)定平衡點(diǎn)，若γ＜0則(，0)為一不穩(wěn)定平衡點(diǎn).而 其穩(wěn)定性與β無關(guān).當(dāng)考慮驅(qū)動(dòng)力時(shí)參數(shù)F對(duì)系統(tǒng)的影響最大，此時(shí)系統(tǒng)會(huì)出現(xiàn)新的平衡點(diǎn)，隨參數(shù)F的增大，系統(tǒng)會(huì)出現(xiàn)分岔現(xiàn)象和混沌現(xiàn)象，因此我們對(duì)于有驅(qū)動(dòng)力的系統(tǒng)應(yīng)該相應(yīng)增大阻尼系數(shù)來保證系統(tǒng)的穩(wěn)定性.

［1］CARR J.Applications of centre manifold theory［M］.New York:Springer－ Verlag，1981.

［2］HOLMES P，MARSDEN J.A partial differential equation with infinitely many periodic orbits:chaotic oscillations of a forced beam［J］.Archive for Rational Mechanics and A-nalysis，1981，76(2):135 －165.

［3］韓強(qiáng)，張善元，楊桂通.一類非線性動(dòng)力系統(tǒng)混沌運(yùn)動(dòng)的研究［J］.應(yīng)用數(shù)學(xué)和力學(xué)，1999，20(8).

［4］韓強(qiáng).幾種結(jié)構(gòu)的動(dòng)力屈曲、分叉和混沌運(yùn)動(dòng)研究［D］.太原:太原工業(yè)大學(xué)，1996.

［5］馬知恩，周義倉.常微分方程定性與穩(wěn)定性方法［M］.北京:科學(xué)出版社.2001:211 －212..

［6］GUILLEMIN V，POLLACK A.Differential topology［M］.American Mathematical Soc.，2010.

［7］HASSAN K.非線性系統(tǒng)［M］.朱義勝，董輝，李作洲，等，譯.北京:電子工業(yè)出版社，2005.

［8］BARREIRA L，VALLS C.Center manifolds for difference equations—smooth parameter dependence［J］.Nonlinear Analysis:Theory，Methods ＆ Applications，2010，73(3):725－749.