分数阶偏微分方程模型在模拟许多复杂的物理现象中是非常有效的模型,在金融经济、半导体、生物、水文、控制理论等领域都有着重要的应用^{[1, 2, 3, 4, 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, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18]}.

注意到很多分数阶微分方程差分方法的文献处理的都是一维问题空间分数阶导数或者时间分数阶导数,对高维问题以及空间时间分数阶导数很少有人涉及^{[3, 9]}.受到以上研究工作的启发,本文研究二维空间时间分数阶色散方程的有限差分格式.由于分数阶问题的隐式差分格式在每一时间步都需要处理非稀疏线性系统,于是隐式差分格式是计算密集型的.通过把偏微分方程差分方法中经常用到的交替方向隐式差分格式推广到分数阶偏微分方程,从而克服了这一问题.

考虑二维空间时间分数阶色散方程:

从物理意义方面考虑,要求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].

令t_n=nτ,n=0,1,…,N;x_i=ih_x,i=0,1,…,N_x;y_j=jh_y,j=0,1,…,N_y,其中$\tau =\frac{T}{N}$为时间步长,h_x=$\frac{{{L}_{x}}}{{{N}_{x}}}$和h_y=$\frac{{{L}_{y}}}{{{N}_{y}}}$为空间步长.令u_i,jⁿ为差分方法的数值解.引入记号Δ_tu(x_i,y_j,t_n)=u(x_i,y_j,t_n+1)－u(x_i,y_j,t_n),按如下方式对时间分数阶导数进行离散:

对于空间分数阶导数D_x^βu(x_i,y_j,t_n+1)和D_y^γu(x_i,y_j,t_n+1),采用转移Grünwald公式^[2]:

由式(4)~式(6),可以得到如下的隐式差分格式:

初边值条件分别离散为

引理1^[2] 系数b_k、g_β^k、g_γ^k(k=0,1,…)满足b_k>0,b_k>b_k+1(k=0,1，…)并且

从式(4)～式(6)知,二维隐式差分格式式(8)与分数阶初边值问题式(1)、式(2)是相容的，并且有局部截断误差O(Δt)+O(Δx)+O(Δy).下面证明差分格式式(8)是稳定的,从而由Lax等价定理^[12]在分数阶方程中的应用^[2],式(8)是收敛的.分别令

引理2 ||Eⁿ||_∞≤||E⁰||_∞,n=0,1,….

证明 当n=0时,令

根据引理1,有$\sum\limits_{k=0}^{l+1}{g_{\beta }^{k}}<0,\sum\limits_{k=0}^{m+1}{g_{\gamma }^{k}}<0$于是

假定||E^k||_∞≤||E⁰||_∞,k=1,2,…,n,令n+1层上的误差|ε_l,mⁿ⁺¹|=$\underset{1\le i\le {{N}_{x}}-1,1\le j\le {{N}_{y}}-1}{\mathop{\max }}\,$|ε_i,jⁿ⁺¹|.由引理1：

证毕

以上分析表明分数阶隐式差分格式式(8)是相容的并且无条件收敛的,因此由Lax等价定理,有以下结论.

定理1 分数阶隐式差分格式式(8)是以O(Δt)+O(Δx)+O(Δy)无条件收敛的.

二维隐式差分格式式(8)是收敛的并且有局部误差O(Δt)+O(Δx)+O(Δy).然而在每一时间步,需要处理(N_x－1)×(N_y－1)个未知量的非稀疏线性系统,于是隐式差分格式是计算密集型的.随着网格的加细和空间维数的提高,计算复杂度变得非常庞大.为了解决这个问题,下面把偏微分方程差分方法中常用的交替方向隐式差分格式推广到分数阶偏微分方程,即在每一时间步只求解一个方向.

定义分数阶偏差分算子:

对于交替方向隐式差分格式,上述算子形式相应变为方向分离乘积的形式:

在计算上,上述形式的交替方向隐式差分格式可表示为如下的迭代形式.在时间层t_n+1:

1) 首先对每个固定的y_j求解x方向,产生中间解u_i,j^*,n+1:

2) 对每个固定的x_i求解y方向:

下面证明隐式差分格式的交替方向方法式(11)、式(12)是相容的和稳定的,因此由Lax等价定理,也是收敛的.交替方向隐式差分格式的额外扰动误差式(10)的阶为O(Δt^2α),因此由式(11)、式(12)定义的交替方向隐式差分格式是相容的,并且有局部截断误差O(Δt^α)+O(Δx)+O(Δy).

假定v_i,jⁿ(0≤i≤N_x,0≤j≤N_y,0≤n≤N)是式(11)和式(12)的逼近解,则误差ε_i,jⁿ=v_i,jⁿ－u_i,jⁿ满足：

引理3 ||Eⁿ⁺¹||_∞≤||E⁰||_∞,n=0,1,….

证明 当n=0时,令

根据引理1和式(11),有

令|ε_r,s¹|=$\underset{1\le i\le {{N}_{x}}-1,1\le j\le {{N}_{y}}-1}{\mathop{\max }}\,$|ε_i,j¹|,根据引理1和式(12),有

假定||E^k||_∞≤||E⁰||_∞,k=1,2,…,n,令|ε_l,m^*,n+1|=$\underset{1\le i\le {{N}_{x}}-1,1\le j\le {{N}_{y}}-1}{\mathop{\max }}\,$|ε_i,j^*,n+1|,根据引理1,得

证毕

通过以上分析,由式(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,所以为了减少计算复杂度,牺牲了一些收敛精度.

考虑二维分数阶色散方程:

验证可知,方程式(13)的精确解为

表 1是交替方向隐式差分格式计算结果与精确解之间的最大误差.最后一列误差率表示上一行最大绝对误差与本行最大绝对误差的比率,即数值方法的收敛速度.可以看到最大绝对误差是以速度2几乎线性下降的,这与理论分析结果是一致的.

通过对空间分数阶导数的转移Grünwald差分近似,分别构造了有界区域上二维空间时间分数阶色散方程的隐式差分格式和交替方向隐式差分格式,应用数学归纳法证明了两种格式的稳定性和收敛性:

1) 隐式差分格式是以O(Δt)+O(Δx)+O(Δy)无条件收敛的.

2) 交替方向隐式差分格式的收敛阶为O(Δt^α)+O(Δx)+O(Δy).

3) 虽然交替方向隐式差分格式在时间方向上的收敛精度低于隐式差分格式,但在计算复杂度上考虑,交替方向隐式差分格式优于隐式差分格式.

4) 数值模拟结果与理论分析是一致的.

对于更广泛的分数阶偏微分方程数值逼近问题,比如无界区域上分数阶偏微分方程的数值方法、分数阶偏微分方程的高阶近似差分格式等问题,将在后续的工作中进行研究.