变阶分数阶阻尼方程的高阶时空离散及高效求解算法

Higher-Order Spatio-Temporal Discretization and Efficient Solution Algorithms for Variable-Order Fractional-Order Damping Equations

摘要: 针对时间分数阶变阻尼方程的数值求解问题，本研究构建了高效离散格式并系统分析了其理论特性。时间方向采用L2_1_σ格式进行离散，空间方向应用二阶中心差分方法，建立了完整的数值离散格式，并证明了格式的收敛性与稳定性。为提升计算效率，重点研究了预处理共轭梯度法与多重网格法两类快速算法，详细对比了二者的计算性能。数值实验不仅验证了理论分析，更通过三种求解方法的耗时对比揭示：预处理共轭梯度法因其优化的算法复杂度和高效的预处理技术，在计算效率方面展现出显著优势。

Abstract: For the numerical solution of the time-fractional order variable damping equation, this study constructs an efficient discretization format and systematically analyzes its theoretical properties. The time direction is discretized by the L2_1_σ formula, and the spatial direction is discretized by the second-order center difference method, which establishes a complete numerical discretization format and proves the convergence and stability of the format. In order to improve the computational efficiency, two types of fast algorithms, namely the preprocessing conjugate gradient method and the multiple grid method, are emphasized, and the computational performances of the two are compared in detail. The numerical experiments not only verify the theoretical analysis, but also through the time-consuming comparison of the three solution methods show that the preprocessing conjugate gradient method has a significant advantage in computational efficiency due to its optimized algorithmic complexity and efficient preprocessing technology.

1. 引言

最近几年，分数阶微分方程的研究与发展日益受到各界关注，目前已广泛渗透到工程学、生物学、化学、电化学以及扩散理论等[1]-[4]多个学科领域之中。随着研究逐步深入，众多研究人员发现，复杂物理系统的记忆性与非局域性，或许会因时间、空间或其他因素而改变[5]。因而，急需更有效的工具来描述和处理更为复杂的模型问题。近年来的诸多实例表明，借助变分数阶算子，研究者能够解决更复杂的实际问题，并取得了一系列有价值的成果[6] [7]。由此，变分数阶微分方程成为了描述复杂动态现象更优的模型选择。

由于变分数阶方程的高度复杂性，使得变阶时间分数阶微分方程解析解求解困难，有必要建立有效的数值求解方法。变分数阶微分方程的数值解研究包括时间和空间，本文重点关注时间变分数阶。Shen等人[8]基于L1公式推导了变阶时间分数阶扩散方程中Caputo分数阶导数的有限差分方法，并通过傅里叶分析验证了该方法的稳定性、收敛性和可解性。Ahmadinia等人[9]推出具有次扩散和超扩散的时变阶分数阶微分方程的近似解。在时空上分别使用有限差分方法、局部不连续伽辽金方法进行近似。Du等人[10]利用每个时间间隔上的特殊点，提出了两种在时间方向上具有二阶精度的差分格式，以求解多维变阶时间分数阶扩散方程，并运用能量方法分析了其收敛性和稳定性。Jia等[11]提出了一种变阶时间分数扩散方程的预条件快速有限元逼近方法。Du等[12]推广了Alikhanov的L2_1_σ公式，并将约简方法应用于变阶分数阶微分算子，得到了变阶时间分数阶波动方程的时间二阶有限差分格式。本文研究了具有阻尼项的非线性变分数阶扩散方程的二阶隐式有限差分方法，并讨论了该方法的稳定性和收敛性。更具体地说，我们考虑下列方程

${}_0^C{D}_t^{\alpha \left(t\right)}u+\frac{\partial u}{\partial t}=\frac{{\partial }^{2}u}{\partial x^2}+f\left(x,t\right),0<x<L,0<t<T$ (1)

初始和边界条件为：

$u\left(x,0\right)=\phi \left(x,0\right),0\le x\le L$ (2)

$u\left(0,t\right)=\varphi \left(t\right),u\left(L,t\right)=\psi \left(t\right),0<t\le T.$ (3)

$0<\alpha \le \overline{\alpha }<1$ ，在Caputo意义下引入时间分数阶导数，即

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

在科学技术日新月异、迅猛发展的当下，在科学与工程计算所涵盖的众多领域里，求解如下大型线性方程组$AX=b$ ，其中$A\in n\times n$ 为正定矩阵，$b\in R^n$ 成为了普遍且关键的需求。多重网格法(MG)作为求解偏微分方程离散系统的高效方法，其核心思想在于通过粗–细网格间的误差校正实现最优计算复杂度O(N)，其中N为未知量数量。在具体实现中，Deuflhard等[13]提出的V循环框架通过在粗网格层采用直接法(如反斜杠)精确求解残差方程，有效避免了迭代误差累积问题。针对误差平滑过程，Brandt [14]指出Gauss-Seidel方法对高频误差的高效消除特性，本文据此设计光滑子，并通过参数动态调整迭代次数以平衡计算效率与局部误差抑制效果。投影算子的构造则遵循Briggs等[15]提出的Galerkin准则，确保粗细网格间传递算子的稀疏性，显著降低大规模问题的内存与计算成本。针对传统MG方法在非结构化网格或复杂插值策略中的编程复杂性，传统MG方法(如V/W循环)的编程复杂性常受限于非结构化网格或复杂插值策略，而本文通过结构化网格与固定模板设计(如有限差分离散下的均匀粗化)简化实现流程，同时继承多重网格的线性复杂度优势(Trottenberg等[16]对结构化MG的高效性分析)。和其他相关的著作可以在[17] [18]中找到。

为了提高计算效率，相关研究持续深入。Wang等[19]发现离散系统呈现Toeplitz-like结构，借助快速傅里叶变换(FFT)，可以在O(NlogN)的时间内完成Toeplitz矩阵的矩阵–向量乘法。因此，结合FFT技术和Toeplitz类结构，Wang等[20]进一步提出，Krylov子空间方法每次迭代的计算成本为O(NlogN)。为了实现更快的收敛速度，Gu等[21] [22]采用了预条件Krylov子空间方法来求解Toeplitz系统。基于此，结合提出的离散系统具有对称的Toeplitz结构，并采用PCG方法作为预条件Krylov子空间方法中的快速算法，从而进一步优化了计算性能。

本文针对含阻尼项的非线性变分数阶扩散方程，构造二阶隐式有限差分格式，分析其稳定性与收敛性，并提出两种快速求解策略：1) 基于有限差分的多重网格法，利用网格嵌套迭代加速求解，克服传统迭代法的收敛缺陷；2) 设计PCG方法中预条件矩阵的优化选取方案，通过降低条件数显著提升计算效率。两种方法均通过结构化离散特性与算法创新实现高效求解。

2. 变分数阶扩散方程的离散化

将区间$\left[0,L\right]$ 分成M等分，$\left[0,T\right]$ 分成N等分，取空间步长$h=\frac{L}{M}$ ，时间步长$\tau =\frac{T}{N}$ ，记$x_j=jh\left(0\le j\le M\right)$ 和$t_n=n\tau \left(\text{0}\le n\le N\right)$ ，${\Omega }_h=\left\{x_j|0\le j\le M\right\}$ ，${\Omega }_{\tau }=\left\{t_k|\text{0}\le n\le N\right\}$ ，$N,M$ 都是正整数。

定义${\Omega }_h$ 上的网格函数$V_h=\left\{v|v=\left(v_0,v_1,\cdots ,v_M\right)\right\},v\in V_h$ ，对于任意网格函数$u,v\in V_h$ ，我们定义以下符号

${\delta }_x^2{v}_j=\frac{1}{h^2}\left(v_{j+1}-2v_j+v_{j-1}\right),{\delta }_x{u}_{i+\frac{1}{2}}=\frac{1}{h}\left(u_{i+1}-u_i\right),$

$\left(u,v\right)=h{\displaystyle \sum _{i=1}^{M-1}{u}_i}{v}_i,\left({\delta }_{x}u,{\delta }_{x}v\right)=h{\displaystyle \sum _{i=0}^{M-1}\left({\delta }_x{u}_{i+\frac{1}{2}}\right)}\left({\delta }_x{v}_{i+\frac{1}{2}}\right)$

$\left(u,v\right)=h{\displaystyle \sum _{j=1}^{M-1}{u}_j}{v}_j,\Vert u\Vert =\sqrt{\left(u,u\right)}$

然后使用[23]中推导的L2_1_σ公式用于离散Caputo分数导数，下面给出以下公式：

$\begin{array}{c}{\Delta }_t^{{\alpha }_n}f\left(t_{n-1+{\sigma }_n}\right)=\frac{{\tau }^{1-{\alpha }_n}}{\Gamma \left(2-{\alpha }_n\right)}{\displaystyle \sum _{k=1}^n{c}_{n-k}^{\left(n,{\alpha }_n\right)}}{\delta }_t{f}^{k-\frac{1}{2}}\\ =\frac{{\tau }^{1-{\alpha }_n}}{\Gamma \left(2-{\alpha }_n\right)}{\displaystyle \sum _{k=0}^{n-1}{c}_k^{\left(n,{\alpha }_n\right)}}{\delta }_t{f}^{n-k-\frac{1}{2}}\\ =\frac{{\tau }^{-{\alpha }_n}}{\Gamma \left(2-{\alpha }_n\right)}{\displaystyle \sum _{k=0}^{n-1}{c}_k^{\left(n,{\alpha }_n\right)}}\left[f\left(t_{n-k}\right)-f\left(t_{n-k-1}\right)\right],1\le n\le N.\end{array}$ (4)

对于$c_0^{\left(1,{\alpha }_n\right)}={\sigma }^{1-{\alpha }_n}$ ，并且当$n\ge 2$ 时有如下定义：

$\{\begin{array}{l}{c}_0^{\left(n,\alpha \right)}=\frac{{\left(1+{\sigma }_n\right)}^{2-{\alpha }_n}-{\sigma }_n^{2-\alpha }}{2-{\alpha }_n}-\frac{{\left(1+{\sigma }_n\right)}^{1-\alpha }-{\sigma }_n^{2-{\alpha }_n}}{2-{\alpha }_n}\\ c_k^{\left(n,\alpha \right)}=\frac{1}{2-{\alpha }_n}\left[{\left(k+1+{\sigma }_n\right)}^{2-{\alpha }_n}-2{\left(k+{\sigma }_n\right)}^{2-{\alpha }_n}+{\left(k-1+{\sigma }_n\right)}^{2-{\alpha }_n}\right]\\ -\frac{1}{2}\left[{\left(k+1+{\sigma }_n\right)}^{1-{\alpha }_n}-2{\left(k+{\sigma }_n\right)}^{1-{\alpha }_n}+{\left(k-1+{\sigma }_n\right)}^{1-{\alpha }_n}\right]\\ c_{n-1}^{\left(n,\alpha \right)}=\frac{1}{2}\left[3{\left(n-1+{\sigma }_n\right)}^{1-{\alpha }_n}-{\left(n-2+{\sigma }_n\right)}^{1-{\alpha }_n}\right]\\ -\frac{1}{2-{\alpha }_n}\left[{\left(n-1+{\sigma }_n\right)}^{2-{\alpha }_n}-{\left(n-2+{\sigma }_n\right)}^{2-{\alpha }_n}\right].\end{array}$

引理1 [23] 设$0<\alpha \left(t\right)<1$ ，$f\left(t\right)\in C^3\left[0,T\right]$ ，则有

${{}_0^C{D}_t^{\alpha \left(t\right)}f\left(t\right)|}_{t=t_{n+{\sigma }_n}}-{\Delta }_t^{\alpha }f\left(t_n+{\sigma }_n\right)=O\left({\tau }^{3-{\alpha }_n}\right),$

其中，${\sigma }_n=1-\frac{{\alpha }_n}{2}$ 。

对于阻尼项的高阶近似，利用泰勒公式进行近似处理，证明过程在[24]中给出。

引理2 [23] 对$0<\alpha \left(t\right)<1$ ，有

$c_0^{\left(n\right)}>c_1^{\left(n\right)}>\cdots >c_{n-1}^{\left(n\right)}>c_n^{\left(n\right)}>\frac{1-{\alpha }_{n+{\sigma }_n}}{2{\left(n+{\sigma }_n\right)}^{{\alpha }_{n+{\sigma }_n}}}>0.$

引理3 [24] 设$u\in C^3\left[0,T\right]$ ，则有

$\begin{array}{c}{\delta }_{\widehat{t}}u\left(t_n\right)\equiv \frac{1}{2\tau }\left[\left(2{\sigma }_n+1\right)u\left(t_{n+1}\right)-4{\sigma }_{n}u\left(t_n\right)+\left(2{\sigma }_n-1\right)u\left(t_{n-1}\right)\right]\\ =\frac{\text{d}u\left(t_{n+{\sigma }_n}\right)}{\text{d}t}+O\left({\tau }^2\right),n\ge 1.\end{array}$

引理4 [25] 设$f\in C^4\left[x_{j-1},x_{j+1}\right]$ ，则有

$f^{″}\left(x_j\right)=\frac{f\left(x_{j-1}\right)-2f\left(x_j\right)+f\left(x_{j+1}\right)}{h^2}-\frac{h^2}{6}{\displaystyle {\int }_0^1\left[f^{\left(4\right)}\left(x_j+\lambda h\right)+f^{\left(4\right)}\left(x_j-\lambda h\right)\right]}{\left(1-\lambda \right)}^3\text{d}\lambda .$

利用引理1~4、式子(4)，对方程(1)在点$\left(x,t\right)=\left(x_j,t_{n+{\sigma }_n}\right)$ 处进行近似，便会得到下面的隐式格式：

${\delta }_{\widehat{t}}u\left(x_j,t_n\right)+{\Delta }_{\tau }^{{\alpha }_{n+{\sigma }_n}}u\left(x_j,t_{n+{\sigma }_n}\right)+O\left({\tau }^2\right)={\delta }_x^{2}u\left(x_j,t_{n+{\sigma }_n}\right)+f\left(x_j,t_{n+{\sigma }_n}\right)+O\left(h^2\right).$ (5)

在网格点$\left(x_j,t_{n+{\sigma }_n}\right)$ 处考虑泰勒公式展开

$u\left(x_j,t_{n+{\sigma }_n}\right)={\sigma }_{n}u\left(x_j,t_{n+1}\right)+\left(1-{\sigma }_n\right)u\left(x_j,t_n\right)+O\left({\tau }^2\right),$

将上式代入(5)中可得：

${\delta }_{\widehat{t}}{U}_j^n+{\Delta }_{\tau }^{{\alpha }_{n+{\sigma }_n}}{U}_j^{n+{\sigma }_n}={\sigma }_n{\delta }_x^2{U}_j^{n+1}+\left(1-{\sigma }_n\right){\delta }_x^2{U}_j^n+f\left(x_j,t_{n+{\sigma }_n}\right)+R_j^n,$ (6)

其中，定义$R_j^n=O\left({\tau }^2+h^2\right)$ $1\le j\le M-1,0\le n\le N-1$ 。

注意到初边值条件(2)~(3)，有

$U_j^0=\phi \left(x_j,0\right),0\le j\le M,$

$U_0^n=\varphi \left(t_n\right),U_M^n=\psi \left(t_n\right).$

在方程(6)中略去小量项$R_j^n$ ，并用数值解$u_i^n$ 代替精确解$U_i^n$ ，得到求解问题(1)~(3)的如下差分格式：

$\{\begin{array}{ll}{\delta }_{\widehat{t}}{u}_j^n+{\Delta }_{\tau }^{{\alpha }_{n+{\sigma }_n}}{u}_j^{n+{\sigma }_n}={\sigma }_n{\delta }_x^2{u}_j^{n+1}+\left(1-{\sigma }_n\right){\delta }_x^2{u}_j^n+f\left(x_j,t_{n+{\sigma }_n}\right),1\le j\le M-1,0\le n\le N-1,\hfill & \text{(7)}\hfill \\ u_j^0=\phi \left(x_j,0\right),0\le j\le M,\hfill & \text{(8)}\hfill \\ u_0^n=\varphi \left(t_n\right),u_M^n=\psi \left(t_n\right),1\le n\le N.\hfill & \text{(9)}\hfill \end{array}$

该格式的无条件稳定性和收敛性将被证明。

3. 差分格式稳定性和收敛性分析

下面将会给出几个引理在之后的证明中将会用到。

引理5 [25] 对任意的网格函数$u^0,u^1,\cdots ,u^{n+1}\in {\mathcal{U}}_h$ ，有

$\left({\delta }_{\widehat{t}}{u}^k,{\sigma }_n{u}^{k+1}+\left(1-{\sigma }_n\right)u^k\right)\ge \frac{1}{4\tau }\left(E^{k+1}-E^k\right),k\ge 1$

其中，$E^{k+1}=\left(2{\sigma }_n+1\right){\Vert u^{k+1}\Vert }^2-\left(2{\sigma }_n-1\right){\Vert u^k\Vert }^2+\left(2{\sigma }_n^2+{\sigma }_n-1\right){\Vert u^{k+1}-u^k\Vert }^2$ ，和$E^{k+1}\ge \frac{1}{{\sigma }_n}{\Vert u^{k+1}\Vert }^2,k\ge 0$ ，并且$E^0=2{\Vert u^0\Vert }^2$ 。

引理6 对任意的网格函数$u^0,u^1,\cdots ,u^{n+1}\in {\mathcal{U}}_h$ ，下面不等式成立

${\displaystyle \sum _{k=0}^n{c}_{n-k}^{\left(n\right)}}\left(u^{k+1}-u^k,{\sigma }_n{u}^{n+1}+\left(1-{\sigma }_n\right)u^n\right)\ge \frac{1}{2}{\displaystyle \sum _{k=0}^n{c}_{n-k}^{\left(n\right)}}\left({\Vert u^{k+1}\Vert }^2-{\Vert u^k\Vert }^2\right).$

引理7 任意的网格函数$u\in {\mathcal{U}}_h$ ，下面不等式成立

${\Vert u\Vert }^2\le \frac{L^2}{6}{\Vert {\delta }_{x}u\Vert }^2.$

3.1. 稳定性分析

定理1 差分格式(7)是稳定的，且有

${\Vert u^1\Vert }^2\le C_1{\Vert u^0\Vert }^2+4{\sigma }^{{\alpha }_{n+{\sigma }_n}}{T}^{1+{\alpha }_{n+{\sigma }_n}}\Gamma \left(2-{\alpha }_{n+{\sigma }_n}\right){\Vert f^{{\sigma }_n}\Vert }^2,$

${\Vert u^m\Vert }^2\le C_2{\Vert u^1\Vert }^2+C_3{\Vert u^0\Vert }^2+8{\sigma }_n{T}^{1+{\alpha }_{n+{\sigma }_n}}\Gamma \left(1-{\alpha }_{n+{\sigma }_n}\right){\displaystyle \sum _{n=1}^{m-1}{\Vert f^{n+{\sigma }_n}\Vert }^2},m\ge 2$

其中，$C_{1-5}$ 是与$h,\tau$ 无关的正数。

证明 在(7)式左右两边同时与$u^{n+{\sigma }_n}={\sigma }_n{u}^{n+1}+\left(1-{\sigma }_n\right)u^n$ 作内积，可得

$\left({\delta }_{\widehat{t}}{u}^{n+{\sigma }_n},u^{n+{\sigma }_n}\right)+\left({\Delta }_{\tau }^{{\alpha }_{n+{\sigma }_n}}{u}^{n+{\sigma }_n},u^{n+{\sigma }_n}\right)=\left({\delta }_x^2{u}^{n+{\sigma }_n},u^{n+{\sigma }_n}\right)+\left(f^{n+{\sigma }_n},u^{n+{\sigma }_n}\right).$

由引理7有

$\begin{array}{c}\left({\delta }_x^2{u}^{n+{\sigma }_n},u^{n+{\sigma }_n}\right)=-\Vert \left({\delta }_x{u}^{n+{\sigma }_n},u^{n+{\sigma }_n}\right)\Vert \\ \le -\frac{6}{L^2}{\Vert u^{n+{\sigma }_n}\Vert }^2\end{array}$

故结合引理5~6以及$\left({\delta }_x^2{u}^{n+{\sigma }_n},u^{n+{\sigma }_n}\right)\le 0$ 有

$\frac{1}{4\tau }\left(E^{n+1}-E^n\right)+\frac{1}{2}{\Delta }_{\tau }^{{\alpha }_{n+{\sigma }_n}}{\Vert u^{n+{\sigma }_n}\Vert }^2\le \left(f^{n+{\sigma }_n},u^{n+{\sigma }_n}\right)\text{.}$ (10)

下面证明过程分两步进行：

第一步，当n = 0时不等式(10)可以简化为：

$\frac{1}{4\tau }\left(E^1-E^0\right)+\frac{1}{2{\tau }^{{\alpha }_{n+{\sigma }_n}}\Gamma \left(2-{\alpha }_{n+{\sigma }_n}\right)}{c}_0^{\left({\alpha }_{n+{\sigma }_n},{\sigma }_n\right)}\left({\Vert u^1\Vert }^2-{\Vert u^0\Vert }^2\right)\le \left(f^{{\sigma }_n},u^{{\sigma }_n}\right)\text{.}$

并使用Cauchy-Schwarz不等式，可以得出

$\begin{array}{l}\left(1+\frac{2\tau {\sigma }_n}{T^{{\alpha }_{n+{\sigma }_n}}\Gamma \left(2-{\alpha }_{n+{\sigma }_n}\right)}{c}_0^{\left(n,{\alpha }_{n+{\sigma }_n}\right)}\right){\Vert u^1\Vert }^2\\ \le 8\tau {\sigma }_n{\epsilon }_1\left({\Vert u^1\Vert }^2+{\Vert u^0\Vert }^2\right)+\frac{\tau {\sigma }_n}{{\epsilon }_1}{\Vert f^{{\sigma }_n}\Vert }^2+\left(\frac{2\tau {\sigma }_n}{T^{{\alpha }_{n+{\sigma }_n}}\Gamma \left(2-{\alpha }_{n+{\sigma }_n}\right)}{c}_0^{\left(n,{\alpha }_{n+{\sigma }_n}\right)}+2{\sigma }_n\right){\Vert u^0\Vert }^2,{\epsilon }_1>0\end{array}$

其中，令${\epsilon }_1=\frac{c_0^{\left(n,{\alpha }_{n+{\sigma }_n}\right)}}{4T^{{\alpha }_{n+{\sigma }_n}}\Gamma \left(2-{\alpha }_{n+{\sigma }_n}\right)}$ 。

则可得

$\begin{array}{c}{\Vert u^1\Vert }^2\le \left(\frac{4\tau {\sigma }_n}{T^{{\alpha }_{n+{\sigma }_n}}\Gamma \left(2-{\alpha }_{n+{\sigma }_n}\right)c_0^{\left(n,{\alpha }_{n+{\sigma }_n}\right)}}+2\sigma \right){\Vert u^0\Vert }^2+\frac{4\tau {\sigma }_n{T}^{{\alpha }_{n+{\sigma }_n}}\Gamma \left(2-{\alpha }_{n+{\sigma }_n}\right)}{c_0^{\left(n,{\alpha }_{n+{\sigma }_n}\right)}}{\Vert f^{{\sigma }_n}\Vert }^2\\ \le \left(\frac{4T^{1-{\alpha }_{n+{\sigma }_n}}{\sigma }^{{\alpha }_{n+{\sigma }_n}}}{\Gamma \left(2-{\alpha }_{n+{\sigma }_n}\right)}+2{\sigma }_n\right){\Vert u^0\Vert }^2+4{\sigma }^{{\alpha }_{n+{\sigma }_n}}{T}^{1+{\alpha }_{n+{\sigma }_n}}\Gamma \left(2-{\alpha }_{n+{\sigma }_n}\right){\Vert f^{{\sigma }_n}\Vert }^2.\end{array}$

令$C_1=\left(\frac{4T^{1-{\alpha }_{n+{\sigma }_n}}{\sigma }^{{\alpha }_{n+{\sigma }_n}}}{\Gamma \left(2-{\alpha }_{n+{\sigma }_n}\right)}+2{\sigma }_n\right)$ 。

第二步，当$n\ge 1$ 时，将(10)从1到$m-1$ 求和，有

$\begin{array}{l}\frac{1}{4\tau }\left(E^m-E^1\right)+\frac{1}{2{\tau }^{{\alpha }_{n+{\sigma }_n}}\Gamma \left(2-{\alpha }_{n+{\sigma }_n}\right)}\left[c_0^{\left(n,{\alpha }_{n+{\sigma }_n}\right)}{\displaystyle \sum _{n=1}^{m-1}{\Vert u^{n+1}\Vert }^2}-{\displaystyle \sum _{j=1}^{m-1}{\displaystyle \sum _{k=2}^n\left(c_{n-k}^{\left(n,{\alpha }_{n+{\sigma }_n}\right)}-c_{n-k+1}^{\left(n,{\alpha }_{n+{\sigma }_n}\right)}\right)}}{\Vert u^m\Vert }^2\right]\\ \le \frac{1}{2{\tau }^{{\alpha }_{n+{\sigma }_n}}\Gamma \left(2-{\alpha }_{n+{\sigma }_n}\right)}{\Vert u^0\Vert }^2{\displaystyle \sum _{n=1}^{m-1}{c}_n^{\left(n,{\alpha }_{n+{\sigma }_n}\right)}}+\frac{1}{2{\tau }^{{\alpha }_{n+{\sigma }_n}}\Gamma \left(2-{\alpha }_{n+{\sigma }_n}\right)}{\Vert u^1\Vert }^2{\displaystyle \sum _{n=1}^{m-1}\left(c_{n-1}^{\left(n,{\alpha }_{n+{\sigma }_n}\right)}-c_n^{(n,{\alpha }_{n+{\sigma }_n})}\right)}\\ +{\displaystyle \sum _{n=1}^{m-1}\Vert f^{n+{\sigma }_n}\Vert }\left({\sigma }_n\Vert u^{n+1}\Vert +\left(1-{\sigma }_n\right)\Vert u^n\Vert \right).\end{array}$ (11)

由引理2来估计不等式左侧第二项，可得

$\begin{array}{l}\frac{1}{2{\tau }^{{\alpha }_{n+{\sigma }_n}}\Gamma \left(2-{\alpha }_{n+{\sigma }_n}\right)}\left[c_0^{\left(n,{\alpha }_{n+{\sigma }_n}\right)}{\displaystyle \sum _{n=1}^{m-1}{\Vert u^{n+1}\Vert }^2}-{\displaystyle \sum _{n=1}^{m-1}{\displaystyle \sum _{k=2}^n\left(c_{n-k}^{\left(n,{\alpha }_{n+{\sigma }_n}\right)}-c_{n-k+1}^{\left(n,{\alpha }_{n+{\sigma }_n}\right)}\right)}}{\Vert u^k\Vert }^2\right]\\ =\frac{1}{2{\tau }^{{\alpha }_{n+{\sigma }_n}}\Gamma \left(2-{\alpha }_{n+{\sigma }_n}\right)}{\displaystyle \sum _{n=2}^m{c}_{m-n}^{\left(n,{\alpha }_{n+{\sigma }_n}\right)}}{\Vert u^n\Vert }^2\\ \ge \frac{1}{2{\tau }^{{\alpha }_{n+{\sigma }_n}}\Gamma \left(2-{\alpha }_{n+{\sigma }_n}\right)}\frac{1-{\alpha }_{n+{\sigma }_n}}{2}{\left(n-1+{\sigma }_n\right)}^{-{\alpha }_{n+{\sigma }_n}}{\displaystyle \sum _{n=2}^m{\Vert u^n\Vert }^2}\\ \ge \frac{1}{4T^{{\alpha }_{n+{\sigma }_n}}\Gamma \left(1-{\alpha }_{n+{\sigma }_n}\right)}{\displaystyle \sum _{n=2}^m{\Vert u^n\Vert }^2}.\end{array}$ (12)

把(12)代入(11)式中，可得

$\begin{array}{l}{\Vert u^m\Vert }^2+\frac{\tau {\sigma }_n}{T^{{\alpha }_{n+{\sigma }_n}}\Gamma \left(1-{\alpha }_{n+{\sigma }_n}\right)}{\displaystyle \sum _{n=2}^m{\Vert u^n\Vert }^2}\\ \le {\sigma }_n{E}^1+\frac{2{\tau }^{1-{\alpha }_{n+{\sigma }_n}}{\sigma }_n\left(4-3{\alpha }_{n+{\sigma }_n}\right){\left(m-1+{\sigma }_n\right)}^{1-{\alpha }_{n+{\sigma }_n}}}{\Gamma \left(3-{\alpha }_{n+{\sigma }_n}\right)}{\Vert u^1\Vert }^2+\frac{2{\tau }^{1-{\alpha }_{n+{\sigma }_n}}{\sigma }_n{(m-1+{\sigma }_n)}^{1-{\alpha }_{n+{\sigma }_n}}}{\Gamma \left(2-{\alpha }_{n+{\sigma }_n}\right)}{\Vert u^0\Vert }^2\\ +4\tau {\sigma }_n{\epsilon }_2{\displaystyle \sum _{n=1}^{m-1}{\left[{\sigma }_n\Vert u^{n+1}\Vert +\left(1-{\sigma }_n\right)\Vert u^n\Vert \right]}^2}+\frac{\tau {\sigma }_n}{{\epsilon }_2}{\displaystyle \sum _{n=1}^{m-1}{\Vert f^{n+{\sigma }_n}\Vert }^2}\\ \le {\sigma }_n\left[\left(4{\sigma }_n^2+4{\sigma }_n-1\right){\Vert u^1\Vert }^2+4{\sigma }_n^2{\Vert u^0\Vert }^2\right]+\frac{2{\sigma }_n\left(4-3{\alpha }_{n+{\sigma }_n}\right)T^{1-{\alpha }_{n+{\sigma }_n}}}{\Gamma (3-{\alpha }_{n+{\sigma }_n})}{\Vert u^1\Vert }^2\\ +\frac{2{\sigma }_n{T}^{1-{\alpha }_{n+{\sigma }_n}}}{\Gamma \left(2-{\alpha }_{n+{\sigma }_n}\right)}{\Vert u^0\Vert }^2+8\tau {\sigma }_n{\epsilon }_2{\displaystyle \sum _{n=1}^m{\Vert u^n\Vert }^2}+\frac{\tau {\sigma }_n}{{\epsilon }_2}{\displaystyle \sum _{n=1}^{m-1}{\Vert f^{n+{\sigma }_n}\Vert }^2},{\epsilon }_2>0\end{array}$ (13)

其中，$E^1\le \left(4{\sigma }_n^2+4{\sigma }_n-1\right){\Vert u^1\Vert }^2+4{\sigma }_n^2{\Vert u^0\Vert }^2$ ，使${\epsilon }_2=\frac{1}{8T^{{\alpha }_{n+{\sigma }_n}}\Gamma \left(1-{\alpha }_{n+{\sigma }_n}\right)}$ 。

可得

$\begin{array}{l}{\Vert u^m\Vert }^2\le \left[\frac{T^{1-{\alpha }_{n+{\sigma }_n}}{\sigma }_n}{\Gamma \left(1-{\alpha }_{n+{\sigma }_n}\right)}+\frac{2{\sigma }_n\left(4-3{\alpha }_{n+{\sigma }_n}\right)T^{1-{\alpha }_{n+{\sigma }_n}}}{\Gamma \left(3-{\alpha }_{n+{\sigma }_n}\right)}+4{\sigma }_n^3+4{\sigma }_n^2-{\sigma }_n\right]{\Vert u^1\Vert }^2\\ +\left(\frac{2{\sigma }_n{T}^{1-{\alpha }_{n+{\sigma }_n}}}{\Gamma \left(2-{\alpha }_{n+{\sigma }_n}\right)}+4{\sigma }_n^3\right){\Vert u^0\Vert }^2+8{\sigma }_n{T}^{1+{\alpha }_{n+{\sigma }_n}}\Gamma \left(1-{\alpha }_{n+{\sigma }_n}\right){\displaystyle \sum _{n=1}^{m-1}{\Vert f^{n+{\sigma }_n}\Vert }^2}.\end{array}$

令$C_3=\frac{2{\sigma }_n{T}^{1-{\alpha }_{n+{\sigma }_n}}}{\Gamma \left(2-{\alpha }_{n+{\sigma }_n}\right)}+4{\sigma }_n^3$ ，$C_3=\frac{2{\sigma }_n{T}^{1-{\alpha }_{n+{\sigma }_n}}}{\Gamma \left(2-{\alpha }_{n+{\sigma }_n}\right)}+4{\sigma }_n^3$ ，由上述证明可得，定理1得证。

3.2. 收敛性分析

定理2 $U_j^n$ 是方程(1)式的精确解，$u_j^n$ 是方程(7)的数值解，且$\frac{\partial U\left(x,0\right)}{\partial t}=0$ ，定义$e_j^n=U_j^n-u_j^n\left(0\le j\le M,0\le n\le N\right)$ ，$e^n={\left[e_1^n,e_2^n,\cdots ,e_{M-1}^n\right]}^{\text{T}}$ ，有

$\Vert e^n\Vert \le C\left({\tau }^2+h^2\right),0\le n\le N,$

其中，C是一个与$h$ 、$\tau$ 无关的正常数。

证明 将

${\delta }_{\widehat{t}}{e}_j^n+{\Delta }_{\tau }^{{\alpha }_{n+{\sigma }_n}}{e}_j^{n+{\sigma }_n}={\sigma }_n{\delta }_x^2{e}_j^{n+1}+\left(1-{\sigma }_n\right){\delta }_x^2{e}_j^n+R_j^{n+{\sigma }_n},$

$e_j^0=0,0\le j\le M,$

$e_0^n=e_M^n=0,1\le n\le N.$

其中，$R_j^{n+{\sigma }_n}=O\left({\tau }^2+h^2\right)$ 。然后，由定理1的证明过程，有$\Vert e^n\Vert \le C\left({\tau }^2+h^2\right)$ ，$0\le n\le N$ 。故定理2得证。

4. 快速数值求解方法的设计与实现

4.1. 多重网格算法

算法1 多重网格法伪代码

1) 初始化：检查输入矩阵与右端项维度匹配，初始化解向量为零，计算初始残差，设定收敛容差与最大迭代次数。

2) 构建多重网格层次：递归生成多层网格结构，通过限制算子生成粗网格矩阵，存储插值与限制算子。

3) 主迭代循环：循环执行V-cycle迭代，更新解向量，计算残差，判断收敛条件。若未收敛则继续迭代。

4) V-cycle递归过程： a. 前光滑：使用Gauss-Seidel迭代平滑误差； b. 残差限制：将细网格残差投影到粗网格； c. 粗网格修正：递归调用V-cycle求解粗网格误差； d. 插值延拓：将粗网格修正插值回细网格； e. 后光滑：再次平滑修正后的解。

5) 限制与插值算子生成：限制算子(R)：通过加权平均将细网格节点映射到粗网格；插值算子(P)：通过线性插值将粗网格值延拓到细网格。

6) 粗网格直接求解：在最粗层使用直接解法(如反斜杠\)精确求解。

7) 光滑迭代：使用Gauss-Seidel方法进行局部误差平滑，迭代次数由参数控制。

8) 收敛判断：计算残差范数，若小于容差则标记收敛，否则继续迭代。

具体来讲，多重网格算法在以下几个关键方面展现出了显著的优势：其一，通过利用稀疏化矩阵结构以及快速矩阵–向量乘技术，大幅降低了计算复杂度，从原本的$O\left(N^2\right)$ 优化至$O\left(N\log N\right)$ ，显著提升了计算效率；其二，定制化的限制–插值算子能够有效地保持分数阶系统的非局部关联特性，从而更好地处理分数阶方程的非局部性问题。

4.2. 预处理共轭梯度法

算法2 PCG算法伪代码

1：初始化r = b - A * x0 , z = solve(M, r), p = z,k = 0 2：while k < max_iter: 3：alpha = dot(r, z) / dot(p, A * p) 4：x = x + alpha * p 5：r_new = r - alpha * A * p 6：if norm(r_new) <ε:break 7：z_new = solve(M, r_new) 8：beta = dot(r_new, z_new) / dot(r, z) 9：p = z_new + beta * p 10：r = r_new,z = z_new,k = k + 1 11：return x

在这一节中，我们考虑用带有预条件子$c_F\left(A\right)$ 的预条件共轭梯度法求解系统$Ax=b$ 的代价。我们知道，在预条件共轭梯度法的每次迭代中，我们必须计算某个向量v的矩阵–向量乘法$Av$ ，并求解系统

$c_F\left(A\right)y=d$

对于某个向量d，设A是一般的$n\times n$ 矩阵，其中每个块$c_F\left(A\right)$ 是A的Chan’s循环预条件子，我们主要考虑Toeplitz系统。

5. 数值算例

为了验证算法的有效性，我们给出了数值算例。考虑变阶分数阶阻尼方程

$\{\begin{array}{l}{}_0^C{D}_t^{\alpha \left(t\right)}u+\frac{\partial u}{\partial t}=\frac{{\partial }^{2}u}{\partial x^2}+f\left(x,t\right),0<x<1,0<t<1\\ u\left(x,0\right)=0,0\le x\le 1\\ u\left(0,t\right)=0,u\left(1,t\right)=0,0<t\le 1.\end{array}$

它的真解为$u\left(x,t\right)=t^{2}x\left(1-x\right)$ ，源项为：

$f\left(x,t\right)=\frac{2x\left(1-x\right)}{\Gamma \left(3-\alpha \left(t\right)\right)}{t}^{2-\alpha \left(t\right)}+2tx\left(1-x\right)+2t^2$ .

定义误差及收敛阶如下：

$Error\left(\tau ,h\right)=\underset{1\le n\le N}{\max }\Vert u^n-U^n\Vert$

$rate={\log }_2\frac{Error\left(2\tau ,2h\right)}{Error\left(\tau ,h\right)},rat{e}_{\tau }={\log }_2\frac{Error\left(2\tau ,h\right)}{Error\left(\tau ,h\right)}.$

本节通过数值实验验证非线性时间分数阶延迟模型数值格式(7)~(9)的可靠性，并对比快速算法与直接法的计算效率。由收敛阶的定义，其次比较不同方法的总CPU时间，得到表1~4。

Table 1. Errors and convergence orders for different α(t) and τ at M = 1/h = 50

表1. 在M = 1/h = 500时，对不同的α(t)和τ的误差和收敛阶

	$\alpha \left(t\right)=\cos \left(t\right)$		$\alpha \left(t\right)=1-t$
$\tau$	误差	收敛阶	误差	收敛阶
1/100	1.2007e−05	—	2.4789e−05	—
1/200	3.0000e−06	2.0001	6.2123e−06	1.9965
1/400	7.4978e−07	2.0000	1.5549e−06	1.9983
1/800	1.8742e−07	2.0002	3.8897e−07	1.9991
1/1600	4.6851e−07	2.0001	9.7272e−08	1.9996

Table 2. Errors and convergence orders when h and τ vary with α(t)

表2. 当h和τ随α(t)变化时的误差和收敛阶

		$\alpha \left(t\right)=\cos \left(t\right)$		$\alpha \left(t\right)=1-t$
h	$\tau$	误差	收敛阶	误差	收敛阶
1/64	1/40	1.1262e−04	—	1.5381e−04	—
1/128	1/80	2.8328e−05	1.9911	3.8687e−05	1.9913
1/256	1/160	7.1035e−06	1.9956	9.7009e−06	1.9957
1/512	1/320	1.7785e−06	1.9978	2.4289e−06	1.9978
1/1024	1/640	4.4497e−07	1.9989	9.7272e−07	1.9989

Table 3. Comparison of CPU(s) of different algorithms at α(t) = exp(−t)

表3. 在α(t) = exp(−t)时不同算法的CPU(s)比较

$h=\tau$	DCPU	PCPU	MCPU
1/128	0.5068	0.1689	0.7808
1/256	0.9980	0.3741	1.3059
1/512	5.9302	1.0706	5.7400
1/1024	44.7123	4.1714	28.6763
1/2048	2402.4512	23.8801	175.2571

Table 4. Performance comparison of different preconditioned subs at α(t) = exp(−t)

表4. 不同预条件子在α(t) = exp(−t)时的性能比较

$h=\tau$	无预条件子-CPU	T. Chan预条件-CPU	对角预条件-CPU
1/128	0.0851	0.1689	0.0748
1/256	0.4589	0.3741	0.3654
1/512	5.6234	1.0706	2.1640
1/1024	50.6234	4.1714	15.2364
1/2048	124.5213	23.8801	51.2486

变阶分数阶导数$\alpha \left(t\right)=\cos \left(t\right)$ 和$\alpha \left(t\right)=1-t$ 两种情形下，不同时间步长对应的误差与收敛阶。数值实验表明，该格式在时间方向可达到二阶精度，见表1。所得数据验证了理论推导的时间收敛性。时间步长$\tau$ 与空间步长h同步变化。实验结果表明，全局收敛阶为二阶，与理论预期一致，见表2。在$\alpha \left(t\right)=\exp \left(-t\right)$ 时，不同算法(直接法DCPU、Chan’s循环预条件子预处理共轭梯度法PCPU、多重网格法MCPU)的CPU时间随网格加密($h=\tau$ 减小)的变化趋势。当网格逐渐加密(即ℎ从1/128减小至1/2048)时，三种算法的计算时间呈现显著差异。直接法的计算耗时随网格加密呈爆炸式增长，从ℎ = 1/128时的0.5068秒急剧上升至ℎ = 1/2048时的2402.45秒。这一现象源于直接法的立方级时间复杂度O(n³)相比之下，预处理共轭梯度法(PCPU)和多重网格法(MCPU)的计算时间增长较为平缓，见表3。这一结果表明，直接法因计算复杂度高，难以适应高精度或大规模问题的需求，而预处理共轭梯度法等快速迭代算法凭借其低复杂度和预处理优势，在计算效率上更具竞争力。为了体现预条件矩阵对PCG算法效率的影响，设计了三组不同预条件子，预条件技术显著提升了PCG算法的效率，见表4。对比不同预处理子的CPU，最终得出Chan循环预条件在预处理时间取得了良好平衡。因此，在实际工程与科学计算中，需综合考虑精度需求与计算成本，优先选择快速迭代算法以提升求解效率，尤其在网格高度加密的场景下。

6. 结论

本文给出了一个数值方案，该方案结合了L2_1_σ公式和二阶中心差分方法，以求解具有阻尼项的非线性变分数阶扩散方程。使用数学归纳法和离散能量法的结合，我们严格证明了该方案的稳定性和收敛性。数值实验验证了理论分析的时间和空间收敛阶接近二阶精度。此外，我们比较了三种求解方法的收敛速度：直接求解、带有T. Chan的循环预调节器的PCG方法、多重网格法。结果表明，PCG方法在实现更快的速度方面优于其他方法。

基金项目

长沙市自然科学基金(kq2502074)、湖南省自然科学基金(2025JJ50034)和湖南省教育厅科学研究基金(23A0266)。

NOTES

^*第一作者。