求解三维空间分数阶对流扩散方程的Douglas-Gunn格式

， ,，

(1.西北工业大学 计算科学研究中心 陕西 西安 710129； 2. 河南工业大学 理学院 河南 郑州 450001)

0 引言

近年来，自然界中的反常扩散现象受到科研人员的广泛关注，为研究其独特的物理过程，常常用分数阶偏微分方程建立相应的数学模型. 其中，在包含对流和超扩散两个物理过程的散布现象中，粒子束的传播与经典的布朗运动模型不再一致，此时把经典对流扩散方程中的空间二阶导数替换成分数阶导数构建的空间分数阶对流扩散方程(SFADE)能更准确地模拟这一现象. 分数阶导数或积分具有非局部性，这给相应方程的求解带来了很大困难，在大多数情况下很难得到解析解，因此研究可靠而有效的数值方法就显得尤为重要. 目前常用的数值解法包括有限差分法[1-2]、有限元法[3-4]、有限体积法[5]、配点法[6]以及谱方法[7]等.

由于三维模型在科学研究中有广泛应用，本文考虑有限区域上带有零Dirichlet边界条件的三维SFADE的数值求解. 分数阶导数是一个非局部算子，这就使得离散SFADE得到的线性系统的刚度矩阵不再是稀疏矩阵，导致计算工作量和存储量都非常大，尤其在多维空间情形下. 目前求解三维SFADE的数值方法还比较少. 文献[8]采用了一种交替方向稳定法(alternating direction implicit method, ADI)差分格式求解三维SFADE，并提高了精度.文献[9]提出了一种求解三维分数阶扩散方程的ADI差分格式.文献[10]研究了一种求解三维空间分数阶扩散方程的快速迭代ADI有限差分方法. 在本文中，我们将提出一种求解三维SFADE的有效的ADI有限差分格式，这种方法是将经典的Douglas-Gunn格式中的二阶中心差分算子推广为包含左、右分数阶导数离散算子及一阶中心差分算子在内的复杂算子得到的，同时给出了该格式的稳定性和收敛性的必要证明. 最后用数值实验验证了理论分析的结果.

1 三维SFADE及其Douglas-Gunn格式

本文考虑三维SFADE及它的初边值条件为

(1)

u(x,y,z,t)=0, (x,y,z,t)∈∂Ω×[0,T],

(2)

u(x,y,z,0)=u0(x,y,z), (x,y,z)∈Ω,

(3)

1.1 Riemann-Liouville分数阶导数的离散

其中系数为

记

(4)

(5)

同时，用中心差分近似对流项的一阶导数，记

/2h1.

(6)

1.2 三维SFADE的有限差分近似

用公式(4)～(6)离散空间导数，时间方向采用Crank-Nicolson格式. 记

然后方程(1)就可以表示为

(7)

(8)

(9)

采用经典的Douglas-Gunn格式分解式(9)得到ADI格式，即

(10)

(11)

(12)

(13)

与经典的Douglas-Gunn格式不同，3个方向的二阶中心差分算子在此处分别被替换为δα,x、δβ,y、δγ,z，它们是包含了左、右分数阶导数离散算子等在内的复杂算子，可以认为是经典Douglas-Gunn格式在求解分数阶方程中的推广. 接下来，我们将给出收敛性和稳定性的必要证明.

2 收敛性和稳定性分析

显然，如果在格式(10)～(13)中消去中间解变量，则得到格式(9)，即格式(10)～(13)和(9)是等价的. 下面用矩阵法证明格式(9)是无条件稳定和收敛的. 首先把方程(9)表示成矩阵形式， 令

(14)

(15)

(16)

(17)

其中Aα，Aβ，Aγ是Toeplitz矩阵，Aα和B分别表示为

其中：I是单位矩阵；符号⊗表示Kronecker积[12].Aβ、Aγ与Aα类似. 利用上述记号，式(9)可以写为

τFn+1/2.

(18)

为了证明式(18)的稳定性和收敛性，下面列出一些相关的引理和定理.

引理1[13]一个n阶实矩阵A是正定的，当且仅当矩阵H=(A+AT)/2是正定的；H是正定的，当且仅当H的特征值都是正的.

引理2[13]设A是一个n阶复矩阵，AH表示A的共轭转置，记H=(A+AH)/2，则对于A的任意特征值λ，它的实部满足不等式λmin(H)≤R(λ(A))≤λmax(H)，这里λmin(H)和λmax(H)分别表示H的最小和最大特征值.

煤炭资源影响着我国社会经济的发展，是社会发展的重要基础资源。随着我国工业技术的迅猛发展，煤炭资源的需求量逐渐增加，为新时代背景下煤炭资源的开发提出了更高的要求。随着自动化技术在各领域中的广泛应用，煤矿机电设备的自动化程度也越来越高，有效满足了煤炭资源供给需求。

引理3设A∈Rm×n,B∈Rr×s,C∈Rp×q，则(A⊗B)⊗C=A⊗(B⊗C).

证明此结论可以由Kronecker积的定义直接得到.

引理4设A,B∈Rm×n,C∈Rs×t，则有(A+B)⊗C=A⊗C+B⊗C，C⊗(A+B)=C⊗A+C⊗B.

证明此结论可以由Kronecker积的定义直接得到.

引理5[12]设A∈Rm×n,B∈Rr×s,C∈Rn×p，D∈Rs×t，则(A⊗B)(C⊗D)=AC⊗BD(∈Rmr×pt).

引理6[12]对于任意的矩阵A和B，有(A⊗B)T=AT⊗BT.

为了叙述并证明下述引理和定理，记‖·‖表示矩阵的2-范数.

引理8设A∈Rn×n是正定矩阵，则对任意的σ∈R且σ>0，有‖(I+σA)-1‖<1.

引理9设A∈Rn×n是正定矩阵，则对任意的σ∈R且σ>0，有‖(I+σA)-1(I-σA)‖<1.

证明由矩阵2-范数的定义并记y=(I+σA)-1x，可得

引理10设A,B,I∈Rn×n，A和B乘积可交换，且(I-A)-1、(I-B)-1存在，则(I+A)与(I-B)-1，(I-A)-1与(I-B)-1也是乘积可交换的.

证明首先，由AB=BA，不难验证(I±A)(I-B)=(I-B)(I±A). 所以有

(I+A)(I-B)-1=(I-B)-1(I-B)(I+A)(I-B)-1=(I-B)-1(I+A)(I-B)(I-B)-1=

(I-B)-1(I+A)，(I-A)-1(I-B)-1=((I-B)(I-A))-1=((I-A)(I-B))-1=(I-B)-1(I-A)-1.

定理3差分格式(9)是无条件稳定的.

(19)

(20)

3 数值结果

下面，我们通过两个数值算例验证本文所提出的数值格式的稳定性和收敛阶，也就是说格式是有效的，并在时间和空间方向都具有较高的二阶收敛精度. 设U和Uh分别表示问题(1)～(3)的解析解和采用格式(10)～(13)得到的数值解，用离散的L∞和L2范数计算全局截断误差，即

算例1在问题(1)～(3)中，取Ω=(0,1)×(0,1)×(0,1)，T=1；对流和扩散系数分别为k1=k2=k3=-1，a1=a2=b1=b2=c1=c2=1；初值取u0(x,y,z)=x3(1-x)3y3(1-y)3z3(1-z)3. 已知的解析解为u(x,y,z,t)=e-tx3(1-x)3y3(1-y)3z3(1-z)3，由以上条件容易算出f(x,y,z,t).

对常系数算例1取优化的步长比例M=N1=N2=N3进行测试. 表1列出的数值结果表明，用格式(10)～(13)计算常系数问题(1)～(3)时，算法是无条件稳定的，而且在时间及空间方向都是二阶收敛的，这和理论分析的结果一致.

算例2在问题(1)～(3)中，取Ω=(0,1)×(0,1)×(0,1)，T=1；对流和扩散系数分别为k1=0.25x，k2=0.25y，k3=0.25z，a1=xα，a2=(1-x)α，b1=yβ，b2=(1-y)β，c1=zγ，c2=(1-z)γ；初值取u0(x,y,z)=x3(1-x)3y3(1-y)3z3(1-z)3．已知的解析解为u(x,y,z,t)=e-tx3(1-x)3y3(1-y)3z3·(1-z)3，由以上条件f(x,y,z,t)容易算出.

表2列出了变系数算例2的数值结果，这里也取优化的步长比例M=N1=N2=N3进行测试，数值结果表明用格式(10)～(13)计算变系数问题(1)～(3)时，算法是无条件稳定的，而且在时间及空间方向也都具有二阶收敛率，这和理论分析的结果是非常吻合的.

表1 算例1在时刻t=1，取M=N1=N2=N3的数值误差和收敛阶Tab.1 The errors and convergence orders for example 1 at t=1 with M=N1=N2=N3

表2 算例2在时刻t=1，取M=N1=N2=N3的数值误差和收敛阶Tab.2 The errors and convergence orders for example 2 at t=1 with M=N1=N2=N3

4 结论

本文将求解三维整数阶抛物方程的经典Douglas-Gunn格式推广到分数阶，提出了一种求解三维SFADE的有效的数值方法，并证明该格式具有无条件稳定性和较高的二阶收敛精度，必要而充足的数值实验验证了理论结果. 最后，由于分数阶导数是非局部算子，对于多维空间问题的求解需要耗费较大的计算工作量和空间存储量，在今后的工作中，我们将考虑开展适当的快速算法，以减少计算花费和加快计算速度.