Caputo型對流方程的間斷伽遼金有限元方法

王 震

(江蘇大學 數學科學學院， 江蘇 鎮江 212013)

0 引言

分數階微積分算子由于其特有的非局部性，能夠對具有記憶或者遺傳特性的材料和過程進行更為準確的描述和刻畫，而受到了越來越多學者的關注[1-3]。但也正是這個特性使得分數階微分方程的解析解很難求得，因此利用數值方法求解分數階微分方程成為一個有重要意義的研究課題。考慮如下的時間分數階對流方程：

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

(1)

初始條件為：

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

邊界條件為：

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

0<α<1

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

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

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

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

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

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

(2)

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

1 預備知識

本節將給出數值分析中一些常用的記號、定義和投影，并且假設C>0表示一個與時間步長和空間步長均無關的有界常數，不同處含義可能不同。首先，對于任意非負整數s，記Hs(Ω)為定義在區間Ω上的標準Sobolev空間，其范數為：

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

令

函數u和v在區間Ω上的L2內積以及相應的L2范數定義如下：

N表示網格單元總數。單元中心和單元長度分別記為：

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)表示定義在網格單元Ii上的k次多項式函數空間。由于函數u在相鄰的2個單元交界處會有左右逼近的2個值，因此定義：

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

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

(3)

式中：C>0是一個與空間步長h無關的常數，且：

2 主要結果

本節構造Caputo型對流方程(1)的全離散DG有限元格式，并分析格式的穩定性、收斂性和誤差估計。對于給定的有限時間T和正整數M，記：

tn=T(n/M)r

式中：r≥1表示分級網格參數，n=0,1,2,…,M表示網格節點。將有限時間區間[0,T]劃分為：

0=t0

非均勻時間步長τn=tn-tn-1,n=1,2,…,M。容易看出，當r=1時，網格是均勻剖分。為書寫方便，引入如下記號：

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

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

為了對時間分數階導數進行離散，需要介紹Caputo分數階導數在非均勻網格上的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

對于n=0,1,…,M-1,有：

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

(5)

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

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

(6)

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

(7)

(8)

2.1 穩定性分析

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

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

如果最大時間步長τM滿足：

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

那么有如下不等式成立：

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

n=0,1,…,M-1

(9)

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

再利用引理1進一步得到：

(10)

由式(5)可知：

(11)

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

證明完畢。

2.2 誤差估計

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

引理3[13]設σ=1-α/2，則對任意的函數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對u(s)的二次插值多項式。

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

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

最大時間步長τM滿足：

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

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

證明令

(12)

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

若記

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

再利用投影算子的定義進一步得到[15]：

(13)

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

(14)

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

此處用到了不等式：

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

因此，有：

結合上式以及n=0的情形可得：

(15)

進一步地，由引理3和引理4可以推出：

(16)

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

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

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

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

3 數值實驗

本節通過一個數值例子來驗證格式(7)的精度和有效性。

例1考慮下面的Caputo型對流方程：

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

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

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

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

4 結論

1) 構造了一個新的求解Caputo型時間分數階對流方程的有限元逼近格式，此格式在時間方向上比已有算法的收斂精度高。

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

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