喬 智，趙維加，黃健飛

(1. 揚(yáng)州大學(xué)數(shù)學(xué)科學(xué)學(xué)院，江蘇揚(yáng)州225002； 2. 青島大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，山東青島266071)

分?jǐn)?shù)階微積分[1]是數(shù)學(xué)的一個(gè)重要的分支，距今已經(jīng)有300多年的歷史。近年來，分?jǐn)?shù)階微積分算子的非局部性使得分?jǐn)?shù)階微分方程在建模復(fù)雜動(dòng)力系統(tǒng)和材料方面的應(yīng)用越來越廣泛，比如各種黏彈性材料[2]、反常擴(kuò)散[3-4]、分形理論[5]、化學(xué)和生物化學(xué)[6]。與整數(shù)階微分方程相比，分?jǐn)?shù)階微分方程的物理意義更清晰，模型更簡單；但是，精確求解分?jǐn)?shù)階微分方程通常比較困難。對(duì)于一些非線性分?jǐn)?shù)階問題，甚至不可能獲得其精確解，因此，研究者廣泛設(shè)計(jì)和使用數(shù)值方法來求解分?jǐn)?shù)階微分方程。

目前，許多論文研究了分?jǐn)?shù)階常微分方程數(shù)值方法。早在1997年，Diethelm[7]對(duì)一類分?jǐn)?shù)階常微分方程構(gòu)造了基于求積公式的隱格式，并給出了嚴(yán)格的誤差分析。這是公認(rèn)的分?jǐn)?shù)階微分方程數(shù)值方法領(lǐng)域比較早的工作之一。隨后，Edwards等[8]研究了線性多項(xiàng)分?jǐn)?shù)階微分方程的數(shù)值解，但是沒有考慮非線性多項(xiàng)分?jǐn)?shù)階微分方程的情形。El-Mesiry等[9]給出了求解非線性多項(xiàng)分?jǐn)?shù)(任意)階微分方程的數(shù)值格式，只通過數(shù)值實(shí)驗(yàn)來驗(yàn)證格式的計(jì)算效果。Kumar等[10]采用積分方程的思想，針對(duì)一類分?jǐn)?shù)階常微分方程設(shè)計(jì)了一個(gè)高階格式，數(shù)值實(shí)驗(yàn)表明該方法有較高的精度。Lin等[11]采用生成函數(shù)法推導(dǎo)出求解分?jǐn)?shù)階常微分方程的高階近似，證明了這些方法的可解性、收斂性和穩(wěn)定性，但是這些方法只能用來求解含有一個(gè)分?jǐn)?shù)階導(dǎo)數(shù)項(xiàng)的方程。Ford等[12]考慮了多項(xiàng)分?jǐn)?shù)階微分方程近似解的基于系統(tǒng)的分解方法，但是只考慮了各分?jǐn)?shù)階指標(biāo)單調(diào)遞增的情況。Huang等[13]詳細(xì)地證明了文獻(xiàn)[10]中所給出的數(shù)值方法至少具有3階精度。此外，關(guān)于分?jǐn)?shù)階偏微分方程數(shù)值方法的研究，可以參考文獻(xiàn)[14-18]。

本文中研究一類多項(xiàng)分?jǐn)?shù)階非線性常微分方程初值問題的數(shù)值方法；為了降低對(duì)解函數(shù)正則性的要求[19]，利用Riemann-Liouville積分算子，將該類多項(xiàng)分?jǐn)?shù)階非線性常微分方程化為與之等價(jià)的積分方程；為了避免求解非線性方程和保證數(shù)值方法的穩(wěn)定性，分別利用復(fù)化右矩形公式和左矩形公式來離散積分方程左邊和右邊的Riemann-Liouville積分，獲得一個(gè)顯式的數(shù)值方法；分別考慮2種不同形式的多項(xiàng)分?jǐn)?shù)階非線性微分方程，并證明2種不同問題數(shù)值方法的收斂性和穩(wěn)定性。

1 情況1下數(shù)值方法的構(gòu)造

考慮情況1，即0<α1<α2<…<αM<1條件下，帶有初值條件u(0)=u0的多項(xiàng)分?jǐn)?shù)階非線性微分方程

(1)

n-1≤αm≤n,

αm(m=1,2,…,M)為階數(shù)，n為正整數(shù)；u(0)為0時(shí)刻的函數(shù)值，u0為常數(shù)；f(t,u)在區(qū)間Ω上連續(xù)，其中Ω=[0,T]×，f關(guān)于第2個(gè)自變量滿足Lipschitz條件，即

|f(t,u1)-f(t,u2)|≤L|u1-u2| ，

其中Lipschitz常數(shù)L>0。

則有

(2)

因此式(2)可化為

(3)

對(duì)于T≤1,假設(shè)f[t,u(t)]在[0,T]上一階連續(xù),分別運(yùn)用復(fù)化乘積右矩形公式和左矩形公式對(duì)式(3)左邊和右邊的Riemann-Liouville積分進(jìn)行逼近，可得

(4)

成立，其中

在式(4)中，移去截?cái)嗾`差O(h),用數(shù)值解un代替精確解u(tn),則有

(5)

其中

由此，式(3)的數(shù)值方法為

(6)

對(duì)于T≥1,對(duì)式(1)作變量替換，令sT=t,則0

(7)

與T≤1的情形相似，式(7)可離散為

(8)

在式(8)中，移去截?cái)嗾`差O(h),則數(shù)值方法變?yōu)?/p>

2 情況1下數(shù)值方法的收斂性與穩(wěn)定性

討論數(shù)值方法式(6)、(9)在情況1下的收斂性與穩(wěn)定性。在證明過程中用到Gronwall不等式[6]。

引理1設(shè)C1>0且獨(dú)立于h>0,C2為正常數(shù)，且滿足不等式

j=0,1,…,n-1,nh≤T,

其中0<α<1,則有

|zn|≤C2Eα[C1Γ(α)Tα],nh≤T,

其中Eα為Mittag-Leffler函數(shù)，定義為

下面給出收斂性分析。

定理1假設(shè)f[t,u(t)]在區(qū)間[0,T]上一階連續(xù),對(duì)充分小的時(shí)間步長h,式(1)的數(shù)值方法式(6)、(9)在情況1下的收斂階數(shù)為1,即|en|=|u(tn)-un|≤Ch,n=1,2,…,N,其中C為與離散化參數(shù)無關(guān)的常數(shù)，并且C可能有不同的值。

證明： 由于式(6)、(9)的收斂性分析類似，因此只證明式(6)的收斂性。f[t,u(t)]在區(qū)間[0,T]上一階連續(xù),根據(jù)中值定理，存在Lk,使得

f[tk,u(tk)]-f(tk,uk)=Lk[u(tk)-uk]=Lkek,

k=1,2,…,N。

成立，其中ek=u(tk)-uk。f(t,u)滿足Lipschitz條件，則|Lk|≤L。對(duì)于0<β<1,1≤k≤n-1,有

(10)

和

(11)

式(10)減式(11)，并進(jìn)行簡單放縮，有

由于tn-k=h(n-k),因此

ChαM-α1|en|+Ch。

考慮到tn-k∈ (0,1],有

CMhαM-αM-1|en|+Ch≤

CMhαM-αM-1|en|+Ch,

因此

(1-CMhαM-αM-1)|en|≤

當(dāng)h充分小時(shí)，由引理1，可得

定理1證畢。

利用定理1的證明思想，可以建立數(shù)值方法式(6)、(9)的穩(wěn)定性結(jié)果。

定理2對(duì)于充分小的h,式(1)的數(shù)值方法式(6)、(9)是無條件穩(wěn)定的，即

由于穩(wěn)定性分析與收斂性分析類似，因此定理2證明略。

3 情況2下數(shù)值方法的構(gòu)造以及收斂性與穩(wěn)定性

考慮式(1),初值條件為u(0)=u0,u′(0)=a且0<α1<α2<…<αM-1<1<αM<2,αM-αM-1>1。應(yīng)用情況1中的離散化方法，可以得到

由此可知，在情況2條件下的數(shù)值方法為

(12)

下面給出數(shù)值方法式(12)的收斂性和穩(wěn)定性分析。

定理3在式(1)中，假設(shè)f[t,u(t)]在區(qū)間[0,T]上一階連續(xù)，且對(duì)充分小的h,在情況2下，數(shù)值方法式(12)的收斂階數(shù)為1，即

|en|=|u(tn)-un|≤Ch,n=1,2,…,N。

證明: 與情況1相似，有

由于αM-αM-1-1>0,…,αM-α1-1>0,因此

對(duì)于充分小的h,根據(jù)Gronwall不等式，有

|en|≤Ch。

定理3證畢。

定理4對(duì)于充分小的h,求解式(1)的數(shù)值方法式(10)是無條件穩(wěn)定的，即

證明: 與情況1的穩(wěn)定性分析相似，等式

由于αM-αM -1-1>0,…,αM-α1-1>0,因此

對(duì)于充分小的h,根據(jù)Gronwall不等式，有

|εn|≤Cη。

定理4證畢。

4 數(shù)值實(shí)驗(yàn)

例1分別在T=1 s和T=2 s時(shí)考慮初值條件為u(0)=0的帶有2項(xiàng)Caputo導(dǎo)數(shù)的方程

其中0<α1<α2<1。該方程的精確解為u(t)=tα2+1。注意到方程右端的光滑性滿足定理1的條件。圖1所示為T=1 s且h=0.025 s時(shí)，精確解和數(shù)值解的比較。從圖中可以看出，數(shù)值解與精確解非常一致。表1所示為當(dāng)T=1 s，α1=0.1，α2=0.3，α1=0.8，α2=0.9且不同時(shí)間步長時(shí)式(6)的誤差和收斂階。表2所示為當(dāng)T=2 s，α1=0.2，α2=0.4，α1=0.6，α2=0.8且不同時(shí)間步長時(shí)式(9)的誤差和收斂階。從表1、2中可以看出，隨著步長h的減小，數(shù)值方法的誤差也減小，收斂階約數(shù)為1。

α1、α2—收斂階數(shù)； t—時(shí)間； u(t)—關(guān)于t的函數(shù)。圖1 時(shí)間T=1 s, α1=0.1,α2=0.3，α1=0.8, α2=0.9時(shí)方程的數(shù)值解與精確解

表1 時(shí)間T=1 s, α1=0.1, α2=0.3，α1=0.8, α2=0.9且不同時(shí)間步長時(shí)式(6)的誤差和收斂階數(shù)

表2 時(shí)間T=2 s,α1=0.2,α2=0.4，α1=0.6,α2=0.8且不同時(shí)間步長時(shí)式(9)的誤差和收斂階數(shù)

例2分別在T=1 s和T=2 s時(shí)，考慮初值條件為u(0)=0的帶有5項(xiàng)Caputo導(dǎo)數(shù)的方程

其中0<α1<α2<α3<α4<α5<1。該方程的精確解為u(t)=tα5+1。圖2所示為當(dāng)T=2 s,h=0.05 s,α1=0.1，α2=0.3，α3=0.5，α4=0.7，α5=0.9時(shí)方程的精確解和數(shù)值解。從圖中可以看出，數(shù)值解與精確解非常吻合。表3所示為T=1 s,α1=0.1，α2=0.3，α3=0.4，α4=0.8，α5=0.9且不同步長時(shí)式(6)的誤差和收斂階數(shù)。表4所示為T=2 s,α1=0.1，α2=0.3，α3=0.5，α4=0.7，α5=0.9且不同步長時(shí)式(9)的誤差和收斂階數(shù)。從表3、4中可以看出，當(dāng)h減小到之前的1/2時(shí)，誤差也減小到之前的1/2，收斂階約數(shù)為1。

α1、α2、α3、α4、α5—收斂階數(shù)； t—時(shí)間； u(t)—關(guān)于t的函數(shù)。圖2 時(shí)間T=2 s,α1=0.1,α2=0.3,α3=0.5,α4=0.7,α5=0.9時(shí)方程的數(shù)值解與精確解

表3 時(shí)間T=1 s,α1=0.1,α2=0.3,α3=0.4,α4=0.8,α5=0.9且不同時(shí)間步長時(shí)式(6)的誤差和收斂階數(shù)

表4 時(shí)間T=2 s,α1=0.1,α2=0.3,α3=0.5,α4=0.7,α5=0.9且不同時(shí)間步長時(shí)式(9)的誤差和收斂階數(shù)

例3 當(dāng)T=2 s時(shí)，考慮初值條件為u(0)=0,u′(0)=0的帶有3項(xiàng)Caputo導(dǎo)數(shù)的方程

其中0<α1<α2<1<α3<2。該問題的精確解為u(t)=tα1+2。表5所示為當(dāng)T=2 s,α1=0.3,α2=0.7 ,α3=1.5時(shí)式(12)的誤差和收斂階數(shù)。從表中可以看出，當(dāng)h減小時(shí)，誤差也隨之減小，并且收斂階數(shù)趨近于1。

表5 時(shí)間T=2 s,α1=0.3,α2=0.7,α3=1.5且不同步長時(shí)式(10)的誤差和收斂階數(shù)

5 結(jié)論

本文中基于等價(jià)的積分方程，對(duì)具有不同初值條件的多項(xiàng)分?jǐn)?shù)階非線性微分方程，構(gòu)造了具有一階精度的顯式數(shù)值方法。針對(duì)2個(gè)不同的初值問題，證明了數(shù)值方法的收斂性和穩(wěn)定性。最后，通過計(jì)算帶有2、5、3項(xiàng)Caputo分?jǐn)?shù)階導(dǎo)數(shù)的數(shù)值算例，驗(yàn)證了理論分析結(jié)果的正確性。在以后的研究中，將嘗試把本文中的數(shù)值方法擴(kuò)展到多項(xiàng)時(shí)間分?jǐn)?shù)階偏微分方程的數(shù)值求解。