， *，

(1.天津大學 機械工程學院 力學系，天津 300072； 2.天津大學 天津市非線性動力學與混沌控制重點實驗室，天津 300072)

1 引 言

分數(shù)階微積分是對整數(shù)階微積分的一種推廣，將微積分運算的階次從傳統(tǒng)的整數(shù)階擴展到了分數(shù)階和復(fù)數(shù)階，至今已有三百多年的發(fā)展歷史，已成為研究反常擴散、多孔介質(zhì)力學、非牛頓流體力學、粘彈性力學以及軟物質(zhì)物理等學科領(lǐng)域的有力數(shù)學工具。

關(guā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ù)階微分方程的初始條件類似，更符合工程實際，因而在有明確物理和工程背景的研究工作中應(yīng)用較多[1-3]。對于Caputo導(dǎo)數(shù)下的微分方程，目前導(dǎo)數(shù)階數(shù)相同時相關(guān)的算法已經(jīng)成熟[4-10]，階數(shù)不同時，文獻[11]考慮了分段線性Chen系統(tǒng)的混沌和多穩(wěn)態(tài)行為，但未給出具體的算法。為計算階數(shù)不同時系統(tǒng)的響應(yīng)，文獻[2，3,17]采用了變通的做法，即用分數(shù)階的Grunwald-Letnikov導(dǎo)數(shù)的算法來計算Caputo導(dǎo)數(shù)，這種處理是否可行，尚未見報道。

為此，本文研究階數(shù)不同時Caputo導(dǎo)數(shù)意義下微分方程的算法。首先，給出一些預(yù)備知識；然后，推導(dǎo)含分數(shù)階導(dǎo)數(shù)項非線性振動系統(tǒng)狀態(tài)方程的標準形式及用于計算Caputo導(dǎo)數(shù)的迭代格式；最后，通過算例說明本文理論分析的正確性和所得計算方法的可靠性，以及用分數(shù)階的Grunwald-Letnikov導(dǎo)數(shù)算法來計算Caputo導(dǎo)數(shù)時可能存在的問題。

2 預(yù)備知識

為敘述方便起見，先介紹本文用到的三種分數(shù)階導(dǎo)數(shù)的定義及其聯(lián)系。

定義2.1定義在區(qū)間[a,b]上的函數(shù)x(t)的p階Caputo導(dǎo)數(shù)為[12]

(1)

式中n-1

定義2.2定義在區(qū)間[a,b]上的函數(shù)x(t)的p階Riemann-Liouville導(dǎo)數(shù)為[13]

(2)

式中n-1

定義2.3函數(shù)x(t)的p階Grunwald-Letnikov導(dǎo)數(shù)定義為[14]

Dp[x(t)]=dp[x(t)]/dtp|t=kh=

(3)

設(shè)x(t)定義在[a,b]上，若x(t)具有n-1階連續(xù)導(dǎo)函數(shù)，且x(n)(t)∈L1(a,b)，則Caputo導(dǎo)數(shù)和Riemann-Liouville導(dǎo)數(shù)滿足關(guān)系[15]

(4)

式中n-1

設(shè)x(t)定義在[a,b]上，若x(t)具有n階連續(xù)導(dǎo)函數(shù)，則Riemann-Liouville導(dǎo)數(shù)與Grunwald-Letnikov導(dǎo)數(shù)等價[16]。

(5)

式中n-1

3 基于Caputo導(dǎo)數(shù)迭代格式的構(gòu)造

在非線性振動系統(tǒng)的研究中，通常取初始時刻t=0，對應(yīng)的Caputo導(dǎo)數(shù)為

(6)

式中n-1

3.1 含分數(shù)階導(dǎo)數(shù)項非線性振動系統(tǒng)狀態(tài)方程的標準形式

利用Caputo導(dǎo)數(shù)算子的疊加關(guān)系，可將含m個分數(shù)階導(dǎo)數(shù)項的單自由度非線性振動系統(tǒng)化為標準形式的狀態(tài)方程：

(7)

式中q=[q1,…,qm +2]T

X(t)=[x1(t),…,xm +2(t)]T

f[X(t),t]={f1[x(t),t],…,fm+2[x(t),t]}T

階數(shù)qj滿足0

先考慮周期激勵下的含單個分數(shù)階導(dǎo)數(shù)項的系統(tǒng)：

(8)

式中 0

在x(t)連續(xù)且分數(shù)階導(dǎo)數(shù)存在的條件下，運用Caputo分數(shù)階算子的疊加關(guān)系[17]：

(9)

進行增維處理[13,18-20]，可將系統(tǒng)(8)化為等價系統(tǒng)：

(10)

式中q=[q1,q2,q3]T

X(t)=[x1(t),x2(t),x3(t)]T

f[X(t),t]= {f1[x(t),t],f2[x(t),t],

f3[x(t),t]}T

當0

(11)

當1

(12)

考慮周期激勵下含多個分數(shù)階導(dǎo)數(shù)項的系統(tǒng)：

(13)

式中 0

同理，可將式(13)化為等價形式：

(14)

式中q1=p1,qj=pj-pj -1(j=2,…,l)

ql +1=1-pl，ql +2=pl +1-1

qj=pj -1-pj -2(j=l+3,…,m+1)

qm +2=2-pm,x(t)=[x1(t),…,xm +2(t)]T

且

(15)

上述分析表明，單自由度含分數(shù)階導(dǎo)數(shù)項(單個或多個)的非線性振動系統(tǒng)，均可轉(zhuǎn)化為標準形式(7)。多自由度非線性振動系統(tǒng)也可做同樣的轉(zhuǎn)化。因此，只要給出分數(shù)階標準微分方程組(7)的算法，就可以解決含分數(shù)階導(dǎo)數(shù)的非線性振動系統(tǒng)的數(shù)值求解問題。

3.2 一般迭代格式的構(gòu)造

由第2節(jié)給出的三種導(dǎo)數(shù)之間的關(guān)系可知

(16)

式中 0

fj[x(tk -1),tk -1]

(17)

式中j=1,…,m+2，從而有

(18)

選取適當?shù)臅r間步長h，通過以上迭代格式，可以得到x(t)的數(shù)值解，即Caputo導(dǎo)數(shù)意義下微分方程(7)的解。為方便起見，稱該算法為Caputo算法。不要求狀態(tài)方程中各分數(shù)階導(dǎo)數(shù)階數(shù)相等，是該算法的主要優(yōu)點。

若去掉方程(18)右端第二項，即得到求GL導(dǎo)數(shù)意義下標準形式方程(7)解的迭代格式，以下簡稱GL算法。不難發(fā)現(xiàn)，只有在零初值條件下，兩種算法的迭代格式才完全一致，對應(yīng)的解才完全相同。而在非零初始條件下，兩種算法的迭代格式形式不同，二者能否得到相同的結(jié)果是一個值得討論的問題，關(guān)系到用GL算法計算Caputo分數(shù)階導(dǎo)數(shù)的可行性。

4 數(shù)值仿真

通過若干算例對本文提出的算法進行驗證。首先，通過有解析解的算例1～3說明本文計算方法的正確性，所選算例均來自文獻[21]。然后，通過算例4說明用基于Grunwald-Letnikov導(dǎo)數(shù)算法計算Caputo導(dǎo)數(shù)下的非線性振動系統(tǒng)存在的問題。

算例1分數(shù)階導(dǎo)數(shù)階數(shù)小于1的情況：

(19)

其精確解為x(t)=1+t。

由Caputo分數(shù)階導(dǎo)數(shù)算子的疊加關(guān)系[17]，對方程(19)進行相關(guān)處理，可將其轉(zhuǎn)化為等價形式:

(20)

取時間步長h=0.001，用3.2節(jié)中的Caputo算法，可得方程(19)的數(shù)值解。數(shù)值解與解析解對比如圖1所示，二者完全重合。

圖1 算例1中數(shù)值解與解析解對比

Fig.1 Comparison of numerical and analytical solutions

of example 1

算例2分數(shù)階導(dǎo)數(shù)階數(shù)大于1的情形：

(21)

同理，可將方程(21)轉(zhuǎn)化為

(22)

從圖2可以看出,數(shù)值解(時間步長取h=0.001)與解析解同樣吻合很好。

算例3含兩個分數(shù)階項的情形：

(23)

同上，可將方程(23)轉(zhuǎn)化為

(24)

從圖3可以看出，數(shù)值解(時間步長取h=0.0001)和解析解吻合很好。

上述算例表明，本文計算方法正確且具有很高的精度。

算例4分數(shù)階Duffing振子系統(tǒng)[2]：

(25)

式中K1>0，0≤p≤2，m，k和c分別為系統(tǒng)的質(zhì)量、線性剛度系數(shù)和線性阻尼系數(shù)，K1為分數(shù)階導(dǎo)數(shù)項的系數(shù)，F(xiàn)1和w分別為外激勵的幅值和頻率，各參數(shù)的取值分別為m=5，k=45，c=0.2，α1=0.5，K1=0.1，p=0.1，F(xiàn)1=100和w=1.0612。

方程(25)可轉(zhuǎn)化為狀態(tài)方程：

對系統(tǒng)初值分三種情況討論，分別將Caputo算法及GL算法所得結(jié)果進行比較(時間步長均取h=0.002)。

可以看出，初始條件均取(0, 0)時，兩種算法等價，計算結(jié)果完全重合。

可以看出，兩種算法所得結(jié)果經(jīng)過不同過渡過程后，達到相同的穩(wěn)態(tài)。

圖2 算例2中數(shù)值解與解析解對比

Fig.2 Comparison of numerical and analytical solutions of example 2

圖3 算例3中數(shù)值解與解析解對比

Fig.3 Comparison of numerical and analytical solutions

of example 3

可以看出，兩種算法所得結(jié)果經(jīng)過不同過渡過程后，分別達到兩個不同的穩(wěn)態(tài)解。此時GL算法不能得到正確的結(jié)果。

從本例可見，用GL算法來計算Caputo導(dǎo)數(shù)下的微分方程，除零初始條件下沒有問題外，因初始條件不同，計算結(jié)果存在如下可能。(1) 僅過渡過程有差異，穩(wěn)態(tài)結(jié)果正確。(2) 過渡過程不同，穩(wěn)態(tài)結(jié)果也不同。這種現(xiàn)象可能與初始條件是否靠近吸引域的邊緣有關(guān)，后續(xù)工作將對此展開詳細討論。

圖4 零初始條件下兩種方法的比較

Fig.4 Comparison of the two methods under the zero initial condition

圖5 初始條件遠離吸引域邊界

Fig.5 Initial condition is far from the boundary of the basin of attraction

圖6 初始條件靠近吸引域邊界

Fig.6 Initial condition is close to the boundary of the basin of attraction

5 結(jié) 論

本文基于Caputo導(dǎo)數(shù)與Riemann-Liouville導(dǎo)數(shù)和Grunwald-Letnikov導(dǎo)數(shù)間的關(guān)系，建立了基于Caputo導(dǎo)數(shù)的數(shù)值離散方法，并通過數(shù)值算例討論了所得算法的正確性及基于Grunwald-Letnikov導(dǎo)數(shù)和Caputo導(dǎo)數(shù)離散方法的差異，比較兩種算法所得結(jié)果可得，用GL算法替代Caputo算法，不僅可能影響正確結(jié)果的計算，也可能影響吸引域的計算。在后續(xù)研究中，還要進一步考慮構(gòu)造該算法的矩陣格式，從而利用Matlab矩陣運算高效的特點，提高其計算效率。