非結構網格下2D Riesz分數階方程的Galerkin有限元方法

卜瑋平

(湘潭大學數學與計算科學學院,湖南 湘潭 411105)

1 引言

最近,分數階微分方程數值方法的研究越來越受關注.有限差分方法(FDM)是求解分數階微分方程最常見的數值方法之一.文獻[1]是研究FDM求解分數階微分方程數值方法的先驅.隨后,利用有限差分方法人們對各種分數階微分方程進行了求解[2-5].譜方法(SM)和譜元方法(SEM)是求解分數階微分方程的重要方法,它們通常具有很高的收斂精度.目前,SM和SEM求解分數階微分方程的工作有文獻[6-8]等.相比于FDM方法和SM(或SEM)方法,有限元方法(FEM)的主要特征是它能夠較易處理復雜區域且對解的光滑性要求較低.文獻[9]首次嘗試用有限元方法求解分數階微分方程,并給出了有限元逼近的理論框架.此后,文獻[10]建立了有限元求解一維和二維分數階對流彌散方程的理論框架.隨后,在有限元求解分數階微分方程方面涌現了越來越多的工作,其中包括文獻[11-15]等.

近年來,關于二維分數階微分方程的有限元方法也有一些研究成果.文獻[10]考慮了基于分數階方向導數的二維分數階對流-彌散方程的有限元方法.通過使用矩陣轉換技術,文獻[16]討論了基于分數階拉普拉斯算子的二維分數階擴散方程有限元方法.文獻[17]發展了基于分數拉普拉斯算子的分數階擴散方程自適應有限元方法.為了求解二維Risez/Riemann-Liouville分數階擴散方程(2DRFDE/2DRLFDE),基于一致三角網格剖分,文獻[18-19]考慮了Galerkin有限元方法.然而,對于不規則區域上2DRFDE/2DRLFDE的有限元方法,通常需要采用非結構網格.基于非結構網格,文獻[20]利用三角形上定義的Lagrange多項式,建立了2DRLFDE的間斷Galerkin方法.文獻[21]考慮了非線性2DRFDE的有限元方法,并描述了有限元方法的實現.文獻[22]討論了二維時空分數階波方程在非規則凸域上的有限元方法.

盡管在文獻[20-22]中已有非結構網格下有限元求解分數階微分方程的工作,但是這些工作使用的是相同的有限元實現技巧.值得注意,上述文中提到的有限元實現方法有不足之處.因此,在這篇文章中將提出一些新的技巧,以改進現有的有限元實現方法.本文的主要貢獻如下:首先,對于剛度矩陣的計算降低了計算花銷,原來的計算花銷為O(Ne3),現在的花銷為O(Ne2),Ne為總剖分單元數;其次,提高了剛度矩陣的元素的精度,對于三角形單元的高斯積分,不同于現有的計算方法,高斯積分區域為被積函數的非零的區域,這比由以往方法得到的剛度矩陣元素更精確;第三,將Riemann-Liouville導數轉化為Caputo導數,簡化了內積的計算.

本文的結構安排如下:在第2節,給出了分數階導數的定義、模型問題及有限元全離散格式的推導;在第3節,首先介紹了現有的剛度矩陣計算方法,然后詳細描述了有限元的實現方法,并與現有方法進行了比較;在第4節,給出了數值實驗來證明方法的有效性;最后,對本文進行了總結.

2 準備工作

令??R2,則x,y方向的左右Riemann-Liouville分數階導數定義如下:

其中γ,n?1<γ≤n,n∈N,a(y),b(y),c(x),d(x)定義如圖1.

圖1 關于?上a(y),b(y),c(x),d(x)的定義

進一步,定義x,y方向的γ(1,3,···)階 Riesz分數階導數如下

類似地,x,y方向γ階右Caputo分數階導數定義如下

在文獻[23]中已經討論了Riemann-Liouville分數階導數和Caputo分數階導數的等價關系,如果u(a(y),y)=0,u(x,c(x))=0,0<γ<1,則

本文考慮如下分數階擴散方程

對上述方程,首先考慮其變分形式和有限元全離散格式.對于任意γ≥0,記γ(?)為Hγ(?)的子集且其元素在 ? 外的零擴張屬于Hγ(R2).令V=α(?)∩β(?).于是,由文獻[18]中引理5可得如下變分問題:尋找u∈V使得

其中 (·,·):=(·,·)L2(?),F(v):=(f,v),且

為了得到上述變分問題的全離散格式,先將?進行剖分.令{Th}是?的一個正則的三角剖分,h為所有三角形單元的最大直徑.定義如下有限元空間

3 有限元方法的實現

為了計算全離散格式(14),需要計算剛度矩陣和荷載向量.由于分數階導數為非局部算子,與整數階微分方程相比,其主要區別在于分數階微分方程剛度矩陣的計算非常復雜.因此,方法實現的重點將放在剛度矩陣的計算上.

這里?e為三角形單元e的面積,(xi,yi),(xj,yj),(xk,yk)是對應頂點i,j,k的坐標.

接下來,考慮B(uh,vh)的計算.令

這里S={Sij}N×N為剛度矩陣,其元素為

3.1 已有的實現方法

目前,文獻[20-22]已經考慮了Sij的計算.然而,應該注意這些文章關于Sij的計算方法是類似的.這里簡單對文獻[20-22]中Sij的計算方法進行描述.由于(17)中四個內積的計算具有相似性,因此僅僅以為例進行討論.考慮分數階算子的非局部性并利用高斯積分可得

這里Ge為三角形單元e上的高斯點,ωl是高斯點(xl,yl)對應的權重.

上面(18)式中的方法可以用來計算Sij,然而該方法有三點不足之處.令Ne為剖分三角形單元的總數.首先,對于Sij的計算需要考慮如下積分在每個三角形單元的高斯積分

因此,為了得到剛度矩陣的一個元素Sij,上述積分需要進行Ne次,而為了獲得整個剛度矩陣,需要計算次,這將使得剛度矩陣的計算花銷隨Ne的增漲而快速增漲.其次,對于下列積分

考慮被積函數在三角形單元上的非零區域.由于分數階算子具有非局部性,顯然存在滿足下列情形的三角形單元:被積函數在三角形的某些部分為零,在其他部分非零.因此,被積函數在滿足上述條件的三角形單元上為間斷函數.如果對這些單元,在整個三角形上運用高斯積分勢必達不到數值積分的相應精度.圖2給了描述上述情形的例子(圖形(a)描述了的非零區域,其中?i為節點i的基函數,?k為節點k的基函數;圖形(b)描述了被積函數的非零區域,它為支集的交集).

圖2 的非零區域和 的非零區域

圖3 在(xp,yp)的計算結果

計算?i(x,y)的左Riemann-Liouville分數階導數在(xp,yp)點的值

3.2 實現方法

針對已有方法的不足,設計一種新的求解二維Riesz分數階擴散方程的有限元實現方法.因為,顯然

因此,為了計算(16)中剛度矩陣S,根據對稱性僅需要計算Sx和Sy.

圖4 Lagrange線性基函數?i(x,y),?j(x,y)的支集

為了計算積分Ii,已有的方法是在整個三角形單元ei上使用高斯積分.然而,上面已經提到這樣做將降低高斯積分的效果,因為分數階算子具有非局部性,的非零區域在某些單元可能只占有一部分(即在該單元為間斷函數),如圖5.

圖5 的支集

圖6 的支集

注意到,對于?i(x,y),由 (9)式可得

于是利用高斯積分有

這里e為的非零區域,ai,?e定義在(15)式中.

假設?j(x,y)及其積分路徑被定義在圖3中,對于高斯點(xp,yp),有

顯然(30)式比(21)式更容易計算.

4 數值實驗

本節給出一個數值例子來驗證方法的有效性.

例 4.1在模型方程(10)中,取P=Q=2,A=4.考慮如下兩種情形:(a)假設考慮問題區域為[0,1]×[0,1],其精確解為u=100x2(1?x)2y2(1?y)2;(b)假設考慮問題區域為,其精確解為,相應的右端函數分別定義如下

當選擇不同的α和β時,表1-表4分別列出了情形(a)和情形(b)所得的數值結果.

可以看出,所得誤差的收斂率是最優收斂率.

表1 基于Lagrange線性多項式與α=0.6,β=0.6計算情形(a)的數值誤差與收斂率

表2 基于Lagrange線性多項式與α=0.7,β=0.8計算情形(a)的數值誤差與收斂率

表3 基于Lagrange線性多項式與α=0.7,β=0.6計算情形(b)的數值誤差與收斂率

表4 基于Lagrange線性多項式與α=0.6,β=0.8計算情形(b)的數值誤差與收斂率

5 總結

本文研究了非結構網格下利用Lagrange線性基函數求解2D Riesz分數階擴散方程的有限元方法實現.首先,描述了現有有限元全離散格式的實現方法,并指出了現有方法的不足之處.隨后,針對這些缺點設計了一種新的實現方法,提高了有限元方法的計算效率和剛度矩陣的精度.最后,給出了數值算例,數值結果驗證了本文所提方法的有效性.