高维时间分数阶扩散方程的快速高效数值格式

A Fast and Efficient Numerical Scheme for High-Dimensional Time-Fractional Order Diffusion Equations

摘要: 本文研究了关于高维时间分数阶扩散方程的高效数值算法问题。提出了两种数值格式，一种为改进的L2数值格式，该格式有效地解决了普通L2格式在初始时刻会出现收敛精度下降的问题，并在空间方向上运用紧致差分算子进行离散，使得该格式的收敛精度达到时间上$3-\alpha$ 阶，空间上4阶。另一种在改进L2格式的基础上，运用指数和方法去逼近时间分数阶导数中的幂函数，从而极大地减少了计算成本，并且保持同样的高收敛精度。针对具有初值奇异性的问题，我们采用梯度网格方法去解决。最后通过多个数值结果证明了这两种算法的有效性和准确性，并展示所提出快速化算法的计算效率。

Abstract: In this paper, the problem of efficient numerical algorithms on high-dimensional time fractional order diffusion equations is investigated. Two numerical schemes are proposed, one is the improved L2 numerical scheme, which efficiently solves the problem of decreasing convergence accuracy at the initial moment of the ordinary L2 scheme, and discretizes it in the spatial direction by applying the compact difference operator, so that the convergence accuracy of this scheme reaches the order in time and the order in space. The other is to improve the L2 scheme by applying the sum of exponentials method to approximate the power function in the fractional order derivatives in time, which greatly reduces the computational cost and maintains the same high convergence accuracy. For problems with initial value singularities, we use the graded meshes method to solve them. Finally, the effectiveness and accuracy of both algorithms are demonstrated through several numerical results and the computational efficiency of the proposed fastening algorithm is shown.

1. 引言

近年来，时间分数阶扩散方程因其能够刻画非局部性和记忆性等复杂特性，在科学研究和工程应用中受到了广泛关注。保持高维时间分数阶扩散方程在数值求解过程中，数值解稳定且具有高精度是相当困难的，特别是解析解的初值具有奇异性的情况。

时间分数阶扩散方程描述了与$0\le \alpha \le 1$ 对应的反常亚扩散现象，由于分数阶导数的数值格式的计算复杂度一般是巨大的。因此，构建分数阶导数的高精度数值格式是求解时间分数阶扩散方程的常用方法，由于紧差分格式具有高精度收敛性，在结合高阶时间数值方法求解时间分数阶扩散方程问题中得到了广泛的应用。在[1]利用时间L1公式建立了一种高效的时间分数阶扩散方程的数值格式。在[2]中，采用时间L1和L2格式构建了时间分数阶扩散方程的高阶稳定性与收敛性数值方法，并进行了严格分析。在[3]中，采用时间L2格式求解时间分数阶扩散方程，但在第一层精度较低。在[4]中，采用时间L2格式，提出两种高阶紧致差分方案去求解Sobolev型一维与二维多项式时间分数阶对流扩散方程，针对解具有初始奇异性问题，采用了梯度网格方案去解决，并通过算例验证了所提方案的精度与效率。在[5]中，采用时间L1格式，构造了一个四阶紧致ADI格式来求解二维时间分数反应–亚扩散方程。在[6]中，采用空间弱Galerkin有限元法，采用梯度网格上的L1方法求解具有弱奇点的多项时间分数阶扩散方程。在[7]中，提出快速改进的L1格式来处理时间分数阶扩散方程的解初始奇异性问题。在[8]中，运用正交分解技术提出了一种时间分数扩散方程的降阶有限差分外推算法，有效地降低了计算成本。在[9]中，提出了运用指数和(SOE)去近似时间分数阶导数中的幂函数，有效地降低了计算成本，并在时间上的精度达到高阶收敛。

根据以上的研究理论基础，本文研究了关于高维时间分数阶扩散方程的高效数值算法，主要构造数值格式，使数值解能具有较高的精度，可保持稳定，并使用快速化算法降低计算成本。本文的主要贡献可分为以下几点：(1) 提出一种改进的L2格式，去处理普通L2格式在第一层精度降低的情况，使得时间收敛阶在整个时间区域上都能达到$3-\alpha$ 阶。同时，空间上结合紧致差分格式，使得该格式的空间收敛精度达到4阶。(2) 在改进的L2格式基础上，运用指数和(SOE)方法去加快计算的速度，降低计算所需要的存储空间，并且能保持同样的收敛精度。(3) 针对具有初值奇异性的问题，结合梯度网格方法去解决，得出的结果仍然能达到理论的收敛精度。

2. 问题描述

本文我们将对以下高维时间分数阶扩散方程进行研究，

$\{\begin{array}{l}{}_0^C{D}_t^{\alpha }u\left(x,t\right)-{\displaystyle \sum _{i=1}^d{\kappa }_i{\partial }_{x_i}^{2}u\left(x,t\right)=f\left(x,t\right)},x\in \Omega ,t\in \left(0,T\right]\\ u\left(x,t\right)=0,x\in \partial \Omega ,t\in \left(0,T\right],u\left(x,0\right)=u_0\left(x\right),x\in \Omega \end{array}$ (1)

其中$x=\left(x_1,x_2,\cdots ,x_d\right)\in {\Re }^d$ ，且$d\ge 1$ 是空间维数，$\Omega ={\displaystyle \prod _{i=1}^d\left(a_i,b_i\right)}\subset {\Re }^d$ 是有界的，$T$ 是正的有界常数，$\Omega$ 的边界记为$\partial \Omega$ ，$\overline{\Omega }$ 表示$\Omega$ 的闭包，${\kappa }_i$ 是扩散常数，且${\kappa }_i>0\left(i=1,2,\cdots ,d\right)$ 。${}_0^C{D}_t^{\alpha }$ 是Caputo分数阶导数，且$\alpha \in \left(0,1\right)$ ，其定义如下[10]：

${}_0^C{D}_t^{\alpha }u\left(x,t\right)=\frac{1}{\Gamma \left(1-\alpha \right)}{\displaystyle {\int }_0^t{\left(t-s\right)}^{-\alpha }}\frac{\partial u\left(x,s\right)}{\partial s}\text{d}s$

其中$\Gamma (\cdot )$ 为欧拉伽马函数。

为了构造全离散数值格式，我们设$M_i\left(i=1,2,\cdots ,d\right)$ 和$N$ 均为正整数，时空网格的大小定义为$\Delta x_i=\left(b_i-a_i\right)/M_i$ ，$\tau =T/N$ ，定义时空离散点为$x_{i,l}=a_i+l\Delta x_i\left(0\le l\le M_i,1\le i\le d\right)$ ，$t_n=n\tau \left(0\le n\le N\right)$ 。对于$1\le i\le d$ ，我们定义${\Theta }_i=\left\{l\in N^{+}|1\le l\le M_i-1\right\}$ 和${\ddot{\Theta }}_i=\left\{l\in N^{+}|0\le l\le M_i\right\}$ ，$\Psi ={\displaystyle \prod _{i=1}^d{\Theta }_i}$ ，$\ddot{\Psi }={\displaystyle \prod _{i=1}^d{\ddot{\Theta }}_i}$ 以及$\partial \ddot{\Psi }=\ddot{\Psi }\backslash \Psi$ 。为简单起见，我们设置一个多索引$L=\left(l_1,l_2,\cdots ,l_d\right)$ ，并记$x_L=\left(x_{1,l_1},x_{2,l_2},\cdots ,x_{d,l_d}\right)$ ，网格函数记为$u_L^n=u\left(x_L,t_n\right)$ ，$f_L^n=f\left(x_L,t_n\right)$ ，其中$L\in \ddot{\Psi }$ ，$0\le n\le N$ 。

我们先对空间方向上进行离散，基于文献[11]的思想，我们定义$x_i\left(i=1,2,\cdots ,d\right)$ 方向上的紧致差分算子如下：

$H_i{u}_L=\{\begin{array}{l}\left(I+\frac{\Delta x_i^2}{12}{\delta }_{x_i}^2\right)u_L,1\le i\le d,L\in \Psi \\ u_L,L\in \partial \ddot{\Psi }\end{array}$ (2)

其中$I$ 为单位矩阵。

3. 改进的L2格式及误差估计

我们现在构造一个高效的数值格式去离散时间上的分数阶导数。由于L2格式在$t_1$ 时刻使用拉格朗日一次插值去近似时间分数阶导数，导致在$t_1$ 上的收敛精度只有$2-\alpha$ 阶，我们提出的改进的L2格式将时间的前两层进行耦合求解，即在$t_1$ 时刻利用拉格朗日二次插值，结合$t_2$ 时刻建立耦合求解，使得其收敛阶能够保持为$3-\alpha$ 阶。首先，我们确定$u\left(x,t\right)$ 在$t_1$ 和$t_2$ 上的值，在区间$\left[t_0,t_2\right]$ 上，$u\left(x,t\right)$ 的近似公式由拉格朗日二次插值给出：

$J_{\left[t_0,t_2\right]}u\left(x,t\right)=u\left(x,t_0\right){\overline{\eta }}_{0,0}\left(t\right)+u\left(x,t_1\right){\overline{\eta }}_{1,0}\left(t\right)+u\left(x,t_2\right){\overline{\eta }}_{2,0}\left(t\right)$

其中${\overline{\eta }}_{i,0}\left(t\right)$ 定义为

${\overline{\eta }}_{0,0}\left(t\right)=\frac{\left(t-t_1\right)\left(t-t_2\right)}{\left(t_0-t_1\right)\left(t_0-t_2\right)},{\overline{\eta }}_{1,0}\left(t\right)=\frac{\left(t-t_0\right)\left(t-t_2\right)}{\left(t_1-t_0\right)\left(t_1-t_2\right)},{\overline{\eta }}_{2,0}\left(t\right)=\frac{\left(t-t_0\right)\left(t-t_1\right)}{\left(t_2-t_0\right)\left(t_2-t_1\right)}$

那么Caputo分数阶导数${}_0^C{D}_t^{\alpha }$ 在$t_1$ 和$t_2$ 时刻的${}_0^C{D}_t^{\alpha }u\left(x,t_1\right)$ 和${}_0^C{D}_t^{\alpha }u\left(x,t_2\right)$ ，可以被近似为

$\begin{array}{l}{}_0^C{D}_t^{\alpha }u\left(x,t_1\right)\approx {}_0^C{D}_{\tau }^{\alpha }u\left(x,t_1\right)=\frac{1}{\Gamma \left(1-\alpha \right)}{\displaystyle {\int }_0^{t_1}{\left(t_1-s\right)}^{-\alpha }}\frac{J_{\left[t_0,t_2\right]}u\left(x,s\right)}{\partial s}\text{d}s\\ =P_1^{0,0}u\left(x,t_0\right)+P_1^{1,0}u\left(x,t_1\right)+P_1^{2,0}u\left(x,t_2\right)\end{array}$ (3)

$\begin{array}{l}{}_0^C{D}_t^{\alpha }u\left(x,t_2\right)\approx {}_0^C{D}_{\tau }^{\alpha }u\left(x,t_2\right)=\frac{1}{\Gamma \left(1-\alpha \right)}{\displaystyle {\int }_0^{t_2}{\left(t_2-s\right)}^{-\alpha }}\frac{J_{\left[t_0,t_2\right]}u\left(x,s\right)}{\partial s}\text{d}s\\ =P_2^{0,0}u\left(x,t_0\right)+P_2^{1,0}u\left(x,t_1\right)+P_2^{2,0}u\left(x,t_2\right)\end{array}$ (4)

其中

$P_1^{j,0}=\frac{1}{\Gamma \left(1-\alpha \right)}{\displaystyle {\int }_0^{t_1}{\left(t_1-s\right)}^{-\alpha }}\frac{{\overline{\eta }}_{j,0}\left(s\right)}{\partial s}\text{d}s,j=0,1,2$

$P_2^{j,0}=\frac{1}{\Gamma \left(1-\alpha \right)}{\displaystyle {\int }_0^{t_2}{\left(t_2-s\right)}^{-\alpha }}\frac{{\overline{\eta }}_{j,0}\left(s\right)}{\partial s}\text{d}s,j=0,1,2$

对于$t_n\left(n\ge 3\right)$ 时刻的分数阶导数离散，我们将积分区间分为$n$ 等分，记${\lambda }_k=\left[t_k,t_{k+1}\right]$ 。$\left(k=1,2,\cdots ,n-1\right)$ 类似地，我们利用拉格朗日插值去近似函数$u\left(x,t_n\right)$ ，则有

$\begin{array}{l}{}_0^C{D}_t^{\alpha }u\left(x,t_n\right)=\frac{1}{\Gamma \left(1-\alpha \right)}{\displaystyle {\int }_0^{t_n}{\left(t_n-s\right)}^{-\alpha }}\frac{\partial u\left(x,s\right)}{\partial s}\text{d}s\\ =\frac{1}{\Gamma \left(1-\alpha \right)}\left[{\displaystyle {\int }_{t_{n-1}}^{t_n}{\left(t_n-s\right)}^{-\alpha }}\frac{\partial u\left(x,s\right)}{\partial s}\text{d}s+{\displaystyle \sum _{k=1}^{n-1}{\displaystyle {\int }_{t_{k-1}}^{t_k}{\left(t_n-s\right)}^{-\alpha }}\frac{\partial u\left(x,s\right)}{\partial s}\text{d}s}\right]\\ \approx {}_0^C{D}_{\tau }^{\alpha }u\left(x,t_n\right)=\frac{1}{\Gamma \left(1-\alpha \right)}\left[{\displaystyle {\int }_{t_{n-1}}^{t_n}{\left(t_n-s\right)}^{-\alpha }}\frac{J_{\left[t_{n-1},t_n\right]}u\left(x,s\right)}{\partial s}\text{d}s+{\displaystyle \sum _{k=1}^{n-1}{\displaystyle {\int }_{t_{k-1}}^{t_k}{\left(t_n-s\right)}^{-\alpha }}\frac{J_{\left[t_{k-1},t_k\right]}u\left(x,s\right)}{\partial s}\text{d}s}\right]\\ =\frac{1}{\Gamma \left(1-\alpha \right)}\left[P_n^{0,n}u\left(x,t_{n-2}\right)+P_n^{1,n}u\left(x,t_{n-1}\right)+P_n^{2,n}u\left(x,t_n\right)+{\displaystyle \sum _{k=1}^{n-1}{P}_n^{0,k}u\left(x,t_k\right)+P_n^{1,k}u\left(x,t_k\right)+P_n^{2,k}u\left(x,t_k\right)}\right]\end{array}$ (5)

其中

$J_{\left[t_{k-1},t_k\right]}u\left(x,t\right)=u\left(x,t_{k-1}\right){\overline{\eta }}_{0,k}\left(t\right)+u\left(x,t_k\right){\overline{\eta }}_{1,k}\left(t\right)+u\left(x,t_{k+1}\right){\overline{\eta }}_{2,k}\left(t\right)$

这里的${\overline{\eta }}_{j,k}\left(j=0,1,2,k=1,2,\cdots ,n-1\right)$ 表示为

${\overline{\eta }}_{0,k}\left(t\right)=\frac{\left(t-t_k\right)\left(t-t_{k+1}\right)}{\left(t_{k-1}-t_k\right)\left(t_{k-1}-t_{k+1}\right)},{\overline{\eta }}_{1,k}\left(t\right)=\frac{\left(t-t_{k-1}\right)\left(t-t_{k+1}\right)}{\left(t_k-t_{k-1}\right)\left(t_k-t_{k+1}\right)},{\overline{\eta }}_{2,k}\left(t\right)=\frac{\left(t-t_{k-1}\right)\left(t-t_k\right)}{\left(t_{k+1}-t_{k-1}\right)\left(t_{k+1}-t_k\right)}$

$P_n^{j,k}\left(j=0,1,2,k=1,2,\cdots ,n-1\right)$ 表示为

$P_n^{j,n}=\frac{1}{\Gamma \left(1-\alpha \right)}{\displaystyle {\int }_{t_{n-1}}^{t_n}{\left(t_n-s\right)}^{-\alpha }}\frac{{\overline{\eta }}_{j,n}\left(s\right)}{\partial s}\text{d}s$

$P_n^{j,k}=\frac{1}{\Gamma \left(1-\alpha \right)}{\displaystyle {\int }_{t_{k-1}}^{t_k}{\left(t_n-s\right)}^{-\alpha }}\frac{{\overline{\eta }}_{j,k}\left(s\right)}{\partial s}\text{d}s$

结合式(3)，(4)和(5)，为简化我们的数值格式，并通过积分求出$P_n^{j,k}\left(j=0,1,2,k=1,2,\cdots ,n-1\right)$ 的值，同时令$\sigma =\Gamma \left(3-\alpha \right){\tau }^{\alpha }$ ，我们有

${}_0^C{D}_{\tau }^{\alpha }u\left(x,t_n\right)=\{\begin{array}{l}{\sigma }^{-1}\left(z_{0}u\left(x,t_0\right)+z_{1}u\left(x,t_1\right)+z_{2}u\left(x,t_2\right)\right),n=1\\ {\sigma }^{-1}\left({\widehat{z}}_{0}u\left(x,t_0\right)+{\widehat{z}}_{1}u\left(x,t_1\right)+{\widehat{z}}_{2}u\left(x,t_2\right)\right),n=2\\ {\sigma }^{-1}[a_{n}u\left(x,t_{n-2}\right)+b_{n}u\left(x,t_{n-1}\right)+c_{n}u\left(x,t_n\right)\begin{array}{c}\\ \end{array}\\ +{\displaystyle \sum _{k=1}^{n-1}{a}_{k}u\left(x,t_{n-k-1}\right)+b_{k}u\left(x,t_{n-k}\right)+c_{k}u\left(x,t_{n-k+1}\right)}],n\ge 3\end{array}$ (6)

其中

$\begin{array}{l}{z}_0=\left(3\alpha -4\right)/2,z_1=2\left(1-\alpha \right),z_2=\alpha /2\\ {\widehat{z}}_0=3\alpha -2/2^{\alpha },{\widehat{z}}_1=-4\alpha /2^{\alpha },{\widehat{z}}_2=\left(\alpha +2\right)/2^{\alpha }\\ a_n=\alpha /2,b_n=-2,c_n=\left(4-\alpha \right)/2\\ a_k=-3\left(2-\alpha \right){\left(k+1\right)}^{1-\alpha }/2+\left(2-\alpha \right)k^{1-\alpha }/2+{\left(k+1\right)}^{2-\alpha }-k^{2-\alpha }\\ b_k=2\left(2-\alpha \right){\left(k+1\right)}^{1-\alpha }-2{\left(k+1\right)}^{2-\alpha }+2k^{2-\alpha }\\ c_k=-\left(2-\alpha \right)\left({\left(k+1\right)}^{1-\alpha }+k^{1-\alpha }\right)/2+{\left(k+1\right)}^{2-\alpha }-k^{2-\alpha }\end{array}$ (7)

引理1 [12] 假设$u\left(x,t\right)\in C^3\left[0,T\right]$ ，且令$r_n\left(\tau \right)={}_0^C{D}_t^{\alpha }u\left(x,t_n\right)-{}_0^C{D}_{\tau }^{\alpha }u\left(x,t_n\right)$ ，$\forall n\ge 1$ ，则有

$\left|r_n\left(\tau \right)\right|\le C_u{\tau }^{3-\alpha }$ (8)

其中$C_u$ 是一个仅关于函数$u$ 的常数。

现在，我们将利用紧致差分算子离散空间上的导数项。我们先介绍以下引理。

引理2 [13] 假设$u\left(x\right)\in C^6\left[a,b\right]$ ，且有$\zeta \left(s\right)=5{\left(1-s\right)}^3-3{\left(1-s\right)}^5$ ，则有

$H_{\Delta x}{f}^{″}\left(x_j\right)={\delta }_x^{2}f\left(x_j\right)\frac{\Delta x^4}{360}{\displaystyle {\int }_0^1\left[f^{\left(6\right)}\left(x_j-s\Delta x\right)+f^{\left(6\right)}\left(x_j+s\Delta x\right)\right]\zeta \left(s\right)}$ (9)

其中$1\le j\le M-1$ 。

我们将点$\left(x_L,t_n\right)$ 带入方程(1)，可得

${}_0^C{D}_t^{\alpha }u\left(x_L,t_n\right)-{\displaystyle \sum _{i=1}^d{\kappa }_i{\partial }_{x_i}^{2}u\left(x_L,t_n\right)}=f\left(x_L,t_n\right)$ (10)

我们在式(10)两边都作用上每个空间方向$x_i$ 的紧致差分算子，记这个高维紧致差分算子为${\displaystyle \prod _{m=1}^d{H}_m}$ ，那么有

${\displaystyle \prod _{m=1}^d{H}_m}\left({}_0^C{D}_t^{\alpha }u\left(x_L,t_n\right)\right)-{\displaystyle \prod _{m=1}^d{H}_m}\left({\displaystyle \sum _{i=1}^d{\kappa }_i{\partial }_{x_i}^{2}u\left(x_L,t_n\right)}\right)={\displaystyle \prod _{m=1}^d{H}_m}f\left(x_L,t_n\right)$ (11)

假设函数$u\left(x,t\right)\in C_{6,3}^{x,t}\left(\Omega \times \left[0,T\right]\right)$ ，根据引理1，引理2和式(2)可以得到方程(1)的全离散格式如下：

${\displaystyle \prod _{m=1}^d{H}_m}\left({}_0^C{D}_{\tau }^{\alpha }{u}_L^n\right)-{\displaystyle \sum _{i=1}^d{\kappa }_i}{\displaystyle \prod _{m=1,m\ne i}^d{H}_m}{\delta }_{x_i}^2{u}_L^n={\displaystyle \prod _{m=1}^d{H}_m}{f}_L^n+R_L^n$ (12)

其中$R_L^n$ 满足以下表达式：

$R_L^n={\displaystyle \sum _{i=1}^d{\kappa }_i}{\displaystyle \prod _{m=1,m\ne i}^d{H}_m}{\displaystyle {\int }_0^1\left[{\partial }_{x_i}^{6}u\left(x_L-{\widehat{s}}_i,t_n\right)+{\partial }_{x_i}^{6}u\left(x_L+{\widehat{s}}_i,t_n\right)\right]\times {\left(1-s_i\right)}^3\left[5-3{\left(1-s_i\right)}^2\right]\text{d}{s}_i}-{\displaystyle \prod _{m=1}^d{H}_m}\left(r_n\left(x_L,t_n\right)\right)$

这里${\widehat{s}}_i=s_i\Delta x_i{e}_i$ ，$e_i$ 为第$i$ 个空间方向的标准正交基，并且$R_L^n$ 满足以下条件

$\left|R_L^n\right|\le C_R\left({\tau }^{3-\alpha }+{\displaystyle \sum _{i=1}^d\Delta }{x}_i^4\right)$ (13)

这里的$C_R$ 是一个不依赖于$\tau$ 和$\Delta x_i\left(i=1,2,\cdots ,d\right)$ 的常数。

为了简化，我们记高维紧致差分格式为

${\Lambda }_{\Delta x}{u}_L={\displaystyle \sum _{i=1}^d{\kappa }_i}{\displaystyle \prod _{m=1,m\ne i}^d{H}_m}{\delta }_{x_i}^2{u}_L,x_L\in \ddot{\Psi }$ (14)

考虑到，当我们的时空步长取得足够小时，即$\tau \to 0$ ，$\Delta x_i\to 0$ 时，有$R_L^n\to 0$ 。再结合式(13)的简化，我们得到方程(1)的一个改进的L2格式如下：

$\{\begin{array}{l}{\displaystyle \prod _{m=1}^d{H}_m}\left(z_0{u}_L^0+z_1{u}_L^1+z_2{u}_L^2\right)-\sigma {\Lambda }_{\Delta x}{u}_L^1=\sigma {\displaystyle \prod _{m=1}^d{H}_m{f}_L^1},n=1\\ {\displaystyle \prod _{m=1}^d{H}_m}\left({\widehat{z}}_0{u}_L^0+{\widehat{z}}_1{u}_L^1+{\widehat{z}}_2{u}_L^2\right)-\sigma {\Lambda }_{\Delta x}{u}_L^2=\sigma {\displaystyle \prod _{m=1}^d{H}_m{f}_L^2},n=2\\ {\displaystyle \prod _{m=1}^d{H}_m}[a_n{u}_L^{n-2}+b_n{u}_L^{n-1}+c_n{u}_L^n\begin{array}{c}\\ \end{array}\\ +{\displaystyle \sum _{k=1}^{n-1}{a}_k{u}_L^{n-k-1}+b_k{u}_L^{n-k}+c_k{u}_L^{n-k+1}}]-\sigma {\Lambda }_{\Delta x}{u}_L^n=\sigma {\displaystyle \prod _{m=1}^d{H}_m{f}_L^n},n\ge 3\end{array}$ (15)

4. 快速化L2格式及误差估计

在本节中，我们主要是在以上提出的改进的L2数值格式基础上，运用指数和法(SOE)去近似分数阶导数中的幂函数${\left(t-s\right)}^{-\alpha -1}$ ，将原来L2格式离散分数阶导数的储存成本和总计算成本分别从$\mathcal{O}\left(N_T\right)$ 和$\mathcal{O}\left(N_T^2\right)$ 减少到$\mathcal{O}\left(N_{\epsilon }\right)$ 和$\mathcal{O}\left(N_T{N}_{\epsilon }\right)$ ，能够极大地减少储存成本和计算时间。

接下来，我们来简单介绍一下，SOE快速化算法是如何作用在L2格式上的。基于数值格式的局限性，我们对于求解当$n=1$ 和$n=2$ 时，不使用快速化算法，仅在$n\ge 3$ 的情况下使用快速化算法。在式(14)中$n\ge 3$ 的情况涉及到前面所有解的求和，这反映了非局部的时间分数阶导数具有记忆性。参考文献[9]，我们将Caputo导数分为当前项$L\left(t_n\right)$ 和历史项$H\left(t_n\right)$ ，其定义如下：

$L\left(t_n\right):=\frac{1}{\Gamma \left(1-\alpha \right)}{\displaystyle {\int }_{t_{n-1}}^{t_n}\frac{{\partial }_{s}u\left(s\right)}{{\left(t_n-s\right)}^{\alpha }}\text{d}s},H\left(t_n\right):=\frac{1}{\Gamma \left(1-\alpha \right)}{\displaystyle {\int }_0^{t_{n-1}}\frac{{\partial }_{s}u\left(s\right)}{{\left(t_n-s\right)}^{\alpha }}\text{d}s}$

当前项$\overline{L}\left(t_n\right)$ 利用L2格式离散，并由式(7)计算可得

$\overline{L}\left(t_n\right)={\sigma }^{-1}\left[\frac{\alpha }{2}u\left(x,t_{n-2}\right)-2u\left(x,t_{n-1}\right)+\frac{4-\alpha }{2}u\left(x,t_n\right)\right]$ (16)

对于历史项$\overline{H}\left(t_n\right)$ ，利用分部积分公式可以得到

$\overline{H}\left(t_n\right)=\frac{1}{\Gamma \left(1-\alpha \right)}\left[\frac{u\left(x,t_{n-1}\right)}{{\tau }^{\alpha }}-\frac{u\left(x,t_0\right)}{t_n^{\alpha }}-\alpha {\displaystyle {\int }_0^{t_{n-1}}\frac{u\left(x,s\right)\text{d}s}{{\left(t_n-s\right)}^{1+\alpha }}}\right]$ (17)

为了近似历史项$\overline{H}\left(t_n\right)$ 最后一项中的幂函数${\left(t-s\right)}^{-\alpha -1},\forall \alpha \in \left(0,1\right)$ 。我们介绍以下引理

引理3 [14] 在区间$\left[\tau ,T\right]$ 上，给定任意规定的误差$\epsilon$ ，则存在正实数$s_i$ 和${\omega }_i\left(i=1,2,\cdots ,N_{\epsilon }\right)$ ，使得对于$\forall \alpha \in \left(0,1\right)$ ，有

$\left|\frac{1}{t^{1+\alpha }}-{\displaystyle \sum _{i=1}^{N_{\epsilon }}{\omega }_i}{\text{e}}^{-s_{i}t}\right|\le \epsilon ,t\in \left[\tau ,T\right]$ (18)

其中

$N_{\epsilon }=\mathcal{O}\left(\log \frac{1}{\epsilon }\left(\log \log \frac{1}{\epsilon }+\log \frac{T}{\tau }\right)+\log \frac{1}{\tau }\left(\log \log \frac{1}{\epsilon }+\log \frac{1}{\tau }\right)\right)$ (19)

对于固定的$\epsilon$ ，如果$T\gg 1$ ，则$N_{\epsilon }=\mathcal{O}\left(\log \left(N_T\right)\right)$ ，当$T\approx 1$ 时，$N_{\epsilon }=\mathcal{O}\left({\log }^2\left(N_T\right)\right)$ 。

利用引理3的结论去近似式(17)中最后一项的幂函数${\left(t-s\right)}^{-\alpha -1}$ ，那么式(17)可以变为

$\overline{H}\left(t_n\right)\approx \frac{1}{\Gamma \left(1-\alpha \right)}\left[\frac{u\left(x,t_{n-1}\right)}{{\tau }^{\alpha }}-\frac{u\left(x,t_0\right)}{t_n^{\alpha }}-\alpha {\displaystyle \sum _{i=1}^{N_{\epsilon }}{\omega }_i{U}_{h,i}^n\left(x,t_n\right)}\right]$ (20)

其中

$U_{h,i}^n\left(x,t_n\right)={\displaystyle {\int }_0^{t_{n-1}}{\text{e}}^{-s_i\left(t_n-s\right)}u\left(x,s\right)\text{d}s}$ (21)

这个方法的关键在于，存在以下递推关系式

$U_{h,i}^n\left(x,t_n\right)={\text{e}}^{-s_i\tau }{U}_{h,i}^n\left(x,t_{n-1}\right)+{\displaystyle {\int }_{t_{n-2}}^{t_{n-1}}{\text{e}}^{-s_i\left(t_n-s\right)}}u\left(x,s\right)\text{d}s$ (22)

从而使得，在计算$U_{h,i}^n\left(x,t_n\right)$ 时，不用算整个区间$\left(0,t_{n-1}\right)$ 上的积分，而只要算${\displaystyle {\int }_{t_{n-2}}^{t_{n-1}}{\text{e}}^{-s_i\left(t_n-s\right)}}u\left(x,s\right)\text{d}s$ 。这个式子，我们可以利用拉格朗日二次插值$J_{\left[t_{n-2},t_{n-1}\right]}u\left(x,t\right)$ 去近似计算，那么有

$\begin{array}{l}{\displaystyle {\int }_{t_{n-2}}^{t_{n-1}}{\text{e}}^{-s_i\left(t_n-s\right)}}u\left(s\right)\text{d}s\approx {\displaystyle {\int }_{t_{n-2}}^{t_{n-1}}{\text{e}}^{-s_i\left(t_n-s\right)}}{J}_{\left[t_{n-2},t_{n-1}\right]}u\left(x,s\right)\text{d}s\\ =A_{i}u\left(t_{n-2}\right)-B_{i}u\left(t_{n-1}\right)+C_{i}u\left(t_n\right)\end{array}$ (23)

其中

$\begin{array}{l}{A}_i=\frac{{\text{e}}^{-2s_i\tau }}{2s_i{\left(s_i\tau \right)}^2}\left[2{\text{e}}^{s_i\tau }-2{\left(s_i\tau \right)}^2-s_i\tau \left(3-{\text{e}}^{s_i\tau }\right)-2\right]\\ B_i=\frac{{\text{e}}^{-2s_i\tau }}{s_i{\left(s_i\tau \right)}^2}\left[2{\text{e}}^{s_i\tau }-2s_i\tau -{\left(s_i\tau \right)}^2{\text{e}}^{s_i\tau }-2\right]\\ C_i=\frac{{\text{e}}^{-2s_i\tau }}{2s_i{\left(s_i\tau \right)}^2}\left[2{\text{e}}^{s_i\tau }-s_i\tau \left(1+{\text{e}}^{s_i\tau }\right)-2\right]\end{array}$

由上可得，类似于式(6)时间分数阶导数可以离散为以下形式

$F_{\tau }^{\alpha }u\left(x,t_n\right)=\{\begin{array}{l}{\sigma }^{-1}\left(z_{0}u\left(x,t_0\right)+z_{1}u\left(x,t_1\right)+z_{2}u\left(x,t_2\right)\right),n=1\\ {\sigma }^{-1}\left({\widehat{z}}_{0}u\left(x,t_0\right)+{\widehat{z}}_{1}u\left(x,t_1\right)+{\widehat{z}}_{2}u\left(x,t_2\right)\right),n=2\stackrel{}{}\\ {\sigma }^{-1}\left[\frac{\alpha }{2}u\left(x,t_{n-1}\right)-2u\left(x,t_{n-1}\right)+\frac{4-\alpha }{2}u\left(x,t_n\right)\right]\\ +\frac{1}{\Gamma \left(1-\alpha \right)}\left[\frac{u\left(x,t_{n-1}\right)}{{\tau }^{\alpha }}-\frac{u\left(x,t_0\right)}{t_n^{\alpha }}-\alpha {\displaystyle \sum _{i=1}^{N_{\epsilon }}{\omega }_i{U}_{h,i}^n\left(x,t_n\right)}\right],n\ge 3\end{array}$ (24)

其中

$U_{h,i}^1=0,U_{h,i}^n={\text{e}}^{-s_i\tau }{U}_{h,i}^{n-1}+A_{i}u\left(x,t_{n-2}\right)-B_{i}u\left(x,t_{n-1}\right)+C_{i}u\left(x,t_n\right),i=1,2,\cdots ,N_{\epsilon }$

接下来，我们时间分数阶导数离散形式(24)的误差，式(24)与式(6)的区别仅在于$n\ge 3$ 时的情况，针对$n\ge 3$ 时，我们有

$\begin{array}{c}{F}_{\tau }^{\alpha }u\left(x,t_n\right)-{}_0^C{D}_{\tau }^{\alpha }u\left(x,t_n\right)=\frac{\alpha }{\Gamma \left(1-\alpha \right)}{\displaystyle \sum _{k=1}^{k-1}{\displaystyle {\int }_{t_{k-1}}^{t_k}\left[{\left(t_n-s\right)}^{-1-\alpha }-{\displaystyle \sum _{i=1}^{N_{\epsilon }}{\omega }_i}{e}^{-s_i\left(t_n-s\right)}\right]J_{\left[t_{n-1},t_n\right]}\text{d}s}}\\ =\frac{\alpha \left(1-\alpha \right)\left(2-\alpha \right)}{\sigma }{\displaystyle \sum _{k=1}^{n-1}\left[{\tilde{A}}_{n-k}u\left(x,t_{k-1}\right)+{\tilde{B}}_{n-k}u\left(x,t_k\right)+{\tilde{C}}_{n-k}u\left(x,t_{k+1}\right)\right]}\end{array}$

其中

$\begin{array}{l}{\tilde{A}}_{n-k}={\tau }^{\alpha }{\displaystyle {\int }_{t_{k-1}}^{t_k}\left[{\left(t_n-s\right)}^{-1}-{\displaystyle \sum _{i=1}^{N_{\epsilon }}{\omega }_i}{\text{e}}^{-s_i\left(t_n-s\right)}\right]}\frac{\left(s-t_k\right)\left(s-t_{k+1}\right)}{2{\tau }^2}\text{d}s\\ {\tilde{B}}_{n-k}=-\tau {\displaystyle {\int }_{t_{k-1}}^{t_k}\left[{\left(t_n-s\right)}^{-1}-{\displaystyle \sum _{i=1}^{N_{\epsilon }}{\omega }_i}{\text{e}}^{-s_i\left(t_n-s\right)}\right]}\frac{\left(s-t_{k-1}\right)\left(s-t_{k+1}\right)}{{\tau }^2}\text{d}s\\ {\tilde{C}}_{n-k}=\tau {\displaystyle {\int }_{t_{k-1}}^{t_k}\left[{\left(t_n-s\right)}^{-1}-{\displaystyle \sum _{i=1}^{N_{\epsilon }}{\omega }_i}{\text{e}}^{-s_i\left(t_n-s\right)}\right]}\frac{\left(s-t_{k-1}\right)\left(s-t_k\right)}{2{\tau }^2}\text{d}s\end{array}$ (25)

通过计算可知，式(25)中的系数满足以下条件

$\left|{\tilde{A}}_{n-k}\right|\le \frac{5}{12}\epsilon {\tau }^{\alpha +1},\left|{\tilde{B}}_{n-k}\right|\le \frac{2}{3}\epsilon {\tau }^{\alpha +1},\left|{\tilde{C}}_{n-k}\right|\le \frac{1}{12}\epsilon {\tau }^{\alpha +1}$ (26)

定理1 假设$u\left(x,t\right)\in C^3\left[0,T\right]$ ，且令${\overline{r}}_n\left(\tau \right)={}_0^C{D}_t^{\alpha }u\left(x,t_n\right)-F_{\tau }^{\alpha }u\left(x,t_n\right)$，$\forall n\ge 1$ ，则有

$\left|{\overline{r}}_n\left(\tau \right)\right|\le C_u{\tau }^{3-\alpha }+\frac{7\alpha \epsilon t_{n-1}}{6\Gamma \left(1-\alpha \right)}\underset{t\in \left[0,t_n\right]}{\max }\left|u\left(t\right)\right|$ (27)

其中$C_u$ 是一个仅关于函数$u$ 的常数。

证明：当$n=1,2$ 时，${\overline{r}}_n\left(\tau \right)=r_n\left(\tau \right)$ ，由引理1可容易证明式(27)结论成立。

当$n\ge 3$ 时，有

$\begin{array}{c}{}_0^C{D}_t^{\alpha }u\left(x,t_n\right)-F_{\tau }^{\alpha }u\left(x,t_n\right)={}_0^C{D}_t^{\alpha }u\left(x,t_n\right)-{}_0^C{D}_{\tau }^{\alpha }u\left(x,t_n\right)-\left(F_{\tau }^{\alpha }u\left(x,t_n\right)-{}_0^C{D}_{\tau }^{\alpha }u\left(x,t_n\right)\right)\\ =r_n\left(\tau \right)-\frac{\alpha \left(1-\alpha \right)\left(2-\alpha \right)}{\sigma }{\displaystyle \sum _{k=1}^{n-1}\left[{\tilde{A}}_{n-k}u\left(x,t_{k-1}\right)+{\tilde{B}}_{n-k}u\left(x,t_k\right)+{\tilde{C}}_{n-k}u\left(x,t_{k+1}\right)\right]}\end{array}$

那么，根据引理1和式(26)，可得

$\begin{array}{c}\left|{\overline{r}}_n\left(\tau \right)\right|\le C_u{\tau }^{3-\alpha }+\frac{\alpha \left(1-\alpha \right)\left(2-\alpha \right)\epsilon \Delta t^{\alpha +1}}{{\kappa }_{\alpha }}{\displaystyle \sum _{k=1}^{n-1}\left[\frac{5}{12}\left|u\left(x,t_{k-1}\right)\right|+\frac{2}{3}\left|u\left(x,t_k\right)\right|+\frac{1}{12}\left|u\left(x,t_{k+1}\right)\right|\right]}\\ \le C_u{\tau }^{3-\alpha }+\frac{7\alpha \epsilon t_{n-1}}{6\Gamma \left(1-\alpha \right)}\underset{t\in \left[0,t_n\right]}{\max }\left|u\left(t\right)\right|\end{array}$

通过证明可得对于$\forall n\ge 1$ ，使得式(27)成立。定理证毕。

同理，我们将点$\left(x_L,t_n\right)$ 带入方程(1)，并作用高维紧致差分算子${\displaystyle \prod _{m=1}^d{H}_m}$ ，可以得到快速化全离散格式为：

${\displaystyle \prod _{m=1}^d{H}_m}\left({}_0^C{D}_{\tau }^{\alpha }{u}_L^n\right)-{\displaystyle \sum _{i=1}^d{\kappa }_i}{\displaystyle \prod _{m=1,m\ne i}^d{H}_m}{\delta }_{x_i}^2{u}_L^n={\displaystyle \prod _{m=1}^d{H}_m}{f}_L^n+{\overline{R}}_L^n$

这里${\overline{R}}_L^n$ 满足以下表达式：

$\begin{array}{c}{\overline{R}}_L^n={\displaystyle \sum _{i=1}^d{\kappa }_i}{\displaystyle \prod _{m=1,m\ne i}^d{H}_m}{\displaystyle {\int }_0^1\left[{\partial }_{x_i}^{6}u\left(x_L-{\widehat{s}}_i,t_n\right)+{\partial }_{x_i}^{6}u\left(x_L+{\widehat{s}}_i,t_n\right)\right]\times {\left(1-s_i\right)}^3\left[5-3{\left(1-s_i\right)}^2\right]\text{d}{s}_i}\\ -{\displaystyle \prod _{m=1}^d{H}_m}\left({\overline{r}}_n\left(x_L,t_n\right)\right)\end{array}$

其中${\widehat{s}}_i=s_i\Delta x_i{e}_i$ ，$e_i$ 为第$i$ 个空间方向的标准正交基，并且$R_L^n$ 满足以下条件

$\left|{\overline{R}}_L^n\right|\le {\overline{C}}_R\left({\tau }^{3-\alpha }+{\displaystyle \sum _{i=1}^d\Delta }{x}_i^4\right)$ (28)

这里的${\overline{C}}_R$ 是一个不依赖于$\tau$ 和$\Delta x_i\left(i=1,2,\cdots ,d\right)$ 的常数。

通过以上分析，当我们的时空步长取得足够小时，即$\tau \to 0$ ，$\Delta x_i\to 0$ 时，有$R_L^n\to 0$ 。再结合式(13)的简化，我们得到方程(1)的一个快速化的L2格式如下：

$\{\begin{array}{l}{\displaystyle \prod _{m=1}^d{H}_m}\left(z_0{u}_L^0+z_1{u}_L^1+z_2{u}_L^2\right)-\sigma {\Lambda }_{\Delta x}{u}_L^1=\sigma {\displaystyle \prod _{m=1}^d{H}_m{f}_L^1},n=1\\ {\displaystyle \prod _{m=1}^d{H}_m}\left({\widehat{z}}_0{u}_L^0+{\widehat{z}}_1{u}_L^1+{\widehat{z}}_2{u}_L^2\right)-\sigma {\Lambda }_{\Delta x}{u}_L^2=\sigma {\displaystyle \prod _{m=1}^d{H}_m{f}_L^2},n=2\\ {\displaystyle \prod _{m=1}^d{H}_m}\left[\frac{\alpha }{2}u\left(x,t_{n-2}\right)-2u\left(x,t_{n-1}\right)+\frac{4-\alpha }{2}u\left(x,t_n\right)\right]\\ +\frac{1}{\Gamma \left(1-\alpha \right)}\left[\frac{u\left(x,t_{n-1}\right)}{{\tau }^{\alpha }}-\frac{u\left(x,t_0\right)}{t_n^{\alpha }}-\alpha {\displaystyle \sum _{i=1}^{N_{\epsilon }}{\omega }_i{U}_{h,i}^n\left(x,t_n\right)}\right]-\sigma {\Lambda }_{\Delta x}{u}_L^n=\sigma {\displaystyle \prod _{m=1}^d{H}_m{f}_L^n},n\ge 3\end{array}$ (29)

5. 数值算例

本节将通过下面的数值算例来比较所提出的快速化L2格式与改进的L2格式在收敛速度和CPU计算时间上的差别，故所有计算均在电脑CPU为Intel(R) Core (TM) i5-9400F CPU @ 2.90 GHz，使用的编译程序为MATLAB。为了避免由MATLAB精度带来的机器误差，我们在计算积分系数时统一使用MATLAB内置函数quadgk进行计算。

为了验证所提出格式的有效性和准确性，我们定义误差$\overline{e}\left(h,\tau \right),h=\Delta x_i,i=1,2,\cdots ,d$ ，其表达式如下：

$\overline{e}\left(h,\tau \right)=\underset{1\le k\le N}{\max }\Vert u\left(x,t_k\right)-u^k\Vert$

当$h$ 足够小时，时间收敛阶$Rate\left(\tau \right)$ 表示如下：

$Rate\left(\tau \right)={\log }_2\left(\frac{\overline{e}\left(h,2\tau \right)}{\overline{e}\left(h,\tau \right)}\right)$

类似地，空间收敛阶$Rate\left(h\right)$ 表示如下：

$Rate\left(h\right)={\log }_2\left(\frac{\overline{e}\left(2h,\tau \right)}{\overline{e}\left(h,\tau \right)}\right)$

例1 假设方程(1)具有光滑的解析解$u\left(x,t\right)=t^4\sin \left(2\pi x\right)$ ，且扩散系数$\kappa =1$ ，那么我们可以得到方程(1)中的$f\left(x,t\right)$ 和$u_0\left(x\right)$ 为

$f\left(x,t\right)=\left(\frac{24}{\Gamma \left(5-\alpha \right)}{t}^{4-\alpha }+4{\pi }^2{t}^4\right)\sin \left(2\pi x\right),u_0\left(x\right)=0$

我们选择空间区域$\Omega =\left[0,1\right]$ ，$M=2^{10},N=2^6:2^9$ 去测试该算例的时间收敛阶，即空间步长$h=\Delta x=2^{-10}$ ，时间步长$\tau =2^{-6}:2^{-9}$ ，对于时间快速化格式，我们选取SOE近似误差为$\epsilon ={10}^{-7}$ 。下表1分别展示了在不同分数阶导数$\alpha =0.3,0.5,0.7$ 的时间收敛阶$Rate\left(\tau \right)$ ，以及改进的L2格式和快速化L2格式的CPU计算时间对比。

Table 1. Comparison of time convergence order calculations between the improved L2 scheme and the fastened L2 scheme

表1. 改进的L2格式和快速化L2格式的时间收敛阶计算结果对比

	$\tau$	改进L2误差	收敛阶	CPU时间	快速化L2 误差	收敛阶	CPU时间
$\alpha =0.3$	1/64	5.8698e−7	-	3.02 s	4.2146e−8	-	2.22 s
	1/128	9.3723e−8	2.6468	11.08 s	6.1203e−9	2.7837	6.14 s
	1/256	1.4838e−8	2.6590	42.46 s	8.4560e−10	2.7837	9.51 s
	1/512	2.3398e−9	2.6648	157.98 s	9.1019e−11	2.8555	20.16 s
$\alpha =0.5$	1/64	2.3328e−6	-	3.02 s	1.6940e−6	-	2.30 s
	1/128	4.1980e−7	2.4742	11.47 s	2.9693e−7	2.5122	5.42 s
	1/256	7.5068e−8	2.4834	44.24 s	5.2010e−8	2.5132	9.77 s
	1/512	1.3377e−8	2.4883	156.55 s	9.0747e−9	2.5188	19.75 s
$\alpha =0.7$	1/64	8.1421e−6	-	3.58 s	7.1767e−6	-	2.18 s
	1/128	1.6696e−6	2.2858	11.88 s	1.4613e−6	2.2961	4.47 s
	1/256	3.4088e−7	2.2921	41.74 s	2.9688e−7	2.2992	9.21 s
	1/512	6.9437e−8	2.2954	151.76 s	6.0232e−8	2.3013	19.16 s

由上表，我们可以观察到，两个数值格式的时间收敛阶都接近为$3-\alpha$ ，这与我们的误差估计理论相符，并且计算得出的误差都差不多，但是随着$N$ 的增大，时间上的计算步骤增加时，明显地可以观察到快速化L2格式的CPU计算时间要小得多。下图1展示了两种格式下的计算时间对比。

Figure 1. Comparison of CPU computation time for two numerical schemes

图1. 两种数值格式的CPU计算时间对比

同时，我们选取$\alpha =0.5$ ，$N=M=2^7$ ，在图2中表示了快速化L2格式中选定的误差$\epsilon$ 从10⁻⁵到10⁻¹⁰时，所对应的$N_{\epsilon }$ 值变化情况，可以发现我们选定的误差越小，所需要计算的计算量是越大的，但相对于整体的计算时间而言是可以忽略不计的。

Figure 2. Effect of different error ε on CPU computation time

图2. 不同误差ε对CPU计算时间的影响

现在，我们选择$M=2^4:2^7,N={\left[M\right]}^{4/\left(3-\alpha \right)}$ 来观察空间方向上的收敛阶，对于时间快速化格式，我们选取SOE近似误差为$\epsilon ={10}^{-7}$ 。下表2分别展示了在不同分数阶导数$\alpha =0.3,0.5,0.7$ 的空间收敛阶$Rate\left(h\right)$ ，以及改进的L2格式和快速化L2格式的CPU计算时间对比。

Table 2. Comparison of spatial convergence order calculations between the improved L2 scheme and the fastened L2 scheme

表2. 改进的L2格式和快速化L2格式的空间收敛阶计算结果对比

	$h$	改进L2误差	收敛阶	CPU时间	快速化L2 误差	收敛阶	CPU时间
$\alpha =0.3$	1/16	9.6583e−5	-	1.66 s	9.6248e−5	-	1.39 s
	1/32	6.0114e−6	4.0060	12.09 s	5.9865e−6	4.0070	4.17 s
	1/64	3.7539e−7	4.0012	92.19 s	3.7360e−7	4.0022	12.79 s
	1/128	2.3459e−8	4.0001	746.59 s	2.3281e−8	4.0043	38.78 s
$\alpha =0.5$	1/16	9.5924e−5	-	2.93 s	9.5589e−5	-	1.95 s
	1/32	5.9687e−6	4.0064	25.95 s	5.9456e−6	4.0069	6.57 s
	1/64	3.7267e−7	4.0014	241.20 s	3.7107e−7	4.0021	21.99 s
	1/128	2.3286e−8	4.0003	2353.45 s	2.3111e−8	4.0050	73.50 s
$\alpha =0.7$	1/16	9.5043e−5	-	7.49 s	9.4819e−6	-	3.36 s
	1/32	5.9137e−6	4.0064	79.00 s	5.8989e−6	4.0067	11.39 s
	1/64	3.6923e−7	4.0015	799.00 s	3.6821e−7	4.0019	41.52 s
	1/128	2.3071e−8	4.0003	9322.15 s	2.2944e−8	4.0043	149.30 s

由上表，我们可以观察到，两个数值格式的空间收敛阶都接近为4，这与我们的误差估计理论相符合，充分验证了我们所提出数值格式的有效性和准确性。

例2 假设方程(1)的解析解不光滑时，即解析解在初始时刻具有奇异性，当$u\left(x,t\right)=t^{\alpha }\sin \left(x\right)$ ，扩散系数$\kappa =1$ ，那么我们可以得到方程(1)中的$f\left(x,t\right)$ 和$u_0\left(x\right)$ 为

$f\left(x,t\right)=\left(\Gamma \left(1+\alpha \right)+t^{\alpha }\right)\sin \left(x\right),u_0\left(x\right)=0.$

我们选择空间区域$\Omega =\left[0,2\pi \right]$ ，$M=2^{10},N=2^6:2^9$ 去测试该算例的时间收敛阶，对于时间快速化格式，我们选取SOE近似误差为$\epsilon ={10}^{-7}$ 。下表3分别展示了在不同分数阶导数$\alpha =0.3,0.5,0.7$ 的时间收敛阶$Rate\left(\tau \right)$ ，以及改进的L2格式和快速化L2格式的CPU计算时间对比。

Table 3. Comparison of time convergence order computational results for solving problems with initial value singularities in the improved L2 scheme and the fastened L2 scheme

表3. 改进的L2格式和快速化L2格式求解具有初值奇异性问题的时间收敛阶计算结果对比

	$\tau$	改进L2误差	收敛阶	CPU时间	快速化L2 误差	收敛阶	CPU时间
$\alpha =0.3$	1/64	3.0470e−2	-	3.06 s	3.0469e−2	-	2.22 s
	1/128	2.5865e−2	0.2363	10.19 s	2.5865e−2	0.2364	4.65 s
	1/256	2.1805e−2	0.2463	39.38 s	2.1806e−2	0.2463	10.01 s
	1/512	1.8273e−2	0.2549	156.43 s	1.8273e−2	0.2550	20.70 s

续表

$\alpha =0.5$	1/64	1.7363e−2	-	5.49 s	1.7363e−2	-	2.95 s
	1/128	1.2712e−2	0.4498	13.38 s	1.2712e−2	0.4498	6.09 s
	1/256	9.2184e−3	0.4637	52.45 s	9.2185e−3	0.4636	10.83 s
	1/512	6.6377e−3	0.4783	199.93 s	6.6377e−3	0.4784	25.37 s
$\alpha =0.7$	1/64	6.7364e−3	-	3.76 s	6.7368e−3	-	2.13 s
	1/128	4.2359e−3	0.6932	11.84 s	4.2361e−3	0.6932	4.51 s
	1/256	2.6422e−3	0.6809	40.47 s	2.6423e−3	0.6809	9.09 s
	1/512	1.6398e−3	0.6881	144.24 s	1.6399e−3	0.6882	18.62 s

由表3我们可以看出，对于具有初值奇异性的解，所计算出的结果并没有达到我们想要的结果，时间收敛阶仅为$\alpha$ 阶，所以我们通过使用梯度网格，加密时间在初始时刻的分布，去解决这类初值奇异性问题。梯度网格将时间剖分重新定义为$t_n={\left(\frac{n}{N}\right)}^r$ ，其中$r$ 为梯度网格参数，当梯度参数取$r=\left(3-\alpha \right)/\alpha$ 时，时间收敛阶能提高到$3-\alpha$ 阶。在文[15]中提出了L2格式下的梯度网格定义以及误差分析，通过数值算例，充分验证了我们提出的改进L2格式和快速化L2格式，可以同样可以使用梯度网格去解决具有初值奇异性的问题。表4验证了我们的结果。

Table 4. Comparison of time-convergence order computational results for solving problems with initial value singularities on gradient grids in improved L2 scheme and fastened L2 scheme

表4. 改进的L2格式和快速化L2格式梯度网格求解具有初值奇异性问题的时间收敛阶计算结果对比

	$\tau$	改进L2误差	收敛阶	CPU时间	快速化L2 误差	收敛阶	CPU时间
$\alpha =0.3$	1/64	9.8170e−3	-	2.84 s	9.7474e−3	-	6.40 s
	1/128	1.5465e−3	2.6662	9.59 s	1.5356e−3	2.6661	14.61 s
	1/256	2.3886e−4	2.6947	36.15 s	2.3718e−4	2.6947	32.34 s
	1/512	3.6780e−5	2.6991	139.90 s	3.6521e−5	2.6991	71.56 s
$\alpha =0.5$	1/64	5.4803e−6	-	2.76 s	5.3762e−6	-	4.31 s
	1/128	9.9887e−7	2.4558	9.51 s	9.8010e−7	2.4555	9.41 s
	1/256	1.7754e−8	2.4921	35.74 s	1.7421e−8	2.4920	20.74 s
	1/512	3.1416e−8	2.4986	139.42 s	3.0826e−9	2.4985	46.98 s
$\alpha =0.7$	1/64	2.7165e−4	-	3.74 s	1.0122e−4	-	3.46 s
	1/128	5.5344e−5	2.2952	12.15 s	2.1427e−5	2.2400	7.57 s
	1/256	1.1245e−5	2.2990	41.83 s	4.4315e−6	2.2735	16.45 s
	1/512	2.2839e−6	2.2998	157.09 s	9.0661e−7	2.2892	35.44 s

通过观察表4，我们发现快速化格式的CPU时间会随着$\alpha$ 的大小变化而产生较大波动，并且会出现在步长较小时，比普通的L2格式计算速度更慢。这是因为我们采取梯度网格去处理初值奇异性问题时，我们的SOE近似方法在计算$N_{\epsilon }$ 的过程中是依据最小步长${\tau }_{\min }$ 决定的，由梯度网格的定义可知，当$\alpha$ 越小时，梯度参数$r$ 越大，导致在初始时刻的$t_n$ 越小，这使得其最小步长${\tau }_{\min }$ 也越小，从而有$N_{\epsilon }$ 越大。下图3能清晰表述我们结论。

Figure 3. $N_{\epsilon }$ for different $\alpha$ under the graded meshes with $\epsilon$

图3. 梯度网格下不同$\alpha$ 的$N_{\epsilon }$ 随着$\epsilon$ 的变化情况

通过观察图3，我们清楚的发现，相对于在均匀网格的情况下，梯度网格的$N_{\epsilon }$ 的值会大得许多，这会导致在$N$ 和$\alpha$ 较小时，没有显示出这么明显的加速效果。总之，该算例验证了我们所提出的数值格式可以结合梯度网格技巧，去解决一些具有初值奇异性的问题。

例3 假设方程(1)具有解析解为$u\left(x,t\right)=t^{3+\alpha }\sin \left(x_1\right)\sin \left(x_2\right)$ ，且扩散系数${\kappa }_1={\kappa }_2=1$ ，那么我们可以得到方程(1)中的$f\left(x,t\right)$ 和$u_0\left(x\right)$ 为

$f\left(x,t\right)=\left(\frac{\Gamma \left(4+\alpha \right)}{6}{t}^3+2t^{3+\alpha }\right)\sin \left(x_1\right)\sin \left(x_2\right),u_0\left(x\right)=0.$

我们选择空间区域${\left[0,2\pi \right]}^2$ ，$M=N=2^5:2^8$ 去测试该算例的时间收敛阶，对于时间快速化格式，我们选取SOE近似误差为$\epsilon ={10}^{-7}$ 。下表5分别展示了在不同分数阶导数$\alpha =0.3,0.5,0.7$ 的时间收敛阶$Rate\left(\tau \right)$ ，以及改进的L2格式和快速化L2格式的CPU计算时间对比。

Table 5. Comparison of time convergence order calculations for two numerical schemes with \alpha = 0.3, 0.5, 0.7

表5. 在\alpha=0.3, 0.5, 0.7时两种数值格式的时间收敛阶计算结果比较

	$\tau$	改进L2误差	收敛阶	CPU时间	快速化L2 误差	收敛阶	CPU时间
$\alpha =0.3$	1/32	2.5050e−5	-	0.64 s	1.4449e−5	-	0.75 s
	1/64	3.6683e−6	2.7716	2.31 s	1.8040e−6	3.0016	2.06 s
	1/128	5.6779e−7	2.6917	16.52 s	2.4388e−7	2.8870	7.98 s
	1/256	8.6113e−8	2.7210	96.86 s	3.3896e−8	2.8469	49.85 s
$\alpha =0.5$	1/32	9.9805e−5	-	0.48 s	7.4814e−5	-	0.76 s
	1/64	1.7699e−5	2.4953	2.24 s	1.2771e−5	2.5504	2.03 s
	1/128	3.1530e−6	2.4889	12.68 s	2.2153e−6	2.5273	7.54 s
	1/256	5.6136e−7	2.4897	90.67 s	3.8648e−7	2.5190	49.04 s
$\alpha =0.7$	1/32	3.6625e−3	-	0.83 s	3.2640e−4	-	0.75 s
	1/64	7.5372e−5	2.2807	3.02 s	6.6443e−5	2.2964	1.97 s
	1/128	1.5432e−6	2.2880	15.42 s	1.3506e−5	2.2984	7.57 s
	1/256	3.1491e−7	2.2929	99.57 s	2.7426e−6	2.3000	48.87 s

同样，我们选择$M=2^3:2^6,N={\left[M\right]}^{4/\left(3-\alpha \right)}$ 来观察空间方向上的收敛阶，对于时间快速化格式，我们选取SOE近似误差为$\epsilon ={10}^{-7}$ 。下表6分别展示了在不同分数阶导数$\alpha =0.3,0.5,0.7$ 的空间收敛阶$Rate\left(h\right)$ ，以及改进的L2格式和快速化L2格式的CPU计算时间对比。

Table 6. Comparison of spatial convergence order calculations for two numerical schemes with \alpha = 0.3, 0.5, 0.7

表6. 在\alpha=0.3,0.5,0.7时两种数值格式的空间收敛阶计算结果比较

	$h$	改进L2误差	收敛阶	CPU时间	快速化L2 误差	收敛阶	CPU时间
$\alpha =0.3$	1/8	9.8786e−4	-	0.27 s	9.6122e−4	-	0.46 s
	1/16	6.1015e−5	4.0170	1.64 s	5.8916e−5	4.0281	1.37 s
	1/32	3.8088e−6	4.0017	12.37 s	3.6572e−6	4.0098	4.40 s
	1/64	2.3832e−7	3.9983	100.29 s	2.2716e−7	4.0089	16.03 s
$\alpha =0.5$	1/8	9.4555e−4	-	0.33 s	9.1152e−4	-	0.58 s
	1/16	5.8728e−5	4.0090	2.92 s	5.6148e−5	4.0209	2.00 s
	1/32	3.6588e−6	4.0046	27.20 s	3.4839e−6	4.0104	6.96 s
	1/64	2.2878e−7	3.9993	237.15 s	2.1649e−7	4.0082	29.10 s
$\alpha =0.7$	1/8	9.4353e−4	-	0.90 s	9.1421e−4	-	0.83 s
	1/16	5.8425e−5	4.0134	8.33 s	5.6357e−5	4.0198	3.26 s
	1/32	3.6418e−6	4.0038	83.81 s	3.5052e−6	4.0070	12.21 s
	1/64	2.2777e−7	3.9990	895.46 s	2.1838e−7	4.0045	54.19 s

通过观察，我们发现所提出的数值格式对于二维算例也具有很好的数值结果，这进一步证明了这两种数值格式的有效性和准确性。

6. 讨论与展望

首先，本文所提出的改进L2格式和L2格式的快速化算法能有效地解决许多时间分数阶微分方程问题，具有较强的实用性，并且得到的数值效果较好，能够达到高阶精度收敛。但我们的解析解需要在空间上达到六阶连续，在解决实际问题中会有一定的局限性。

其次，在针对具有奇异性的初值问题时，我们仅在均匀网格上做了误差估计，可在数值算例中，我们可以发现这两种格式均可以结合梯度网格去处理这类问题，并且能达到与均匀网格理论的相同精度，这说明我们所提出格式的准确性和有效性。

最后，我们接下来的工作将沿着以下几个方向进行扩展，从而提高所提出格式的高效性。我们将以更严谨的理论证明在梯度网格下的误差估计以及稳定性、收敛性，运用更好的办法去解决MATLAB带来的机器误差去提高计算的效率。

7. 总结

本文提出了两种L2数值格式去解决时间分数阶扩散方程，在空间上结合紧致差分格式，使得这两种格式在时间上的收敛阶达到$3-\alpha$ 阶，在空间上的收敛阶达到4阶，并给出误差估计的证明。同时，通过数值算例对比了两种格式的差异性，体现了SOE快速化算法的有效性，在保证收敛阶稳定的同时，能够在很大程度上降低计算成本，并且验证了该数值格式可以结合梯度网格技巧去解决具有初值奇异性的问题。同理推出这两种格式在任意维度上都具有很好的计算效果。

基金项目

国家自然科学基金(12361083，12461077和11961009)，贵州省高等学校大数据分析与智能计算重点实验室(黔教技[2023] 012号)，贵州省科技厅基金(QHKJC-ZK [2024] YB497)，贵州民族大学科研基金资助项目(GZMUZK [2023] CXTD05)，贵州省研究生科研基金项目(2024YJSKYJJ218)。

NOTES

^*通讯作者。