閆羽媛,梁宗旗

(集美大學(xué)理學(xué)院，福建 廈門 361021)

0 引言

分?jǐn)?shù)階微積分理論及其數(shù)值逼近是數(shù)學(xué)的一個重要分支,它是傳統(tǒng)的整數(shù)階微積分理論的推廣。近年來，分?jǐn)?shù)階微分方程及其應(yīng)用得到了廣泛的關(guān)注,其主要歸因于分?jǐn)?shù)階微積分理論自身的迅速發(fā)展及其在數(shù)學(xué)、物理、信號和圖像處理等各學(xué)科中的廣泛應(yīng)用。分?jǐn)?shù)階微分方程是廣義非整數(shù)階的微分方程,它能獲取在時間和空間上具有冪律內(nèi)存內(nèi)核的非局部關(guān)系,為描述不同物質(zhì)的記憶和繼承性質(zhì)提供了強有力的工具。研究這類方程具有明確的物理背景,同時也開啟了分?jǐn)?shù)階微分和積分方程理論方面的研究。分?jǐn)?shù)階微分方程適用于描述現(xiàn)實世界中具有記憶以及遺傳性質(zhì)的物理行為或材料問題,因而更廣泛地應(yīng)用于反常擴(kuò)散、粘性力學(xué)、流體力學(xué)、信號處理等領(lǐng)域。在實際求解中,分?jǐn)?shù)階微分方程的解析解很難顯式給出,即使是線性分?jǐn)?shù)階方程的解析解，也大多含有一些特殊函數(shù),而這些函數(shù)的計算本身相當(dāng)困難。所以，如何構(gòu)造高效的數(shù)值方法來模擬分?jǐn)?shù)階微分方程是當(dāng)今最熱門的研究課題。關(guān)于分?jǐn)?shù)階導(dǎo)數(shù),現(xiàn)在有幾種不同的定義,如空間分?jǐn)?shù)階Laplace算子[1-2]、Riesz空間分?jǐn)?shù)階導(dǎo)數(shù)[3]等。大多數(shù)文獻(xiàn)中涉及的Riesz分?jǐn)?shù)階導(dǎo)數(shù)均為Riemann-Liouville框架下的分?jǐn)?shù)階導(dǎo)數(shù)。但Pandey等[4]指出,Caputo定義能夠避免一些分?jǐn)?shù)階導(dǎo)數(shù)存在的質(zhì)量守恒誤差、超奇異反常積分、常數(shù)的導(dǎo)數(shù)非零、初始條件中的分?jǐn)?shù)階導(dǎo)數(shù)無直觀的物理意義等問題。分?jǐn)?shù)階導(dǎo)數(shù)定義為在Caputo意義下的Riesz分?jǐn)?shù)階導(dǎo)數(shù),也稱為Riesz-Caputo分?jǐn)?shù)階導(dǎo)數(shù)[5]。

Caputo分?jǐn)?shù)階導(dǎo)數(shù)分為左側(cè)導(dǎo)數(shù)和右側(cè)導(dǎo)數(shù)。左側(cè)Caputo分?jǐn)?shù)階導(dǎo)數(shù)的離散主要是L1和L2離散：Murio[6]對α階(0<α<1)Caputo型時間分?jǐn)?shù)階低擴(kuò)散方程作L1離散建立差分方程,其截斷誤差為O(Δt);Lin等[7]也類似地利用L1離散,并將精度提高到O(Δt2-α);Sun等[8]構(gòu)造了L1-2格式進(jìn)行離散,并證明其收斂階O(Δt3-α);Du等[9]改進(jìn)L1-2格式為一致的L2格式,對α階(1<α<2)左側(cè)Caputo分?jǐn)?shù)階導(dǎo)數(shù)進(jìn)行離散,得到了一致的O(Δt4-α)精度。杜瑞蓮等[10]對α階(0<α<1)右側(cè)Caputo分?jǐn)?shù)階導(dǎo)數(shù)作L2-1插值逼近,并證明其收斂階為O(Δt3-α)。為了得到Riesz-Caputo分?jǐn)?shù)階導(dǎo)數(shù)的高階格式,需要研究α階(1<α<2)右側(cè)Caputo分?jǐn)?shù)階導(dǎo)數(shù)的數(shù)值逼近格式。本文主要構(gòu)造了α階(1<α<2)右側(cè)Caputo分?jǐn)?shù)階導(dǎo)數(shù)的一種高階數(shù)值逼近格式,通過引入新變量,先用L2-1格式構(gòu)造數(shù)值離散,進(jìn)一步改善得到一致的L2逼近格式,并給出其嚴(yán)格的誤差估計,證明了收斂階為O(Δt4-α)。

1 右側(cè)Caputo分?jǐn)?shù)階導(dǎo)數(shù)的L2-1格式

根據(jù)右側(cè)Caputo分?jǐn)?shù)階導(dǎo)數(shù)的定義知,當(dāng)α=n-1時，該導(dǎo)數(shù)退化為一般形式的整數(shù)階導(dǎo)數(shù)。為方便起見,本文只討論1<α<2的情況。

假設(shè)g(t)∈C3[t,b],記t=tk

(1)

將式(1)中區(qū)間[tj-1,tj](j∈[k+1,N-1])上用L2插值P2,jg(t)來近似代替g(t),區(qū)間[tN-1,tN]上用L1插值P1,Ng(t)來近似代替g(t),P2,jg(t)=g(tj-1)(tj-t)(tj+1-t)/(2Δt2)+g(tj)(t-tj-1)(tj+1-t)/Δt2+g(tj+1)(t-tj)(tj-1-t)/(2Δt2),P1,Ng(t)=g(tN-1)(tN-t)/Δt+g(tN)(t-tN-1)/Δt,對應(yīng)的截斷誤差分別為:g(t)-P2,jg(t)=g?(ζi)(t-tj-1)(t-tj)(t-tj+1)/6,t∈[tj-1,tj],ζi∈(tj-1,tj+1),g(t)-P1,Ng(t)=g″(ξ)(t-tN-1)(t-tN)/2,t∈[tN-1,tN],ζ∈(tN-1,tN)。

式(1)中右側(cè)α階(1<α<2)Caputo分?jǐn)?shù)階導(dǎo)數(shù)可以寫為:

(2)

其中:

(3)

注意到,

(4)

(5)

(6)

其中:bl=l2-α-(l-1)2-α，1≤l≤N-k;el=((l-1)3-α-l3-α)/(3-α)+(3(l-1)2-α-l2-α)/2;l≥1。

將式(4)～式(6)代入式(2)中得，

(7)

引理1當(dāng)N≥k+3時,對于任意的α(1<α<2),式(7)中系數(shù)wl有下列性質(zhì):1)wN-1>|wN-2|;2)wl>0，l≠N-2;3)wN-1>wN-3;4)wN-3≥wN-4≥…≥wk。

證明根據(jù)wl的定義及ej、bj的性質(zhì),易得1)～3)。下面證明性質(zhì)4)。

當(dāng)k+1≤l≤N-3時,wl=eN-l-1+bN-l-1-eN-l=[(N-l-2)3-α-2(N-l-1)3-α+(N-l)3-α]/(3-α)+[(N-l-2)2-α-2(N-l-1)2-α+(N-l)2-α]/2。

令w(x)=[(x-1)3-α-2x3-α+(x+1)3-α]/(3-α)+[(x-1)2-α-2x2-α+(x+1)2-α]/2,則有w′(x)=[(x-1)2-α-2x2-α+(x+1)2-α]+(2-α)[(x-1)1-α-2x1-α+(x+1)1-α]/2。

令h(x)=x2-α+(2-α)x1-α/2,則w′(x)=h(x-1)-2h(x)+h(x+1)。由于h″(x)=(2-α)(1-α)x-α+(2-α)(1-α)x-α-1/2=(2-α)(1-α)(x-α-(α/2)x-α-1)<0,故h(x)是一個凸函數(shù)。根據(jù)凸函數(shù)的性質(zhì)得h(x-1)-2h(x)+h(x+1)>0。特別地,當(dāng)x=N-l-1時,w′(N-l-1)=h(N-l-2)-2h(N-l-1)+h(N-l)>0。證畢。

為了研究式(7)中的系數(shù),引入引理2。

引理2假設(shè)f(t)∈C4[t，b],記fj=f(tj)。當(dāng)N≥k+3時,則有:

(8)

(9)

(10)

證明先討論k+1≤j≤N-2的情況。應(yīng)用Taylor展開式,

(11)

(12)

(13)

(14)

注意到,

(15)

(16)

(17)

(18)

(19)

(20)

利用引理2,當(dāng)N≥k+3時,式(7)變?yōu)?/p>

(21)

其中,

Rk+1/2=O(Δt4-α)+(Rk+Rk+1)/2。

(22)

在式(21)中,系數(shù)dk的具體形式如下:

|Rk+1/2|≤(1/Γ(2-α))[((6-α)(α-1)/(6(2-α)(3-α)(4-α))-1/12)

(23)

證明由式(3)有,

(24)

其中,

(6-α)/(6(2-α)(3-α)(4-α))f(4)(ζk)Δt4-α,ζk∈(tk,tk+2)。

(25)

(26)

(27)

2 α階(1<α<2)右側(cè)Caputo分?jǐn)?shù)階導(dǎo)數(shù)的L2插值逼近

定理1中的式(23),由于在區(qū)間[tN-1,tN]上采用L1插值,其整體誤差不是一致的O(Δt4-α)格式。實際計算中該區(qū)間上的誤差對整體誤差不會產(chǎn)生較大的影響,為此本節(jié)將利用L2在區(qū)間[tN-1,tN]上作插值逼近,使其整體誤差得到一致的O(Δt4-α)精度。

假設(shè)函數(shù)f(t)在區(qū)間[tk,tN+1]上有定義,可以在每一個小區(qū)間[tj-1,tj](j≥k+1)上均利用二次插值多項式P2,jg(t)逼近函數(shù)g(t)。

(28)

其中，

(29)

(30)

s2-αds+[(N-k-1)3-α-(N-k)3-α]/2|≤C(N-k)2-α，C是一個正常數(shù)?？梢缘玫?/p>

(31)

其中,

(6-α)/(6(2-α)(3-α)(4-α))f(4)(ζk)Δt4-α,ζk∈(tk,tk+2)，

(32)

(33)

3 結(jié)論

本文先利用L2-1格式對α階(1<α<2)右側(cè)Caputo分?jǐn)?shù)階導(dǎo)數(shù)進(jìn)行插值逼近,并討論了其系數(shù)性質(zhì)。為了改善由L1插值在區(qū)間[tN-1,b]上帶來的缺陷,進(jìn)一步提出L2插值逼近格式,通過增加一定的限制條件,得到了一致的O(Δt4-α)階精度,并分別對二者的截斷誤差進(jìn)行估計?；谏鲜鐾茖?dǎo),本文提供的是一種對α階(1<α<2)右側(cè)Caputo分?jǐn)?shù)階導(dǎo)數(shù)的高階差值逼近格式。下一步的研究可以將這種高階差分格式應(yīng)用到求解其他含有Caputo分?jǐn)?shù)階導(dǎo)數(shù)的差分方程并進(jìn)行誤差估計,除此之外，還可以構(gòu)造α階(1<α<2)Riesz-Caputo分?jǐn)?shù)階導(dǎo)數(shù)的高精度數(shù)值逼近格式,以更高效地解決一些時間或時空分?jǐn)?shù)階微分方程。