秦 浩，溫少芳，申永軍，邢海軍，王 軍

(1.石家莊鐵道大學 交通運輸學院 省部共建交通工程結構力學行為與系統(tǒng)安全國家重點實驗室，石家莊 050043；2.石家莊鐵道大學 機械工程學院，石家莊 050043)

分數階微積分已經存在了300多年，但是目前還是一個比較新鮮的概念。早在整數階微積分創(chuàng)立初期，法國學家Hospital和德國數學家Leibniz就提出了分數階的概念[1-3]，但是由于缺乏實際應用背景的支撐等諸多方面原因使得它長期以來沒有得到研究和發(fā)展。隨著自然科學、社會科學的發(fā)展和工程的實際需要，各種復雜的系統(tǒng)和計算機技術引入到工程領域，分數階微積分的定義、特性和計算等才得以在工程領域應用。目前，分數階微積分在工程領域中的應用主要分為兩類：一類是分數階微積分的控制，能夠提高控制的效果和系統(tǒng)的魯棒性；另一類是利用分數階微積分建模，很多工程材料介于理想固體與理想流體之間，具有記憶特性，用分數階微積分來模擬更加準確，更能反應真實的本構關系[4-8]。

早期學者在解釋生活中的經典現(xiàn)象時提出了許多假說，其中大部分是基于整數階模型模擬自然現(xiàn)象，并且模擬后的結果能夠滿足當時的需求，因此人們忽視了分數階模型。隨著研究的深入，研究者們發(fā)現(xiàn)對于一些具有黏彈性特性的材料單純用整數階的模型無法滿足精度的需要，工程系統(tǒng)中的非光滑性、不連續(xù)性、參數激勵、間隙等，使得系統(tǒng)響應變得更加復雜，并且會降低系統(tǒng)的整體性能，甚至會導致系統(tǒng)出現(xiàn)分岔和混沌的動力學現(xiàn)象[9-14]。因此人們將整數階模型進一步優(yōu)化為分數階模型，用分數階微分方程來模擬非經典現(xiàn)象，能夠更精確地反應材料的本構關系，更好地滿足實際的需求。在控制方面，分數階相較于整數階的優(yōu)點在于能夠控制或改變系統(tǒng)的閉環(huán)控制特性，增強控制效果，提高系統(tǒng)的魯棒性。所以大量的學者開始研究分數階微分方程的典型力學特性和分數階微積分對整個動力系統(tǒng)力學特性的影響[15]。

隨著分數階微積分的發(fā)展，分數階微積分方程的求解成為學者們研究的一個熱點方向。分數階微積分的計算較為復雜，大部分的文獻中對分數階的處理采用數值計算的方法，除此之外還有定性和解析的方法。數值計算方法是處理分數階微積分的一個主流分支，應用較為廣泛，如Diethelm等[16]討論了一種可以應用于分數階微分方程的求數值解的Adams型預測-矯正法；Firdous等[17]提出了一種基于哈爾小波求解分數階微分方程的新的運算矩陣方法用于求解線性和非線性分數階微分方程的數值結果；Javidi等[18]提出了一種基于同倫微擾法和拉普拉斯變換求解分數階偏微分方程的數值方法。定性研究則側重分析方程解的數目和穩(wěn)定性的變化，比如Trigeassou等[19]利用李雅普諾夫方法通過頻率分布分數階微積分模型研究了分數階微分方程的穩(wěn)定性；Qi等[20]提出了分數階的階次大于1小于2時Caputo分數階導數的一個新定理和Caputo分數階系統(tǒng)的兩個新定理，方便了使用李雅普諾夫方法證明分數階非線性時變系統(tǒng)的穩(wěn)定性。解析研究是通過求得近似解析解對方程進行定量分析，比如Shen等[21-23]利用平均法分別研究了分數階Duffing振子和分數階van der Pol振子的近似解析解，提出了等效線性剛度和等效線性阻尼的概念，分析了分數階微分的參數對不同振子動力學特性的影響；姜源等[24]采用平均法研究了含分數階微分項的Duffing振子的超諧與亞諧聯(lián)合共振，得到了一階近似解析解，分析了分數階微分的參數對系統(tǒng)動力學特性的影響，并且提出了超、亞諧聯(lián)合共振時等效線性阻尼和等效線性剛度的概念；顧曉輝等[25]利用多尺度法得到了兩個諧波激勵作用下分數階Duffing振子的一類組合共振的一階近似解析解，并且研究分析了定常解的穩(wěn)定性；Tabejieu等[26]通過平均法研究了分數階階次對非線性Rayleigh梁振幅幅值的影響，證明了隨著分數階階次的增加，梁的共振振幅減小，同時利用Melnikov方法分析了混沌產生的必要條件。

關于Duffing振子產生混沌必要條件的問題，很多學者利用Melnikov方法對其進行了研究。如：Shen等[27]基于Melnikov方法分析了帶有位移延遲反饋和速度延遲反饋的整數階Duffing振子產生混沌的必要條件，分析了位移延遲反饋和速度延遲反饋分別影響混沌產生必要條件的規(guī)律。Sun等[28]利用Melnikov方法研究了具有位移延遲反饋和速度延遲反饋的整數階Duffing系統(tǒng)在諧波激勵下的混沌行為，導出了同宿分岔產生混沌的必要條件的解析解。Wen等[29]用Melnikov方法研究了強迫激勵下的整數階Duffing振子在位移延遲反饋和速度延遲反饋控制下的異宿分岔和混沌現(xiàn)象。張思進等[30]研究了具有立方非線性項和外部激勵項的二自由度非線性碰振系統(tǒng)的動力學特性，分別運用改進后的局部亞諧Melnikov方法和數值方法得到了二自由度碰振系統(tǒng)穩(wěn)定周期運動的存在條件。Abtahi[31]采用Melnikov方法研究了陀螺衛(wèi)星復雜自旋軌道動力學中的混沌現(xiàn)象。多數文獻主要針對整數階Duffing振子產生混沌的問題進行研究，關于分數階Duffing振子產生混沌必要條件多采用數值解，解析結果還未見相關文獻。下面重點關注應用Melnikov方法研究含有分數階微分項的Duffing振子的分岔和混沌現(xiàn)象的問題。Melnikov方法是一種預測系統(tǒng)混沌的一階近似方法，本文基于Melnikov方法提出一種分數階Duffing振子產生混沌必要條件的新的研究思路。首先將分數階微分項進行了處理，利用等效阻尼和等效剛度的概念，將分數階微分項等效成三角函數與指數函數的形式。接著利用Melnikov方法研究分數階Duffing振子的同宿軌分岔與混沌的產生條件。最后將解析結果與數值結果進行了比較，驗證了解析結果的準確度，并分析了分數階的階次和系數對分數階Duffing振子產生混沌的必要條件的影響。

1 分數階Duffing振子產生混沌的必要條件

研究如下的含分數階微分項的Duffing振子

(1)

(2)

式中，Γ(n)為Gamma函數，Γ(y+1)=yΓ(y)。

設式(1)的解滿足

x(t)=acos(ωt+θ)

(3a)

(3b)

利用Caputo定義求解分數階微分項，將式(3)代入式(2)，得到

(4)

引入s=t-u，則式(4)變?yōu)?/p>

(5)

引入兩個基本公式[32]

(6a)

(6b)

代入式(5)，得到

(7)

代入方程原參數，得到

(8)

本文主要研究式(1)經過長時間振動穩(wěn)定以后的動力學現(xiàn)象，所以在后面的研究中將高階項略掉，只取一階近似項，即

(9)

式(9)得到的結果與很多文獻中得到的結論相同，分數階微分項不僅起著阻尼的作用，同時也有剛度的作用，可以等價為等效線性阻尼和等效線性剛度。下面利用等效線性剛度與等效線性阻尼的概念，將式(9)代入式(1)變?yōu)?/p>

(10)

整理式(10)可得

(11)

引入c=εc1,λ=ελ1,f=εf1。

其中，ε為小參數，則式(11)變?yōu)?/p>

(12)

當ε=0時，則系統(tǒng)為未擾動系統(tǒng)

(13)

在式(13)中，存在一個同宿軌道滿足

(14)

(15)

對式(14)進行積分計算，得到

(16)

計算式(16)的定積分，整理得到

(17)

將式(12)轉化為

(18a)

其中，

(18b)

(18c)

(18d)

同宿軌道的位移為

(19a)

(19b)

建立Melnikov函數[33]

(20)

將式(20)中的定積分分為兩部分M1(t0)和M2(t0)，分別經過定積分計算，得到

(21a)

根據函數的奇偶性質，式(21a)變?yōu)?/p>

(21b)

(22)

聯(lián)立式(20)～式(22)，可以得到small意義下混沌的必要條件為

(23a)

代入原系統(tǒng)參數，式(23a)可變?yōu)?/p>

(23b)

當分數階微分項產生的等效剛度小于系統(tǒng)的固有剛度時，將式(23b)中的絕對值去掉，得到分數階Duffing系統(tǒng)產生混沌的必要條件為

(24)

2 數值仿真

為了驗證本文解析結果的正確性，下面將數值結果與本文解析結果進行對比。數值仿真利用Diethelm等研究中的預測矯正法對分數階Duffing系統(tǒng)式(1)進行計算，將計算結果中不穩(wěn)定的響應去掉，選取穩(wěn)態(tài)響應的部分，分析系統(tǒng)的穩(wěn)態(tài)響應情況。

本文中選取的基本參數c=0.2，k=1，α=1，ω=1.1，λ=0.1，分數階微分項的階次以p=0.5為例進行仿真驗證。系統(tǒng)初值首先取為“1，0，0，0”，計算時間取為10 000 s，舍去前5 000 s，確保數值的穩(wěn)定性。并使振幅f由0～0.5變化，選擇步長為0.005，得到系統(tǒng)式(1)產生分岔到混沌的過程示于圖1(a)中。從圖1(a)中可以看出系統(tǒng)產生混沌的臨界外激勵振幅fmin約為0.23，放大f=0.23周圍的分岔圖見圖1(b)，橫軸為外激勵振幅f從0.20～0.26，步長選取0.001。

從圖1(b)中可以看出混沌產生的臨界振幅fmin約為0.234。為了證明f=0.234為產生混沌的臨界點，下面f分別取為0.21，0.22，0.227，0.233，0.234時，分析其系統(tǒng)響應的特性如下：

(1) 當f=0.21時，系統(tǒng)的相圖和時間歷程圖如圖2(a)所示，此時的系統(tǒng)運動狀態(tài)為單周期運動。

(2) 當f=0.22時，系統(tǒng)的相圖和時間歷程圖如圖2(b)所示，此時的系統(tǒng)運動狀態(tài)為周期2的運動。

(3) 當f=0.227時，系統(tǒng)的相圖和時間歷程圖如圖2(c)所示，此時的系統(tǒng)運動狀態(tài)為周期4的運動。

(4) 當f=0.233時，系統(tǒng)的相圖和時間歷程圖如圖2(d)所示，此時的系統(tǒng)運動狀態(tài)為擬周期的運動。

圖2 相圖與時間歷程圖Fig.2 Phase portrait and time history

(5) 當f=0.234時，Duffing系統(tǒng)的相圖和時間歷程圖如圖3所示，從圖中可以看出系統(tǒng)此時處于混沌狀態(tài)。

由圖1～圖3的仿真結果，可以看出：當p=0.5時，外激勵振幅f=0.234是利用數值方法得到的分數階Duffing振子產生混沌的最小值。再根據本文解得的產生混沌必要條件的解析式(24)計算，得到系統(tǒng)產生混沌的必要條件的解析結果為fmin=0.22。對比數值解與解析解，差值為0.014。Melnikov方法本身是一種求解產生混沌必要條件的一階近似方法，所以數值解與解析解必然會存在著誤差。為了判斷差值是否在合理范圍內，將分數階微分項設為零，其他參數取值同上。使用龍格庫塔法求解此整數階Duffing系統(tǒng)產生混沌的fmin數值解為0.192，而利用Melnikov方法計算的整數階Duffing振子產生混沌的必要條件fmin為0.173，整數階情況下數值解與解析解的差值為0.019。兩個差值屬于同一量級，說明這種誤差主要源自Melnikov方法本身。而本文中提出的利用等效阻尼和等效剛度將分數階系統(tǒng)等效替換成整數階系統(tǒng)進行分析，精確度較高，方法可行。

圖1 分岔圖Fig.1 Bifurcation diagram

圖3 相圖與時間歷程圖Fig.3 Phase portrait and time history

在數值仿真過程中，還發(fā)現(xiàn)改變系統(tǒng)的初值時，出現(xiàn)了另外一條通過倍周期分岔通往混沌的路徑，如：當初值取為“-1，0，0，0”時，其他系統(tǒng)參數的取值同上，得到的系統(tǒng)分岔圖如圖4所示。當f分別為0.21，0.22，0.227，0.233和0.234時，系統(tǒng)的相圖和時間歷程圖如下：

(1) 當f=0.21時，系統(tǒng)的相圖和時間歷程圖如圖5(a)所示，此時的系統(tǒng)運動狀態(tài)為單周期運動。

(2) 當f=0.22時，系統(tǒng)的相圖和時間歷程圖如圖5(b)所示，此時的系統(tǒng)運動狀態(tài)為周期2的運動。

(3) 當f=0.227時，系統(tǒng)的相圖和時間歷程圖如圖5(c)所示，此時的系統(tǒng)運動狀態(tài)為周期4的運動。

(4) 當f=0.233時，系統(tǒng)的相圖和時間歷程圖如圖5(d)所示，此時的系統(tǒng)運動狀態(tài)為擬周期的運動。

(5) 當f=0.234時，分數階Duffing系統(tǒng)的相圖和時間歷程圖如圖6所示,與圖3(a)和圖3(b)運動狀態(tài)相同，均為混沌狀態(tài)。

通過圖1和圖4、圖2和圖5、圖3和圖6的兩兩對比分析可知：當f取值小于0.234時，選取上述不同初值時，分數階Duffing振子的運動特性都是相同的，相圖方向相反，分別為圖3和圖6中的混沌相圖的左右半支。但是不同初值時，產生混沌的臨界位置并沒有改變，fmin仍為0.234。所以，可以證明改變系統(tǒng)的初值不會改變系統(tǒng)產生混沌的必要條件,但是取不同初值分數階Duffing振子會通過不同的倍周期分岔路徑到達混沌。主要原因是該分數階Duffing振子具有雙穩(wěn)態(tài)特性，在文中參數條件下，外激勵f取值小于圖4所示的分岔點f1處時，系統(tǒng)處于單周期運動，存在兩個穩(wěn)態(tài)解。取不同初值時，將進入不同的吸引域從而形成不同的穩(wěn)態(tài)解。從兩個穩(wěn)態(tài)解出發(fā)，隨著外激勵參數的變化都能通過倍周期分岔到達混沌的狀態(tài)。

圖4 p=0.5的分岔圖Fig.4 Bifurcation diagram for p=0.5

圖5 相圖與時間歷程圖Fig.5 Phase portrait and time history

圖6 相圖與時間歷程圖Fig.6 Phase portrait and time history

3 分數階參數的影響分析

下面主要研究分數階微分項的階次和系數對分數階Duffing振子產生分岔和混沌的必要條件的影響。

首先分析分數階微分項的階次對系統(tǒng)產生混沌的必要條件的影響。選取的基本參數為c=0.2，k=1，α=1，ω=1.1，λ=0.1，分數階階次p在[0 1]取值，利用本文解得的解析結果，即式(24)，得到系統(tǒng)產生混沌必要條件的激勵振幅fmin與分數階階次p的關系曲線，示于圖7中用實線表示。從圖7可以看出，當分數階階次p從0增大到1的過程中，系統(tǒng)產生混沌的必要條件fmin是逐漸增大的。為了驗證本文的解析結果的精確度，在不同p值時，利用數值仿真的方法也進行了驗證，圖中圓圈分別代表不同分數階階次下的數值仿真結果，從兩者之間的對比可以看出有一定的差值，本文得到的解析值在定量上小于數值結果，但整體趨勢相同，定性分析是一致的?？梢姺謹惦A階次的增加能抑制系統(tǒng)的分岔和混沌的產生。

圖7 數值解和解析解比較Fig.7 The comparisons between the numerical and analytical solutions

接著分析分數階微分項的系數對系統(tǒng)產生混沌的必要條件的影響。選取的基本參數為c=0.2，k=1，α=1，ω=1.1，分數階階次p在[0 1]取值，如圖8(a)所示。經過對比發(fā)現(xiàn)改變分數階項系數λ變大時，系統(tǒng)產生混沌的必要條件要求fmin也要不斷增大，可見分數階微分項的系數的增大同樣也會抑制分數階Duffing振子產生分岔和混沌的現(xiàn)象。

選取分數階階次p為0.5時，系統(tǒng)產生混沌的必要條件fmin隨著分數階系數改變的變化規(guī)律見圖8(b)所示，從圖中可以看出，隨著分數階微分項系數λ的增加，fmin也會逐漸增大，當分數階微分項系數λ遠遠小于線性剛度k時，解析解和數值解之間的差別較小。當分數階微分項系數λ繼續(xù)增大，分數階微分項的等效線性剛度接近于系統(tǒng)的線性剛度k時，fmin出現(xiàn)增大較快的現(xiàn)象，此時解析解和數值解之間的差值會加大。這是因為利用等效線性阻尼和等效線性剛度來替代分數階微分項的方法省略了式(8)的高階項，當分數階微分項的等效線性剛度與線性剛度k越接近時，省略的兩項高階項所占的比值就會越大，所以此時解析解和數值解差值會加大。因此本文的解析方法要求分數階微分項的系數要遠小于系統(tǒng)的線性剛度k。

圖8 分數階系數的影響Fig.8 The influence of fractional order coefficients

4 結 論

本文利用Melnikov方法研究了含有分數階微分項的Duffing振子產生混沌的必要條件，將分數階Duffing振子中分數階微分項替換成等效線性阻尼與等效線性剛度，用Melnikov方法分析等效后的整數階系統(tǒng)的同宿軌分岔和產生混沌的必要條件，并得到了其解析結果，并通過數值迭代的方法驗證了其精確度。通過對比不同階次情況下數值解與解析解得到的產生混沌的閾值，發(fā)現(xiàn)解析解的趨勢與數值迭代模擬一致，臨界激勵振幅fmin會隨著分數階階次的增加而增大，定性結果相同，Melnikov方法分析系統(tǒng)混沌產生的必要條件是一階近似結果，因此數值結果與解析解之間的定量差異是可以接受的。最后分析了分數階微分項的系數對系統(tǒng)產生混沌的必要條件的影響。在數值模擬過程中,還發(fā)現(xiàn)在分數階Duffing振子中，當選取不同的初值時，存在兩條通過倍周期分岔通往混沌的路徑，并通過分析其動力學響應證實了這一現(xiàn)象。本文提出的將分數階微分項替換為三角函數和指數函數形式的等效線性阻尼和線性剛度的方法可以推廣應用于其他類似的分數階系統(tǒng)中，為分數階系統(tǒng)的動力學研究提供了新的研究思路。