时空分数扩散方程的非均匀L1/谱元方法：数值分析及实现

《Applied Numerical Mathematics》：Nonuniform L1/spectral element scheme for the time-space fractional diffusion equation: Numerical analysis and implementation

　　闵采|李长平|童薇薇
上海大学数学系，中国上海200444

**摘要**
本文提出了一种高效算法，用于求解具有时间方向Caputo导数和空间方向Riesz导数的时空分数扩散方程。为解决初始时刻的奇异性问题，在分级网格上采用L1近似来处理Caputo导数。同时，利用谱元法进行空间离散化以提高计算精度。本文严格证明了所提算法的数值稳定性和误差估计，并详细展示了相应的数值实现方法。通过数值算例验证了该算法的可行性和高效性。

**引言**
分数微分方程是描述空间中的长程相互作用以及时间上的历史依赖性的有效工具。这类方程被广泛用于物理、化学、力学和工程学等多个领域中材料与现象的建模。由于其在实际应用中的重要性，求解分数微分方程已成为当前研究热点之一。

分数微分算子的非局部性和可积奇异性使得求解分数微分方程面临挑战，尽管有些方程的解可以与特殊函数有紧密联系[1]。近几十年来，许多学者致力于设计高效的分数微分方程求解算法，并开展相应的数值分析研究。例如基于有限差分法[2][3]、有限元法[4][5][6]以及谱方法[7][8]的算法。与现有的针对时间或空间分数微分方程的算法相比，针对时空分数微分方程的数值方案仍存在较大不足。因此，开发高效的时空分数微分方程数值方法仍有广阔的研究空间。

设计数值方案的关键在于离散化技术的选择。众所周知，微分方程中的分数导数可能导致初始时刻或边界处的弱奇异性。在这种情况下，非均匀网格常被用作离散化的基础结构。实现这一目标的一种常见方法是通过基于非均匀网格的有限差分公式来近似分数导数，相关文献可见[9][10][11][12]及其中引用的内容。实际上，在使用有限元法[13]、多域谱方法[14]或谱元法[15][16]对分数微分方程进行离散化时，也可以应用非均匀网格。谱元法是一种高阶数值技术，主要用于求解偏微分方程，尤其在流体力学、波动传播和电磁学等领域。它结合了有限元法的几何灵活性与谱方法的高精度和指数级收敛性。尽管已有部分算法在非均匀网格上使用有限差分法或谱元法求解时间或空间分数微分方程，但目前尚无算法同时结合这两种技术来数值求解时空分数微分方程，更没有详细的数值分析及逐步实现指导。

本文研究以下时空分数扩散方程的数值解法：
$$
\{CD_0^{t_\gamma}u(x,t) - RZD_{x^\nu}u(x,t) = f(x,t),\ (x,t) \in \Omega \times (0,T],\ u(x,0) = u_0(x),\ x \in \Omega,\ u(a,t) = u(b,t) = 0,\ t \in (0,T]
$$
其中 $\Omega = (a,b)$，$0 < \gamma < 1 < \nu < 2$，$f(x,t)$ 和 $u_0(x)$ 为已知函数。这里
$$
CD_0^{t_\gamma}u(x,t) = \frac{1}{\Gamma(1-\gamma)}\int_0^t (t-s)^{-\gamma}\frac{\partial^2 u(x,s)}{\partial s^2}\,ds,\ t>0,\ 0<\gamma<1,
$$
且通常假定对于任意 $x \in [a,b]$，$u(x,t) \in AC[0,T]$ 关于 $t$ 连续，这只是保证上述积分存在的充分条件。符号 $RZD_{x^\nu}$ 表示目前常用的Riesz导数[17]：
$$
RZD_{x^\nu}u(x,t) = -\frac{1}{2}\cos\left(\frac{\pi\nu}{2}\right)\left(RLD_{a,x^\nu}u(x,t) + RLD_{x,b}^\nu u(x,t)\right),\ x \in (a,b),\ 1<\nu<2,
$$
其中
$$
RLD_{a,x^\nu}u(x,t) = \frac{1}{\Gamma(2-\nu)}\frac{\partial^2 u(x,t)}{\partial x^2}\int_a^x (x-\zeta)^{1-\nu}u(\zeta,t)\,d\zeta,
$$
$$
RLD_{x,b}^\nu u(x,t) = (-1)^{2-\nu}\frac{1}{\Gamma(2-\nu)}\frac{\partial^2 u(x,t)}{\partial x^2}\int_x^b (\zeta-x)^{1-\nu}u(\zeta,t)\,d\zeta,
$$
分别为左右Riemann-Liouville导数。同样，为保证上述两个积分的存在，还假定对于任意 $t \in [0,T]$，$u(x,t) \in AC^2[a,b]$ 关于 $x$ 二次连续。需要说明的是，此处所用的Riesz导数定义并非最原始的定义，详细解释可见[18]，但目前人们已习惯使用该定义。综合考虑解的弱奇异性与计算效率，本文选择非均匀L1公式来近似Caputo导数，并采用谱元法进行空间离散化。

本文结构如下：第2节介绍分数微积分的基础知识、与分数导数相关的函数空间以及正交多项式；第3节基于分级网格上的L1公式和谱元法，推导并分析方程(1.1)的完全离散化方案；第4节展示所提算法的数值实现方法；第5节通过数值实验验证理论分析结果；第6节总结全文主要结论。

**章节节选**

**基础知识**
首先介绍分数微积分所需的定义、基本性质以及与分数导数相关的函数空间。
**定义1[19]**：设 $\phi(x) \in L^1(\Omega)$，其中 $\Omega = (a,b)$。左右Riemann-Liouville积分定义为：
$$
RLD_{a,x}^{-\mu}\phi(x) = \frac{1}{\Gamma(\mu)}\int_a^x (x-\zeta)^{-\mu}\phi(\zeta)\,d\zeta,\ \mu>0,
$$
$$
RLD_{x,b}^{-\mu}\phi(x) = \frac{1}{\Gamma(\mu)}\int_x^b (\zeta-x)^{-\mu}\phi(\zeta)\,d\zeta,\ \mu>0.
$$
特别地，分数积分具有半群性质：
$$
RLD_{a,x}^{-\mu}RLD_{a,x}^{-\nu}\phi(x) = RLD_{a,x}^{-(\mu+\nu)}\phi(x),\quad RLD_{x,b}^{-\mu}RLD_{x,b}^{-\nu}\phi(x) = RLD_{x,b}^{-(\mu+\nu)}\phi(x),
$$
其中 $\mu...$

**方程(1.1)的非均匀L1/谱元方案**
本节首先介绍非均匀L1公式和谱元法的基本原理，随后推导方程(1.1)的非均匀L1/谱元方案，最后证明该方案的无条件稳定性与收敛性。

**数值实现**
下面详细展示方程(3.8)的数值实现方法。对于 $i=1,2,\dots,M$，模态基函数定义为：
$$
\phi_j^i(x) = \begin{cases}
L_j(2x - x_i - 1 - \xi_{hi}), & x \in I_i,\\
0, & x \in \Omega \setminus I_i,
\end{cases}
$$
其中 $L_j(x)$ 表示 $x \in I=[-1,1]$ 上的第 $j$ 阶勒让德多项式。节点基函数定义为：
$$
\phi^i(x) = \begin{cases}
\frac{1}{2} + 2x - x_i - 1 - \xi_{hi}, & x \in I_i,\\
\frac{1}{2} - 2x - x_i + 1 - \xi_{hi} + 1, & x \in I_{i+1},\\
0, & \text{其他情况},
\end{cases}
$$
其中 $i=1,2,\dots,M-1$，且 $\phi^0(x) = \phi^M(x) = 0$。因此，方程(3.3)中的谱元空间可表示为：
$$
V_{N_0,\nu}^2(\Omega_h) = \text{span}\{\phi^0(1),\dots,\phi^{N_1-2}(1),\dots,\phi^0(M)\}.
$$

**数值算例**
本节通过数值算例展示所提方案(3.8)的适用性，包括存在精确解的情况和不存在精确解的情况。所有数值实验均在Core i5处理器上使用Matlab 2022b完成。此处数值误差指最终时刻的L2范数误差，即 $\|u(x,T) - u_{hn}^T\|_2$。为简化计算，所有子区间均设置 $N_i = N$，几何网格则取 $q=0.5$。
**算例1**：设 $\gamma \in (0,1)$，$\nu \in (1,2)$，考虑以下问题...

**结论**
本文提出了一种适用于具有时间方向Caputo导数和空间方向Riesz导数的时空分数扩散方程的可靠数值方法。通过非均匀L1公式处理解在初始时刻的弱奇异性，同时采用可与均匀网格和几何网格结合使用的谱元法进行空间离散化，从而提升计算效率。本文还给出了该算法的详细数值实现方法。

**CRediT作者贡献声明**
闵采：写作——审阅与编辑，写作——初稿撰写，验证，监督，方法设计。
李长平：写作——审阅与编辑，验证，监督，概念构思。
童薇薇：写作——审阅与编辑，写作——初稿撰写，可视化。

**利益冲突声明**
作者声明不存在任何利益冲突。