卜瑋平

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

1 引言

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

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

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

本文的結(jié)構(gòu)安排如下:在第2節(jié),給出了分?jǐn)?shù)階導(dǎo)數(shù)的定義、模型問題及有限元全離散格式的推導(dǎo);在第3節(jié),首先介紹了現(xiàn)有的剛度矩陣計算方法,然后詳細(xì)描述了有限元的實現(xiàn)方法,并與現(xiàn)有方法進(jìn)行了比較;在第4節(jié),給出了數(shù)值實驗來證明方法的有效性;最后,對本文進(jìn)行了總結(jié).

2 準(zhǔn)備工作

令??R2,則x,y方向的左右Riemann-Liouville分?jǐn)?shù)階導(dǎo)數(shù)定義如下:

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

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

進(jìn)一步,定義x,y方向的γ(1,3,···)階 Riesz分?jǐn)?shù)階導(dǎo)數(shù)如下

類似地,x,y方向γ階右Caputo分?jǐn)?shù)階導(dǎo)數(shù)定義如下

在文獻(xiàn)[23]中已經(jīng)討論了Riemann-Liouville分?jǐn)?shù)階導(dǎo)數(shù)和Caputo分?jǐn)?shù)階導(dǎo)數(shù)的等價關(guān)系,如果u(a(y),y)=0,u(x,c(x))=0,0<γ<1,則

本文考慮如下分?jǐn)?shù)階擴(kuò)散方程

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

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

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

3 有限元方法的實現(xiàn)

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

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

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

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

3.1 已有的實現(xiàn)方法

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

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

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

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

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

圖2 的非零區(qū)域和 的非零區(qū)域

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

計算?i(x,y)的左Riemann-Liouville分?jǐn)?shù)階導(dǎo)數(shù)在(xp,yp)點的值

3.2 實現(xiàn)方法

針對已有方法的不足,設(shè)計一種新的求解二維Riesz分?jǐn)?shù)階擴(kuò)散方程的有限元實現(xiàn)方法.因為,顯然

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

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

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

圖5 的支集

圖6 的支集

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

于是利用高斯積分有

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

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

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

4 數(shù)值實驗

本節(jié)給出一個數(shù)值例子來驗證方法的有效性.

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

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

可以看出,所得誤差的收斂率是最優(yōu)收斂率.

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

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

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

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

5 總結(jié)

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