李亞杰，吳志強，蘭奇遜，郝 穎，張祥云

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

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

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

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

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

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

(1)

式中：m-1

(2)

式中，m-1

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

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

(3)

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

(4)

由文獻[20-23]，分數階導數項可以等效為阻尼力和回復力的線性組合，故引入以下的等效系統(tǒng)

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

(5)

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

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

(6)

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

(7)

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

(8)

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

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

(9)

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

則有

(10)

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

(11)

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

(12)

其中

(13)

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

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

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

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

圖2 相圖Fig.2 Phase portraits

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

(14)

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

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

(15)

其中，

(16)

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

(17)

式中，Bk(t)(k=1,2)為標準的維納過程，且

(18)

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

(19)

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

(20)

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

(21)

其對應的邊界條件為

(22)

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

其中，C為歸一化常數，滿足

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

(25)

其中，

(26)

3 隨機P分岔

隨機P分岔是指概率密度函數曲線峰值數目的變化，為得到P分岔的臨界參數條件，以下從奇異性分析的角度來分析參數變化對系統(tǒng)隨機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)

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

(29)

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

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

(30)

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

由于參數關系在三維曲面中不容易刻畫和顯示，這里我們只給出轉遷集的二維截面來表示噪聲強度和分數階導數階次的影響。

4 分數階導數階次p對系統(tǒng)隨機P分岔行為的影響

根據圖1中確定性吸引子的分布，不失一般性，我們分別在系統(tǒng)單穩(wěn)態(tài)、雙穩(wěn)態(tài)及三穩(wěn)態(tài)區(qū)間內取定分數階導數階次p的值，并根據式(28)及式(30)計算聯(lián)合噪聲激勵下系統(tǒng)對應的轉遷集。由于當分數階導數階次p在區(qū)間[0，0.063 6]中取值時，系統(tǒng)轉遷集為空集，因此我們僅在區(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，具體轉遷集分別如圖3所示。

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

圖3 不同分數階導數階次p下的轉遷集(以D2和D1為開折參數)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ū)域中幅值的概率密度函數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)可見：概率密度曲線出現多峰的參數區(qū)域由2個近似三角形區(qū)域圍成，特別地，二者重合的區(qū)域4構成了系統(tǒng)幅值穩(wěn)態(tài)概率密度曲線的三峰區(qū)域。當參數在(D2,D1)參數平面區(qū)域1中取值時，概率密度曲線在離原點較遠處有一明顯峰值，如圖4(a)所示；當參數在區(qū)域2中取值時，概率密度曲線在離原點較遠處有2個區(qū)分不明顯的峰值，系統(tǒng)同時存在大、小極限環(huán)，如圖4(b)所示；當參數在區(qū)域3中取值時，概率密度曲線在離原點較遠處仍有一明顯峰值，但在原點附近概率明顯不為0，系統(tǒng)同時存在平衡點與大極限環(huán)，如圖4(c)所示；當參數在區(qū)域4中取值時，概率密度曲線在原點以外還存在2個峰值，系統(tǒng)表現為平衡點與大、小極限環(huán)共存，如圖4(d)所示；當參數在區(qū)域5中取值時，概率密度曲線與區(qū)域3中定性相同，概率密度曲線偏離原點的峰值所對應的幅值a要小于圖4(c)中峰值所對應的a值，系統(tǒng)同時存在平衡點與小極限環(huán)，如圖4(e)所示。

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

5 結 論

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