王 震

(江蘇大學(xué) 數(shù)學(xué)科學(xué)學(xué)院， 江蘇 鎮(zhèn)江 212013)

0 引言

分?jǐn)?shù)階微積分算子由于其特有的非局部性，能夠?qū)哂杏洃浕蛘哌z傳特性的材料和過程進(jìn)行更為準(zhǔn)確的描述和刻畫，而受到了越來越多學(xué)者的關(guān)注[1-3]。但也正是這個(gè)特性使得分?jǐn)?shù)階微分方程的解析解很難求得，因此利用數(shù)值方法求解分?jǐn)?shù)階微分方程成為一個(gè)有重要意義的研究課題?？紤]如下的時(shí)間分?jǐn)?shù)階對(duì)流方程：

x∈Ω=(a,b),t∈(0,T]

(1)

初始條件為：

u(x,0)=u0(x),x∈Ω

邊界條件為：

u(a,t)=0,t∈(0,T]

0<α<1

初值u0(x)和源項(xiàng)f(x,t)是已知函數(shù)，T是給定的有限時(shí)間，Γ(·)是Gamma函數(shù)。容易看出，問題 (1) 是一個(gè)時(shí)間分?jǐn)?shù)階雙曲型方程。為了數(shù)值計(jì)算的方便，假設(shè)：

u(b,t)=0,t∈(0,T]

到目前為止，關(guān)于整數(shù)階對(duì)流方程數(shù)值方法的研究已經(jīng)取得了很多成果[4]，而這類方程在分?jǐn)?shù)階導(dǎo)數(shù)意義下的發(fā)展卻很緩慢，直到近幾年才有零星結(jié)果。Li等[5]考慮了3類方程(即分?jǐn)?shù)階對(duì)流擴(kuò)散方程、分?jǐn)?shù)階對(duì)流占優(yōu)擴(kuò)散方程和分?jǐn)?shù)階對(duì)流占優(yōu)分?jǐn)?shù)階擴(kuò)散方程)的建模。隨后，他們?cè)谖墨I(xiàn)[6]中建立了求解空間分?jǐn)?shù)階對(duì)流方程的有限差分方法。最近，Li和Wang[7-8]研究了一維和二維Caputo型分?jǐn)?shù)階對(duì)流方程的間斷伽遼金(discontinuous Galerkin,DG)有限元方法。注意到時(shí)間分?jǐn)?shù)階偏微分方程的解在初始時(shí)刻可能會(huì)表現(xiàn)出一些弱正則性，為了保證格式在時(shí)間方向上的一致高精度，他們采用非均勻網(wǎng)格上的L1差分方法離散方程的時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)，空間方向使用DG有限元方法近似，從而構(gòu)造了全離散數(shù)值格式，并對(duì)格式的穩(wěn)定性、收斂性和誤差估計(jì)進(jìn)行了詳細(xì)討論。從結(jié)果可以看出，文獻(xiàn)[7-8]中構(gòu)造的數(shù)值格式在時(shí)間方向上的精度最多只能達(dá)到2-α階。

將繼續(xù)研究解在初始時(shí)刻具有弱正則性的Caputo型對(duì)流方程 (1) 的有效數(shù)值解法，即假設(shè)問題(1)的解滿足：

‖?lu(x,t)/?tl‖H2(Ω)≤Cu(1+tα-l),

l=0,1,2,3,t∈(0,T]

(2)

式中：Cu>0是一個(gè)與u的正則性有關(guān)的常數(shù)。首先使用非均勻網(wǎng)格上的Alikhanov公式離散Caputo時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)，空間方向使用DG有限元方法逼近，進(jìn)而給出全離散數(shù)值格式。然后對(duì)格式的穩(wěn)定性、收斂性和誤差估計(jì)進(jìn)行嚴(yán)格的證明，并通過數(shù)值算例進(jìn)一步驗(yàn)證理論分析的正確性和算法的有效性。結(jié)果表明，通過選取合適的分級(jí)網(wǎng)格參數(shù)，格式在時(shí)間方向的誤差可以達(dá)到最優(yōu)收斂階2，比之前的2-α階要高，空間方向是k+1階收斂的，其中k是分段多項(xiàng)式的次數(shù)。

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

本節(jié)將給出數(shù)值分析中一些常用的記號(hào)、定義和投影，并且假設(shè)C>0表示一個(gè)與時(shí)間步長和空間步長均無關(guān)的有界常數(shù)，不同處含義可能不同。首先，對(duì)于任意非負(fù)整數(shù)s，記Hs(Ω)為定義在區(qū)間Ω上的標(biāo)準(zhǔn)Sobolev空間，其范數(shù)為：

‖·‖s=‖·‖Hs(Ω)

令

函數(shù)u和v在區(qū)間Ω上的L2內(nèi)積以及相應(yīng)的L2范數(shù)定義如下：

N表示網(wǎng)格單元總數(shù)。單元中心和單元長度分別記為：

xi=(xi-1/2+xi+1/2)/2,hi=xi+1/2-xi-1/2

Vh={v∈L2(Ω)∶v|Ii∈Pk(Ii),i=1,2,…,N}

式中：k≥0,Pk(Ii)表示定義在網(wǎng)格單元Ii上的k次多項(xiàng)式函數(shù)空間。由于函數(shù)u在相鄰的2個(gè)單元交界處會(huì)有左右逼近的2個(gè)值，因此定義：

記u在單元交界處的跳躍值為：

[|u|]i+1/2=u+|xi+1/2-u-|xi+1/2

(3)

式中：C>0是一個(gè)與空間步長h無關(guān)的常數(shù)，且：

2 主要結(jié)果

本節(jié)構(gòu)造Caputo型對(duì)流方程(1)的全離散DG有限元格式，并分析格式的穩(wěn)定性、收斂性和誤差估計(jì)。對(duì)于給定的有限時(shí)間T和正整數(shù)M，記：

tn=T(n/M)r

式中：r≥1表示分級(jí)網(wǎng)格參數(shù)，n=0,1,2,…,M表示網(wǎng)格節(jié)點(diǎn)。將有限時(shí)間區(qū)間[0,T]劃分為：

0=t0

非均勻時(shí)間步長τn=tn-tn-1,n=1,2,…,M。容易看出，當(dāng)r=1時(shí)，網(wǎng)格是均勻剖分。為書寫方便，引入如下記號(hào)：

tn+σ=tn+στn+1,un+σ=u(x,tn+σ)

un,σ=σun+1+(1-σ)un,n=0,1,…,M-1

為了對(duì)時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)進(jìn)行離散，需要介紹Caputo分?jǐn)?shù)階導(dǎo)數(shù)在非均勻網(wǎng)格上的Alikhanov逼近公式為[10]：

(4)

n≥1, 0≤j≤n-1

n≥1, 0≤j≤n-1

令

▽tuj+1=uj+1-uj, 0≤j≤n, 0≤n≤M-1

對(duì)于n=0,1,…,M-1,有：

通過參閱文獻(xiàn)[11]，可以得到下面的不等式：

(5)

式中：ωβ(t)=tβ-1/Γ(β),πA>0是一個(gè)常數(shù)。

方程(1)在t=tn+σ時(shí)刻的變分形式為：尋找u∈H1(Ω)，使得：

(6)

式中：v是空間H1(Ω)中任意的檢測(cè)函數(shù)。

(7)

(8)

2.1 穩(wěn)定性分析

本小節(jié)考慮全離散DG格式(7)的穩(wěn)定性，首先介紹如下2個(gè)引理。

式中：{φn+1,ψn+1|0≤n≤M-1}是2個(gè)非負(fù)序列。

如果最大時(shí)間步長τM滿足：

τM≤(2πAΓ(2-α)Λ)-1/α

那么有如下不等式成立：

引理2[12]如果σ=1-α/2，那么對(duì)任意的函數(shù)u成立不等式：

n=0,1,…,M-1

(9)

根據(jù)引理2和Cauchy-Schwarz不等式可得：

再利用引理1進(jìn)一步得到：

(10)

由式(5)可知：

(11)

將式(11)代入到式(10)中，可得：

證明完畢。

2.2 誤差估計(jì)

本小節(jié)給出方程(1)的全離散DG格式的最優(yōu)誤差估計(jì)。下面先介紹2個(gè)非常重要的引理。

引理3[13]設(shè)σ=1-α/2，則對(duì)任意的函數(shù)u(t)∈C3(0,T]，以下不等式成立：

0≤n≤M-1

式中：

n=1,2,…,M-1

n=1,2,…,M-1

2≤i≤n≤M-1

此處I2,1u(s)表示在ts-1、ts和ts+1對(duì)u(s)的二次插值多項(xiàng)式。

引理4[13-14]如果u∈C[0,T]∩C3(0,T]且滿足條件(2)，那么有：

u(x,t)∈L∞(0,T;Hk+1(Ω))

最大時(shí)間步長τM滿足：

令σ=1-α/2，則有如下估計(jì)：

式中：有界常數(shù)C>0不依賴于時(shí)間步長M和空間步長h。

證明令

(12)

由方程(6)和(7)可得誤差方程：

若記

然后將式(12)代入上述誤差方程可得：

再利用投影算子的定義進(jìn)一步得到[15]：

(13)

由引理2和Cauchy-Schwarz不等式可得：

(14)

根據(jù)式(3)，有下面的估計(jì)式成立：

此處用到了不等式：

τn+1≤CTM-rnr-1,n=0,1,…,M-1

因此，有：

結(jié)合上式以及n=0的情形可得：

(15)

進(jìn)一步地，由引理3和引理4可以推出：

(16)

將式(16)代入式(14)，然后利用Cauchy-Schwarz不等式可得：

再根據(jù)式(5)以及引理1便可推出：

Chk+1≤CM-min{rα,2}+Chk+1

最后利用插值性質(zhì)(3)和三角不等式即可完成定理的證明。證明完畢。

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

本節(jié)通過一個(gè)數(shù)值例子來驗(yàn)證格式(7)的精度和有效性。

例1考慮下面的Caputo型對(duì)流方程：

(x,t)∈(0,1)×(0,1]

(17)

初始條件為：

u(x,0)=0,x∈(0,1)

邊界條件為：

u(0,t)=0,t∈(0,1]

式中：

2π(tα+t3)cos(2πx)

精確解為u(x,t)=(tα+t3)sin(2πx)。

容易看出，方程(17)的精確解在初始時(shí)刻有一定的弱正則性，因此可以利用全離散格式 (7) 來求解上述方程。首先，驗(yàn)證格式在時(shí)間方向上的精度。選取N(N=1 000)足夠大，使得空間方向的誤差不影響時(shí)間方向的精度。表1和表2分別列出了當(dāng)分級(jí)網(wǎng)格參數(shù)r取1和(3-α)/α?xí)r，在不同的分?jǐn)?shù)階導(dǎo)數(shù)α下的時(shí)間方向L2模誤差和收斂階，結(jié)果表明全離散DG格式(7)在時(shí)間方向的收斂階為O(M-min{2-α,2})，這與定理2中的結(jié)果是一致的。如果選取r=(3-α)/α，那么格式在時(shí)間方向可以達(dá)到最優(yōu)的收斂階2。其次，驗(yàn)證格式在空間方向上的精度。固定時(shí)間方向的剖分為M=2 000，選取分級(jí)網(wǎng)格參數(shù)r=(3-α)/α，以k=1為例，從表3可以看出格式(7)在空間方向是k+1階收斂的。

表1 當(dāng)r=1,N=1 000,T=1時(shí)，格式(7)時(shí)間方向L2模誤差及收斂階

表2 當(dāng)r=(3-α)/α,N=1 000,T=1時(shí)，格式(7)時(shí)間方向L2模誤差及收斂階

表3 當(dāng)k=1,M=2 000,T=1時(shí)，格式(7)空間方向L2模誤差及收斂階

4 結(jié)論

1) 構(gòu)造了一個(gè)新的求解Caputo型時(shí)間分?jǐn)?shù)階對(duì)流方程的有限元逼近格式，此格式在時(shí)間方向上比已有算法的收斂精度高。

2) 證明了全離散數(shù)值格式的穩(wěn)定性、收斂性和誤差估計(jì)。通過數(shù)值算例檢驗(yàn)了當(dāng)前算法的有效性以及數(shù)值格式的收斂精度，為以后進(jìn)一步研究高維情形歸納出一種基本的分析框架。

3) 所提出的數(shù)值算法主要適用于方程是線性的情形，將其推廣到非線性方程并且證明格式的收斂性是一個(gè)難點(diǎn)。接下來將考慮構(gòu)造合適的離散分?jǐn)?shù)階Gronwall不等式來克服這一困難。