代 躍，周曉軍

(貴州師范大學(xué) 數(shù)學(xué)科學(xué)學(xué)院，貴州 貴陽 550025)

0 引言

分?jǐn)?shù)階微積分和經(jīng)典微積分研究幾乎同時(shí)開始[1]，但由于分?jǐn)?shù)階微積分的實(shí)際應(yīng)用受限，及缺乏物理背景的支持，發(fā)展緩慢。但分?jǐn)?shù)階微積分具有描述物質(zhì)記憶功能和遺傳效應(yīng)的特征，使得它比整數(shù)階導(dǎo)數(shù)的描述更加精確。近年來，分?jǐn)?shù)階微積分被廣泛的應(yīng)用在自然科學(xué)和工程計(jì)算[2-3]的諸多領(lǐng)域內(nèi)。如被用來描述粘彈性材料、材料的電性質(zhì)、電磁波、電化學(xué)過程，動(dòng)力系統(tǒng)中的控制理論，生物系統(tǒng)中的電導(dǎo)、混沌、神經(jīng)細(xì)胞中離子的反常擴(kuò)散、生化、水文、金融及其他應(yīng)用領(lǐng)域中的問題。分?jǐn)?shù)階微積分理論也受到越來越多的國內(nèi)外學(xué)者的廣泛關(guān)注，特別是從實(shí)際問題抽象出來的分?jǐn)?shù)階微分方程成為很多數(shù)學(xué)工作者的研究熱點(diǎn)。

由于分?jǐn)?shù)階微分的解析解很難求得，所以像整數(shù)階一樣借助數(shù)值方法求解是一個(gè)好的選擇。因此，研究分?jǐn)?shù)階微分方程的數(shù)值算法成為一個(gè)熱點(diǎn)。至今，已提出了許多經(jīng)典的數(shù)值方法[4-7]：Diethelm等[4]提出一種Adams型的預(yù)估-校正方法，并給出了這種算法的誤差分析；Lin等[5]在Lubich方法的基礎(chǔ)上，提出了一種線性多步法來求解分?jǐn)?shù)階微分方程初值問題，并證明了該方法的相容性、收斂性和穩(wěn)定性；Kumar等[6]利用求解Volterra方程的block-by-block方法設(shè)計(jì)了一類新的分?jǐn)?shù)階block-by-block方法，數(shù)值算例也證明了此方法的有效性和穩(wěn)定性,但并沒有作出證明；李[7]推導(dǎo)出了分?jǐn)?shù)階積分的Chebyshev小波矩陣，并利用它來求解了非線性分?jǐn)?shù)階微分方程。

本文考慮非線性分?jǐn)?shù)階微分方程組初值問題

(1)

已有一些文章對(duì)該方程組進(jìn)行研究，Varsha Daftardar-Gejji等[8]對(duì)方程組解的存在性、唯一性以及穩(wěn)定性進(jìn)行了論證。林永華等[9]根據(jù)分?jǐn)?shù)階導(dǎo)數(shù)的高階近似提出了分?jǐn)?shù)階微分方程組的高階近似法。童啟秀等[10]利用求解普通積分方程的Adams技巧，建立了分?jǐn)?shù)階微分方程組的一種顯式數(shù)值算法。代群等[11]使用變分迭代法求解分?jǐn)?shù)階微分方程組，并改進(jìn)了校正函數(shù)。欒新等[12]根據(jù)延遲校正法的思想來設(shè)計(jì)求解分?jǐn)?shù)階微分方程組初值問題的高精度格式。

田獻(xiàn)珍等[13]將求解整數(shù)解方程的Euler法應(yīng)用于非線性分?jǐn)?shù)階常微分方程模型。本文對(duì)Euler法進(jìn)行推廣，將方程組(1)轉(zhuǎn)化為積分方程組后對(duì)每一個(gè)積分方程應(yīng)用Euler法，得到了方程組(1)的隱式算法，并證明了算法的收斂性及穩(wěn)定性，最后給出了數(shù)值實(shí)例。

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

引理1[14]若連續(xù)函數(shù)y(t)是初值問題

(2)

的解，則y(t)是第二類Volterra積分方程：

的解，反之亦然。

由引理1可將方程組(1)轉(zhuǎn)化為如下的Volterra積分方程組：

(3)

引理2[15]設(shè)β,γ≥0，h>0，Nh≤(b-a)，0

則必有

ηm≤βEα(γΓ(α)(mh))α,m=0,1,2,…，N

其中

是單參數(shù)Mittag-Leffler函數(shù)。

2 Euler法

對(duì)方程組(3)的每一個(gè)積分方程采用文獻(xiàn)[13]中所用的積分方法。首先進(jìn)行網(wǎng)格剖分：取均勻網(wǎng)格節(jié)點(diǎn)tj=jh,j=0,1,2,…,N,h=T/N是積分步長。(2)式中取t=tm得:

類似于數(shù)值積分的處理方法:在子區(qū)間[tj,tj+1],j=0,1,2,…,m上，若用fi(tj+1,y1(tj+1),y2(tj+1),…,yz(tj+1))近似代替fi(t,y1(t),y2(t),…,yz(t))，得：

令

yi,j≈yi(tj),gi,j=gi(tj),fi,j=fi(tj,y1(tj),y2(tj),…,yz(tj))。

得

(4)

(5)

我們稱式(4)和(5)為分?jǐn)?shù)階隱式Euler法。

3 收斂性分析

引理3[13]當(dāng)z=1時(shí)，分?jǐn)?shù)階微分方程組(1)退化為分?jǐn)?shù)階微分方程(2)。設(shè)α>0，y(t)充分光滑，且CDαy(t)∈C1[0,T]，函數(shù)f(t,y)關(guān)于第二個(gè)變量滿足Lipschitz條件，則分?jǐn)?shù)階隱式Euler法的誤差滿足：

定理1 設(shè)fi(t,y1(t),y2(t),…,yz(t))∈C3(0,T),i=1,2,…,z,z≥2。則有分?jǐn)?shù)階常微分方程組隱式Euler法(4)的誤差滿足：

為證明這個(gè)定理，需要用到引理3。

證明

令Y(t)=(y1(t),y2(t),…,yz(t))，則：fi(t,y1(t),y2(t),…,yz(t))=fi(t,Y(t))。

當(dāng)t=t0時(shí)，|yi(t0)-yi,0|=0，等式成立。

當(dāng)t=t1時(shí)，

|yi(t1)-yi,1|=

由引理3得

|yi(t1)-yi,1|

其中

令

得

當(dāng)t=tm時(shí)，

|yi(tm)-yi,m|=

由引理3：

又有：

4 算法穩(wěn)定性分析

證明

由

兩式相減，得：

因?yàn)?/p>

所以有

則由Gronwall不等式可得存在常數(shù)C，滿足：

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

例1 考慮非線性分?jǐn)?shù)階微分方程組

初始條件為：

y1(0)=y2(0)=0,0≤t≤1，

y1′(0)=y2′(0)=0,(αi>1)。

該非線性分?jǐn)?shù)階微分方程組的一組精確解為：

表1 α1=0.5,α2=0.8和α1=1.3,α2=1.6，y1的最大誤差隨時(shí)間步長h的變化與收斂階Tab.1 α1=0.5,α2=0.8 and α1=1.3,α2=1.6 change of the maximum error of y1 with time step h and order of convergence

表2 α1=0.5,α2=0.8和α1=1.3,α2=1.6，y2的最大誤差隨時(shí)間步長h的變化與收斂階Tab.2 α1=0.5,α2=0.8 and α1=1.3,α2=1.6 Change of the maximum error of y2 with time step h and order of convergence

圖1 α1=0.5,α2=0.8和α1=1.3,α2=1.6在h=0.01時(shí)y1(t)的數(shù)值解和精確解曲線Fig.1 Analytical solution and numerical solution for y1(t)with α1=0.5,α2=0.8 and α1=1.3,α2=1.6 when h=0.01

圖2 α1=0.5,α2=0.8和α1=1.3,α2=1.6在h=0.01時(shí)y2(t)的數(shù)值解和精確解曲線Fig.2 Analytical solution and numerical solution for y2(t)with α1=0.5,α2=0.8 and α1=1.3,α2=1.6 when h=0.01

例2 考慮非線性分?jǐn)?shù)階微分方程組

初始條件為：y1(0)=y2(0)=y3(0)=0,0≤t≤1,

y1′(0)=y2′(0)=y3′(0)=0,(αi>1)。

該非線性分?jǐn)?shù)階微分方程組的一組精確解為：

表3 α1=0.5,α2=0.6,α3=0.7和α1=1.1,α2=1.2,α3=1.3，y1的最大誤差隨時(shí)間步長h的變化與收斂階Tab.3 α1=0.5,α2=0.6,α3=0.7 and α1=1.1,α2=1.2,α3=1.3Change of the maximum error of y1 with time step h and order of convergence

表4 α1=0.5,α2=0.6,α3=0.7和α1=1.1,α2=1.2,α3=1.3，y2的最大誤差隨時(shí)間步長h的變化與收斂階Tab.4 α1=0.5,α2=0.6,α3=0.7 and α1=1.1,α2=1.2,α3=1.3Change of the maximum error of y2 with time step h and order of convergence

表5 α1=0.5,α2=0.6,α3=0.7和α1=1.1,α2=1.2,α3=1.3，y3的最大誤差隨時(shí)間步長h的變化與收斂階Tab.5 α1=0.5,α2=0.6,α3=0.7 and α1=1.1,α2=1.2,α3=1.3Change of the maximum error of y3 with time step h and order of convergence

圖3 α1=0.5,α2=0.6,α3=0.7在h=0.01時(shí)y1(t)的數(shù)值解和精確解曲線Fig.3 Analytical solution and numerical solution for y1(t) with α1=0.5,α2=0.6,α3=0.7 whenh=0.01

圖4 α1=0.5,α2=0.6,α3=0.7在h=0.01時(shí)y2(t)的數(shù)值解和精確解曲線Fig.4 Analytical solution and numerical solution for y2(t)with α1=0.5,α2=0.6,α3=0.7 whenh=0.01

圖5 α1=0.5,α2=0.6,α3=0.7在h=0.01時(shí)y3(t)的數(shù)值解和精確解曲線Fig.5 Analytical solution and numerical solution for y3(t)with α1=0.5,α2=0.6,α3=0.7 whenh=0.01

6 結(jié)論

本文提出用分?jǐn)?shù)階Euler法來解一類非線性分?jǐn)?shù)階常微分方程組,并證明了該數(shù)值方法的收斂性和穩(wěn)定性,所舉的數(shù)值例子也證實(shí)了分?jǐn)?shù)階Euler法是解這類非線性分?jǐn)?shù)階常微分方程組的一個(gè)有效方法。