李亞杰，吳志強(qiáng)，蘭奇遜，郝 穎，張祥云

(1.河南城建學(xué)院 數(shù)理學(xué)院，河南 平頂山 467036；2.天津大學(xué) 機(jī)械學(xué)院，天津 300072；3.天津市非線性動力學(xué)與混沌控制重點(diǎn)實(shí)驗室，天津 300072；4.燕山大學(xué) 建筑工程與力學(xué)學(xué)院，河北 秦皇島 066004;5.天津中德應(yīng)用技術(shù)大學(xué) 基礎(chǔ)課部，天津 300350)

分?jǐn)?shù)階導(dǎo)數(shù)是對整數(shù)階導(dǎo)數(shù)的推廣，將函數(shù)的求導(dǎo)階次從傳統(tǒng)的整數(shù)階擴(kuò)展到了非整數(shù)階的情形，至今已有三百多年的發(fā)展歷史。由于整數(shù)階導(dǎo)數(shù)定義的局限性，其不能表達(dá)黏性物質(zhì)的記憶特性，而分?jǐn)?shù)階導(dǎo)數(shù)定義含有卷積部分，能很好地表達(dá)記憶效果，表示出隨時間累積的效應(yīng)，因而與傳統(tǒng)的整數(shù)階微積分系統(tǒng)相比，含分?jǐn)?shù)階微積分的系統(tǒng)更具有優(yōu)越性，是描述記憶特征的恰當(dāng)數(shù)學(xué)工具[1-4]。近年來，關(guān)于分?jǐn)?shù)階微積分及其應(yīng)用的研究引起了不同領(lǐng)域科研人員的廣泛關(guān)注，已成為研究反常擴(kuò)散、多孔介質(zhì)力學(xué)、非牛頓流體力學(xué)、黏彈性力學(xué)、軟物質(zhì)物理等學(xué)科領(lǐng)域的有力數(shù)學(xué)工具。由于分?jǐn)?shù)階導(dǎo)數(shù)能更準(zhǔn)確的描述各種反應(yīng)變化過程，很多問題可以用分?jǐn)?shù)階微分方程來更好地描述[5-9]，因而研究分?jǐn)?shù)階微分方程中的典型力學(xué)特性和分?jǐn)?shù)階參數(shù)對系統(tǒng)的影響也是十分必要和有著重要意義的。

近年來相關(guān)學(xué)者研究了在不同噪聲激勵下非線性多穩(wěn)態(tài)系統(tǒng)的動力學(xué)行為，取得了豐碩的成果:整數(shù)階系統(tǒng)方面，文獻(xiàn)[10-14]研究了Duffing-van der Pol振子在Lévy噪聲、色噪聲、諧和與隨機(jī)噪聲聯(lián)合激勵下系統(tǒng)響應(yīng)的穩(wěn)態(tài)概率密度問題，通過分析系統(tǒng)穩(wěn)態(tài)概率密度函數(shù)性質(zhì)的變化，討論了噪聲振子的隨機(jī)P分岔現(xiàn)象，得到了系統(tǒng)雙峰概率密度函數(shù)的解析表達(dá)式，結(jié)果表明系統(tǒng)參數(shù)與噪聲強(qiáng)度均能誘導(dǎo)系統(tǒng)發(fā)生隨機(jī)P分岔行為。Wu等[15]利用奇異性理論研究了一、三穩(wěn)態(tài)van der Pol-Duffing振子在乘性色噪聲激勵下系統(tǒng)穩(wěn)態(tài)響應(yīng)的隨機(jī)P分岔行為，得到了系統(tǒng)幅值穩(wěn)態(tài)概率密度函數(shù)的表達(dá)式，分析了系統(tǒng)參數(shù)及噪聲強(qiáng)度對系統(tǒng)隨機(jī)P分岔的影響；分?jǐn)?shù)階系統(tǒng)方面，Chen等[16]研究了諧和與寬帶噪聲聯(lián)合激勵下含分?jǐn)?shù)階導(dǎo)數(shù)阻尼Duffing振子的響應(yīng)，發(fā)現(xiàn)分?jǐn)?shù)階導(dǎo)數(shù)階次的變化可以導(dǎo)致系統(tǒng)發(fā)生隨機(jī)P分岔。Huang等[17]討論了分?jǐn)?shù)階單自由度非線性系統(tǒng)在高斯白噪聲激勵下的響應(yīng)和穩(wěn)態(tài)概率密度函數(shù)。Yang等[18]應(yīng)用Zhuravlev非光滑變換和隨機(jī)平均法研究了高斯白噪聲激勵下含Caputo型分?jǐn)?shù)階阻尼的非線性碰撞振動系統(tǒng)的隨機(jī)分岔。Li等[19]研究了在加性及乘性色噪聲激勵下含分?jǐn)?shù)階導(dǎo)數(shù)的廣義Duffing-van der Pol系統(tǒng)的雙穩(wěn)態(tài)隨機(jī)分岔現(xiàn)象，發(fā)現(xiàn)線性阻尼系數(shù)、分?jǐn)?shù)階導(dǎo)數(shù)階次及噪聲強(qiáng)度的變化均可導(dǎo)致系統(tǒng)發(fā)生隨機(jī)P分岔行為。

由于分?jǐn)?shù)階導(dǎo)數(shù)的復(fù)雜性，其分析方法也變得更加困難，系統(tǒng)參數(shù)對振動特性的研究多限于定性分析，未能找出參數(shù)影響的臨界條件，影響此類系統(tǒng)的分析和設(shè)計，且目前關(guān)于三穩(wěn)態(tài)分?jǐn)?shù)階系統(tǒng)的隨機(jī)P分岔問題，還未見報道。針對以上情況，本文以聯(lián)合噪聲激勵下廣義van der Pol方程的非線性振動為例，用奇異性方法求得分?jǐn)?shù)階系統(tǒng)的轉(zhuǎn)遷集曲線，得到了系統(tǒng)發(fā)生隨機(jī)P分岔的臨界參數(shù)條件，并分析了轉(zhuǎn)遷集參數(shù)平面內(nèi)各區(qū)域中系統(tǒng)幅值穩(wěn)態(tài)概率密度的曲線類型。通過Monte Carlo模擬的方法，將所得數(shù)值結(jié)果與文中求得的解析結(jié)果進(jìn)行了比對，可以看出數(shù)值解與解析解的吻合情況較好，從而驗證了本文理論分析的正確性。

1 等價系統(tǒng)的推導(dǎo)

關(guān)于分?jǐn)?shù)階導(dǎo)數(shù)的定義較多，常用的有Riemann-Liouville導(dǎo)數(shù)和Caputo導(dǎo)數(shù)，但Riemann-Liouville導(dǎo)數(shù)對應(yīng)的初始條件沒有物理意義，而Caputo導(dǎo)數(shù)所描述系統(tǒng)的初始條件有著清楚的物理意義且形式上和整數(shù)階微分方程的初始條件類似，故在本文中我們采用Caputo導(dǎo)數(shù)。定義在區(qū)間[a,b]上的函數(shù)x(t)的p階Caputo導(dǎo)數(shù)定義為

(1)

式中：m-1

(2)

式中，m-1

本文研究聯(lián)合噪聲激勵下含分?jǐn)?shù)階導(dǎo)數(shù)項的廣義van der Pol系統(tǒng)

ξ1(t)+x(t)ξ2(t)

(3)

E[ξk(t)]=0,E[ξk(t)ξl(t-τ)]=2Dklδ(τ)

(4)

由文獻(xiàn)[20-23]，分?jǐn)?shù)階導(dǎo)數(shù)項可以等效為阻尼力和回復(fù)力的線性組合，故引入以下的等效系統(tǒng)

[K(p,w)+w2]x=ξ1(t)+x(t)ξ2(t)

(5)

式中，C(p,w),K(p,w)分別為等效阻尼力和回復(fù)力的系數(shù)。

式(3)和式(5)的誤差為

(6)

由均方誤差最小的必要條件[24]

(7)

將式(6)代入式(7)中可得

(8)

式(3)的解可設(shè)為如下形式

x(t)=a(t)cosφ(t)

(9)

式中，φ(t)=wt+θ。

則有

(10)

將式(9)、式(10)代入式(8)，并關(guān)于φ進(jìn)行積分平均可得

(11)

從而，等價式(5)可以表示為以下形式

(12)

其中

(13)

2 幅值的穩(wěn)態(tài)概率密度

對于式(12)，取定線性及非線性阻尼的系數(shù)分別為ε=-0.1,α1=1.51,α2=2.85,α3=1.693,α4=0.312,w=1，為方便討論參數(shù)影響，圖1給出了當(dāng)D1=D2=0時，對應(yīng)的確定性系統(tǒng)極限環(huán)幅值隨分?jǐn)?shù)階導(dǎo)數(shù)階次p變化的分岔曲線。

圖1 確定性系統(tǒng)分岔圖Fig.1 Bifurcation diagram of the deterministic system

可以看出，當(dāng)p在區(qū)間[0.063 7,0.137 7]變化時，系統(tǒng)有2個吸引子：平衡點(diǎn)和大極限環(huán)，如圖2(a)所示；當(dāng)p在區(qū)間[0.137 7,0.140 6]變化時,系統(tǒng)有3個吸引子：平衡點(diǎn)、小極限環(huán)和大極限環(huán)，如圖2(b)所示；當(dāng)p在區(qū)間[0.140 6,0.152 2]變化時，系統(tǒng)有2個吸引子：平衡點(diǎn)和小極限環(huán)，如圖2(c)所示。

圖2 相圖Fig.2 Phase portraits

為求解式(12)響應(yīng)幅值的穩(wěn)態(tài)概率密度函數(shù)，引入如下變換[25]

(14)

式中,a(t),θ(t)分別為系統(tǒng)的幅值和初始相位，均為隨機(jī)過程。

將式(14)代入式(12)，由確定性平均法可得

(15)

其中，

(16)

式(15)可以被看作Stratonovich隨機(jī)微分方程，通過加入相應(yīng)的Wong-Zakai修正項，可以將其轉(zhuǎn)化為如下的It隨機(jī)微分方程

(17)

式中，Bk(t)(k=1,2)為標(biāo)準(zhǔn)的維納過程，且

(18)

對式(17)關(guān)于Φ進(jìn)行隨機(jī)平均[26]，可得如下的平均It方程

(19)

式中，B1(t)與B2(t)為2個相互獨(dú)立的單位Wiener過程，且

(20)

式(20)表明，振動幅值a(t)的平均It微分方程與θ(t)是相互獨(dú)立的，故a(t)為一維的隨機(jī)過程。振動幅值a(t)的FPK方程可表示為

(21)

其對應(yīng)的邊界條件為

(22)

基于以上邊界條件，式(21)的穩(wěn)態(tài)解即為系統(tǒng)幅值的穩(wěn)態(tài)概率密度函數(shù)

其中，C為歸一化常數(shù)，滿足

將式(20)代入式(23)，可得振動幅值a穩(wěn)態(tài)概率密度函數(shù)的具體表達(dá)式

(25)

其中，

(26)

3 隨機(jī)P分岔

隨機(jī)P分岔是指概率密度函數(shù)曲線峰值數(shù)目的變化，為得到P分岔的臨界參數(shù)條件，以下從奇異性分析的角度來分析參數(shù)變化對系統(tǒng)隨機(jī)P分岔的影響。

為方便起見，p(a)可表示為

p(a)=4CR(a,D1,D2,ε,w,p,α1,α2,α3,α4)exp[Q(a,D1,D2,ε,w,p,α1,α2,α3,α4)]

(27)

其中,

(28)

根據(jù)奇異性理論[27],概率密度函數(shù)需滿足如下的2個條件

(29)

將式(27)代入式(29)，可得如下條件

H={R′+RQ′=0,R″+2R′Q′+RQ″+RQ′2=0}

(30)

式中，H為概率密度函數(shù)曲線峰值數(shù)目變化的條件。

由于參數(shù)關(guān)系在三維曲面中不容易刻畫和顯示，這里我們只給出轉(zhuǎn)遷集的二維截面來表示噪聲強(qiáng)度和分?jǐn)?shù)階導(dǎo)數(shù)階次的影響。

4 分?jǐn)?shù)階導(dǎo)數(shù)階次p對系統(tǒng)隨機(jī)P分岔行為的影響

根據(jù)圖1中確定性吸引子的分布，不失一般性，我們分別在系統(tǒng)單穩(wěn)態(tài)、雙穩(wěn)態(tài)及三穩(wěn)態(tài)區(qū)間內(nèi)取定分?jǐn)?shù)階導(dǎo)數(shù)階次p的值，并根據(jù)式(28)及式(30)計算聯(lián)合噪聲激勵下系統(tǒng)對應(yīng)的轉(zhuǎn)遷集。由于當(dāng)分?jǐn)?shù)階導(dǎo)數(shù)階次p在區(qū)間[0，0.063 6]中取值時，系統(tǒng)轉(zhuǎn)遷集為空集，因此我們僅在區(qū)間[0.063 7，0.137 7]中取(a)p=0.1，區(qū)間[0.137 8，0.140 5]中取(b)p=0.14，區(qū)間[0.140 6，0.152 2]中取(c)p=0.15及區(qū)間[0.152 3，0.18]中取(d)p=0.17，具體轉(zhuǎn)遷集分別如圖3所示。

為便于比較分析，以上在轉(zhuǎn)遷集各區(qū)域中均用數(shù)字進(jìn)行了標(biāo)注，數(shù)字標(biāo)號相同的區(qū)域表示系統(tǒng)幅值穩(wěn)態(tài)概率密度曲線是定性相同的，方便起見，我們僅對轉(zhuǎn)遷集圖3(d)各子區(qū)域中系統(tǒng)幅值的穩(wěn)態(tài)概率密度曲線進(jìn)行分析，將所得解析結(jié)果與式(1)的Monte Carlo數(shù)值模擬結(jié)果進(jìn)行對比，具體結(jié)果如圖4所示。

圖3 不同分?jǐn)?shù)階導(dǎo)數(shù)階次p下的轉(zhuǎn)遷集(以D2和D1為開折參數(shù))Fig.3 Transition sets under different values of fractional derivative’s order p(taking D2 and D1 as the unfolding parameters)

圖4 圖3(d)不同子區(qū)域中幅值的概率密度函數(shù)p(a)Fig.4 PDF of amplitude p(a)in different sub-areas of Fig.3(d)(taking D2 and D1 as the unfolding parameters)

由圖3(d)可見：概率密度曲線出現(xiàn)多峰的參數(shù)區(qū)域由2個近似三角形區(qū)域圍成，特別地，二者重合的區(qū)域4構(gòu)成了系統(tǒng)幅值穩(wěn)態(tài)概率密度曲線的三峰區(qū)域。當(dāng)參數(shù)在(D2,D1)參數(shù)平面區(qū)域1中取值時，概率密度曲線在離原點(diǎn)較遠(yuǎn)處有一明顯峰值，如圖4(a)所示；當(dāng)參數(shù)在區(qū)域2中取值時，概率密度曲線在離原點(diǎn)較遠(yuǎn)處有2個區(qū)分不明顯的峰值，系統(tǒng)同時存在大、小極限環(huán)，如圖4(b)所示；當(dāng)參數(shù)在區(qū)域3中取值時，概率密度曲線在離原點(diǎn)較遠(yuǎn)處仍有一明顯峰值，但在原點(diǎn)附近概率明顯不為0，系統(tǒng)同時存在平衡點(diǎn)與大極限環(huán)，如圖4(c)所示；當(dāng)參數(shù)在區(qū)域4中取值時，概率密度曲線在原點(diǎn)以外還存在2個峰值，系統(tǒng)表現(xiàn)為平衡點(diǎn)與大、小極限環(huán)共存，如圖4(d)所示；當(dāng)參數(shù)在區(qū)域5中取值時，概率密度曲線與區(qū)域3中定性相同，概率密度曲線偏離原點(diǎn)的峰值所對應(yīng)的幅值a要小于圖4(c)中峰值所對應(yīng)的a值，系統(tǒng)同時存在平衡點(diǎn)與小極限環(huán)，如圖4(e)所示。

以上分析結(jié)果表明，分?jǐn)?shù)階導(dǎo)數(shù)階次p的取值不同，其對應(yīng)的聯(lián)合噪聲下的轉(zhuǎn)遷集圖形亦不相同，這意味著，分?jǐn)?shù)階導(dǎo)數(shù)項的加入能對系統(tǒng)幅值穩(wěn)態(tài)概率密度起到調(diào)控的作用，可以從分?jǐn)?shù)階導(dǎo)數(shù)階次p、加性噪聲強(qiáng)度D1及乘性噪聲強(qiáng)度D23個方面來控制聯(lián)合噪聲作用下系統(tǒng)的運(yùn)動形式。圖3中任意相鄰的2個區(qū)域所對應(yīng)系統(tǒng)幅值的穩(wěn)態(tài)概率密度曲線是定性不同的，參數(shù)(D2,D1)的取值越過圖3中任意一條線，系統(tǒng)均將發(fā)生隨機(jī)P分岔行為，故轉(zhuǎn)遷集曲線即為系統(tǒng)發(fā)生隨機(jī)P分岔的臨界參數(shù)條件，且圖4中解析結(jié)果與Monte Carlo數(shù)值模擬結(jié)果吻合較好，這也進(jìn)一步驗證了本文理論分析的正確性。

5 結(jié) 論

本文對聯(lián)合噪聲激勵下含分?jǐn)?shù)階導(dǎo)數(shù)項的廣義van der Pol系統(tǒng)的三穩(wěn)態(tài)隨機(jī)P分岔現(xiàn)象進(jìn)行了研究。根據(jù)均方誤差最小原則，將原系統(tǒng)轉(zhuǎn)化為了等價的整數(shù)階系統(tǒng)，運(yùn)用隨機(jī)平均法得到了系統(tǒng)幅值的穩(wěn)態(tài)概率密度函數(shù)。并利用奇異性理論，得到了系統(tǒng)發(fā)生隨機(jī)P分岔的臨界參數(shù)條件，通過選擇合適的開折參數(shù)，可使系統(tǒng)響應(yīng)維持在單穩(wěn)態(tài)或平衡點(diǎn)附近小幅振動，可避免系統(tǒng)發(fā)生大幅振蕩或非線性跳躍現(xiàn)象造成失穩(wěn)，可為相關(guān)系統(tǒng)設(shè)計提供理論指導(dǎo)。對原系統(tǒng)進(jìn)行Monte Carlo數(shù)值模擬驗證了所得理論結(jié)果的正確性，分析得出：分?jǐn)?shù)階導(dǎo)數(shù)階次p及噪聲強(qiáng)度D2,D1均可引起系統(tǒng)發(fā)生隨機(jī)P分岔行為，通過選取相應(yīng)的參數(shù)p及(D2,D1)可以實(shí)現(xiàn)系統(tǒng)幅值穩(wěn)態(tài)概率密度曲線峰值1～3的變化。