姜宏杰，劉祚秋，呂中榮，劉廣

中山大學(xué)航空航天學(xué)院，廣東 廣州 510006

與整數(shù)階導(dǎo)數(shù)相比，分?jǐn)?shù)階導(dǎo)數(shù)具有方程簡(jiǎn)潔和參數(shù)意義明確等優(yōu)點(diǎn)，通過(guò)分?jǐn)?shù)階微分系統(tǒng)描述的現(xiàn)象更加接近于真實(shí)的自然現(xiàn)象，因此分?jǐn)?shù)階導(dǎo)數(shù)在非線性振動(dòng)［1］、控制理論［2］、參數(shù)識(shí)別［3-5］等領(lǐng)域有廣泛的應(yīng)用。vdP（van der Pol）系統(tǒng)［6］最初是荷蘭工程師Balthazar van der Pol在研究三極管的振蕩效應(yīng)中提出的。作為一種典型的自激系統(tǒng)，其在不同的參數(shù)條件下能夠展現(xiàn)出多種豐富的動(dòng)力學(xué)現(xiàn)象，因此基于vdP系統(tǒng)的研究產(chǎn)生了眾多的學(xué)術(shù)成果。

目前，學(xué)者們也針對(duì)含分?jǐn)?shù)階導(dǎo)數(shù)的vdP 系統(tǒng)進(jìn)行了研究。Barbosa 等［7］通過(guò)分析分?jǐn)?shù)階vdP 系統(tǒng)的相圖、分岔圖以及借譜分析等，較為詳細(xì)地研究了分?jǐn)?shù)階vdP 振動(dòng)系統(tǒng)在時(shí)域和頻域的動(dòng)力學(xué)行為；Tavazoei 研究了分?jǐn)?shù)階vdP 系統(tǒng)的極限環(huán)與初始條件的關(guān)系，并給出了存在極限環(huán)的解析判據(jù)；Chen 等［8］用相軌跡圖和龐加萊映射研究了分?jǐn)?shù)階vdP系統(tǒng)的混沌動(dòng)力學(xué)行為，發(fā)現(xiàn)分?jǐn)?shù)階階次會(huì)直接導(dǎo)致系統(tǒng)進(jìn)入混沌狀態(tài)。毛北行等［9］基于Lyapunov穩(wěn)定性理論給出了分?jǐn)?shù)階Duffling-van der Pol系統(tǒng)混沌同步問(wèn)題的兩個(gè)充分性條件。韋鵬等［10］以含分?jǐn)?shù)階vdP 系統(tǒng)為對(duì)象，研究了其超諧共振時(shí)的動(dòng)力學(xué)特性，得到了超諧共振周期響應(yīng)的穩(wěn)定性判斷準(zhǔn)則，同時(shí)提出等效非線性阻尼和非線性穩(wěn)定性條件參數(shù)的概念。

上述是從定性角度來(lái)描述和分析分?jǐn)?shù)階vdP系統(tǒng)的各種性質(zhì)，定量的研究則需要獲取系統(tǒng)較為精確的數(shù)值解或解析/半解析解［11-12］。和整數(shù)階非線性系統(tǒng)不同的是，分?jǐn)?shù)階微分算子的引入使得系統(tǒng)的解析/半解析解的求解難度急劇上升，相關(guān)的成果非常少。所以，現(xiàn)階段我們更加關(guān)注分?jǐn)?shù)階vdP 系統(tǒng)的數(shù)值求解。對(duì)于整數(shù)階非線性系統(tǒng)，數(shù)值計(jì)算法包括有限差分法［13］、Runge-Kutta 法［14］、精細(xì)積分法［15］、預(yù)估校正法［16］等。對(duì)于分?jǐn)?shù)階微分系統(tǒng)，即便是數(shù)值方法也尚在發(fā)展當(dāng)中。

Newmark-β法是一種非常成熟的數(shù)值方法，它是基于ti時(shí)刻的運(yùn)動(dòng)狀態(tài)ui、u˙i、u¨i，然后進(jìn)一步計(jì)算ti+1時(shí)刻狀態(tài)，逐步獲得結(jié)構(gòu)動(dòng)力響應(yīng)的逐步積分法。Newmark-β法通過(guò)簡(jiǎn)單的參數(shù)設(shè)置便可以無(wú)條件穩(wěn)定，且編程簡(jiǎn)單易于實(shí)現(xiàn)，這些特點(diǎn)使其在整數(shù)階微分系統(tǒng)中得到了廣泛應(yīng)用。在本文中，我們將基于Adams離散，提出一種針對(duì)Caputo分?jǐn)?shù)階導(dǎo)數(shù)的離散格式，然后進(jìn)一步基于Newmark-β法構(gòu)造出完整的逐步迭代格式，最后通過(guò)Newton-Raphson迭代［17］來(lái)進(jìn)一步求得分?jǐn)?shù)階vdP系統(tǒng)的響應(yīng)。

1 基于Adams離散和Newmark-β法的迭代格式

1.1 對(duì)分?jǐn)?shù)階算子Dαx進(jìn)行Adams離散

對(duì)于如下形式的分?jǐn)?shù)階vdP方程［18］

其中μ為阻尼系數(shù)，Dαx表示分?jǐn)?shù)階微分算子，α為對(duì)應(yīng)的分?jǐn)?shù)階階次。分?jǐn)?shù)階導(dǎo)數(shù)Dαx通常有Riemann-Liouvill 定義、Caputo 定義和Grünwald-Letnikov 定義［19-21］。本文中，分?jǐn)?shù)階算子Dαx采用Caputo 定義［22］。當(dāng)分?jǐn)?shù)階次0 ＜α＜1時(shí)，有［23］

其中求積系數(shù)di，n為

1.2 結(jié)合Newmark-β法的逐步迭代格式

當(dāng)分?jǐn)?shù)階階次α= 0.5時(shí)，分?jǐn)?shù)階vdP系統(tǒng)可以被完整地表示為

式中β和γ為Newmark-β法的控制參數(shù)，Δt=tn-tn-1為時(shí)間步長(zhǎng)。進(jìn)一步，將方程（10）中的狀態(tài)變量代入到方程（9）中，可以得到系統(tǒng)（8）的完整迭代格式

可以看出，系數(shù)A、B、C、D是通過(guò)β、γ、tn、α等表示的常數(shù)，則方程（12）是關(guān)于tn時(shí)刻未知位移xn的非線性代數(shù)方程。

1.3 通過(guò)Newton-Raphson迭代法求解位移方程

方程（12）可以通過(guò)Newton-Raphson 迭代法來(lái)求解。Newton-Raphson 迭代本質(zhì)上是將非線性方程f(x) = 0 進(jìn)行泰勒展開(kāi)，并保留一階項(xiàng)得到線性化方程來(lái)尋求方程f(x) = 0 的根。對(duì)于方程（12），有f(xn) =Axn+Bxn2+Cxn3+D= 0，則f(xn)的泰勒展開(kāi)為

對(duì)于足夠小的δ，忽略二階及以上的項(xiàng)，可得

需要注意的是，通過(guò)Newton-Raphson 迭代求解xn，需要選取一個(gè)合適的初值［26］。通常選取上一時(shí)刻tn-1的位移xn-1作為下一時(shí)刻的迭代初值。獲得位移xn后，將其代入方程（10），得到系統(tǒng)的速度和加速度。自此，系統(tǒng)在tn時(shí)刻所有的狀態(tài)變量全部獲得。重復(fù)多次，就可獲得分?jǐn)?shù)階vdP 系統(tǒng)完整的時(shí)程響應(yīng)。

2 數(shù)值算例

2.1 0 < α < 1的分?jǐn)?shù)階vdP系統(tǒng)

首先假定分?jǐn)?shù)階階次α= 0.5，即以系統(tǒng)（8）為例，取參數(shù)μ= 0.8.Newmark-β法中的參數(shù)取值為β=0.25、γ= 0.5. 由初始條件x0= 0，x˙0= 1可以計(jì)算出x¨0，進(jìn)一步計(jì)算獲得t1的系數(shù)A、B、C和D，然后逐步獲得系統(tǒng)的數(shù)值響應(yīng)。圖1為α= 0.5時(shí)系統(tǒng)（8）的時(shí)程響應(yīng)，圖2是對(duì)應(yīng)的相圖。

圖1 系統(tǒng)（8）的時(shí)程響應(yīng)Fig.1 The time histories of system(8)

圖2 系統(tǒng)（8）的相圖Fig.2 Phase diagram of system(8)

由圖1可以看出，在最初的15 s內(nèi)，系統(tǒng)出現(xiàn)了瞬態(tài)響應(yīng)，但是逐漸地系統(tǒng)進(jìn)入穩(wěn)態(tài)響應(yīng)，對(duì)應(yīng)的相圖則出現(xiàn)了穩(wěn)定的極限環(huán)。這說(shuō)明本文所提的數(shù)值方法不會(huì)出現(xiàn)失穩(wěn)現(xiàn)象，可以穩(wěn)定獲得系統(tǒng)的數(shù)值響應(yīng)。當(dāng)分?jǐn)?shù)階階次α→1時(shí)，系統(tǒng)（8）將會(huì)退化為對(duì)應(yīng)的整數(shù)階vdP系統(tǒng)，即

通過(guò)整數(shù)階的數(shù)值方法驗(yàn)證本文所提方法的有效性。當(dāng)α= 1 時(shí)，四階Runge-Kutta（RK）法和本文所提方法獲得的時(shí)程曲線和相圖，如圖3所示。從圖中可以看出，無(wú)論是時(shí)程曲線還是相圖，兩種方法獲取的結(jié)果吻合很好。這證明本文的方法對(duì)于整數(shù)階微分系統(tǒng)仍舊是有效的。

圖3 α = 1時(shí)，四階RK法和本文算法的時(shí)程曲線（a）和相圖（b）Fig.3 The time histories(a)and phase(b)obtained by fourth-order RK method and the algorithm of this paper when α = 1

2.2 1 < α < 2的分?jǐn)?shù)階vdP系統(tǒng)

不失一般性，可以假定系統(tǒng)（1）中的分?jǐn)?shù)階階次α= 1.5，μ= 0.8. 則系統(tǒng)（1）有如下形式

與α= 0.5的討論類似，結(jié)合Adams離散和Newmark-β法，可以得到系統(tǒng)（17）的逐步迭代格式

與上一個(gè)算例類似，在tn時(shí)刻系數(shù)A′，B′，C′和D′也是β、γ、tn、α表示的常數(shù)，則方程（19）是關(guān)于未知位移xn的非線性代數(shù)方程。將xn代入方程（10）得到系統(tǒng)的速度與加速度。

本算例中，Newmark-β法中的參數(shù)取值與2.1節(jié)相同，時(shí)間步長(zhǎng)取Δt= 0.1.α= 1.5時(shí)，系統(tǒng)（17）的時(shí)程響應(yīng)和相圖，如圖4和圖5所示。

圖4 系統(tǒng)（17）的時(shí)程曲線Fig.4 The time histories of system(17)

圖5 系統(tǒng)（17）的相圖Fig.5 Phase diagram of system(17)

從圖4所示的時(shí)間歷程可以看出，由于初始狀態(tài)的影響，系統(tǒng)出現(xiàn)了瞬態(tài)響應(yīng)，經(jīng)過(guò)一段時(shí)間后才進(jìn)入穩(wěn)態(tài)。且，由于本算例中α= 1.5，其時(shí)程曲線與α= 0.5時(shí)明顯不同，相圖也有明顯的區(qū)別。此外，α=1.5對(duì)應(yīng)的響應(yīng)頻率比α= 0.5時(shí)要小得多。同樣，當(dāng)分?jǐn)?shù)階階次α→2時(shí)，系統(tǒng)（17）也會(huì)退化為對(duì)應(yīng)的整數(shù)階vdP系統(tǒng)。

α= 2時(shí)，四階RK 法和本文所提方法獲得的時(shí)程曲線和相圖，如圖6所示。從圖6可以看到，本文方法和四階RK 法的結(jié)果吻合非常好。對(duì)比圖4和圖6（a），可以發(fā)現(xiàn)α= 1.5和α= 2時(shí)vdP 系統(tǒng)的響應(yīng)完全不同；對(duì)于α= 1.5 的分?jǐn)?shù)階vdP 系統(tǒng)，系統(tǒng)的速度和加速度都存在明顯的尖峰。對(duì)比圖5 和圖6（b），可以發(fā)現(xiàn)α= 1.5對(duì)應(yīng)的相圖更接近于一個(gè)方形，而α= 2時(shí)的相圖則是一個(gè)橢圓。同樣的，本文所提的離散格式和數(shù)值方法也可以很好地求解1 ＜α＜2的分?jǐn)?shù)階vdP系統(tǒng)。

圖6 α = 2時(shí)，四階RK法和本文算法的時(shí)程曲線（a）和相圖（b）Fig.6 The time histories and phase obtained by fourth-order RK method （a）and the algorithm of this paper （b）when α = 2

3 結(jié) 論

本文基于Adams離散和Newmark-β法，提出了一種求解分?jǐn)?shù)階van der Pol系統(tǒng)的逐步迭代格式，并分別討論了0 ＜α＜1 和1 ＜α＜2 對(duì)應(yīng)的vdP 系統(tǒng)的數(shù)值求解。結(jié)果表明，本文所提方法可以穩(wěn)定獲得分?jǐn)?shù)階微分系統(tǒng)的數(shù)值解。此外，我們還討論了α→1和α→2時(shí)，即分?jǐn)?shù)階vdP 系統(tǒng)退化為整數(shù)階vdP 系統(tǒng)時(shí)的情形。結(jié)果表明，本文所提算法和四階RK法獲得的結(jié)果吻合非常好。原則上，本文所提的迭代格式適用于其他類型的分?jǐn)?shù)階微分系統(tǒng)。