二維分數(shù)階擴散方程的最小二乘混合有限元計算

張洪光 陳煥貞

( 山東師范大學數(shù)學與統(tǒng)計學院, 250358, 濟南 )

1 引 言

眾所周知, 分數(shù)階擴散方程描述了反常擴散、混沌動力學、粘彈性問題等重要物理現(xiàn)象[1-3]. 大量的實驗結果表明, 分數(shù)階擴散模型給出了對這些實際問題較二階擴散方程更為精確的刻畫. 但一般的分數(shù)階擴散方程難以求得其解析解, 只能采用恰當?shù)臄?shù)值方法進行離散求解. 為此, 人們提出了大量的數(shù)值模擬方法. 這些數(shù)值方法主要分為兩類: 一類是基于差分框架下的數(shù)值方法, 如文獻[4]、[5],另一類是基于Galerkin框架下的數(shù)值方法,如文獻[6]、[7]、[8]、[9]等.但由于分數(shù)階算子非局部性質(zhì)所導致的計算復雜性, 大部分數(shù)值方法都是圍繞一維分數(shù)階擴散問題進行討論的. 本文針對二維分數(shù)階擴散問題, 借鑒文獻[10]中對一維問題提出的最小二乘思想與Galerkin框架, 通過引入擴散通量作為中間變量, 提出了相應的最小二乘混合變分形式, 證明了極小問題與變分問題的等價性. 據(jù)此, 建立了數(shù)值求解二維分數(shù)階擴散問題的最小二乘混合有限元離散格式, 證明了離散解的存在性與唯一性, 并通過數(shù)值算例說明離散格式具有較好的逼近性質(zhì).

2 分數(shù)階擴散問題與變分形式

對0<β<1, 我們考慮如下2-β階二維擴散方程

(1)

這里,未知函數(shù)p表示擴散濃度,f(x,y)表示源項, 介質(zhì)的擴散系數(shù)假定為常數(shù).為二維梯度算子.

與一階散度算子類似, 對一給定的向量函數(shù)u=(u1,u2)T, 我們定義1-β階散度算子為

1-β·u=0D1-βxu1+0D1-βyu2,

其中,0D1-βx,0D1-βy分別為Riemann-Liouville二維分數(shù)階導數(shù)算子[11],

在(1)中令u=-p, 則原方程改寫為下列一階系統(tǒng)，

u=-p(x,y), (x,y)∈Ω=[0,1]2，
1-β·u=f(x,y), (x,y)∈Ω=[0,1]2，
p(x,y)=0, (x,y)∈?Ω.

(2)

從而,上述一階系統(tǒng)的求解可轉(zhuǎn)化為下列極小問題: 求(p,u)∈P×U使得

(3)

其中

(4)

這里,P和U為適當選取的Sobolev空間以保證下節(jié)中定義的離散格式解的存在唯一性. 例如，

空間P與U的范數(shù)定義為

與(3)，(4)相應的變分形式可定義為求(p,u)∈P×U滿足

(1-β·u,1-β·v)=(f,1-β·v), ?v∈U,
(p,q)=(-u,q), ?q∈P.

(5)

對以上定義的極小問題(3)和變分形式(5), 我們有下列等價性結論.

定理1極小問題(3)與變分形式(5)是等價的.

證定義二次函數(shù)Φ(s)為

Φ(s)=J(u+sv,p+sq),?s≥0.

利用雙線性形式的定義(4)和范數(shù)的定義將二次函數(shù)Φ(s)展開為

如果(p,u)是極小問題(3)的解, 根據(jù)變分法原理, 我們得到如下結論,

Φ′(0)=(u+p,v+q)+(1-β·u-f,1-β·v)=0.

所以(p,u)是變分形式(5)的解.

反之,若(p,u)是變分形式(5)的解, ?s≥0，

3 離散格式

本節(jié)中, 我們將會給出有限元空間與離散格式, 并且證明離散解的存在性與唯一性.

對區(qū)域[0,1]×[0,1]進行一致矩形剖分, 剖分節(jié)點為0=x0

在該剖分基礎下, 我們定義有限元空間如下:

Uh:=span{φij,ψij}?U.

其中

對于真解p的逼近采用雙線性有限元空間Ph. 若記對應于內(nèi)點(xi,yj)的基函數(shù)為φij(x,y), 則

Ph:=span{φij}?P.

據(jù)此, 我們定義如下的混合有限元離散格式：求(ph,uh)∈Ph×Uh滿足

(1-β·uh,1-β·vh)=(f,1-β·vh), ?vh∈Uh,
(ph,qh)=(-uh,qh), ?qh∈Ph.

(6)

定理2離散格式(6)的解是存在唯一的.

為了將離散格式(6)的代數(shù)結構表達清楚, 我們把基函數(shù)φij,ψij,φij重新進行編號, 使得

φm=φij,m=(N+1)j+i+1;ψn=ψij,n=(N+1)i+j+1;φk=φij,k=(i-1)(N-1)+j;

其中,m,n=1,2,3,…,N(N+1),k=1,2,3,…,(N-1)2.

注意到uh∈Uh,ph∈Ph, 所以uh和ph可以表示為如下形式:

令vh依次取遍空間Uh中的基函數(shù),qh依次取遍空間Ph中的基函數(shù)， 經(jīng)計算后可將有限元離散格式(6)化為以下線性代數(shù)方程組，

A1X1=b,
A2X2=l,

其中

對m,n=1,2,…N(N+1),

C=(cmn),cmn=(1-β·φm,1-β·φn),
D=(dmn),dmn=(1-β·φm,1-β·ψn),
E=(emn),emn=(1-β·ψm,1-β·φn),
F=(fmn),fmn=(1-β·ψm,1-β·ψn);

A2=(amn),amn=(φm,φn);1-β·φm,f),bn=(1-β·ψn,f);φm),ln=(-uh,φn).

4 數(shù)值實驗

算例在方程(1)中, 取β=1/10, 源項f取為：

則方程(1)的真解為

p(x,y)=x2y2(1-x)2(1-y)2,

擴散通量為

u=(u1,u2)=(-(2x-6x2+4x3)y2(1-y)2,-(2y-6y2+4y3)x2(1-x)2).

我們對最小二乘混合離散格式(6)編程，并利用Matlab程序進行計算，計算結果由表1、圖1和圖2給出. 表1給出p-ph的H1與L2誤差及收斂階, 同時給出u-uh的L2誤差及收斂階. 表1說明p-ph的H1與L2誤差階分別達到1階與2階, 與插值逼近的誤差階相符合.u-uh的誤差階略低于插值逼近的誤差階.

表1 p-ph的H1與L2誤差及收斂階，u-uh的L2誤差及收斂階

圖1 p, u1和u2的圖像

圖2 ph, u1h和u2h的圖像

圖1、圖2分別給出了真解和與有限元逼近解的圖像, 從左至右依次為p,u1,u2以及ph,u1h,u2h.由此看出，最小二乘有限元解對真解具有較好的逼近性質(zhì).