二維空間時間分數階色散方程的差分方法

張英晗，楊小遠

(北京航空航天大學 數學與系統科學學院，北京100191)

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

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

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

1 問題描述

考慮二維空間時間分數階色散方程:

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

從物理意義方面考慮，要求0 ＜α≤1、1 ＜β、γ≤2，并且u(0，y，t)=u(x，0，t)=0.這意味著在下邊界沒有示蹤劑泄漏.函數f(x，y，t)可用來表示源和匯.當α =1、β =γ =2 時，方程式(1)還原為經典的二維色散方程.以下假設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］系數滿足bk＞0，bk＞bk+1(k=0，1，…)并且

從式(4)～式(6)知，二維隱式差分格式式(8)與分數階初邊值問題式(1)、式(2)是相容的，并且有局部截斷誤差O(Δt)+ O(Δx)+O(Δy).下面證明差分格式式(8)是穩定的，從而由Lax 等價定理［12］在分數階方程中的應用［2］，式(8)是收斂的.分別令

則式(8)可表示為

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

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

證明 當n=0 時，令

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

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

以上分析表明分數階隱式差分格式式(8)是相容的并且無條件收斂的，因此由Lax 等價定理，有以下結論.

定理1 分數階隱式差分格式式(8)是以O(Δt)+O(Δx)+O(Δy)無條件收斂的.

3 交替方向隱式差分格式

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

定義分數階偏差分算子:

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

對于交替方向隱式差分格式，上述算子形式相應變為方向分離乘積的形式:

這產生了一個額外的誤差擾動

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

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

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

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

證明 當n=0 時，令

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

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

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

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

通過以上分析，由式(11)、式(12)定義的分數階交替方向隱式差分方法是相容的和無條件穩定的，因此由Lax 等價定理在分數階方程中的應用［2］，有如下結論.

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

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

4 數值實驗

考慮二維分數階色散方程:

式中:(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

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

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

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

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

5 結 論

通過對空間分數階導數的轉移Grünwald 差分近似，分別構造了有界區域上二維空間時間分數階色散方程的隱式差分格式和交替方向隱式差分格式，應用數學歸納法證明了兩種格式的穩定性和收斂性:

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

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

3)雖然交替方向隱式差分格式在時間方向上的收斂精度低于隱式差分格式，但在計算復雜度上考慮，交替方向隱式差分格式優于隱式差分格式.

4)數值模擬結果與理論分析是一致的.

對于更廣泛的分數階偏微分方程數值逼近問題，比如無界區域上分數階偏微分方程的數值方法、分數階偏微分方程的高階近似差分格式等問題，將在后續的工作中進行研究.

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.