一种用于二维变温空间分数扩散方程的Crank-Nicolson ADI紧凑差分格式

《Computers & Mathematics with Applications》：A Crank-Nicolson ADI compact difference scheme for the two-dimensional tempered space-fractional diffusion equation

　　邱泽山|胡东东|李建新|江泰强
新疆政法学院基础教学部，中国图木舒克843900

摘要
本文关注二维双边Tempered空间分数阶扩散方程的高效数值方案。为此，我们结合Crank-Nicolson方法和紧凑差分方法，设计了一种时间上二阶、空间上四阶的交替方向隐式方案。所提出的方案被证明是无条件稳定且收敛的。最后，我们提供了一些数值结果来验证我们的理论发现以及所提出方法的有效性。

引言
本文数值研究了二维（2D）Tempered空间分数阶扩散方程：
$$
\partial_t u(x,y,t) = \alpha^\lambda_1 \partial_x u(x,y,t) + \beta^\lambda_2 \partial_y u(x,y,t) + f(x,y,t), \quad (x,y,t) \in \Omega \times (0,T),
u(x,y,0) = \phi(x,y), \quad (x,y) \in \Omega', \quad u(x,y,t) = \psi(x,y,t), \quad (x,y,t) \in \partial_{\Omega} \times [0,T],
$$
其中 $1 < \alpha, \beta < 2, \lambda_1, \lambda_2 > 0$, $\Omega = (a,b) \times (c,d)$。符号 $\partial_x^{\alpha,\lambda_1}$ 和 $\partial_y^{\beta,\lambda_2}$ 分别表示对 $u(x, y, t)$ 的双边Riemann-Liouville（RL）Tempered分数阶导数：
$$
\partial_x^{\alpha,\lambda_1} u(x,y,t) = l_x D_x^{-\alpha,\lambda_1} u(x,y,t) + r_x D_x^{+\alpha,\lambda_1} u(x,y,t),
\partial_y^{\beta,\lambda_2} u(x,y,t) = l_y D_y^{-\beta,\lambda_2} u(x,y,t) + r_y D_y^{+\beta,\lambda_2} u(x,y,t),
$$
其中 $D_z^{\pm\xi,\lambda}u(z)$ 是在[1],[2],[3]中提到的标准化Tempered分数阶导数，定义如下：
$$
D_z^{-\xi,\lambda}u(z) = \frac{D_z^{\xi,\lambda}u(z) - \lambda\xi u(z) - \xi\lambda^{1-\xi}\partial_u(z)\partial_z}{D_z^{\xi,\lambda}},
D_z^{\xi,\lambda}u(z) = \frac{D_z^{\xi,\lambda}u(z) + \lambda\xi u(z) + \xi\lambda^{1-\xi}\partial_u(z)\partial_z}{D_z^{\xi,\lambda}},
$$
这里 $z \in \{x, y\}$, $\xi \in \{\alpha, \beta\}$, $\lambda \in \{\lambda_1, \lambda_2\}$。Tempered分数阶导数 $D_z^{\pm\xi,\lambda}u(z)$ 的表达式为：
$$
D_z^{-\xi,\lambda}u(z) = e^{-\lambda z}\Gamma\left(2-\xi\right) \left(\frac{\partial^2}{\partial z^2}\int_{z-\xi}^{\tau}\int_{z}^{\tau}\lambda^{\tau}\tau u(\tau)(z-\tau)\xi^{-1} d\tau\right), \quad i = \{a\text{ if } z = x, \quad j = \{b\text{ if } z = y\},
$$
其中 $l_x, r_x, l_y, r_y$ 是非负常数，且 $l_x + r_x > 0$, $l_y + r_y > 0$。边界条件为：
$$
u(a,y,t) = 0 \text{ if } l_x \neq 0; \quad u(b,y,t) = 0 \text{ if } r_x \neq 0; \quad u(x,c,t) = 0 \text{ if } l_y \neq 0; \quad u(x,d,t) = 0 \text{ if } r_y \neq 0.
$$
在过去的几十年中，分数阶偏微分方程（FPDEs）已广泛应用于各个领域，包括等离子体物理[4],[5],[6]、地下水水文学[7],[8]、金融建模[9],[10]和生物系统[11]。由于分数阶导数的非局部性，获得FPDEs的解析解存在重大挑战。在某些情况下，这样的解甚至完全无法得到。因此，有必要开发相应的FPDEs的数值解。时间上的分数阶导数通常用Caputo导数表示，而空间上的分数阶导数大多用RL导数描述。Lubich的分数阶线性多步方法[12]及其改进版本在RL导数的近似中占据了主导地位。关于FPDEs的数值方案的研究已经取得了大量成果，读者可以参考[13],[14],[15],[16],[17],[18],[19],[20],[21],[22],[23],[24],[25],[26],[27],[28]及其中的参考文献。

近年来，作为分数阶导数的一种推广，Tempered分数阶导数成为研究的热点。由此产生的Tempered FPDEs在物理和数学上都具有优势[2]，并在金融[29],[30]、统计物理[31],[32]和地球物理学[33],[34]等领域得到广泛应用。通过扩展经典的分数阶导数近似方法，研究人员开发了一些有效的RL Tempered分数阶导数的数值方案，包括一阶Tempered移位Grünwald差分算子[1]、二阶Tempered加权移位Grünwald差分（TWSGD）算子[3]，以及三阶[35]和四阶[36] Tempered准紧凑差分算子。相关研究还涵盖了空间Tempered[37],[38],[39],[40],[41],[42],[43]和时空Tempered FPDEs[44],[45]。

在空间Tempered扩散问题的背景下：Baeumer和Meerschaert[1]提出了一种用于对流-扩散方程的一阶Tempered移位Grünwald差分算子；Li和Deng[3]结合二阶Tempered加权Grünwald差分算子和Crank-Nicolson（CN）方法来求解双边Tempered分数阶扩散方程；Chen和Deng[46]推导出了一种具有时间上一阶和空间上二阶精度的数值方案，用于涉及截断Lévy飞行的分数阶密度扩散方程。受分数阶扩散方程的准紧凑差分方法的启发，Yu等人[35]和Qiu[36]分别开发了Tempered空间分数阶扩散方程的准紧凑差分方案。然而，由于Tempered分数阶算子的兼容性问题，这些研究仅适用于单边方程，仅能达到时间上一阶精度和空间上s阶（$s = 3,4$）精度。Zhao等人[39]研究了非线性变系数Tempered分数阶扩散方程，提出了非线性和线性化隐式方案，并为相应的非线性一次性系统开发了预处理策略。此外，关于Tempered分数阶扩散方程产生的离散方程的快速解算法可以在[42],[43]中找到。

据我们所知，关于双边Tempered分数阶扩散方程的四阶高精度数值方案的研究仍然相对较少。值得注意的是，Guo等人[32]采用了[47]中提出的加权移位Lubich差分（WSLD）算子的思想，构建了一个用于双边空间Tempered分数阶扩散方程的四阶方案。值得注意的是，该方案依赖于计算域外的网格点，其数值性能不如避免使用这些外部点的方案。

本文重点研究二维Tempered分数阶扩散方程，并提出了一种具有时间上二阶精度和空间上四阶精度的数值方案。与现有方法相比，我们提出的方案不仅精度高，而且适用于双边和单边扩散情况。此外，所提出的数值方案在计算上更为简单。

本文的其余部分组织如下。第2节首先介绍一些初步知识，并提出一个四阶Tempered紧凑差分（TCD）近似公式；第3节讨论数值方法，分为两部分，包括数值方案的构建及其稳定性和收敛性的分析；第4节设计了一些数值实验来验证数值方案的有效性；第5节总结了本文。

**初步知识**
定义分数阶Sobolev空间 $F_{\lambda^{n}+\alpha}(R)$：
$$
F_{\lambda^{n}+\alpha}(R) = \{ \nu \mid \nu \in L^1(R), \int_{R} (\vert \lambda\| + \vert w \|)^{n+\alpha} |\nu^w| dw < \infty \},
$$
其中 $\nu^w$ 表示函数 $\nu(x)$ 的傅里叶变换，形式为 $\nu^w = \int_{R} \nu(x)e^{-iwx} dx$。

**引理2.1（[1, 3]）** 设 $1 < \alpha < 2, \lambda \geq 0$，移位数 $p$ 是整数，$h$ 是步长，$\nu(x)$ 定义在有界区间 $[a, b]$ 上，并在 $x \in (-\infty, a) \cup (b, +\infty)$ 上经过零扩展后属于 $F_{\lambda^{n}+\alpha}(R)$。Tempered和移位Grünwald型差分算子定义为：

**数值方案的推导**
为了推导求解方程（1）的数值方案，我们对空间网格和时间网格进行划分：设 $(x_i, y_j) \mid x_i = a + ihx, y_j = c + jhy, \quad h_x = b - a, \quad h_y = d - c, 0 \leq i \leq M_x, 0 \leq j \leq M_y$，$T_n = n\tau, \quad \tau = TN, 0 \leq n \leq N$。为方便后续讨论，我们定义一些符号：
$$
t_{n+1/2} = \frac{t_n + t_{n+1}}{2}, \quad u_{i,j}^n = u(x_i, y_j, t_n), \quad f_{i,j}^{n+1/2} = f(x_i, y_j, t_{n+1/2), \quad \delta_t u_{i,j}^{n+1/2} = u_{i,j}^{n} - u_{i,j}^n \tau.
$$
将Tempered紧凑差分算子 $\Psi_x^{\alpha,\lambda_1}$ 和 $\Psi_y^{\beta,\lambda_2}$ 应用于方程（1）的两边，得到：
$$
\Psi_x^{\alpha,\lambda_1} \Psi_y^{\beta,\lambda_2} \partial_t u(x,y,t) = \Psi_x^{\alpha,\lambda_1} \Psi_y^{\beta,\lambda_2} \partial_t u(x,y,t).
$$

**数值实验**
在本节中，我们数值研究了二维空间Tempered分数阶扩散方程。所有数值实验都在配备Intel Core i5-10400 CPU和8 GB RAM的桌面计算机上使用MATLAB R2016a进行。所有数值示例都采用了ADI方法来降低计算成本。通过一系列数值实验，评估了所提方案的时间和空间精度。此外，我们将我们的方法与CN-TWSGD方法[3]进行了比较。

**结论**
本文开发了一个用于二维Tempered空间分数阶扩散方程的四阶ADI紧凑差分方案。设计此类方案的主要困难在于左右Tempered分数阶导数之间的不兼容性，这是由于Tempered分数阶导数中固有的指数因子造成的。这种不兼容性限制了之前的研究仅适用于单边问题。为了解决这个问题，我们直接构建了一个四阶紧凑差分方案。

**利益声明**
作者没有利益冲突。