张英晗，杨小远

(北京航空航天大学 数学与系统科学学院，北京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.