， *，

(1.天津大學(xué) 機(jī)械工程學(xué)院 力學(xué)系，天津 300072； 2.天津大學(xué) 天津市非線性動(dòng)力學(xué)與混沌控制重點(diǎn)實(shí)驗(yàn)室，天津 300072)

1 引 言

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

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

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

2 預(yù)備知識(shí)

為敘述方便起見，先介紹本文用到的三種分?jǐn)?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ù)等價(jià)[16]。

(5)

式中n-1

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

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

(6)

式中n-1

3.1 含分?jǐn)?shù)階導(dǎo)數(shù)項(xiàng)非線性振動(dòng)系統(tǒng)狀態(tài)方程的標(biāo)準(zhǔn)形式

利用Caputo導(dǎo)數(shù)算子的疊加關(guān)系，可將含m個(gè)分?jǐn)?shù)階導(dǎo)數(shù)項(xiàng)的單自由度非線性振動(dòng)系統(tǒng)化為標(biāo)準(zhǔn)形式的狀態(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

先考慮周期激勵(lì)下的含單個(gè)分?jǐn)?shù)階導(dǎo)數(shù)項(xiàng)的系統(tǒng)：

(8)

式中 0

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

(9)

進(jìn)行增維處理[13,18-20]，可將系統(tǒng)(8)化為等價(jià)系統(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

當(dāng)0

(11)

當(dāng)1

(12)

考慮周期激勵(lì)下含多個(gè)分?jǐn)?shù)階導(dǎo)數(shù)項(xiàng)的系統(tǒng)：

(13)

式中 0

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

(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)

上述分析表明，單自由度含分?jǐn)?shù)階導(dǎo)數(shù)項(xiàng)(單個(gè)或多個(gè))的非線性振動(dòng)系統(tǒng)，均可轉(zhuǎn)化為標(biāo)準(zhǔn)形式(7)。多自由度非線性振動(dòng)系統(tǒng)也可做同樣的轉(zhuǎn)化。因此，只要給出分?jǐn)?shù)階標(biāo)準(zhǔn)微分方程組(7)的算法，就可以解決含分?jǐn)?shù)階導(dǎo)數(shù)的非線性振動(dòng)系統(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)

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

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

4 數(shù)值仿真

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

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

(19)

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

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

(20)

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

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

Fig.1 Comparison of numerical and analytical solutions

of example 1

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

(21)

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

(22)

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

算例3含兩個(gè)分?jǐn)?shù)階項(xiàng)的情形：

(23)

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

(24)

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

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

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

(25)

式中K1>0，0≤p≤2，m，k和c分別為系統(tǒng)的質(zhì)量、線性剛度系數(shù)和線性阻尼系數(shù)，K1為分?jǐn)?shù)階導(dǎo)數(shù)項(xiàng)的系數(shù)，F(xiàn)1和w分別為外激勵(lì)的幅值和頻率，各參數(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)方程：

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

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

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

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

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

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

Fig.3 Comparison of numerical and analytical solutions

of example 3

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

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

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

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

圖5 初始條件遠(yuǎn)離吸引域邊界

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é)果的計(jì)算，也可能影響吸引域的計(jì)算。在后續(xù)研究中，還要進(jìn)一步考慮構(gòu)造該算法的矩陣格式，從而利用Matlab矩陣運(yùn)算高效的特點(diǎn)，提高其計(jì)算效率。