空间分数阶扩散方程数值解中结构化方程组的预处理方法研究

A Research of Preconditioning Method for Structured Systems in Numerical Solutions of Spatial Fractional Diffusion Equations

摘要: 空间分数阶扩散方程能有效地描述众多科学领域中的反常扩散现象。而大多数FDE难以得到解析解，需通过建立离散格式以获得高精度的数值解。数值求解FDE通常归结为线性方程组的求解，而预处理技术则是加速迭代求解的关键。近年来，基于不同离散格式下系数矩阵的结构和性质，学者们研究了高效的预处理方法，显著地降低了计算成本。针对数值求解空间分数阶扩散方程问题，本文整理和分析了方程不同形式下的离散情形和预处理方法，并为预处理进一步的研究提供思路参考。

Abstract: Spatial fractional diffusion equations effectively describe anomalous diffusion phenomena in various scientific fields. However, most FDEs are difficult to solve analytically, necessitating the establishment of discrete schemes to obtain high-precision numerical solutions. Numerical solutions of FDEs typically reduce to solving linear systems, where preconditioning techniques are crucial for accelerating iterative solvers. In recent years, scholars have investigated efficient preconditioning methods based on the structure and properties of coefficient matrices under different discretization schemes, significantly reducing computational costs. For the numerical solution of spatial fractional diffusion equations, this paper organizes and analyzes the discrete scenarios and preconditioning methods for various forms of spatial fractional diffusion equations, providing insights and references for further research in preconditioning.

1. 问题背景

分数阶微分方程基于分数阶微积分的概念，是对传统整数阶微分方程的扩展，包含一个或多个分数阶导数项[1] [2]。分数阶微分算子具有非局部性，适用于描述具有记忆效应和遗传特性的材料[3]-[5]。通过在空间变量中引入分数阶导数，可用于刻画物质在空间上的反常扩散现象。空间分数阶扩散方程在生物医学、材料科学等多个研究领域中得到了广泛的应用[6] [7]。

1.1. 空间分数阶扩散方程

对于一维变系数双侧空间分数阶扩散方程的初边值问题：

$\frac{\partial u\left(x,t\right)}{\partial t}=d_{+}\left(x,t\right)\frac{{\partial }^{\alpha }u\left(x,t\right)}{{\partial }_{+}{x}^{\alpha }}+d_{-}\left(x,t\right)\frac{{\partial }^{\alpha }u\left(x,t\right)}{{\partial }_{-}{x}^{\alpha }}+f\left(x,t\right)$ (1)

$u\left(x_L,t\right)=u\left(x_R,t\right)=0,0\le t\le T$ (2)

$u\left(x,0\right)=u_0\left(x\right),x\in \left(x_L,x_R\right)$ (3)

其中$\alpha$ 为分数，左侧分数阶导数$\frac{{\partial }^{\alpha }u\left(x,t\right)}{{\partial }_{+}{x}^{\alpha }}$ 和右侧分数阶导数$\frac{{\partial }^{\alpha }u\left(x,t\right)}{{\partial }_{-}{x}^{\alpha }}$ 可为G-L分数阶导数或R-L分数阶导数，$d_{\pm }\left(x,t\right)$ 为扩散系数且$d_{\pm }\left(x,t\right)\ge 0$ ，$f\left(x,t\right)$ 为源项，(2)为边值条件，(3)为初值条件。

定义1 Grünwald-Letnikov (G-L)分数阶导数：设$\alpha$ 是一个正实数，令$n-1\le \alpha <n$ ，n为正整数，函数$f\left(t\right)$ 定义在区间$\left[a,b\right]$ 上，称

${}_a{}^{ }D{}_t^{\alpha }f\left(t\right)={\lim }_{h\to 0}{h}^{-\alpha }{\displaystyle {\sum }_{k=0}^{\left[\left(t-a\right)/h\right]}{\left(-1\right)}^k\left(\begin{array}{c}\alpha \\ k\end{array}\right)f\left(t-kh\right)}$ (4)

为函数$f\left(t\right)$ 的$\alpha$ 阶Grünwald-Letnikov分数阶导数，其中$t\in \left[a,b\right]$ ，$\left[z\right]$ 为不超过$z$ 的最大整数，$\left(\begin{array}{c}\alpha \\ k\end{array}\right)$ 表示二项式系数$\left(\begin{array}{c}\alpha \\ k\end{array}\right)=\frac{\alpha \left(\alpha -1\right)\cdots \left(\alpha -k+1\right)}{k!}$ 。

定义2 Riemann-Liouville (R-L)分数阶导数：设$\alpha$ 是一个正实数，令$n-1\le \alpha <n$ ，n为正整数，函数$f\left(t\right)$ 定义在区间$\left[a,b\right]$ 上，称

${}_a{}^{ }D{}_t^{\alpha }f\left(t\right)=\frac{{\text{d}}^n}{\text{d}{t}^n}\left(\frac{1}{\Gamma \left(n-\alpha \right)}{\displaystyle {\int }_a^t\frac{f\left(\tau \right)\text{d}\tau }{{\left(t-\tau \right)}^{\alpha -n+1}}}\right)$ (5)

为函数$f\left(t\right)$ 的$\alpha$ 阶Riemann-Liouville分数阶导数，其中$t\in \left[a,b\right]$ 。

1.2. 离散化

由于空间分数阶扩散方程的解析解往往难以显式给出，故高效的数值模拟已成为当前相关领域研究的前沿问题[7]-[10]。此外，基于加速器计算的混合精度求解方法能够显著提高数值解的精度[11]。求解分数阶偏微分方程的数值算法主要包括有限差分法、有限元法、谱方法和级数逼近法等[9] [12]-[14]。本节将主要介绍利用有限差分法对时间导数和空间分数阶导数进行离散，从而得到方程的离散形式。

通过定义时间步长$\Delta t=T/M$ 和空间网格大小$h=\left(x_R-x_L\right)/N$ (M和N都是正整数)，我们可以得到时间和空间上的离散点$t_m=m\Delta t\left(m=0,1,\cdots ,M\right)$ ，$x_i=x_L+ih\left(i=0,1,\cdots ,N\right)$ 。(1)中的时间导数可以用标准一阶时差商进行离散[15]-[17]，即：

$\frac{\partial u\left(x_i,t_m\right)}{\partial t}=\frac{u\left(x_i,t_m\right)-u\left(x_i,t_{m-1}\right)}{\Delta t}+O\left(\Delta t\right)$ (6)

当对空间分数阶导数利用移位Grünwald近似[18] [19]，可得：

$\frac{{\partial }^{\alpha }u\left(x_i,t_m\right)}{{\partial }_{+}{x}^{\alpha }}=h^{-\alpha }{\displaystyle {\sum }_{k=0}^{i+1}{g}_k^{\left(\alpha \right)}u}\left(x_{i-k+1},t_m\right)+O\left(h\right)$ (7)

$\frac{{\partial }^{\alpha }u\left(x_i,t_m\right)}{{\partial }_{-}{x}^{\alpha }}=h^{-\alpha }{\displaystyle {\sum }_{k=0}^{N-i+1}{g}_k^{\left(\alpha \right)}u}\left(x_{i+k-1},t_m\right)+O\left(h\right)$ (8)

其中$g_k^{\left(\alpha \right)}={\left(-1\right)}^k\left(\begin{array}{c}\alpha \\ k\end{array}\right),k=1,2,\cdots$ 。由此，我们可以得到方程无条件稳定的隐式有限差分格式：

$\frac{u\left(x_i,t_m\right)-u\left(x_i,t_{m-1}\right)}{\Delta t}=\frac{d_{+}\left(x_i,t_m\right)}{h^{\alpha }}{\displaystyle \sum }_{k=0}^{i+1}{g}_k^{\left(\alpha \right)}u\left(x_{i-k+1},t_m\right)+\frac{d_{-}\left(x_i,t_m\right)}{h^{\alpha }}{\displaystyle \sum }_{k=0}^{N-i+1}{g}_k^{\left(\alpha \right)}u\left(x_{i+k-1},t_m\right)+f\left(x_i,t_m\right)$ (9)

其中$i=1,2,\cdots ,N$ ，$m=1,2,\cdots ,M$ 。

1.3. 预处理

利用空间分数阶扩散方程的离散形式，可以进一步转化为线性方程组的求解问题。不同的离散格式对应着具有不同结构特征的方程组。由于分数阶微分算子的非局部性，系数矩阵通常是大型稠密矩阵，因此通常选择Krylov子空间法来进行迭代求解[20]-[22]。然而，方程组的收敛性质受系数矩阵特征值分布的影响，当特征值分布较为发散时，Krylov子空间迭代求解的效率会很低，因此需要对方程组进行预处理。特别地，在约束优化问题中，多重网格预处理使得PCG迭代求解的次数显著减少[23]。近年来，针对不同特点的结构化方程组，学者们开展了多种预处理方法的研究，旨在提高求解效率[18] [24]-[26]。

不同的预处理方法取决于方程组的构造和对系数矩阵的近似，对于数值求解空间分数阶扩散方程中的结构化方程组，学者们通常基于矩阵结构和性质来构造不同的预处理矩阵，并应用到Krylov子空间法中求解该线性方程组。针对数值求解空间分数阶PDE中出现的结构化方程组，本文整理了近些年来提出的不同预处理方法，并进行了对比分析和总结。

本文具体结构如下：在第2部分，介绍了空间分数阶扩散方程不同离散格式下对应的结构化方程组，在第3部分，整理并对比分析了关于这些方程组的预处理方法。最后，本文对数值求解空间分数阶扩散方程中的预处理方法作出总结，并为进一步的研究提供了参考思路。

2. 结构化方程组

2.1. 对角矩阵与Toeplitz矩阵相乘结构的方程组

对于一维变系数双侧空间PDE，对时间导数采用一阶时差商离散格式，对空间导数采用移位Grünwald差分公式，得到数值求解框架[27]：

$A^{\left(m\right)}{u}^{\left(m\right)}=u^{\left(m-1\right)}+f^{\left(m\right)},m=1,2,\cdots ,M$ (10)

其中，系数矩阵$A^{\left(m\right)}$ 的形式为：

$I+D_{+}^{\left(m\right)}T+D_{-}^{\left(m\right)}{T}^{\text{T}}$ (11)

I为单位矩阵，$D_{\pm }^{\left(m\right)}$ 为对角矩阵，T为Toeplitz矩阵，${\delta }^{\left(m\right)}={\left[\delta \left(x_1,t_m\right),\delta \left(x_2,t_m\right),\cdots ,\delta \left(x_N,t_m\right)\right]}^{\text{T}}\left(\delta =u,f\right)$ ，系数矩阵的主要结构为对角矩阵与Toeplitz矩阵的乘积。空间上，$D_{\pm }^{\left(m\right)}$ 的对角元素对应于空间各离散点的扩散值，因此矩阵的对角元素并不相同，从而使得系数矩阵$A^{\left(m\right)}$ 不再具有Toeplitz结构。时间上，$D_{\pm }^{\left(m\right)}$ 的元素受时间步数m的影响，导致系数矩阵在时域上不断变化。系数矩阵在时空结构上的变化使得求解方程组具有较高的计算复杂度。

2.2. 对角矩阵与Toeplitz矩阵相加结构的方程组

对于一维单变系数双侧空间PDE，对时间导数采用一阶时差商离散格式，对空间导数采用移位Grünwald差分公式，得到数值求解框架[28]：

$A^{\left(m\right)}{u}^{\left(m\right)}=\Delta t{f}^{\left(m\right)}+D^{\left(m\right)}{u}^{\left(m-1\right)},m=1,2,\cdots ,M$ (12)

其中，系数矩阵$A^{\left(m\right)}$ 的形式为：

$D^{\left(m\right)}+T$ (13)

$D^{\left(m\right)}$ 为对角矩阵，T为对称Toeplitz矩阵，系数矩阵的结构为对角矩阵与Toeplitz矩阵的和。同样地，$D^{\left(m\right)}$ 的对角元素取决于空间离散点上的扩散值，使得系数矩阵$A^{\left(m\right)}$ 未能保持严格的对称Toeplitz结构。虽然矩阵T在整个时间过程中保持不变，但$D^{\left(m\right)}$ 的动态变化赋予了$A^{\left(m\right)}$ 在时间维度上额外的复杂性。

2.3. Toeplitz类矩阵之和结构的方程组

对于一维变系数双侧空间PDE，对时间导数采用边值公式，空间导数采用四阶quasi-compact离散格式[29]，得到数值求解框架：

$AU=b$ (14)

其中，系数矩阵A的形式为：

$B\otimes K-\Delta tC\otimes J$ (15)

K为三对角Toeplitz矩阵、J为对称Toeplitz矩阵、B和C均为块Toeplitz矩阵，系数矩阵的结构为Toeplitz类矩阵的组合。解矩阵U的定义为：

$U={\left(u_1^{\text{T}},u_2^{\text{T}},\cdots ,u_N^{\text{T}}\right)}^{\text{T}}$ (16)

其中$u_i={\left(u\left(x_1,t_i\right),u\left(x_2,t_i\right),\cdots ,u\left(x_N,t_i\right)\right)}^{\text{T}},i=1,2,\cdots ,M$ 。与2.1和2.2中的数值求解框架不同，此框架中矩阵U包含了所有时间层上各离散点的数值解，因此系数矩阵A是一个具有Kronecker积的大规模矩阵。由于其复杂的结构和庞大的规模，若直接使用Krylov子空间迭代法进行求解将面临较大的计算量。

3. 预处理方法

针对数值求解空间分数阶扩散方程中不同特点的结构化方程组，本节分别整理了对应的预处理方法，并从不同角度进行了对比分析。

3.1. 针对对角矩阵与Toeplitz矩阵相乘结构的方程组的预处理

3.1.1. 预处理形式及性质

关于方程组(10)，忽略时间层面的影响，对每一步的方程组进行预处理，简写为：

$Au=b$ (17)

此时系数矩阵为

$A=I+D_{+}^{}T+D_{-}^{}{T}^{\text{T}}$ (18)

Pan等人[27]通过近似系数矩阵的逆来构造预处理矩阵。记$d_{+}\left(x_i\right)=d_{+,i},d_{-}\left(x_i\right)=d_{-,i},i=1,2,\cdots ,N$ ，则$D_{+}=diag\left(d_{+,1},d_{+,2},\cdots ,d_{+,N}\right)$ 和$D_{-}=diag\left(d_{-,1},d_{-,2},\cdots ,d_{-,N}\right)$ 。定义$K_i\triangleq I+d_{+,i}T+d_{-,i}{T}^{\text{T}}$ ，显然所有的$K_i$ 都是Toeplitz矩阵。根据$e_i^{\text{T}}A=e_i^{\text{T}}{K}_i$ ，那么$e_i^{\text{T}}{A}^{-1}\approx e_i^{\text{T}}{K}_i^{-1}$ ，其中$e_i$ 表示单位矩阵的第i列，由于$K_i$ 是Toeplitz矩阵，因此可用循环矩阵来近似$K_i$ 。设C为矩阵T的循环近似，记$C_i=I+d_{+,i}C+d_{-,i}{C}^{\text{T}}$ 。由此得到了预处理矩阵$B_1$ ，其逆定义为：

$B_1^{-1}={\displaystyle {\sum }_{i=1}^N{e}_i{e}_i^{\text{T}}{C}_i^{-1}}$ (19)

然而，以$B_1$ 为预处理矩阵时，每次迭代的计算量为O(N)次FFT，为了减少计算成本，Pan利用插值方法构造了实用的预处理矩阵。设$\lambda$ 为某实部为正的复数，$\overline{\lambda }$ 表示$\lambda$ 的复共轭。定义函数

$q_{\lambda }\left(x\right)\triangleq \frac{1}{1+\lambda d_{+}\left(x\right)+\overline{\lambda }{d}_{-}\left(x\right)}$

$p_{\lambda }\left(x\right)={\phi }_1\left(x\right)q_{\lambda }\left({\tilde{x}}_1\right)+{\phi }_2\left(x\right)q_{\lambda }\left({\tilde{x}}_2\right)+\cdots +{\phi }_l\left(x\right)q_l\left({\tilde{x}}_l\right)$ 是基于$l$ ($l\ll N$ )个点的$q_{\lambda }\left(x\right)$ 分段线性插值，且${\left\{\left({\tilde{x}}_j,q_{\lambda }\left({\tilde{x}}_j\right)\right)\right\}}_{j=1}^l\subset {\left\{\left(x_i,q_{\lambda }\left(x_i\right)\right)\right\}}_{i=1}^N$ 。

记${\tilde{C}}_j\triangleq I+d_{+}\left({\tilde{x}}_j\right)C+d_{-}\left({\tilde{x}}_j\right)C^{\text{T}}$，我们已知循环矩阵可被傅里叶矩阵对角化，即$C=F\wedge F^{\text{*}}$ ，${\tilde{C}}_j=F{\tilde{\wedge }}_j{F}^{\text{*}}$，其中${\tilde{\wedge }}_j=I+d_{+}\left({\tilde{x}}_j\right)\wedge +d_{-}\left({\tilde{x}}_j\right){\wedge }^{\text{*}}\left(j=1,2,\cdots ,l\right)$。使用插值函数$p_{\lambda }\left(x\right)$ 来近似$C_i^{-1}$ ，可得

$C_i^{-1}\approx F\left({\displaystyle {\sum }_{j=1}^l{\phi }_j\left(x_i\right){\tilde{\wedge }}_j^{-1}}\right)F^{\text{*}}={\displaystyle {\sum }_{j=1}^l{\phi }_j\left(x_i\right){\tilde{C}}_j^{-1}},i=1,2,\cdots ,N$

其中${\phi }_j\left(x_i\right)$ 是插值系数，可以得到预处理矩阵$B_2$ ，其逆定义为：

$B_2^{-1}={\displaystyle {\sum }_{j=1}^l{\Phi }_{j}F{\tilde{\wedge }}_j^{-1}{F}^{\text{*}}}$(20)

其中，${\Phi }_j=diag\left({\phi }_j\left(x_1\right),{\phi }_j\left(x_2\right),\cdots ,{\phi }_j\left(x_N\right)\right)$ ，将$B_2^{-1}$ 作用于任意向量，计算量仅为$O\left(lN\log N\right)$ 。关于矩阵近似的理论分析表明，$B_2^{-1}-A^{-1}$ 可以被表示为小范数矩阵和低秩矩阵之和。

对于方程组(17)，Huang等人[30]构造了含参的预处理矩阵。记$H=\frac{1}{2}\left(T+T^{\text{T}}\right)$ 和$S=\frac{1}{2}\left(T-T^{\text{T}}\right)$ 分别

为矩阵T的Hermitian部分和skew-Hermitian部分，则$A=I+\left(D_{+}+D_{-}\right)H+\left(D_{+}-D_{-}\right)S$ 。当$D_{+}+D_{-}$ 非奇异时，方程组(17)可改写为：

$\left(D_{+}+D_{-}\right)\left[{\left(D_{+}+D_{-}\right)}^{-1}+H\right]u=\left(D_{+}+D_{-}\right)\left[-{\left(D_{+}+D_{-}\right)}^{-1}\left(D_{+}-D_{-}\right)Su+{\left(D_{+}+D_{-}\right)}^{-1}b\right]$ (21)

通过引入参数$\theta \left(\theta >0\right)$ ，得到不平衡缩放对角和Toeplitz分裂(LSDTS)迭代：

$\{\begin{array}{l}H{u}^{\left(k+\frac{1}{2}\right)}=-\left[{\left(D_{+}+D_{-}\right)}^{-1}+{\left(D_{+}+D_{-}\right)}^{-1}\left(D_{+}-D_{-}\right)S\right]u^{\left(k\right)}+{\left(D+W\right)}^{-1}b\\ \left[\theta I+{\left(D_{+}+D_{-}\right)}^{-1}\right]u^{\left(k+1\right)}=\left(\theta I-H\right)u^{\left(k+\frac{1}{2}\right)}-{\left(D_{+}+D_{-}\right)}^{-1}\left(D_{+}-D_{-}\right)S{u}^{\left(k+\frac{1}{2}\right)}+{\left(D_{+}+D_{-}\right)}^{-1}b\end{array}$ (22)

LSDTS迭代可以等效地表述为标准的一步矩阵分裂迭代：

$u^{\left(k+1\right)}=L\left(\theta \right)u^{\left(k\right)}+q,k=0,1,2,\cdots$

其中

$M\left(\theta \right)=\frac{1}{\theta }\left(D_{+}+D_{-}\right)H\left[\theta I+{\left(D_{+}+D_{-}\right)}^{-1}\right]$ (23)

$N\left(\theta \right)=-\frac{1}{\theta }\left(D_{+}+D_{-}\right)H\left(\theta I-H\right)H^{-1}{\left(D_{+}+D_{-}\right)}^{-1}-\left(D_{+}-D_{-}\right)S$ (24)

$L\left(\theta \right)=M{\left(\theta \right)}^{-1}N\left(\theta \right),q=M{\left(\theta \right)}^{-1}b$ (25)

为快速求解$M^{-1}\left(\theta \right)$ ，因H为对称Toeplitz矩阵，可用基于正弦变换的$\tau$ 矩阵来近似[31]-[34]，得到预处理矩阵：

$M_{\tau }\left(\theta \right)=\frac{1}{\theta }\left(D_{+}+D_{-}\right)H_{\tau }\left[\theta I+{\left(D_{+}+D_{-}\right)}^{-1}\right]$ (26)

预处理矩阵$M_{\tau }\left(\theta \right)$ 具有以下近似性质。

定理1 对$M_{\tau }{\left(\theta \right)}^{-1}A$ 的特征值位于以(1, 0)为中心，半径为$\varsigma \left(\theta \right)$ 的复平面圆盘内。 $\varsigma \left(\theta \right)=\max \left|\frac{\theta -\lambda }{\lambda }\right|-\theta \Phi \tan \left(\frac{\alpha \pi }{2}\right)$ ，${\lambda }_{\min }\le \lambda \le {\lambda }_{\max }$ ，其中${\lambda }_{\min }$ 和${\lambda }_{\max }$ 是矩阵H的最小和最大特征值。当$\Phi <-\frac{1}{{\lambda }_{\max }}\cot \left(\frac{\alpha \pi }{2}\right)$ 和$0<\theta <\frac{2{\lambda }_{\min }}{1-\Phi {\lambda }_{\min }\tan \left(\frac{\alpha \pi }{2}\right)}$ 成立时，$\varsigma \left(\theta \right)<1$ 。

3.1.2. 预处理对比与分析

通过对比分析，我们发现这两种预处理方法从不同的角度出发构造预处理矩阵。Pan等人基于近似矩阵逆的思想，利用“若两个矩阵同一行的元素相等，则其逆矩阵的同一行元素近似相等”的原理，构造了由N个Toeplitz矩阵的逆组成的逆预处理矩阵。为了减少计算量，对该逆预处理矩阵进行了逐步优化，最终得到了高效可用的预处理矩阵。而Huang等人则从直接近似系数矩阵的角度出发，基于矩阵分裂构造了具有对称Toeplitz结构的预处理矩阵，并进一步引入了$\tau$ 近似来提升算法效率。这两种方法的相同点在于，构造有效的预处理矩阵的关键都在于参数的合理选取，他们均利用了系数矩阵的特性以提高数值求解的效率。这两种预处理的构造方法和理论基础均不同，各具特色，各有优劣之处。

Pan的逐行逆近似预处理方法中随着插值点数量的增加，收敛所需的迭代次数虽然减少，然而形成和应用预处理矩阵的成本随插值点的数量成比例地增长。因此，在确定插值点的数量时需要进行权衡。Huang的预处理方法的有效性在很大程度上依赖于参数的选择。针对这一问题，我们提出一种直接的无参预处理方法，即利用矩阵分裂$A=P-M$ ，其中$P=\left(D_{+}+D_{-}\right)H$ ，$M=-\left(I+\left(D_{+}-D_{-}\right)S\right)$ ，进而构造预处理矩阵$P_{\tau }=\left(D_{+}+D_{-}\right)H_{\tau }$ 。Huang方法的特点在于对称Toeplitz结构的构造，对此我们的想法是通过求解$\text{min}{\Vert H-A\Vert }_F$ ，得到近似A的对称Toeplitz矩阵H，进一步利用$\tau$ 近似得到预处理矩阵$\tau \left(H\right)$ 。无论采用何种改进的构造方法，如何在近似程度和计算量之间实现平衡是较为关键的。

3.2. 针对对角矩阵与Toeplitz矩阵相加结构方程组的预处理

3.2.1. 预处理形式及性质

对于方程组(12)，对某一时刻的方程组进行预处理，方程组简写形式同(17)，此时系数矩阵为：

$A=D+T$ (27)

Yang等人[35]通过平均对角矩阵保证了系数矩阵的对称Toeplitz结构。即令$\omega =mean\left(D\right)$ ，这里$mean\left(D\right)$ 是指D的对角元素平均值，从而得到$P_{}=\omega I+T$ ，进一步用$\tau$ 矩阵近似，得到预处理矩阵$P_{\tau }=\tau \left(\omega I+T\right)$ 。

预处理矩阵$P_{\tau }$ 具有以下近似性质。

定理2 $\min \left(\frac{d_{\min }}{\omega },\frac{1}{2}\right)\le \lambda \left(P_{\tau }^{-1}A\right)\le \max \left(\frac{d_{\max }}{\omega },\frac{3}{2}\right)$ ，其中$d_{\min }$ 和$d_{\max }$ 是矩阵D对角元素的最小、最大值。

Bai等人[36]提出了对角和Toeplitz分裂的预处理方法，寻求了有效的预处理矩阵。通过引入参数$\theta$ ($\theta >0$ )，我们可以构造出A的两种分裂形式：

$A=\left(\theta I+D\right)-\left(\theta I-T\right)$ 和$A=\left(\theta I+T\right)-\left(\theta I-D\right)$ (28)

因此，可以通过两步迭代方法进行求解，即：

$\{\begin{array}{c}\left(\theta I+D\right)u^{\left(k+\frac{1}{2}\right)}=\left(\theta I-T\right)u^{\left(k\right)}+b \\ \left(\theta I+T\right)u^{\left(k+1\right)}=\left(\theta I-D\right)u^{\left(k+\frac{1}{2}\right)}+b\end{array}$ (29)

将上述的式子写成标准的一步迭代方法：

$P\left(\theta \right)u^{\left(k+1\right)}=Q\left(\theta \right)u^{\left(k\right)}+b$

$P\left(\theta \right)=\frac{1}{2\theta }\left(\theta I+D\right)\left(\theta I+T\right),Q\left(\theta \right)=\frac{1}{2\theta }\left(\theta I-D\right)\left(\theta I-T\right)$ (30)

为快速求解$P^{-1}\left(\theta \right)$ ，用$\tau$ 矩阵近似T，得到预处理矩阵

$P_{\tau }\left(\theta \right)=\frac{1}{2\theta }\left(\theta I+D\right)\left(\theta I+\tau \left(T\right)\right)$ (31)

预处理矩阵$P_{\tau }\left(\theta \right)$ 具有以下近似性质。

定理3 $P_{\tau }{\left(\theta \right)}^{-1}A$ 的所有特征值都聚集在以(1, 0)为中心，半径为$\sigma \left(\theta \right)$ 的复平面圆内。 $\sigma \left(\theta \right)=\max \left\{\frac{\left|\theta -d\right|}{\theta +d}\right\}\max \left\{\frac{\left|\theta -\lambda \right|}{\theta +\lambda }\right\}<1$ ，$d_{\min }\le d\le d_{\max }$ ，$t_{\min }\le \lambda \le t_{\max }$ 。其中$d_{\min }$ 和$d_{\max }$ 为对角矩阵D中的最小和最大元素值，$t_{\min }$ 和$t_{\max }$ 为矩阵T中的最小和最大特征值。

3.2.2. 预处理对比与分析

从近似思想来看，Yang从系数矩阵的整体出发，将对角矩阵近似为数量矩阵，这是一种相对简单的策略；而Bai则采用了基于矩阵分裂的较为复杂方法。从近似效果上看，Yang构造的预处理矩阵的有效性高度依赖于对角矩阵元素，而Bai的预处理方法能够确保预处理后方程组的系数矩阵的特征值聚集在1附近，从而保证了预处理的有效性。因此，相较之下，Bai的预处理方法更优。

Bai所构造的预处理矩阵同样涉及到参数的选取问题，可以研究最优参数的选择。此外，我们的想法是可以借鉴Pan的逐行逆近似预处理方法，亦或利用SMW公式来构造近似逆预处理矩阵。具体而言，使$A=A_0+{\displaystyle {\sum }_{i=1}^n{x}_k{y}_k^{\text{T}}}$ ，那么$A_{k+1}=A_k+x_{k+1}{y}_{k+1}^{\text{T}}$ 。令$x_k$ 、$y_k$ 满足$r_k=1+y_k^{\text{T}}{A}_{k-1}^{-1}{x}_k\ne 0$ ，则

$A^{-1}=A_0^{-1}-{\displaystyle {\sum }_{i=1}^n\frac{1}{r_k}{A}_{k-1}^{-1}{x}_k{y}_k^{\text{T}}{A}_{k-1}^{-1}}=A_0^{-1}-\Phi {\Omega }^{-1}{\Psi }^{\text{T}}$ ，其中$\Phi =\left[A_0^{-1}{x}_1{A}_1^{-1}{x}_2\cdots A_{n-1}^{-1}{x}_n\right]$ ， ${\Omega }^{-1}=diag\left(r_1^{-1},r_2^{-1},\cdots ,r_n^{-1}\right)$ ，${\Psi }^{\text{T}}={\left[y_1^{\text{T}}{A}_0^{-1}{y}_2^{\text{T}}{A}_1^{-1}\cdots y_n^{\text{T}}{A}_{n-1}^{-1}\right]}^{\text{T}}$ 。其中，可令$A_0=D$ ，利用这种矩阵分解方法，我们可以得到系数矩阵逆的较为精确的近似。

3.3. 针对Toeplitz类矩阵之和结构方程组的预处理

3.3.1. 预处理形式及性质

Zhou等[29]对方程组(14)，构造Kronecker积分裂(Kronecker product splitting)预处理。引入参数$\theta >0$ ，根据Kronecker积的性质，矩阵A允许如下的分裂：

$Q=\left(B+\theta C\right)\otimes K-C\otimes \left(\Delta tJ+\theta K\right)$ (32)

和

$Q=-C\otimes \left(\Delta tJ-\theta K\right)+\left(B-\theta C\right)\otimes K$ (33)

通过交替迭代这两个分裂，可以建立KPS迭代：

$\{\begin{array}{l}\left[\left(B+\theta C\right)\otimes K\right]U^{\left(k+\frac{1}{2}\right)}=\left[C\otimes \left(\Delta tJ+\theta K\right)\right]U^{\left(k\right)}+g\\ \left[-C\otimes \left(\Delta tJ-\theta K\right)\right]U^{\left(k+1\right)}=-\left[\left(B-\theta C\right)\otimes K\right]U^{\left(k+\frac{1}{2}\right)}+g\end{array}$

我们可以去掉中间变量$U^{k+\frac{1}{2}}$ ，将KPS迭代写成标准形式$U^{\left(k+1\right)}=G\left(\theta \right)U^{\left(k\right)}+Q\left(\theta \right)$ ，其中 $G\left(\theta \right)=\left[{\left(B+\theta C\right)}^{-1}\left(B-\theta C\right)\right]\otimes \left[{\left(\Delta tJ-\theta K\right)}^{-1}\left(\Delta tJ+\theta K\right)\right]$ ，$Q\left(\theta \right)=-2\theta {\left(B+\theta C\right)}^{-1}\otimes {\left(\Delta tJ-\theta K\right)}^{-1}g$ 。

实际上，此迭代也可由$A=P\left(\theta \right)-R\left(\theta \right)$ 产生，其中

$P\left(\theta \right)=-\frac{1}{2\theta }\left(B+\theta C\right)\otimes \left(\Delta tJ-\theta K\right)$ (34)

$R\left(\theta \right)=-\frac{1}{2\theta }\left(B-\theta C\right)\otimes \left(\Delta tJ+\theta K\right)$ (35)

此外，$G\left(\theta \right)=P{\left(\theta \right)}^{-1}R\left(\theta \right)=I-P{\left(\theta \right)}^{-1}A$ 。

预处理矩阵$P\left(\theta \right)$ 具有以下近似性质。

定理4 当$\rho \left(G\left(\theta \right)\right)=\max \left|\frac{\lambda -\theta }{\lambda +\theta }\right|\left|\frac{\mu +\theta }{\mu -\theta }\right|<\max \left|\frac{\lambda -\theta }{\lambda +\theta }\right|<1$ 时，$\lambda \in \sigma \left(C^{-1}B\right)$ ，$\mu \in \sigma \left(\Delta t{K}^{-1}J\right)$ ，$P{\left(\theta \right)}^{-1}A$ 的特征值完全包含在以(1, 0)圆心且半径小于1的圆中。

3.3.2. 预处理分析

Zhou基于矩阵分裂方法构造了预处理矩阵，从近似程度上看，可知预处理后的系数矩阵具有良好的谱性质。使用预处理Krylov子空间法求解时，主要的计算量在于预处理矩阵的逆和任一向量的乘积，但Zhou的方法中直接计算$P{\left(\theta \right)}^{-1}$ 的代价较高。为此，我们可以对$P\left(\theta \right)$ 做进一步的近似。具体而言，矩阵B和C可以用Strang型块循环矩阵来近似，K和J同样可以用循环矩阵来近似。由此，预处理矩阵可表示为：

$S\left(P\left(\theta \right)\right)=-\frac{1}{2\theta }\left(s\left(B\right)+\theta s\left(C\right)\right)\otimes \left(\Delta ts\left(J\right)-\theta s\left(K\right)\right)$

根据循环矩阵可以被傅里叶矩阵对角化，我们可得：

$S{\left(P\left(\theta \right)\right)}^{-1}=-2\theta \left(F\left({\Lambda }_B^{-1}+\frac{1}{\theta }{\Lambda }_C^{-1}\right)F^{*}\right)\otimes \left(F\left(\Delta t{\Lambda }_J^{-1}-\frac{1}{\theta }{\Lambda }_K^{-1}\right)F^{*}\right)$

其中，${\Lambda }_{\delta }^{}$ 是由矩阵$s\left(\delta \right)$ 的特征值构成的对角矩阵。由于系数矩阵A为Toeplitz类矩阵，我们的一个思路是利用生成函数构造近似逆预处理矩阵。具体而言，我们可以结合扩散系数和源项对矩阵A的生成函数提出近似。基于生成函数，我们可以估计出系数矩阵的特征值分布。进一步地，利用HNP过滤器或Tikhonov过滤器[37]，并结合逆傅里叶变换(IFFT)，可以构造出有效的循环近似逆预处理矩阵。

4. 总结

针对数值求解空间分数阶扩散方程中的结构化方程组，学者们围绕系数矩阵的结构和性质开展了一系列预处理方法的研究，显著提高了Krylov子空间迭代求解的效率。本文对近年来提出的预处理方法进行了整理，并从不同角度进行了对比分析。我们发现，对于任何结构的方程组，都可以通过两种途径来构造预处理矩阵：一是直接近似系数矩阵，二是近似系数矩阵的逆。在直接近似系数矩阵的方法中，典型的有矩阵分裂和平均对角矩阵。矩阵分裂方法通过引入参数，利用矩阵的不同分裂形式进行交替迭代，从而构造出可调控的含参预处理矩阵。平均对角矩阵法主要依赖于对角矩阵对系数矩阵的影响，因而不能确保得到高质量的近似。在近似系数矩阵逆的方法中，逐行逆近似能够有效地构造出系数矩阵的逆形式。然而，这一方法需要通过插值对计算过程进行优化，在此过程中，平衡好插值点的数量与形成预处理矩阵的成本是一个关键问题。结果表明，相比原有的方程组，以上预处理后的方程组的收敛性质明显改善，Krylov子空间迭代求解的速度显著升高。

关于空间分数阶扩散方程数值解中结构化方程组的预处理方法研究，如何使得预处理后系数矩阵的特征值更加聚集并减少计算成本，本文认为还可以从如下方面着手：

一、获得更有效的近似逆预处理矩阵：可以尝试结合扩散系数和源项，对矩阵元素序列进行傅里叶变换，从而得到系数矩阵生成函数的近似。进一步根据特征值的分布进行滤波，并结合IFFT构造循环近似逆预处理矩阵。此外，还可以应用 SMW公式，通过矩阵分解得到系数矩阵逆的精确形式。在构造过程中，需在近似程度与构建预处理矩阵的成本之间进行权衡。

二、构造更精确的对称Toeplitz近似矩阵：对于任意矩阵$A\in C^{n\times n}$ ，其对称Toeplitz近似矩阵H，可以尝试通过求解最小值问题$\text{min}{\Vert H-A\Vert }_F$ 来获得。那么，我们可以从系数矩阵的整体或其关键部分出发，利用此最优问题构造出更精确的对称Toeplitz近似矩阵，进一步基于正弦变换得到预$\tau$ 处理矩阵。

基金项目

北京工商大学学科建设经费资助项目(项目编号：STKY202508)。

北京工商大学大学生创新创业训练计划项目(项目编号：X202410011027)。

NOTES

^*通讯作者。