馬召召 楊慶超 周瑞平

1) (武漢理工大學(xué)能源與動(dòng)力工程學(xué)院，武漢 430063)

2) (海軍工程大學(xué)艦船與海洋學(xué)院，武漢 430033)

Lyapunov 指數(shù)是識(shí)別系統(tǒng)非線性動(dòng)力學(xué)特征的重要標(biāo)志，但是目前的算法通用性不足且計(jì)算流程復(fù)雜.本文在經(jīng)典的Lyapunov 指數(shù)算法的基礎(chǔ)上，基于攝動(dòng)理論提出了一種適用于不連續(xù)系統(tǒng)的Lyapunov 指數(shù)計(jì)算方法.首先，以系統(tǒng)狀態(tài)參數(shù)初始值和沿相空間每個(gè)基本矢量的擾動(dòng)量為初始條件，確定相軌跡.其次，采取差商近似導(dǎo)數(shù)法，獲得Jacobi 矩陣的近似矩陣.然后，對Jacobi 矩陣進(jìn)行特征值提取，得到系統(tǒng)的Lyapunov 指數(shù)譜.最后，將新算法應(yīng)用到二自由度干摩擦沖擊振蕩器系統(tǒng)實(shí)例中，并將計(jì)算結(jié)果與同步方法的計(jì)算結(jié)果進(jìn)行對比，對新算法的有效性進(jìn)行驗(yàn)證.該算法不僅適用于離散系統(tǒng)和連續(xù)時(shí)間系統(tǒng)，而且能夠快速計(jì)算復(fù)雜不連續(xù)系統(tǒng)的Lyapunov 指數(shù)，為確定復(fù)雜不連續(xù)系統(tǒng)的動(dòng)力學(xué)行為提供了新思路.

1 引言

Lyapunov 指數(shù)是在局部和全局穩(wěn)定性準(zhǔn)則下確定系統(tǒng)穩(wěn)定性的數(shù)值，它表征了系統(tǒng)在相空間中相鄰軌道間收斂或發(fā)散的平均指數(shù)率[1].從數(shù)學(xué)角度，Lyapunov 指數(shù)是狀態(tài)轉(zhuǎn)移矩陣的特征值，決定了系統(tǒng)吸引子上軌跡無窮小擾動(dòng)的演化進(jìn)程.因此，Lyapunov 指數(shù)被認(rèn)為是衡量系統(tǒng)對初始條件敏感性的重要指標(biāo).所有Lyapunov 指數(shù)值的和決定了系統(tǒng)在相空間中體積的發(fā)散度，所有Lyapunov指數(shù)值的和為負(fù)值是存在全局穩(wěn)定吸引子的充要條件，并且表明系統(tǒng)是耗散的;若所有Lyapunov指數(shù)均為負(fù)值，則吸引子為穩(wěn)定點(diǎn);對于n維系統(tǒng)，如果前m個(gè)Lyapunov 指數(shù)等于零而另外n—m個(gè)Lyapunov 指數(shù)均為負(fù)值，則吸引子為m維環(huán)面;若最大Lyapunov 指數(shù)為正值，則為混沌運(yùn)動(dòng)，若存在不止1 個(gè)Lyapunov 指數(shù)為正值，則為超混沌運(yùn)動(dòng).對于保守系統(tǒng)，所有Lyapunov 指數(shù)值的和為零，若所有Lyapunov 指數(shù)值的和為正值，則表明系統(tǒng)是全局不穩(wěn)定的，相軌跡呈指數(shù)發(fā)散.

Benettin 等[2]及Shimada 和Nagashima[3]首先提出了1 個(gè)基于Oseledets 定理[4]計(jì)算Lyapunov 指數(shù)譜的有效算法，文獻(xiàn)[5-8]對算法進(jìn)行了改進(jìn).近年來，提出了基于系統(tǒng)擾動(dòng)的數(shù)量積和導(dǎo)數(shù)計(jì)算Lyapunov 指數(shù)的方法[9-13]，該算法可以計(jì)算連續(xù)系統(tǒng)的Lyapunov 指數(shù)譜，但是如果系統(tǒng)中存在不連續(xù)性則會(huì)出現(xiàn)較大誤差或數(shù)值溢出現(xiàn)象.

為了確定具有沖擊或干摩擦的機(jī)械系統(tǒng)、分段線性特性電振蕩器等不連續(xù)系統(tǒng)的Lyapunov 指數(shù)，可以采用等效的連續(xù)函數(shù)近似不連續(xù)分量[14]、最小二乘陰影法[15]和偽軌道法[16]消除狀態(tài)擾動(dòng)矢量的不連續(xù)性等方法，雖然不需要重構(gòu)相空間和Jacobi 矩陣，但容易引起系統(tǒng)動(dòng)力學(xué)特性的質(zhì)變.Takens[17]和Wolf 等[18]根據(jù)時(shí)間序列信號(hào)重構(gòu)系統(tǒng)吸引子，提出了計(jì)算不連續(xù)系統(tǒng)最大Lyapunov 指數(shù)的算法.多年以來，文獻(xiàn)[19-27]針對該類方法的數(shù)據(jù)量、穩(wěn)定性和計(jì)算速率等方面進(jìn)行了評估改進(jìn)，該類方法的優(yōu)點(diǎn)是物理意義明顯，便于理解，計(jì)算結(jié)果不易受到拓?fù)鋸?fù)雜性的影響，有一定的抗噪能力.近年來，部分學(xué)者又提出了映射法[28-30]、補(bǔ)償矩陣法[31]、完全同步法[32-36]、克隆動(dòng)力學(xué)法[37]等，這些方法對一些特性類型的不連續(xù)系統(tǒng)具有較高的準(zhǔn)確性，但算法通用性不強(qiáng)且計(jì)算流程較為復(fù)雜.

本文提出了一種計(jì)算不連續(xù)系統(tǒng)Lyapunov指數(shù)譜的新方法.該方法基于攝動(dòng)理論，使用正交矢量型初始擾動(dòng)計(jì)算Jacobi 矩陣，對Jacobi 矩陣進(jìn)行特征值提取，得到系統(tǒng)的Lyapunov 指數(shù)譜.該方法僅需要選定初始條件，無需確切知道系統(tǒng)的方程式，便能夠確定系統(tǒng)的最終狀態(tài)，而且對不連續(xù)系統(tǒng)和連續(xù)系統(tǒng)均適用.

2 方 法

為了克服現(xiàn)有Lyapunov 指數(shù)算法的不足，以一般性的理論為依據(jù)，使得新算法適用于更多的系統(tǒng)并更具通用性;另一方面，以經(jīng)典的數(shù)值算法為基礎(chǔ)，對廣泛應(yīng)用的Jacobi 法進(jìn)行改進(jìn)，保證新算法計(jì)算的準(zhǔn)確性.

2.1 算法理論

以常微分方程(ODE)形式描述一種n維自治連續(xù)時(shí)間系統(tǒng):

其中，x ∈Rn為狀態(tài)向量，t ∈R為時(shí)間，f為Rn →Rn上的連續(xù)可微函數(shù).假設(shè)St(x0) 為滿足ODE(1)初始條件x0的解(軌跡)，則有

Lyapunov 指數(shù)的定義是基于以下Jacobi 矩陣建立的:

(2)式對x0進(jìn)行微分獲得如下變分方程:

矩陣Ut(x0) 顯示了初始條件x0的無窮微擾對軌跡St(x0)的影響:

其中Δs(t)是由初始擾動(dòng)量Δx0引起St(x0) 的擾動(dòng)，即: Δs(t)St(x0+Δx0)-St(x0) .向量Δs(t)的長度可以從下式中求得

令ui(t)，i1，···，n為Ut(x0)TUt(x0)的特征值.Ut(x0)是1個(gè)實(shí)矩陣，所以Ut(x0)TUt(x0) 是實(shí)對稱的.因此，ui(t)∈R和對應(yīng)的特征向量可以構(gòu)成正交基vi(t)，i1，···，n.由于任意1個(gè)Δx0∈Rn，都有 |Δs(t)|2≥0，可得到矩陣Ut(x0)TUt(x0)為半正定，且ui(t)≥0.因此，如果Δx0與vi(t) 平行，則 |Δs(t)||Δx0|可以得到:若初始擾動(dòng)量Δx0的1個(gè)分量與vi(t) 平行，則其軌跡長度的影響因子為可用以下公式表示:

其中 [a，b]表示向量a，b之間的數(shù)量積.

Lyapunov 指數(shù)的定義為

從(8)式可以得到，在時(shí)間t足夠長的情況下，

因此，對于每個(gè)Lyapunov 指數(shù)值λi，在相空間中都存在1 個(gè)對應(yīng)的方向，使得初始擾動(dòng)量Δx0在該方向上的投影長度乘以 eλit作為系統(tǒng)在時(shí)間t內(nèi)的演化.可用以下公式表示:

這意味著每個(gè)Lyapunov 指數(shù)是初始擾動(dòng)量Δx0的分量沿相空間某個(gè)方向的平均指數(shù)收縮率(如果λi ＜0)或擴(kuò)展率(如果λi ＞0).

Lyapunov 指數(shù)的概念不僅適用于連續(xù)時(shí)間系統(tǒng)，而且適用于離散時(shí)間系統(tǒng).上述理論為連續(xù)時(shí)間和離散時(shí)間系統(tǒng)提供了Lyapunov 指數(shù)計(jì)算的一般方法.首先，對于足夠長的時(shí)間t，使用(4)式計(jì)算Jacobi 矩陣Ut(x0) ;然后，計(jì)算Ut(x0)TUt(x0)的特征值ui;最后，應(yīng)用(8)式計(jì)算每個(gè)Lyapunov指數(shù)值.

2.2 連續(xù)系統(tǒng)經(jīng)典算法

基于以上推導(dǎo)公式建立了一種連續(xù)系統(tǒng)的經(jīng)典算法[8].該算法利用了最大的擾動(dòng)子空間Wi中的每個(gè)初始擾動(dòng)均以 eλit的速度變化這一特點(diǎn)，即平行于vi的初始擾動(dòng)量Δx0分量的長度在時(shí)間t上以因數(shù) eλit進(jìn)行擴(kuò)展或收縮.

假設(shè)n維系統(tǒng)的Lyapunov 指數(shù)λ1≥λ2≥···≥λn，幾乎任意1 個(gè)初始擾動(dòng)量 Δx0都具有與v1平行的分量，則該分量長度的變化因數(shù)為 eλ1t.如果λ1是最大Lyapunov 指數(shù)，則其分量將成為主導(dǎo)量.由于擾動(dòng)量 Δx0的方向與v1對齊，則 Δx0的長度以近似指數(shù)速率λ1變化，即 |Δs(t)|≈|Δx0|eλ1t.這意味著λ1決定了幾乎所有初始擾動(dòng)改變其長度的平均指數(shù)速率.為了計(jì)算n個(gè)Lyapunov 指數(shù)譜，通過Gram-Schmidt 正交歸一化獲得相互正交的初始擾動(dòng)量 Δx，減少擾動(dòng)量之間對齊的影響.為了計(jì)算λ2，當(dāng) Δx0與v1對齊時(shí)，必須計(jì)算另1 個(gè)與 Δx1正交的新擾動(dòng)，且該擾動(dòng)的長度以近似指數(shù)速率λ2變化，并與v2對齊，如圖1 所示.

圖1 兩個(gè)擾動(dòng)量 Δx，Δx0 的正交歸一化Fig.1.Orthonormalization of two perturbations Δx， Δx0 .

以此類推，為了計(jì)算λ3，必須計(jì)算1 個(gè)與前面兩個(gè)擾動(dòng)量正交的新擾動(dòng)，則其與v1和v2正交，以近似指數(shù)速率λ3變化，并與v3對齊.總的來說，若計(jì)算n個(gè)Lyapunov 指數(shù)譜，必須計(jì)算n個(gè)不同的擾動(dòng)，第1 個(gè)是根據(jù)等式(5)自由演化的，其他第i個(gè)新擾動(dòng)都與第 1，···，(i-1) 個(gè)擾動(dòng)量保持正交.

該算法唯一的約束是正交向量必須與原始向量跨越相同的擾動(dòng)子空間.同時(shí)，由于矩陣Ut(x0)在時(shí)間t較長的情況下變得計(jì)算困難，甚至無法計(jì)算，使得該經(jīng)典算法在實(shí)踐中效果不佳.

2.3 不連續(xù)系統(tǒng)Lyapunov 指數(shù)的算法

2.3.1 算法思想

在連續(xù)系統(tǒng)Lyapunov 指數(shù)經(jīng)典算法[8]的基礎(chǔ)上，基于算法理論詳述利用差商近似導(dǎo)數(shù)法獲得Jacobi 矩陣計(jì)算不連續(xù)系統(tǒng)Lyapunov 指數(shù)的具體算法.

算法思想:首先設(shè)定初始條件x0和計(jì)算參數(shù);然后求解系統(tǒng)，得到充分接近實(shí)際系統(tǒng)吸引子的軌跡;再利用差商近似倒數(shù)方法計(jì)算不連續(xù)系統(tǒng)的Jacobi 矩陣;利用Gram-Schmidt 正交歸一化獲得相互正交的初始擾動(dòng)量，并將不同的擾動(dòng)量用于求解相應(yīng)的Lyapunov 指數(shù)，進(jìn)而得到系統(tǒng)的Lyapunov 指數(shù)譜.

2.3.2 算法實(shí)現(xiàn)的具體方法

下面詳述該算法在不連續(xù)系統(tǒng)中實(shí)現(xiàn)的具體方法.

對于連續(xù)系統(tǒng)，方程(1)中的矢量場f在相空間軌跡的任意點(diǎn)x處是連續(xù)可微分的，利用經(jīng)典算法[8]可以得出系統(tǒng)Lyapunov 指數(shù).否則，不能直接從變分方程(4)獲得Jacobi 矩陣Ut(xj) .針對不連續(xù)系統(tǒng)，進(jìn)行適應(yīng)性修改.首先考慮連續(xù)時(shí)間不連續(xù)系統(tǒng)，系統(tǒng)(1)的解St(x) 實(shí)際上是將初始狀態(tài)x轉(zhuǎn)換為時(shí)間t之后系統(tǒng)(1)狀態(tài)的映射.Ut(x)為在點(diǎn)x處求得軌跡St的Jacobi 矩陣:

其中e1，···，en是Rn中的標(biāo)準(zhǔn)基，Ut(x) 的 列 是St(x)關(guān)于后續(xù)狀態(tài)變量的偏導(dǎo)數(shù).根據(jù)導(dǎo)數(shù)定義，(10)式中的偏導(dǎo)數(shù)可以表示為Δ→ 0 時(shí)差商的極限.因此，可通過以下方式估算矩陣Ut(x) :

在選擇Δ值時(shí)，需要在與Δ成比例的截?cái)嗾`差和取決于實(shí)際使用數(shù)據(jù)類型的近似誤差之間取得平衡，矩陣Ut(x) 的近似算法如圖2 所示.

圖2 矩陣 Ut(x) 逼近算法示意圖Fig.2.Illustration of the algorithm of the Ut(x) matrix approximation.

為了估計(jì)Jacobi 矩陣Ut(x)，需要利用未擾動(dòng)的初始條件x以及沿著相空間的每個(gè)受Δ值擾動(dòng)的初始條件x+Δei，i1，···，n，確定映射St(x)的值.總之，估算一次矩陣Ut(x)，需要對映射St(x)進(jìn)行n+1次計(jì)算:St(x)，St(x+Δe1)，···，St(x+Δen)，并代入(11)式得到近似Jacobi 矩陣Ut(x) .新算法是將上述Jacobi 矩陣的估算方法與經(jīng)典Lyapunov 指數(shù)算法[8]相結(jié)合，其中的所有其他步驟都與經(jīng)典Lyapunov 指數(shù)算法相同，主要區(qū)別在于Jacobi 矩陣是從近似公式(11)而不是(4)式獲得.令tj ∈R+，j0，1，···，J為1個(gè)遞增的時(shí)間序列，xjx(tj) .假設(shè)初始條件t00并且x0x(t0)x0，則 Δtjtj -tj-1和x(j-1)代入(11)式中就可以得到矩陣(x(j-1)) 的近似值.這種方法的最大優(yōu)點(diǎn)是，不再要求軌跡上每個(gè)點(diǎn)的向量場f必須連續(xù).應(yīng)用(11)式時(shí)，只要滿足對于每個(gè)tj，軌跡(x) 在點(diǎn)x(j-1)處對于所有狀態(tài)變量是可微的，則系統(tǒng)中不論存在多少不連續(xù)性或不連續(xù)的類型均能滿足要求.同時(shí)，對時(shí)間序列tj僅要求tj ∈R+，j0，1，···，J是遞增的，tj值大小沒有限制.而且，不需要明確知道定義矢量場f的方程式，只要有一種計(jì)算軌跡(x) 的方法就足夠了.(x)的值無論是從軌跡(x) 的顯式定義或隱式定義，還是從數(shù)值過程，甚至是從實(shí)驗(yàn)中獲得都沒有關(guān)系，該算法均有效.

以上提出的方法可以通過調(diào)整使其適用于離散時(shí)間不連續(xù)系統(tǒng).與連續(xù)時(shí)間系統(tǒng)類似，離散時(shí)間系統(tǒng)的U(x) 可以理解為在點(diǎn)x處評估映射F的Jacobi 矩陣.可以表示為

因此，可以通過以下方式估算矩陣U(x) :

估算一次矩陣U(x)，需要對映射F進(jìn)行n+1次計(jì)算:F(x)，F(xiàn)(x+Δe1)，···，F(xiàn)(x+Δen)，并代入(13)式得到近似Jacobi 矩陣U(x) .

基于以上分析，以連續(xù)時(shí)間不連續(xù)系統(tǒng)為例，不連續(xù)系統(tǒng)Lyapunov 指數(shù)算法如圖3 所示.

圖3 不連續(xù)系統(tǒng)Lyapunov 指數(shù)算法的示意圖Fig.3.Illustration of the algorithm of discontinuous system Lyapunov exponent.

不連續(xù)系統(tǒng)Lyapunov 指數(shù)算法實(shí)現(xiàn)的步驟如下.

2)使用初始條件x(j-1)在不連續(xù)系統(tǒng)中獲得新狀態(tài)x(j).

4)計(jì)算新的擾動(dòng):

5)擾動(dòng)正交化:

6)對于j1，2，···，T，重復(fù)以上2)—5)過程.然后，使用以下公式估算不連續(xù)系統(tǒng)的整個(gè)Lyapunov 指數(shù)譜:

3 數(shù)值仿真研究

3.1 二自由度干摩擦沖擊振蕩器系統(tǒng)

二自由度干摩擦沖擊振蕩器系統(tǒng)由連接到剛度為k1的彈簧的質(zhì)量塊M和阻尼系數(shù)為c1的黏滯阻尼器組成，另1 個(gè)質(zhì)量塊m位于質(zhì)量塊M上，其相對于M的位置由緩沖器A和B限制，質(zhì)量塊之間存在摩擦力FT和阻尼系數(shù)為c2的黏性阻尼力.兩自由度沖擊振子是1 個(gè)質(zhì)量-彈簧-阻尼系統(tǒng)，其中兩個(gè)質(zhì)量塊可以相互碰撞，并在相對運(yùn)動(dòng)過程中存在干摩擦，如圖4 所示.

圖4 二自由度干摩擦沖擊振蕩器模型Fig.4.Model of the two-degrees-of-freedom (2-DoF) mechanical oscillator with impacts and friction.

二自由度沖擊振蕩器可由以下方程組描述:

其中M和m分別是外部和內(nèi)部振動(dòng)體的質(zhì)量;x1，x2是它們的位置;k1是彈簧剛度;c1，c2是阻尼系數(shù);FT是摩擦力;A是強(qiáng)迫振幅;Ω是諧波激勵(lì)的角頻率;δ0是內(nèi)部質(zhì)量塊m與外部質(zhì)量塊M之間的間隙;μ是摩擦系數(shù)，R是恢復(fù)系數(shù);tc是發(fā)生碰撞的時(shí)間.方程組(16)定義了振蕩器系統(tǒng)的演化，方程組(18)描述了在系統(tǒng)發(fā)生碰撞時(shí)刻tc處不連續(xù)性的狀態(tài)轉(zhuǎn)換.

將(16)—(18)式轉(zhuǎn)換為無量綱表示得到

在這個(gè)振蕩器系統(tǒng)中，沖擊和干摩擦的兩種不連續(xù)性都存在.假設(shè)對于任何初始條件y(0)，使得|y3-y1|＜y0，可以數(shù)值估計(jì)軌跡Sτ(y(0)) .唯一的限制是兩個(gè)物體在碰撞過程中的速度沒有定義.此外，還假設(shè)當(dāng) |y3(0)-y1(0)|＜y0和y2/y4時(shí)，軌跡Sτ(y(0))相對于初始條件是可微的.

3.2 數(shù)值仿真結(jié)果

為了將數(shù)值仿真結(jié)果與不同算法[35]獲得的結(jié)果進(jìn)行比較驗(yàn)證，控制參數(shù)η在[1.2，1.45]范圍內(nèi)變化，并設(shè)定以下參數(shù)值:γ0.693，σ0.50，μ0.02，β0.05，y00.80，p1.00，R0.60，初始條件為零.采用自適應(yīng)步長的龍格-庫塔(Runge-Kutta)方法對二自由度沖擊振蕩器系統(tǒng)進(jìn)行了數(shù)值求解.在開始計(jì)算Lyapunov 指數(shù)和將軌跡點(diǎn)保存到分岔圖之前，每組參數(shù)針對運(yùn)行該振蕩器系統(tǒng)，以減少瞬態(tài)影響.其中，Jacobi矩陣(x(j-1)) 通過(11)式獲得，為了保持方法的準(zhǔn)確性，Δ值應(yīng)與變分方程(4)的積分絕對誤差處于同一階或更低階，并考慮沖擊振蕩器系統(tǒng)中數(shù)據(jù)類型的近似誤差，從而確定Δ值取10—6.

由于最小擾動(dòng) Δy4的收縮速度過快，導(dǎo)致在計(jì)算自然對數(shù)時(shí)出現(xiàn)了數(shù)值問題.因此，僅給出了該振蕩器系統(tǒng)的3 個(gè)最大的Lyapunov 指數(shù).為了觀察Lyapunov 指數(shù)值的穩(wěn)定性，在估算程序的每次迭代之后，再利用(15)式計(jì)算結(jié)果.在每100 次迭代后，估計(jì)最后100 個(gè)λ1，λ2，λ3值的標(biāo)準(zhǔn)差.當(dāng)所有標(biāo)準(zhǔn)偏差低于閾值10—4時(shí)，終止計(jì)算，并取最后100 次迭代中λ1，λ2，λ3的平均值作為最終值.圖5給出了3 個(gè)最大Lyapunov 指數(shù)圖和系統(tǒng)狀態(tài)分岔圖.可以觀察到，Lyapunov 指數(shù)值恰當(dāng)?shù)胤从沉讼到y(tǒng)在分岔圖中的行為.因此，通過新算法在二自由度干摩擦沖擊振蕩器系統(tǒng)(19)—(21)中的應(yīng)用，證明了該Lyapunov 指數(shù)新算法的有效性.

圖5 二自由度干摩擦沖擊振蕩器系統(tǒng)的分叉圖和3 個(gè)最大Lyapunov 指數(shù)圖 (a) 狀態(tài)變量 y1;(b) 狀態(tài)變量 y2 ;(c) 狀態(tài)變量y3Fig.5.Bifurcation diagram and the corresponding diagram of the three largest Lyapunov exponents of the 2-DoF mechanical oscillator system with impacts and friction:(a) State variable y1;(b) state variable y2 ;(c) state variable y3 .

最后，對使用本文算法獲得的Lyapunov 指數(shù)計(jì)算結(jié)果與同步方法[35]得到的計(jì)算結(jié)果進(jìn)行比較，可以得到，使用兩種算法獲得的最大Lyapunov 指數(shù)圖高度一致.從而，兩種方法也相互驗(yàn)證了算法的正確性.

4 討論

在物理環(huán)境中，具有沖擊、干摩擦等因素的不連續(xù)系統(tǒng)，難以直接從變分方程中獲得Jacobi 矩陣，經(jīng)典的Lyapunov 指數(shù)算法[8]變得難以實(shí)現(xiàn)，因此，本文算法基于攝動(dòng)理論，在經(jīng)典Lyapunov指數(shù)算法基礎(chǔ)上，采用差商近似導(dǎo)數(shù)獲得近似Jacobi 矩陣的方法，以損失一定Jacobi 矩陣計(jì)算精度為代價(jià)提升算法普適性.

本文算法應(yīng)用于二自由度干摩擦沖擊振蕩器系統(tǒng)取得了較好的結(jié)果.計(jì)算得到了3 個(gè)Lyapunov指數(shù)譜圖和系統(tǒng)中3 個(gè)狀態(tài)變量的分叉圖，其中最大Lyapunov 指數(shù)也正確地反映了系統(tǒng)狀態(tài)變量在分叉圖中的軌跡變化，表明可以采用差商近似導(dǎo)數(shù)的方式獲得近似完備的Jacobi 矩陣，通過與連續(xù)系統(tǒng)經(jīng)典Lyapunov 指數(shù)算法相結(jié)合，能夠有效地應(yīng)用于包含沖擊、干摩擦等因素的不連續(xù)系統(tǒng).在本算法中，差商近似導(dǎo)數(shù)法利用n+1 個(gè)不同初始條件向量:x0，x0+Δe1，···，x0+Δen在時(shí)間t上對系統(tǒng)積分求解Ut(x)，因此，求解Ut(x) 相當(dāng)于求解n+1 階非線性系統(tǒng)的計(jì)算量.在經(jīng)典算法中，通常變分方程(4)與系統(tǒng)軌跡St(x) 一起求解，可以將變分方程(4)視為n個(gè)線性時(shí)變系統(tǒng)，因此，求解Ut(x) 相當(dāng)于求解1 個(gè)n階非線性系統(tǒng)和n個(gè)線性時(shí)變系統(tǒng)的計(jì)算量.只要系統(tǒng)狀態(tài)的時(shí)間導(dǎo)數(shù)值的計(jì)算量不比等式(4)中Ut(x) 后續(xù)列的時(shí)間導(dǎo)數(shù)的計(jì)算量大很多，則兩種算法的計(jì)算量是相當(dāng)?shù)?同時(shí)，許多非線性系統(tǒng)的矢量場f是通過狀態(tài)變量的和與積來定義的，都能夠滿足這個(gè)條件，例如著名的Duffing 或Lorenz 方程[11].差商近似導(dǎo)數(shù)法獲得近似Jacobi 矩陣中所需的每個(gè)軌跡St(x)，St(x+Δe1)，···，St(x+Δen)的積分是完全獨(dú)立的，每個(gè)軌跡都可以在使用單獨(dú)處理器的獨(dú)立進(jìn)程中進(jìn)行求解，使得新算法非常適合于多處理器計(jì)算模式，提升運(yùn)算效率.

本文算法目前適應(yīng)于對連續(xù)系統(tǒng)和不連續(xù)系統(tǒng)的Lypunov 指數(shù)譜進(jìn)行提取，而實(shí)際應(yīng)用中不連續(xù)系統(tǒng)的不連續(xù)性是多樣的，這就使得Lyapunov 指數(shù)的計(jì)算過程繁冗復(fù)雜.后續(xù)研究將圍繞如何繞過不連續(xù)系統(tǒng)的不連續(xù)點(diǎn)計(jì)算Lyapunov指數(shù)譜展開.

5 結(jié)論

本文基于經(jīng)典Lyapunov 指數(shù)算法和攝動(dòng)理論提出了一種計(jì)算不連續(xù)系統(tǒng)Lyapunov 指數(shù)的新算法，在經(jīng)典Lyapunov 指數(shù)算法基礎(chǔ)上，利用相軌跡上的正交矢量型初始微擾動(dòng)獲得近似Jacobi 矩陣.該算法不需要知道定義矢量場或者映射的方程式，只需要有一種獲得相軌跡的方法就足夠了，并且適用于連續(xù)或不連續(xù)的離散時(shí)間和連續(xù)時(shí)間系統(tǒng).通過數(shù)值仿真分析，驗(yàn)證了本文Lyapunov 指數(shù)算法簡單、有效和魯棒.本文提出的算法可以應(yīng)用于包括不連續(xù)滑動(dòng)振蕩器和許多其他不連續(xù)系統(tǒng)的大量不連續(xù)系統(tǒng)的研究.