張英晗，楊小遠(yuǎn)

(北京航空航天大學(xué) 數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院，北京100191)

分?jǐn)?shù)階偏微分方程模型在模擬許多復(fù)雜的物理現(xiàn)象中是非常有效的模型，在金融經(jīng)濟(jì)、半導(dǎo)體、生物、水文、控制理論等領(lǐng)域都有著重要的應(yīng)用［1-5］.例如分?jǐn)?shù)階擴(kuò)散方程可描述反常擴(kuò)散系統(tǒng)中的粒子輸運(yùn)現(xiàn)象［2-3］，空間分?jǐn)?shù)階偏微分方程可用來模擬超擴(kuò)散現(xiàn)象(系統(tǒng)中微粒子流的傳播速率比傳統(tǒng)的Brownian 運(yùn)動(dòng)預(yù)測模型預(yù)測的速度快［2］).

隨著分?jǐn)?shù)階微分方程理論研究的不斷完善，分?jǐn)?shù)階偏微分方程數(shù)值方法的研究引起了學(xué)者們越來越多的重視.最近幾年，出現(xiàn)了許多分?jǐn)?shù)階偏微分方程數(shù)值方法的成果.Liu 等［4］提出了一種直線法，把空間分?jǐn)?shù)階Fokker-Planck 方程轉(zhuǎn)化為一階微分系統(tǒng)進(jìn)而求其數(shù)值解.Meerschaert等［2-3］和Tadjeran 等［5-6］分別提出了空間分?jǐn)?shù)階偏微分方程的顯示Euler 方法、隱式Euler 方法和分?jǐn)?shù)階Crank-Nicolson(C-N)方法，進(jìn)行了穩(wěn)定性和收斂性分析.Liu 等［7］討論了Levy-Feller 對流彌散過程的隨機(jī)步和有限差分逼近.Zhuang 等［8］考慮了一類變階的分?jǐn)?shù)階對流擴(kuò)散方程的顯式和隱式Euler 格式，進(jìn)行了收斂分析.更多分?jǐn)?shù)階偏微分方程數(shù)值方法的文獻(xiàn)，可參考文獻(xiàn)［4-18］.

注意到很多分?jǐn)?shù)階微分方程差分方法的文獻(xiàn)處理的都是一維問題空間分?jǐn)?shù)階導(dǎo)數(shù)或者時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)，對高維問題以及空間時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)很少有人涉及［3，9］.受到以上研究工作的啟發(fā)，本文研究二維空間時(shí)間分?jǐn)?shù)階色散方程的有限差分格式.由于分?jǐn)?shù)階問題的隱式差分格式在每一時(shí)間步都需要處理非稀疏線性系統(tǒng)，于是隱式差分格式是計(jì)算密集型的.通過把偏微分方程差分方法中經(jīng)常用到的交替方向隱式差分格式推廣到分?jǐn)?shù)階偏微分方程，從而克服了這一問題.

1 問題描述

考慮二維空間時(shí)間分?jǐn)?shù)階色散方程:

式中:0 ＜t≤T;0 ＜x ＜Lx;0 ＜y ＜Ly，初邊值條件

從物理意義方面考慮，要求0 ＜α≤1、1 ＜β、γ≤2，并且u(0，y，t)=u(x，0，t)=0.這意味著在下邊界沒有示蹤劑泄漏.函數(shù)f(x，y，t)可用來表示源和匯.當(dāng)α =1、β =γ =2 時(shí)，方程式(1)還原為經(jīng)典的二維色散方程.以下假設(shè)d(x，y，t)≥0并且e(x，y，t)≥0，即流體是從下邊界到上邊界.在上述條件下方程式(1)存在唯一的光滑解［11］.

2 隱式差分格式

式中:bk=(k+1)1－α－k1－α，k=0，1，…，N.

式中:

分別令

由式(4)～式(6)，可以得到如下的隱式差分格式:

式中:0≤n≤N －1;1≤i≤Nx－1;1≤j≤Ny－1.式(7)又可表示為

初邊值條件分別離散為

式中:0≤n≤N－1;0≤i≤Nx;0≤j≤Ny.

引理1［2］系數(shù)滿足bk＞0，bk＞bk+1(k=0，1，…)并且

從式(4)～式(6)知，二維隱式差分格式式(8)與分?jǐn)?shù)階初邊值問題式(1)、式(2)是相容的，并且有局部截?cái)嗾`差O(Δt)+ O(Δx)+O(Δy).下面證明差分格式式(8)是穩(wěn)定的，從而由Lax 等價(jià)定理［12］在分?jǐn)?shù)階方程中的應(yīng)用［2］，式(8)是收斂的.分別令

則式(8)可表示為

式中:1≤i≤Nx－1;1≤j≤Ny－1;n=0，1，….令

引理2 || En||∞≤ ||E0||∞，n=0，1，….

證明 當(dāng)n=0 時(shí)，令

因此 ||E1||∞≤ ||E0||∞.

即 ||En+1||∞≤||E0||∞. 證畢

以上分析表明分?jǐn)?shù)階隱式差分格式式(8)是相容的并且無條件收斂的，因此由Lax 等價(jià)定理，有以下結(jié)論.

定理1 分?jǐn)?shù)階隱式差分格式式(8)是以O(shè)(Δt)+O(Δx)+O(Δy)無條件收斂的.

3 交替方向隱式差分格式

二維隱式差分格式式(8)是收斂的并且有局部誤差O(Δt)+O(Δx)+O(Δy).然而在每一時(shí)間步，需要處理(Nx－1)×(Ny－1)個(gè)未知量的非稀疏線性系統(tǒng)，于是隱式差分格式是計(jì)算密集型的.隨著網(wǎng)格的加細(xì)和空間維數(shù)的提高，計(jì)算復(fù)雜度變得非常龐大.為了解決這個(gè)問題，下面把偏微分方程差分方法中常用的交替方向隱式差分格式推廣到分?jǐn)?shù)階偏微分方程，即在每一時(shí)間步只求解一個(gè)方向.

定義分?jǐn)?shù)階偏差分算子:

則隱式差分格式式(8)可表示為

對于交替方向隱式差分格式，上述算子形式相應(yīng)變?yōu)榉较蚍蛛x乘積的形式:

這產(chǎn)生了一個(gè)額外的誤差擾動(dòng)

在計(jì)算上，上述形式的交替方向隱式差分格式可表示為如下的迭代形式.在時(shí)間層tn+1:

2)對每個(gè)固定的xi求解y 方向:

下面證明隱式差分格式的交替方向方法式(11)、式(12)是相容的和穩(wěn)定的，因此由Lax等價(jià)定理，也是收斂的.交替方向隱式差分格式的額外擾動(dòng)誤差式(10)的階為O(Δt2α)，因此由式(11)、式(12)定義的交替方向隱式差分格式是相容的，并且有局部截?cái)嗾`差O(Δtα)+O(Δx)+O(Δy).

引理3 ||En+1||∞≤ ||E0||∞，n=0，1，….

證明 當(dāng)n=0 時(shí)，令

根據(jù)引理1 和式(11)，有

因此 ||E*，1||∞≤|| E0||∞.

因此 ||E1||∞≤ ||E0||∞.

因此 ||En+1||∞≤|| E0||∞. 證畢

通過以上分析，由式(11)、式(12)定義的分?jǐn)?shù)階交替方向隱式差分方法是相容的和無條件穩(wěn)定的，因此由Lax 等價(jià)定理在分?jǐn)?shù)階方程中的應(yīng)用［2］，有如下結(jié)論.

定理2 分?jǐn)?shù)階交替方向隱式差分式(11)、式(12)是以O(shè)(Δtα)+O(Δx)+O(Δy)無條件收斂的.

注意到交替方向隱式差分格式是以階O(Δtα)+O(Δx)+O(Δy)無條件收斂的，隱式差分格式(8)是以O(shè)(Δt)+O(Δx)+O(Δy)無條件收斂的，而0 ＜α≤1，所以為了減少計(jì)算復(fù)雜度，犧牲了一些收斂精度.

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

考慮二維分?jǐn)?shù)階色散方程:

式中:(x，y)∈(0，1)2;0≤t≤T;初值u(x，y，0)=x2y2，Dirichlet 邊值

u(x，0，t)=u(0，y，t)=0

u(1，y，t)=(t2+1)y2

u(x，1，t)=(t2+1)x2

驗(yàn)證可知，方程式(13)的精確解為

u(x，y，t)=(t2+1)x2y2

表1 是交替方向隱式差分格式計(jì)算結(jié)果與精確解之間的最大誤差.最后一列誤差率表示上一行最大絕對誤差與本行最大絕對誤差的比率，即數(shù)值方法的收斂速度.可以看到最大絕對誤差是以速度2 幾乎線性下降的，這與理論分析結(jié)果是一致的.

表1 交替方向隱式差分格式最大誤差與誤差率Table 1 Maximum error and error rate of alternating directions implicit scheme

5 結(jié) 論

通過對空間分?jǐn)?shù)階導(dǎo)數(shù)的轉(zhuǎn)移Grünwald 差分近似，分別構(gòu)造了有界區(qū)域上二維空間時(shí)間分?jǐn)?shù)階色散方程的隱式差分格式和交替方向隱式差分格式，應(yīng)用數(shù)學(xué)歸納法證明了兩種格式的穩(wěn)定性和收斂性:

1)隱式差分格式是以O(shè)(Δt)+ O(Δx)+O(Δy)無條件收斂的.

2)交替方向隱式差分格式的收斂階為O(Δtα)+O(Δx)+O(Δy).

3)雖然交替方向隱式差分格式在時(shí)間方向上的收斂精度低于隱式差分格式，但在計(jì)算復(fù)雜度上考慮，交替方向隱式差分格式優(yōu)于隱式差分格式.

4)數(shù)值模擬結(jié)果與理論分析是一致的.

對于更廣泛的分?jǐn)?shù)階偏微分方程數(shù)值逼近問題，比如無界區(qū)域上分?jǐn)?shù)階偏微分方程的數(shù)值方法、分?jǐn)?shù)階偏微分方程的高階近似差分格式等問題，將在后續(xù)的工作中進(jìn)行研究.

References)

［1］ Arrarás A，Portero L，Jorge J C.Convergence of fractional step mimetic finite difference discretizations for semilinear parabolic problems［J］.Applied Numerical Mathematics，2010，60(4):473-485.

［2］ Meerschaert M，Tadjeran C.Finite difference approximations for fractional advection-dispersion flow equations［J］.Journal of Computational and Applied Mathematics，2004，172(1):65-77.

［3］ Meerschaert M，Scheffler H，Tadjeran C.Finite difference methods for two-dimensional fractional dispersion equation［J］.Journal of Computational Physics，2006，211(2):249-261.

［4］ Liu F，Anh V，Turner I.Numerical solution of the space fractional Fokker-Planck equation［J］.Journal of Computational and Applied Mathematics，2004，166(1):209-219.

［5］ Tadjeran C，Meerschaert M，Scheffler H.A second-order accurate numerical approximation for the fractional diffusion equation［J］.Journal of Computational Physics，2006，213(2):205-213.

［6］ Tadjeran C，Meerschaert M.A second-order accurate numerical method for the two-dimensional fractional diffusion equation［J］.Journal of Computational Physics，2007，220(2):813-823.

［7］ Liu Q，Liu F，Turner I，et al.Approximation for the Levy-Feller advection-dispersion process by random walk and finite difference method［J］.Journal of Computational Physics，2007，222(1):57-70.

［8］ Zhuang P，Liu F，Anh V，et al.Numerical methods for the variable-order fractional advection-diffusion equation with a nonlinear source term［J］.SIAM Journal of Numerical Analysis，2009，47(3):1760-1781.

［9］ Liu F，Zhuang P，Anh V，et al.Stability and convergence of the difference methods for the space-time fractional advection-diffusion equation［J］.Applied Mathematics and Computation，2007，191(1):12-20.

［10］ Podlubny I.Fractional differential equations［M］.New York:Academic Press，1999:51-87.

［11］ Ervin J S，Roop J P.Variational solution of fractional advection dispersion equations on bounded domains in Rd［J］.Numerical Mathematics-Theory Methods and Applications，2007，23(2):256-281.

［12］ Richtmyer R D，Morton K W.Difference methods for initial-value problems［M］.New York:Krieger Publishing Co.，1994:133-200.

［13］ Ghazizadeh H R，Maerefat M，Azimi A.Explicit and implicit finite difference schemes for fractional Cattaneo equation［J］.Journal of Computational Physics，2010，229(19):7042-7057.

［14］ Gao G，Sun Z.A compact finite difference scheme for the fractional sub-diffusion equations［J］.Journal of Computational Physics，2011，230(3):586-595.

［15］ Sousa E.Finite difference approximations for a fractional advection diffusion problem［J］.Journal of Computational Physics，2009，228(11):4038-4054.

［16］ Li C，Zeng F.The finite difference methods for fractional ordinary differential equations［J］.Numerical Functional Analysis and Optimization，2014，34(2):149-179.

［17］ Zhang Y，Liao Z，Lin H.Finite difference methods for the time fractional diffusion equation on non-uniform meshes［J］.Journal of Computational Physics，2014，265(1):195-210.

［18］ Brunner H，Han H，Yin D.Artificial boundary conditions and finite difference approximations for a time-fractional diffusion wave equation on a two-dimensional unbounded spatial domain［J］.Journal of Computational Physics，2014，276 (3):541-562.