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及它的初边值条件为

$ \begin{array}{l} \frac{{\partial u\left( {x,y,z,t} \right)}}{{\partial t}} = {a_{1{x_L}}}D_x^\alpha u\left( {x,y,z,t} \right) + {\alpha _{2x}}D_{{x_R}}^\alpha u\left( {x,y,z,t} \right) + {b_{1{y_L}}}D_y^\beta u\left( {x,y,z,t} \right) + \\ {b_{1y}}D_{{y_R}}^\beta u\left( {x,y,z,t} \right) + {c_{1{z_L}}}D_z^\gamma u\left( {x,y,z,t} \right) + {c_{2z}}D_{{z_R}}^\gamma \mathit{u}\left( {x,y,z,t} \right) + {k_1}\frac{{\partial u\left( {x,y,z,t} \right)}}{{\partial x}} + \\ {k_2}\frac{{\partial u\left( {x,y,z,t} \right)}}{{\partial y}} + {k_3}\frac{{\partial u\left( {x,y,z,t} \right)}}{{\partial z}} + f\left( {x,y,z,t} \right),\left( {x,y,z,t} \right) \in \mathit{\Omega } \times \left( {0,T} \right], \end{array} $

(1)

$ u\left( {x,y,z,t} \right) = 0,\;\;\;\;\;\left( {x,y,z,t} \right) \in \partial \mathit{\Omega } \times \left[ {0,T} \right], $

(2)

$ u\left( {x,y,z,0} \right) = {u_0}\left( {x,y,z} \right),\;\;\;\;\left( {x,y,z} \right) \in \mathit{\Omega }, $

(3)

其中：1 < α, β, γ < 2；Ω=(x_L, x_R)×(y_L, y_R)×(z_L, z_R)⊂R³；a₁, a₂, b₁, b₂, c₁, c₂≥0是3个空间方向的左、右扩散系数；k₁、k₂、k₃分别是3个空间方向的对流系数；f(x, y, z, t)是源项.方程(1)中的分数阶导数是Riemann-Liouville型的，即函数v(x)的α(1 < α < 2)阶Riemann-Liouville左导数和右导数分别定义为$ _{{{x}_{L}}}D_{x}^{\alpha }v\left( x \right)=\frac{1}{\Gamma (2-\alpha )}\frac{{{\partial }^{2}}}{\partial {{x}^{2}}}\int_{{{x}_{L}}}^{x}{{{\left( x-\xi \right)}^{1-\alpha }}v\left( \xi \right)\text{d}\xi } $和$ _{x}D_{{{x}_{R}}}^{\alpha }v\left( x \right)=\frac{1}{\Gamma (2-\alpha )}\frac{{{\partial }^{2}}}{\partial {{x}^{2}}}\int_{x}^{{{x}_{R}}}{{{\left( \xi -x \right)}^{1-\alpha }}v\left( \xi \right)\text{d}\xi }$.

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

设N₁、N₂、N₃和M为正整数，h₁=(x_R－x_L)/N₁，h₂=(y_R－y_L)/N₂，h₃=(z_R－z_L)/N₃，τ=T/M分别是一致的空间步长和时间步长，由此定义的空间和时间的剖分为x_i=x_L+ih，i=0, 1, …, N₁，y_j=y_L+jh，j=0, 1, …, N₂，z_m=z_L+mh，m=0, 1, …, N₃，t_n=nτ，n=0, 1, …, M.设u_{i, j, m}ⁿ表示u(x_i, y_j, z_m, t_n)的近似值，f_{i, j, m}ⁿ=f(x_i, y_j, z_m, t_n).采用文献[11]中的方法离散方程(1)中的分数阶导数.以x方向为例，有

$ _{{x_L}}D_x^\alpha u\left( {x,{y_j},{z_m},{t_n}} \right)\left| {_{x = {x_i}}} \right. = \frac{1}{{{\rm{\Gamma }}\left( {4 - \alpha } \right)h_1^\alpha }}\sum\limits_{l = 0}^{i + 1} {g_l^\alpha u\left( {{x_{i - l + 1}},{y_j},{z_m},{t_n}} \right)} + O\left( {h_1^2} \right), $

$ _xD_{{x_R}}^\alpha u\left( {x,{y_j},{z_m},{t_n}} \right)\left| {_{x = {x_i}}} \right. = \frac{1}{{{\rm{\Gamma }}\left( {4 - \alpha } \right)h_1^\alpha }}\sum\limits_{l = 0}^{{N_1} - i + 1} {g_l^\alpha u\left( {{x_{i + l - 1}},{y_j},{z_m},{t_n}} \right)} + O\left( {h_1^2} \right), $

其中系数为

$ g_l^\alpha = \left\{ {\begin{array}{*{20}{c}} \begin{array}{l} 1,\\ - 4 + {2^{3 - \alpha }},\\ 6 - {2^{5 - \alpha }} + {3^{3 - \alpha }},\\ {\left( {l + 1} \right)^{3 - \alpha }} - 4{l^{3 - \alpha }} + 6{\left( {l - 1} \right)^{3 - \alpha }} - 4{\left( {l - 2} \right)^{3 - \alpha }} + {\left( {l - 3} \right)^{3 - \alpha }}, \end{array}&\begin{array}{l} l = 0,\\ l = 1,\\ l = 2,\\ l \ge 3. \end{array} \end{array}} \right. $

记

$ \delta _{\alpha ,x}^ + u_{i,j,m}^n = \frac{1}{{{\rm{\Gamma }}\left( {4 - \alpha } \right)h_1^\alpha }}\sum\limits_{l = 0}^{i + 1} {g_l^\alpha u_{i - l + 1,j,m}^n} , $

(4)

$ \delta _{\alpha ,x}^ - u_{i,j,m}^n = \frac{1}{{{\rm{\Gamma }}\left( {4 - \alpha } \right)h_1^\alpha }}\sum\limits_{l = 0}^{{N_1} - i + 1} {g_l^\alpha u_{i + l - 1,j,m}^n} . $

(5)

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

$ {D_{\alpha ,x}}u_{i,j,m}^n = \left( {u_{i + 1,j,m}^n - u_{i - 1,j,m}^n} \right)/2{h_1}. $

(6)

为了便于表示，进一步引入记号$ \tilde{D} $_{α, x}u_{i, j, m}ⁿ=k₁D_{α, x}u_{i, j, m}ⁿ, $ \tilde{\delta } $_{α, x}⁺u_{i, j, m}ⁿ=a₁δ_{α, x}⁺u_{i, j, m}ⁿ, $ \tilde{\delta } $_{α, x}^－u_{i, j, m}ⁿ=a₂δ_{α, x}^－u_{i, j, m}ⁿ.y和z方向的记号可以类似表示.

1.2 三维SFADE的有限差分近似

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

$ {\delta _{\alpha ,x}}: = \tilde \delta _{\alpha ,x}^ + + \tilde \delta _{\alpha ,x}^ - + {{\tilde D}_{\alpha ,x}},{\delta _{\beta ,y}}: = \tilde \delta _{\beta ,y}^ + + \tilde \delta _{\beta ,y}^ - + {{\tilde D}_{\beta ,y}},{\delta _{\gamma ,z}}: = \tilde \delta _{\gamma ,z}^ + + \tilde \delta _{\gamma ,z}^ - + {{\tilde D}_{\gamma ,z}}, $

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

$ \begin{array}{l} \left( {1 - \frac{\tau }{2}{\delta _{\alpha ,x}} - \frac{\tau }{2}{\delta _{\beta ,y}} - \frac{\tau }{2}{\delta _{\gamma ,z}}} \right)u\left( {{x_i},{y_j},{z_m},{t_{n + 1}}} \right) = \\ \left( {1 + \frac{\tau }{2}{\delta _{\alpha ,x}} + \frac{\tau }{2}{\delta _{\beta ,y}} + \frac{\tau }{2}{\delta _{\gamma ,z}}} \right)u\left( {{x_i},{y_j},{z_m},{t_n}} \right) + \tau f\left( {{x_i},{y_j},{z_m},{t_{n + 1/2}}} \right) + R_{i,j,m}^{n + 1}. \end{array} $

(7)

存在正常数$ \tilde{C} $，使得|R_{i, j, m}ⁿ⁺¹|≤$ \tilde{C} $τ(τ²+h₁²+h₂²+h₃²).接下来，在方程(7)中用近似值u_{i, j, m}ⁿ代替函数值u(x_i, y_j, z_m, t_n)，并去掉高阶项，得到方程(1)的全离散格式，

$ \left( {1 - \frac{\tau }{2}{\delta _{\alpha ,x}} - \frac{\tau }{2}{\delta _{\beta ,y}} - \frac{\tau }{2}{\delta _{\gamma ,z}}} \right)u_{i,j,m}^{n + 1} = \left( {1 + \frac{\tau }{2}{\delta _{\alpha ,x}} + \frac{\tau }{2}{\delta _{\beta ,y}} + \frac{\tau }{2}{\delta _{\gamma ,z}}} \right)u_{i,j,m}^n + \tau f_{i,j,m}^{n + 1/2}. $

(8)

为便于计算，下面构造ADI差分格式.将方程(8)左边加上高阶项($ \frac{{{\tau }^{2}}}{4} $δ_{α, x}δ_{β, y}+$ \frac{{{\tau }^{2}}}{4} $δ_{α, x}δ_{γ, z}+$ \frac{{{\tau }^{2}}}{4} $δ_{β, y}δ_{γ, z})·(u_{i, j, m}ⁿ⁺¹－u_{i, j, m}ⁿ)－$ \frac{{{\tau }^{3}}}{8} $δ_{α, x}δ_{β, y}δ_{γ, z}(u_{i, j, m}ⁿ⁺¹+u_{i, j, m}ⁿ)，再把适当的部分移项到方程的右边并分解因式，得

$ \left( {1 - \frac{\tau }{2}{\delta _{\alpha ,x}}} \right)\left( {1 - \frac{\tau }{2}{\delta _{\beta ,y}}} \right)\left( {1 - \frac{\tau }{2}{\delta _{\gamma ,z}}} \right)u_{i,j,m}^{n + 1} = \left( {1 + \frac{\tau }{2}{\delta _{\alpha ,x}}} \right)\left( {1 + \frac{\tau }{2}{\delta _{\beta ,y}}} \right) \\ \left( {1 + \frac{\tau }{2}{\delta _{\gamma ,z}}} \right)u_{i,j,m}^n + \tau f_{i,j,m}^{n + 1/2}. $

(9)

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

$ \left( {1 - \frac{\tau }{2}{\delta _{\alpha ,x}}} \right)\Delta {u^ * } = \left( {\tau {\delta _{\alpha ,x}} + \tau {\delta _{\beta ,y}} + \tau {\delta _{\gamma ,z}}} \right)u_{i,j,m}^n + \tau f_{i,j,m}^{n + 1/2}; $

(10)

$ \left( {1 - \frac{\tau }{2}{\delta _{\beta ,y}}} \right)\Delta {u^{ * * }} = \Delta {u^ * }; $

(11)

$ \left( {1 - \frac{\tau }{2}{\delta _{\gamma ,z}}} \right)\Delta u = \Delta {u^{ * * }}; $

(12)

$ \Delta u = u_{i,j,m}^{n + 1} - u_{i,j,m}^n. $

(13)

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

2 收敛性和稳定性分析

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

$ \begin{array}{*{20}{l}} \begin{array}{l} {{\boldsymbol{U}}^n} = \left[ {u_{1,1,1}^n,u_{2,1,1}^n, \cdots ,u_{{N_1} - 1,1,1}^n,u_{1,2,1}^n,u_{2,2,1}^n, \cdots ,\\ u_{{N_1} - 1,2,1}^n, \cdots ,u_{1.{N_2} - 1,1}^n,u_{2.{N_2} - 1,1}^n, \cdots ,} \right.\\ u_{{N_1} - 1.{N_2} - 1,1}^n,u_{1,1,2}^n,u_{2,1,2}^n, \cdots , \end{array}\\ \begin{array}{l} u_{{N_1} - 1,1,2}^n, \cdots ,u_{1,{N_2} - 1,2}^n,u_{2,{N_2} - 1,2}^n, \cdots ,u_{{N_1} - 1,{N_2} - 1,2}^n, \cdots ,u_{1,{N_2} - 1,{N_3} - 1}^n,\\ {\left. {u_{2,{N_2} - 1,{N_3} - 1}^n, \cdots ,u_{{N_1} - 1,{N_2} - 1,{N_3} - 1}^n} \right]^{\rm{T}}}, \end{array} \end{array} $

(14)

为了书写简单起见，我们用记号Uⁿ=(([u_{i, j, m}ⁿ]_{i=1, 2, …，N1－1})_{j=1, 2, …，N2－1})_{m=1, 2, …，N3－1}表示式(14)，类似的表示还有Fⁿ=(([f_{i, j, m}ⁿ]_{i=1, 2, …，N1－1})_{j=1, 2, …，N2－1})_{m=1, 2, …，N3－1}.记

$ {{\mathit{\boldsymbol{\tilde A}}}_x} = \frac{{{a_1}\tau }}{{2\Gamma \left( {4 - \alpha } \right)h_1^\alpha }}\mathit{\boldsymbol{I}} \otimes \mathit{\boldsymbol{I}} \otimes {\mathit{\boldsymbol{A}}_\alpha } + \frac{{{a_2}\tau }}{{2\Gamma \left( {4 - \alpha } \right)h_1^\alpha }}\mathit{\boldsymbol{I}} \otimes \mathit{\boldsymbol{I}} \otimes \mathit{\boldsymbol{A}}_\alpha ^{\rm{T}} + \frac{{{k_1}\tau }}{{4{h_1}}}\mathit{\boldsymbol{I}} \otimes \mathit{\boldsymbol{I}} \otimes \mathit{\boldsymbol{B}}, $

(15)

$ {{\mathit{\boldsymbol{\tilde A}}}_y} = \frac{{{b_1}\tau }}{{2\Gamma \left( {4 - \beta } \right)h_2^\beta }}\mathit{\boldsymbol{I}} \otimes {\mathit{\boldsymbol{A}}_\beta } \otimes \mathit{\boldsymbol{I}} + \frac{{{b_2}\tau }}{{2\Gamma \left( {4 - \beta } \right)h_2^\beta }}\mathit{\boldsymbol{I}} \otimes \mathit{\boldsymbol{A}}_\beta ^{\rm{T}} \otimes \mathit{\boldsymbol{I}} + \frac{{{k_2}\tau }}{{4{h_2}}}\mathit{\boldsymbol{I}} \otimes \mathit{\boldsymbol{B}} \otimes \mathit{\boldsymbol{I}}, $

(16)

$ {{\mathit{\boldsymbol{\tilde A}}}_z} = \frac{{{c_1}\tau }}{{2\Gamma \left( {4 - \gamma } \right)h_3^\gamma }}{\mathit{\boldsymbol{A}}_\gamma } \otimes \mathit{\boldsymbol{I}} \otimes \mathit{\boldsymbol{I}} + \frac{{{c_2}\tau }}{{2\Gamma \left( {4 - \gamma } \right)h_3^\gamma }}\mathit{\boldsymbol{A}}_\gamma ^{\rm{T}} \otimes \mathit{\boldsymbol{I}} \otimes \mathit{\boldsymbol{I}} + \frac{{{k_3}\tau }}{{4{h_3}}}\mathit{\boldsymbol{B}} \otimes \mathit{\boldsymbol{I}} \otimes \mathit{\boldsymbol{I}}, $

(17)

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

$ {\mathit{\boldsymbol{A}}_\alpha } = \left[ {\begin{array}{*{20}{c}} {g_1^\alpha }&{g_0^\alpha }&0& \cdots &0&0\\ {g_2^\alpha }&{g_1^\alpha }&{g_0^\alpha }&0& \cdots &0\\ {g_3^\alpha }&{g_2^\alpha }&{g_1^\alpha }&{g_0^\alpha }&{g_2^\alpha }&0\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots \\ {g_{{N_1} - 2}^\alpha }& \cdots&\cdots&\cdots &{g_1^\alpha }&{g_0^\alpha }\\ {g_{{N_1} - 1}^\alpha }&{g_{{N_1} - 2}^\alpha }&{g_{{N_1} - 3}^\alpha }& \cdots &{g_2^\alpha }&{g_1^\alpha } \end{array}} \right],\mathit{\boldsymbol{B}} = \left[ {\begin{array}{*{20}{c}} 0&1&0& \cdots &0&0\\ { - 1}&0&1&0& \cdots &0\\ 0&{ - 1}&0&1& \cdots &0\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots \\ 0&0&0& \cdots &0&1\\ 0&0&0& \cdots &{ - 1}&0 \end{array}} \right], $

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

$ \left( {\mathit{\boldsymbol{I}} - {{\mathit{\boldsymbol{\tilde A}}}_x}} \right)\left( {\mathit{\boldsymbol{I}} - {{\mathit{\boldsymbol{\tilde A}}}_y}} \right)\left( {\mathit{\boldsymbol{I}} - {{\mathit{\boldsymbol{\tilde A}}}_z}} \right){\mathit{\boldsymbol{U}}^{n + 1}} = \\ \left( {\mathit{\boldsymbol{I}} + {{\mathit{\boldsymbol{\tilde A}}}_x}} \right)\left( {\mathit{\boldsymbol{I}} + {{\mathit{\boldsymbol{\tilde A}}}_y}} \right)\left( {\mathit{\boldsymbol{I}} + {{\mathit{\boldsymbol{\tilde A}}}_z}} \right){\mathit{\boldsymbol{U}}^n} + \tau {\mathit{\boldsymbol{F}}^{n + 1/2}}. $

(18)

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

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

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

定理1^[8]  设A_α是式(15)中的Toeplitz矩阵，则对于A_α的任意特征值λ，有R(λ(A_α)) < 0，并且A_α是负定矩阵.同时，R(λ(d₁A_α+d₂A_α^T)) < 0，d₁, d₂≥0, d₁²+d₂²≠0.

引理3  设A∈R^m×n, B∈R^r×s, C∈R^p×q，则(A⊗B)⊗C=A⊗(B⊗C).

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

引理4  设A, B∈R^m×n, C∈R^s×t，则有(A+B)⊗C=A⊗C+B⊗C，C⊗(A+B)=C⊗A+C⊗B.

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

引理5^[12]  设A∈R^m×n, B∈R^r×s, C∈R^n×p，D∈R^s×t，则(A⊗B)(C⊗D)=AC⊗BD(∈R^mr×pt).

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

引理7^[12]  设矩阵A∈R^n×n有特征值{λ_i}_i=1ⁿ，矩阵B∈R^m×m有特征值{μ_j}_j=1^m.则矩阵A⊗B的mn个特征值为λ₁μ₁, …, λ₁μ_m, λ₂μ₁, …, λ₂μ_m, …, λ_nμ₁, …，λ_nμ_m.

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

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

证明  由矩阵2-范数的定义，有$ {{\left\| {{(\mathit{\boldsymbol{I}}+\sigma \mathit{\boldsymbol{A}})}^{-1}} \right\|}^{2}}=\underset{x\ne 0}{\mathop{\text{max}}}\, \frac{({{(\mathit{\boldsymbol{I}}+\sigma \mathit{\boldsymbol{A}})}^{-1}}x, {{(\mathit{\boldsymbol{I}}+\sigma \mathit{\boldsymbol{A}})}^{-1}}x)}{(x, x)} $.设y=(I+σA)^－1x，则有$ {{\left\| {{(\mathit{\boldsymbol{I}}+\sigma \mathit{\boldsymbol{A}})}^{-1}} \right\|}^{2}}=\underset{y\ne 0}{\mathop{\text{max}}}\, \frac{(y, y)}{((\mathit{\boldsymbol{I}}+\sigma \mathit{\boldsymbol{A}})y, (\mathit{\boldsymbol{I}}+\sigma \mathit{\boldsymbol{A}})y)}=\frac{1}{\underset{y \ne 0}{\mathop{\text{min}}}[1+2\sigma \frac{(\mathit{\boldsymbol{A}}y, y)}{(y, y)}+{{\sigma }^{2}}\frac{\left( \mathit{\boldsymbol{A}}y, \mathit{\boldsymbol{A}}y \right)}{(y, y)}]}<1 $.

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

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

$ \begin{array}{l} {\left\| {\left( {\mathit{\boldsymbol{I}} - \sigma \mathit{\boldsymbol{A}}} \right){{\left( {\mathit{\boldsymbol{I}} + \sigma \mathit{\boldsymbol{A}}} \right)}^{ - 1}}} \right\|^2} = \mathop {\max }\limits_{x \ne 0} \frac{{\left( {\left( {\mathit{\boldsymbol{I}} - \sigma \mathit{\boldsymbol{A}}} \right){{\left( {\mathit{\boldsymbol{I}} + \sigma \mathit{\boldsymbol{A}}} \right)}^{ - 1}}x,\left( {\mathit{\boldsymbol{I}} - \sigma \mathit{\boldsymbol{A}}} \right){{\left( {\mathit{\boldsymbol{I}} + \sigma \mathit{\boldsymbol{A}}} \right)}^{ - 1}}x} \right)}}{{\left( {x,x} \right)}} = \\ \mathop {\max }\limits_{y \ne 0} \frac{{\left( {\left( {\mathit{\boldsymbol{I}} - \sigma \mathit{\boldsymbol{A}}} \right)y,\left( {\mathit{\boldsymbol{I}} - \sigma \mathit{\boldsymbol{A}}} \right)y} \right)}}{{\left( {\left( {\mathit{\boldsymbol{I}} + \sigma \mathit{\boldsymbol{A}}} \right)y,\left( {\mathit{\boldsymbol{I}} + \sigma \mathit{\boldsymbol{A}}} \right)y} \right)}} = \mathop {\max }\limits_{y \ne 0} \frac{{\left( {y,y} \right) - 2\sigma \left( {\mathit{\boldsymbol{A}}y,y} \right) + {\sigma ^2}\left( {\mathit{\boldsymbol{A}}y,\mathit{\boldsymbol{A}}y} \right)}}{{\left( {y,y} \right) + 2\sigma \left( {\mathit{\boldsymbol{A}}y,y} \right) + {\sigma ^2}\left( {\mathit{\boldsymbol{A}}y,\mathit{\boldsymbol{A}}y} \right)}} < 1. \end{array} $

引理10  设A, B, I∈R^n×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).所以有

$ \begin{array}{l} \left( {\mathit{\boldsymbol{I}} + \mathit{\boldsymbol{A}}} \right){\left( {\mathit{\boldsymbol{I}} - \mathit{\boldsymbol{B}}} \right)^{ - 1}} = {\left( {\mathit{\boldsymbol{I}} - \mathit{\boldsymbol{B}}} \right)^{ - 1}}\left( {\mathit{\boldsymbol{I}} - \mathit{\boldsymbol{B}}} \right)\left( {\mathit{\boldsymbol{I}} + \mathit{\boldsymbol{A}}} \right){\left( {\mathit{\boldsymbol{I}} - \mathit{\boldsymbol{B}}} \right)^{ - 1}} = \\{\left( {\mathit{\boldsymbol{I}} - \mathit{\boldsymbol{B}}} \right)^{ - 1}}\left( {\mathit{\boldsymbol{I}} + \mathit{\boldsymbol{A}}} \right)\left( {\mathit{\boldsymbol{I}} - \mathit{\boldsymbol{B}}} \right){\left( {\mathit{\boldsymbol{I}} - \mathit{\boldsymbol{B}}} \right)^{ - 1}} = \\ {\left( {\mathit{\boldsymbol{I}} - \mathit{\boldsymbol{B}}} \right)^{ - 1}}\left( {\mathit{\boldsymbol{I}} + \mathit{\boldsymbol{A}}} \right),{\left( {\mathit{\boldsymbol{I}} - \mathit{\boldsymbol{A}}} \right)^{ - 1}}{\left( {\mathit{\boldsymbol{I}} - \mathit{\boldsymbol{B}}} \right)^{ - 1}} = \\ {\left( {\left( {\mathit{\boldsymbol{I}} - \mathit{\boldsymbol{B}}} \right)\left( {\mathit{\boldsymbol{I}} - \mathit{\boldsymbol{A}}} \right)} \right)^{ - 1}} = {\left( {\left( {\mathit{\boldsymbol{I}} - \mathit{\boldsymbol{A}}} \right)\left( {\mathit{\boldsymbol{I}} - \mathit{\boldsymbol{B}}} \right)} \right)^{ - 1}} = {\left( {\mathit{\boldsymbol{I}} - \mathit{\boldsymbol{B}}} \right)^{ - 1}}{\left( {\mathit{\boldsymbol{I}} - \mathit{\boldsymbol{A}}} \right)^{ - 1}}. \end{array} $

定理2  由式(15)~(17)定义的矩阵$ \mathit{\boldsymbol{\tilde{A}}} $_x、$ \mathit{\boldsymbol{\tilde{A}}} $_y、$ \mathit{\boldsymbol{\tilde{A}}} $_z是负定的.

证明  记$ {{\mathit{\boldsymbol{A}}}_{x}}=\frac{{{a}_{1}}\tau }{2\Gamma (4-\alpha )h_{1}^{\alpha }}{{\mathit{\boldsymbol{A}}}_{\alpha }}+\frac{{{a}_{2}}\tau }{2\Gamma (4-\alpha )h_{1}^{\alpha }}\mathit{\boldsymbol{A}}_{\alpha }^{\text{T}}+\frac{{{k}_{1}}\tau }{4{{h}_{1}}}\mathit{\boldsymbol{B}} $，则有

$ \mathit{\boldsymbol{A}}_x^{\rm{T}} = \frac{{{a_1}\tau }}{{2{\rm{\Gamma }}\left( {4 - \alpha } \right)h_1^\alpha }}\mathit{\boldsymbol{A}}_\alpha ^{\rm{T}} + \frac{{{a_2}\tau }}{{2{\rm{\Gamma }}\left( {4 - \alpha } \right)h_1^\alpha }}{\mathit{\boldsymbol{A}}_\alpha } + \frac{{{k_1}\tau }}{{4{h_1}}}{\mathit{\boldsymbol{B}}^{\rm{T}}},\frac{{{\mathit{\boldsymbol{A}}_x} + \mathit{\boldsymbol{A}}_x^{\rm{T}}}}{2} = \frac{{\left( {{a_1} + {a_2}} \right)\tau }}{{2{\rm{\Gamma }}\left( {4 - \alpha } \right)h_1^\alpha }}\frac{{{\mathit{\boldsymbol{A}}_\alpha } + \mathit{\boldsymbol{A}}_\alpha ^{\rm{T}}}}{2}, $

由定理1，可得$ R(\lambda (\frac{{{\mathit{\boldsymbol{A}}}_{x}}+\mathit{\boldsymbol{A}}_{x}^{\text{T}}}{2}))<0 $，也就是$ \lambda (\frac{{{\mathit{\boldsymbol{A}}}_{x}}+\mathit{\boldsymbol{A}}_{x}^{\text{T}}}{2})<0 $，而且由引理3，4和6，有$\frac{{{{\mathit{\boldsymbol{\tilde{A}}}}}_{x}}+\mathit{\boldsymbol{\tilde{A}}}_{x}^{\text{T}}}{2}=\mathit{\boldsymbol{I}}\otimes \mathit{\boldsymbol{I}}\otimes (\frac{{{\mathit{\boldsymbol{A}}}_{x}}+\mathit{\boldsymbol{A}}_{x}^{\text{T}}}{2}) $，再根据引理7，可以得到$ \lambda (\frac{{{{\mathit{\boldsymbol{\tilde{A}}}}}_{x}}+\mathit{\boldsymbol{\tilde{A}}}_{x}^{\text{T}}}{2})<0$.最后由引理2和引理1可得R(λ($ \mathit{\boldsymbol{\tilde{A}}} $_x)) < 0，且$ \mathit{\boldsymbol{\tilde{A}}} $_x是负定的.类似地，可以证明$ \mathit{\boldsymbol{\tilde{A}}} $_y和$ \mathit{\boldsymbol{\tilde{A}}} $_z是负定的.

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

证明  我们利用差分格式(9)的矩阵形式(18)证明.设Uⁿ和Vⁿ分别是格式(18)对应于初值U⁰和V⁰的解，记εⁿ：=Uⁿ－Vⁿ=(([ε_{i, j, m}ⁿ]_{i=1, 2, …，N1－1})_{j=1, 2, …，N2－1})_{m=1, 2, …，N3－1}.由式(18)可得εⁿ满足方程

$ \left( {\mathit{\boldsymbol{I}} - {{\mathit{\boldsymbol{\tilde A}}}_x}} \right)\left( {\mathit{\boldsymbol{I}} - {{\mathit{\boldsymbol{\tilde A}}}_y}} \right)\left( {\mathit{\boldsymbol{I}} - {{\mathit{\boldsymbol{\tilde A}}}_z}} \right){\mathit{\boldsymbol{\varepsilon }}^{n + 1}} = \left( {\mathit{\boldsymbol{I}} + {{\mathit{\boldsymbol{\tilde A}}}_x}} \right)\left( {\mathit{\boldsymbol{I}} + {{\mathit{\boldsymbol{\tilde A}}}_y}} \right)\left( {\mathit{\boldsymbol{I}} + {{\mathit{\boldsymbol{\tilde A}}}_z}} \right){\mathit{\boldsymbol{\varepsilon }}^n}. $

(19)

矩阵$ \mathit{\boldsymbol{\tilde{A}}} $_x、$ \mathit{\boldsymbol{\tilde{A}}} $_y、$ \mathit{\boldsymbol{\tilde{A}}} $_z是乘积可交换的.事实上，只需验证(15)~(17)中Kronecker积形式的矩阵都是乘积可交换的.由引理3和5，可得(I⊗I⊗A_α)(I⊗A_β⊗I)=(I⊗(I⊗A_α))(I⊗(A_β⊗I))=I⊗((I⊗A_α)(A_β⊗I))=I⊗(A_β⊗A_α)，同时(I⊗A_β⊗I)(I⊗I⊗A_α)=(I⊗(A_β⊗I))(I⊗(I⊗A_α))=I⊗((A_β⊗I)(I⊗A_α))=I⊗(A_β⊗A_α)，等等.由$ \mathit{\boldsymbol{\tilde{A}}} $_x、$ \mathit{\boldsymbol{\tilde{A}}} $_y、$ \mathit{\boldsymbol{\tilde{A}}} $_z乘积可交换并利用引理10，方程(19)可以写为εⁿ=((I－$ \mathit{\boldsymbol{\tilde{A}}} $_z)^－1(I+$ \mathit{\boldsymbol{\tilde{A}}} $_z))ⁿ((I－$ \mathit{\boldsymbol{\tilde{A}}} $_y)^－1(I+$ \mathit{\boldsymbol{\tilde{A}}} $_y))ⁿ((I－$ \mathit{\boldsymbol{\tilde{A}}} $_x)^－1(I+$ \mathit{\boldsymbol{\tilde{A}}} $_x))ⁿε⁰.由定理2和引理9可知‖(I－$ \mathit{\boldsymbol{\tilde{A}}} $_x)^－1(I+$ \mathit{\boldsymbol{\tilde{A}}} $_x)‖ < 1，所以(I－$ \mathit{\boldsymbol{\tilde{A}}} $_x)^－1(I+$ \mathit{\boldsymbol{\tilde{A}}} $_x)的谱半径小于1，因此当n→∞时，((I－$ \mathit{\boldsymbol{\tilde{A}}} $_x)^－1(I+$ \mathit{\boldsymbol{\tilde{A}}} $_x))ⁿ收敛到零矩阵，也就是对任意的n≥0，((I－$ \mathit{\boldsymbol{\tilde{A}}} $_x)^－1(I+$ \mathit{\boldsymbol{\tilde{A}}} $_x))ⁿ有界.同理可证((I－$ \mathit{\boldsymbol{\tilde{A}}} $_y)^－1(I+$ \mathit{\boldsymbol{\tilde{A}}} $_y))ⁿ和((I－$ \mathit{\boldsymbol{\tilde{A}}} $_z)^－1(I+$ \mathit{\boldsymbol{\tilde{A}}} $_z))ⁿ对任意的n≥0有界.这就证明了差分格式(9)是无条件稳定的.

定理4  设u(x_i, y_j, z_m, t_n)是问题(1)~(3)的解，u_{i, j, m}ⁿ是差分格式(9)的解.记e_{i, j, m}ⁿ=u(x_i, y_j, z_m, t_n)－u_{i, j, m}ⁿ，eⁿ=(([e_{i, j, m}ⁿ]_{i=1, 2, …，N1－1})_{j=1, 2, …，N2－1})_{m=1, 2, …，N3－1}，那么存在一个正常数C，使得‖eⁿ‖ < C(τ²+h₁²+h₂²+h₃²).

证明  记Rⁿ=(([R_{i, j, m}ⁿ]_{i=1, 2, …，N1－1})_{j=1, 2, …，N2－1})_{m=1, 2, …，N3－1}，那么方程(7)减去(9)的矩阵形式为

$ \left( {\mathit{\boldsymbol{I}} - {{\mathit{\boldsymbol{\tilde A}}}_x}} \right)\left( {\mathit{\boldsymbol{I}} - {{\mathit{\boldsymbol{\tilde A}}}_y}} \right)\left( {\mathit{\boldsymbol{I}} - {{\mathit{\boldsymbol{\tilde A}}}_z}} \right){\mathit{\boldsymbol{e}}^{n + 1}} = \left( {\mathit{\boldsymbol{I}} + {{\mathit{\boldsymbol{\tilde A}}}_x}} \right)\left( {\mathit{\boldsymbol{I}} + {{\mathit{\boldsymbol{\tilde A}}}_y}} \right)\left( {\mathit{\boldsymbol{I}} + {{\mathit{\boldsymbol{\tilde A}}}_z}} \right){\mathit{\boldsymbol{e}}^n} + {\mathit{\boldsymbol{R}}^{n + 1}}. $

(20)

因为$ \mathit{\boldsymbol{\tilde{A}}} $_x、$ \mathit{\boldsymbol{\tilde{A}}} $_y、$ \mathit{\boldsymbol{\tilde{A}}} $_z乘积可交换，根据引理10，方程(20)可以写为eⁿ⁺¹=(I－$ \mathit{\boldsymbol{\tilde{A}}} $_z)^－1(I+$ \mathit{\boldsymbol{\tilde{A}}} $_z)(I－$ \mathit{\boldsymbol{\tilde{A}}} $_y)^－1(I+$ \mathit{\boldsymbol{\tilde{A}}} $_y)·(I－$ \mathit{\boldsymbol{\tilde{A}}} $_x)^－1(I+$ \mathit{\boldsymbol{\tilde{A}}} $_x)eⁿ+(I－$ \mathit{\boldsymbol{\tilde{A}}} $_z)^－1(I－$ \mathit{\boldsymbol{\tilde{A}}} $_y)^－1(I－$ \mathit{\boldsymbol{\tilde{A}}} $_x)^－1Rⁿ⁺¹.用矩阵的2-范数作用上式两端，并利用引理8和9，可得$ \left\| {{\mathit{\boldsymbol{e}}}^{n}} \right\|<\left\| {{\mathit{\boldsymbol{e}}}^{0}} \right\|+\sum\limits_{k=0}^{n-1}{\left\| {{\mathit{\boldsymbol{R}}}^{k+1}} \right\|}$.因e⁰=0，所以存在正常数C，使得‖eⁿ‖ < C(τ²+h₁²+h₂²+h₃²).

3 数值结果

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

$ {\left\| {{\mathit{\boldsymbol{U}}_h} - \mathit{\boldsymbol{U}}} \right\|_{{L^\infty }}}: = \mathop {\max }\limits_{\begin{array}{*{20}{c}} {1 \le m \le {N_3} - 1}\\ {1 \le j \le {N_2} - 1}\\ {1 \le i \le {N_1} - 1} \end{array}} \left| {u_{i,j,m}^M - u\left( {{x_i},{y_j},{z_m},T} \right)} \right|,{\left\| {{\mathit{\boldsymbol{U}}_h} - \mathit{\boldsymbol{U}}} \right\|_{{L^2}}}: = {\left( {\sum\limits_{m = 1}^{{N_3} - 1}\\ {\sum\limits_{j = 1}^{{N_2} - 1} {\sum\limits_{i = 1}^{{N_1} - 1} {{{\left| {u_{i,j,m}^M - u\left( {{x_i},{y_j},{z_m},T} \right)} \right|}^2}{h_1}{h_2}{h_3}} } } } \right)^{1/2}}. $

算例1  在问题(1)~(3)中，取Ω=(0, 1)×(0, 1)×(0, 1)，T=1；对流和扩散系数分别为k₁=k₂=k₃=－1，a₁=a₂=b₁=b₂=c₁=c₂=1；初值取u₀(x, y, z)=x³(1－x)³y³(1－y)³z³(1－z)³.已知的解析解为u(x, y, z, t)=e^－tx³(1－x)³y³(1－y)³z³(1－z)³，由以上条件容易算出f(x, y, z, t).

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

算例2  在问题(1)~(3)中，取Ω=(0, 1)×(0, 1)×(0, 1)，T=1；对流和扩散系数分别为k₁=0.25x，k₂=0.25y，k₃=0.25z，a₁=x^α，a₂=(1－x)^α，b₁=y^β，b₂=(1－y)^β，c₁=z^γ，c₂=(1－z)^γ；初值取u₀(x, y, z)=x³(1－x)³y³(1－y)³z³(1－z)³.已知的解析解为u(x, y, z, t)=e^－tx³(1－x)³y³(1－y)³z³·(1－z)³，由以上条件f(x, y, z, t)容易算出.

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

4 结论

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