細長剛體塊系統中的混沌現象及其Lyapunov指數譜計算

蔣金凱，杜正東

(四川大學數學學院，成都 610064)

1 引 言

在非光滑動力系統中有多種多樣的分岔現象.近年來，隨著非光滑動力系統在工程應用領域得到越來越多的應用[1，2]，對其動力學行為的研究也日益深入[3].其中，由不連續誘導分岔產生的非光滑動力系統的混沌現象是微分方程與動力系統研究的一個熱點.另一方面，原本應用于光滑動力系統的Lyapunov指數譜計算方法也被推廣應用在各種非光滑動力系統中[4，5].Lyapunov指數譜度量了系統相鄰軌道的平均收斂程度，可以作為一種判定系統混沌的指標.

細長剛體塊模型是一個典型的分段光滑系統，在研究地震對剛體結構的影響等方面有諸多應用.研究表明[6-8]，該模型的動力學行為非常豐富，不僅可以通過鞍結點分岔、音叉分岔、倍周期分岔等產生各種碰撞周期軌道，進而導致混沌，還可能通過異宿分岔產生混沌.本文首先將文獻[4]中的Lyapunov指數譜計算方法推廣應用于類似細長剛體塊的其它非光滑動力系統，然后通過數值模擬驗證：當外部周期激勵的振幅改變時，系統可以通過異宿分岔產生混沌，最后利用Lyapunov指數譜驗證了這種混沌現象.

我們首先引入以下概念.

定義2.1當動力系統的軌道反復穿越同一截面時，在連續運動的軌道上用一個截面(稱Poincaré截面)將其橫截，利用軌道在截面上穿過的情況就可以判斷運動的狀態.記S為Poincaré截面，xn(n=0,1,2,…)為軌道前一次穿過S的點，xn+1為軌道后一次穿過S的點，那么，xn+1可看成xn的映射：

xn+1=f(xn)，n=0,1,2,...,

稱為Poincaré映射.

顯然，由上述映射f可定義一個離散動力系統.

定義2.2對上述的離散動力系統，考慮一個起點在x0處的無限小位移矢量y0，其在映射f下演變過程可用以下方程表示：

yn+1=Df(xn)·yn，n=0,1,2,...，

其中Df表示映射f的Jacobi矩陣.|yn|/|y0|描述了無限小位移的增減比率.因此，系統在滿足初始條件x0及沿初始方向u0=y0/|y0|的Lyapunov指數定義為

Lyapunov指數描述了相空間相鄰軌道的平均指數發散率的數值特征.若映射f是n維的，則存在n個Lyapunov指數，分別對應于n個線性無關的初始方向.將這n個數值按從大到小的順序排列起來，所得結果稱為離散動力系統的Lyapunov指數譜.

圖1為細長剛體塊系統的力學模型.假設該剛體塊是質地均勻的，其重心與幾何中心R重合.同時，在施加外部水平激勵時，若不考慮滑動和彈跳現象，則該剛體塊在激勵作用下反復繞兩個底角A,B轉動.在每次變換轉動角時，剛體塊與地面發生碰撞.當α?1時，其無量綱運動方程由以下常微分方程[6]表示:

圖1 細長剛體塊模型

(1)

其中，x表示剛體塊的角位移，k表示系統的阻尼系數，β,ω分別表示外部水平激勵的振幅和頻率.當x=0時，碰撞發生，tA,tB分別表示發生碰撞的前后時刻.由于是剛性碰撞，因此tA=tB.r表示碰撞的能量恢復系數.進一步，我們將系統(1)改寫為如下自治系統:

y(tA)=ry(tB),x=0

(2)

由于碰撞的存在，系統軌道在經過碰撞面時不連續，這導致系統的向量場在該處的Jacobi矩陣不存在.因此，光滑系統中利用Jacobi矩陣計算Lyapunov指數譜的方法便不再適用.為了解決這一問題，我們采用類似文獻[4]中的方法，結合局部碰撞映射構造Poincaré映射，將非光滑系統轉化為離散動力系統，進而采用離散動力系統計算Lyapunov指數譜的方法.為了達到這一目的，我們首先介紹如何構建系統(2)的Poincaré映射及計算該映射的Jacobi矩陣.

3.1 Poincaré映射的建立

為了敘述方便，我們首先給出以下子空間的記號:

y>0,φ=φa},

y<0,φ=φb},

Σ+={(x,y,φ)∈R2×S1|x=0,y>0}=

{(y,φ)∈R×S1|y>0},

Σ-={(x,y,φ)∈R2×S1|x=0,y<0}=

{(y,φ)∈R×S1|y<0}.

(3)

(4)

(5)

上式中I為3維單位矩陣，·和*分別表示向量的內積和2階張量.

(i)從定相位面到碰撞面的局部子映射

(ii)從碰撞面到碰撞面的局部子映射

利用系統中的碰撞關系式，可得到局部映射P2(P5)的Jacobi矩陣為

(iii)從碰撞面到定相位面的局部子映射

由上述三類局部映射構成的復合映射為

上述映射的Jacobi矩陣可通過鏈式求導法則分別求得

下面接著定義系統(2)在兩個非碰撞階段的子映射，分別記為

因此，由式(3)定義的Poincaré映射P可以表示為以上映射的復合，即

(6)

由復合映射的求導法則，Poincaré映射P的Jacobi矩陣可寫成

(7)

3.2 Lyapunov指數譜的計算方法

根據式(3)(7)定義的Poincaré映射P及其Jacobi矩陣，令

(8)

則離散動力系統Lyapunov指數譜通常定義為[6]

(9)

在直接利用(9)式計算Lyapunov指數譜的過程中，由于要求迭代的次數N→∞，經常出現病態問題，導致M中的元素或者變得很大(對于混沌吸引子)，或者變為0(對于周期吸引子)，使得計算結果與實際情況不符.為了消除該問題，常采用QR分解將矩陣MN轉化為正交矩陣QN和上三角矩陣RN，文獻[4]中給出了其具體的計算過程.因此，Lyapunov指數譜的計算式轉化為

(10)

在得到系統的Lyapunov指數譜后，我們可以根據其特征判定系統此時的狀態.其中，周期解的所有Lyapunov指數都為負，擬周期解至少存在一個Lyapunov指數為0，且其余的Lyapunov指數都為負，而混沌解至少存在一個Lyapunov指數為正.

4 數值模擬

對細長剛體塊系統的研究可追溯到上世紀九十年代.Hogan在文獻[6]中研究了該系統可能存在的(1,n)型穩定周期軌道的參數范圍，其中(m,n)型周期軌道表示該周期軌道在一個周期內經歷了2m次碰撞，且其周期是外部激勵周期的n倍，并解析的給出了其發生鞍結點分岔與倍周期分岔的條件.進一步，Hogan還計算了系統的異宿軌的擾動穩定流形與擾動不穩定流形的解析形式，利用Melnikov方法得到了系統發生異宿分岔的臨界參數條件為[7]

(11)

其中，

并推測系統在滿足參數條件β>βM時可能會產生混沌，但沒有通過數值模擬進行驗證.下面我們選取滿足β>βM的適當參數，通過相應的數值模擬來觀察系統是否會產生混沌.

選取系統固定參數k=1,ω=5,r=0.925.此時βM≈12.4702.當β=15時，選取初始迭代點為(10-8,0)，去除前10 000個暫態點的影響后，系統存在穩定的對稱(1,1)型周期軌道，如圖2所示.其中，實線表示系統的碰撞面，虛線表示此時系統的穩態軌道.圖3展示了該吸引子的Lyapunov指數譜序列的收斂情況.從該圖中可見，序列的收斂速度和收斂性都很好，且系統在該參數條件下的兩個Lyapunov指數(分別以實線和虛線表示)的收斂值都小于0，其中最大Lyapunov指數λ1≈-0.039，表明系統此時處于周期運動狀態，與數值模擬結果相符.

圖2 (1,1)型周期吸引子相平面圖

圖3 周期吸引子Lyapunov指數譜收斂序列

當β=23.6時，去除暫態點影響后，系統吸引子的相平面圖如圖4所示.為了進一步確認系統是否處于混沌狀態，圖5展示了該吸引子的Lyapunov指數譜序列的收斂情況.如圖所示，系統在該參數條件下具有一個正的Lyapunov指數λ1≈0.260，由此可說明系統此時正處于混沌狀態.

圖4 混沌吸引子相平面圖

圖5 混沌吸引子Lyapunov指數譜收斂序列

通過以上的數值模擬，我們驗證了系統發生異宿分岔后的確會導致混沌.

5 結 論

通過對細長剛體塊系統進行數值模擬，我們證實了文獻[7]中關于異宿分岔會導致混沌的猜測.值得注意的是，這可能只是該系統產生混沌的必要條件，因為在相同的參數條件下，當改變系統初始迭代點時，其產生的穩態現象可能不同.此外，我們推廣了文獻[4]中一類碰撞振動系統利用Poincaré映射計算Lyapunov指數譜的方法.與文獻[4]不同的是，本文所研究的剛體塊系統在碰撞面兩端的向量場是不同且不連續的，而這導致了碰撞階段子映射的Jacobi矩陣計算上的差異，且構建的Poincaré映射P中包含了兩次碰撞過程.雖然在計算上略顯繁雜，但由于碰撞的特性相同，因而可將該方法推廣至具有多個碰撞面的Filippov系統.同時，本文利用Lyapunov指數譜驗證了剛體塊系統發生的混沌現象.數值模擬結果與Lyapunov指數計算結果的相容性也證明了我們推廣的Lyapunov指數計算方法的正確性和有效性.