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

Preconditioning Methods for Structured Systems in the Numerical Solution of Space-Fractional Cahn-Hilliard Equations

摘要: 空间分数阶Cahn-Hilliard方程能精确描述具有反常扩散、复杂界面动力学或长程关联系统中的相演化过程。由于方程的非局部性与刚性导致解析解难以获得，需建立能量稳定的离散格式。数值求解空间分数阶Cahn-Hilliard方程通常归结为求解线性方程组，构造有效预处理器可加速迭代求解。近年来，学者们基于不同离散格式下系数矩阵的结构与性质提出了不同的预处理方法，显著提高了求解效率。针对空间分数阶Cahn-Hilliard方程数值求解问题，本文分别整理并分析了一维和高维空间分数阶Cahn-Hilliard方程在不同离散情形下的相应预处理方法，并为后续研究提供了建议。

Abstract: The space-fractional Cahn-Hilliard equation provides an accurate description of phase evolution processes in systems exhibiting anomalous diffusion, complex interfacial dynamics, or long-range interactions. Due to the nonlocal nature and stiffness of the equation, analytical solutions are difficult to obtain, necessitating the development of energy-stable discretization schemes. The numerical solution of the space-fractional Cahn-Hilliard equation typically reduces to solving linear systems, and constructing effective preconditioners can accelerate iterative solutions. In recent years, researchers have proposed various preconditioning methods based on the structure and properties of coefficient matrices arising from different discretization schemes, significantly improving solution efficiency. Addressing the numerical solution of space-fractional Cahn-Hilliard equations, this paper reviews and analyzes corresponding preconditioning methods for one-dimensional and higher-dimensional space-fractional Cahn-Hilliard equations under different discretization scenarios, and provides recommendations for future research.

1. 问题背景

空间分数阶Cahn-Hilliard方程是经典Cahn-Hilliard方程的扩展，以空间分数阶导数代替了经典Cahn-Hilliard方程中的局部二阶拉普拉斯算子，通过其非局部积分特性刻画长程空间相互作用，描述反常扩散与复杂界面动力学。空间分数阶Cahn-Hilliard方程被广泛应用于生物学[1]、图像修复[2]-[4]和材料科学[5]等领域。

1.1. 空间分数阶Cahn-Hilliard方程

对于常系数双侧空间分数阶Cahn-Hilliard方程的初边值问题：

$\{\begin{array}{ll}\frac{\partial u\left(x,t\right)}{\partial t}={ℒ}_x^{\alpha }v\left(x,t\right),\hfill & \left(x,t\right)\in \Omega \times \left(0,T\right],\hfill \\ v\left(x,t\right)=-{ϵ}^2{ℒ}_x^{\alpha }u\left(x,t\right)+f\left(u\left(x,t\right)\right),\hfill & \left(x,t\right)\in \Omega \times \left(0,T\right],\hfill \\ u\left(x,t\right)=v\left(x,t\right)=0,\hfill & \left(x,t\right)\in \partial \Omega \times \left(0,T\right],\hfill \\ u\left(x,0\right)=u_0\left(x\right),\hfill & x\in \overline{\Omega },\hfill \end{array}$ (1)

其中$\Omega ={\displaystyle {\prod }_{i=1}^d\left(a_i,b_i\right)}\subset {ℝ}^d$ 为空间域，$\overline{\Omega }=\Omega \cup \partial \Omega$ ，$x=\left(x_1,x_2,\cdots ,x_d\right)\in {ℝ}^d$ ，$\alpha =\left({\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_d\right)$ 为分数且$1<{\alpha }_i<2$ ($i=1,\cdots ,d$ )，$u\left(x,t\right)$ 用于描述相场浓度，$f\left(u\right)=F^{\prime }\left(u\right)$ 是具有$F\left(u\right)=\frac{1}{4}{\left(u^2-1\right)}^2$ 的非线性双阱势，$ϵ$ 为表示界面宽度的正常数。${ℒ}_x^{\alpha }$ 通常表示为Riesz分数阶算子。

定义1 d维Riesz分数阶导数[6]：d维Riesz分数阶算子${ℒ}_x^{\alpha }$ 定义为

${ℒ}_x^{\alpha }u\left(x,t\right)={\displaystyle {\sum }_{i=1}^d{ℒ}_{x_i}^{{\alpha }_i}u\left(x,t\right)}={\displaystyle {\sum }_{i=1}^d{c}_{{\alpha }_i}}\left({}_{a_i}{D}_{x_i}^{{\alpha }_i}u\left(x,t\right)+{}_{x_i}{D}_{b_i}^{{\alpha }_i}u\left(x,t\right)\right)$ (2)

其中$c_{{\alpha }_i}=-\frac{1}{2\cos \left(\frac{\pi {\alpha }_i}{2}\right)}>0$ ，左、右Riemann-Liouville分数阶导数定义为

${}_{a_i}{D}_{x_i}^{{\alpha }_i}u\left(x,t\right)=\frac{1}{\Gamma \left(2-{\alpha }_i\right)}\frac{{\partial }^2}{\partial x_i^2}{\displaystyle {\int }_{a_i}^{x_i}\frac{u\left(x_1,x_2,\cdots ,x_{i-1},\eta ,x_{i+1},\cdots ,x_d,t\right)}{{\left(x_i-\eta \right)}^{{\alpha }_i-1}}\text{d}\eta },$

${}_{x_i}{D}_{b_i}^{{\alpha }_i}u\left(x,t\right)=\frac{1}{\Gamma \left(2-{\alpha }_i\right)}\frac{{\partial }^2}{\partial x_i^2}{\displaystyle {\int }_{x_i}^{b_i}\frac{u\left(x_1,x_2,\cdots ,x_{i-1},\eta ,x_{i+1},\cdots ,x_d,t\right)}{{\left(\eta -x_i\right)}^{{\alpha }_i-1}}\text{d}\eta }.$

1.2. 预处理

分数阶Cahn-Hilliard方程通常包含非线性项与非局部分数阶导数，这使得几乎不可能求解到严格的解析解，需要通过离散化过程构造近似的数值格式，在离散点上逼近原方程的解，从而获得近似的数值解。常见的数值算法：有限差分法[7] [8]、谱方法[9] [10]、有限元法[11]被广泛应用于众多文献中。

在得到空间分数阶Cahn-Hilliard方程的全离散格式后，可以将问题转化为求解一个线性方程组。由于空间分数阶Cahn-Hilliard方程将两个非局部积分–微分方程耦合在一起，所产生的线性系统的规模通常是大而稠密的，因此我们一般选择使用Krylov子空间法进行迭代求解[12]-[15]。然而，Krylov子空间法的收敛速度会依赖于线性方程组的系数矩阵的谱性质，特征值分布越聚集，Krylov子空间法的收敛速度越快，反之，收敛速度越慢甚至发散，因此对线性方程组进行预处理是关键。

不同结构的线性方程组取决于对原始空间分数阶Cahn-Hilliard方程采用不同的离散化方法，而不同的预处理方法取决于对方程组的构造与系数矩阵的近似，因此，学者们通常会基于此来构造不同的预处理器，从而提高Krylov子空间法的迭代求解效率[16]。针对数值求解一维和高维空间分数阶Cahn-Hilliard方程中产生的结构化方程组，本文介绍了近年来提出的不同预处理方法，并进行了对比分析与总结。

本文的具体结构如下：在第2部分，分别介绍了一维和高维空间分数阶Cahn-Hilliard方程在不同离散格式下对应的结构化方程组。在第3部分，整理并对比分析了关于这些方程组的预处理方法。最后，本文总结了数值求解一维和高维空间分数阶Cahn-Hilliard方程的预处理方法，并为进一步研究提供了建议。

2. 结构化方程组

2.1. 一维情况

2.1.1. 对称不定块系统

对于一维常系数双侧空间分数阶Cahn-Hilliard方程，结合SAV保结构技术，对时间导数采用二阶蛙跳格式进行离散化，对空间导数采用有限差分法进行离散化，数值求解框架中包含求解以下方程组[7]：

$\left[\begin{array}{cc}{T}_{\beta }& I\\ I& T_{\alpha }\end{array}\right]\left[\begin{array}{c}{u}^{k+1}\\ v^{k+1}\end{array}\right]=\left[\begin{array}{c}{y}_1^k\\ y_2^k\end{array}\right]$ (3)

其中$I$ 为单位矩阵，$T_{\alpha }$ 为对称正定的Toeplitz矩阵，$T_{\beta }$ 为对称负定的Toeplitz矩阵，从而使得整体系数矩阵在性质上是一个病态且不定的矩阵，在结构上是一个对称的块矩阵。对于不定系统，通常使用MINRES方法来处理。系数矩阵的病态性说明其条件数极大，从而会导致迭代收敛缓慢，且由于存在着稠密的Toeplitz矩阵，会导致求解方程组具有较高的计算复杂度。

2.1.2. 单位矩阵与“Toeplitz-Like”结构矩阵相加的系统

对于一维常系数双侧空间分数阶Cahn-Hilliard方程，结合稳定化技术[8]，对时间导数采用Crank-Nicolson格式进行离散化，对空间导数采用有限差分法进行离散化，数值求解框架中包含求解以下方程组[17]：

$M{u}^{k+1}=b^k$ (4)

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

$I+T$ (5)

$I$ 为单位矩阵，$T$ 是由Toeplitz矩阵通过加法、乘法等运算组合而成的矩阵，即具有“Toeplitz-like”结构的矩阵。若系数矩阵$M$ 非对称，我们通常会选择将其对称化以便于后续使用共轭梯度法(CG)来高效求解，因此如何构造对称正定系统和解决Toeplitz矩阵的稠密性以降低计算复杂度成为了关键问题。

2.2. 高维情况

对称不定块系统：对于高维常系数双侧空间分数阶Cahn-Hilliard方程，采用SAV保结构技术，对时间导数采用二阶蛙跳格式进行离散化，对空间导数采用有限差分法进行离散化，数值求解框架中包含求解以下方程组[6]：

$\left[\begin{array}{cc}-G_{\alpha }& I\\ I& G_{\alpha }\end{array}\right]\left[\begin{array}{c}{u}^{k+1}\\ v^{k+1}\end{array}\right]=\left[\begin{array}{c}{y}_1^k\\ y_2^k\end{array}\right]$ (6)

其中$I$ 为单位矩阵，$G_{\alpha }$ 是由一维Toeplitz矩阵通过Kronecker积组合而成的对称正定块矩阵，使得整体系数矩阵在性质上是一个病态且不定的矩阵，在结构上是一个对称的块矩阵。其病态性与Toeplitz矩阵的稠密性会降低MINRES方法的迭代求解效率。

3. 预处理方法

针对上述数值求解一维和高维空间分数阶Cahn-Hilliard方程中不同特点的结构化方程组，本节分别整理了相应的预处理方法，并从不同方面进行了对比分析。

3.1. 针对一维情况的预处理

3.1.1. 针对对称不定块系统的预处理形式及性质

关于方程组(3)，Huang等人[7]基于离散正弦变换与舒尔补精确分裂构造了一个对称正定预处理器。针对线性系统的系数矩阵

$A=\left[\begin{array}{cc}{ϵ}^2{A}_{\beta }& I\\ I& -\Delta t{A}_{\alpha }\end{array}\right],$

其中$A_{\gamma }$ ($\gamma =\alpha$ 或$\beta$ )是将空间分数阶导数离散化所产生的对称负定的Toeplitz矩阵，通过基于离散正弦变换的$\tau$ 矩阵[18]来近似$A_{\gamma }$ ，使得求其逆只需要$O\left(M\log M\right)$ 的计算复杂度，定义矩阵

$P_1=\left[\begin{array}{cc}{ϵ}^2\tau \left(A_{\beta }\right)& I\\ I& -\Delta t\tau \left(A_{\alpha }\right)\end{array}\right].$

由于对于任何$\gamma$ ，$\tau \left(A_{\gamma }\right)$ 都是负定的，因此矩阵$P_1$ 是不定的，还不能直接用作预处理器。为了找寻一个正定预处理器，通过正弦变换将$P_1$ 对角化，记${\Lambda }_1$ 和${\Lambda }_2$ 为包含${ϵ}^2\tau \left(A_{\beta }\right)$ 和$-\Delta t\tau \left(A_{\alpha }\right)$ 所有特征值的对角矩阵，S为正弦变换矩阵，即

${ϵ}^2\tau \left(A_{\beta }\right)=S{\Lambda }_{1}S,-\Delta t\tau \left(A_{\alpha }\right)=S{\Lambda }_{2}S.$

于是可以得到$P_1$ 的分裂形式

$P_1=\left[\begin{array}{cc}S& 0\\ 0& S\end{array}\right]\left[\begin{array}{cc}{\Lambda }_1& I\\ I& {\Lambda }_2\end{array}\right]\left[\begin{array}{cc}S& 0\\ 0& S\end{array}\right],$

其中0是与S大小相同的零矩阵。再利用舒尔补，可以获得块对角矩阵的一种分解

$\left[\begin{array}{cc}{\Lambda }_1& I\\ I& {\Lambda }_2\end{array}\right]=\left[\begin{array}{cc}I& 0\\ {\Lambda }_1^{-1}& I\end{array}\right]\left[\begin{array}{cc}{\Lambda }_1& 0\\ 0& -{\Lambda }_1^{-1}+{\Lambda }_2\end{array}\right]\left[\begin{array}{cc}I& {\Lambda }_1^{-1}\\ 0& I\end{array}\right].$

经验证${\Lambda }_1$ 为负，${\Lambda }_2$ 为正，为得到一个正定的预处理器，以$-{\Lambda }_1$ 替换${\Lambda }_1$ ，记为矩阵P，被描述为如下对角化形式：

$P=\left[\begin{array}{cc}S& 0\\ 0& S\end{array}\right]\left[\begin{array}{cc}I& 0\\ {\Lambda }_1^{-1}& I\end{array}\right]\left[\begin{array}{cc}{\Lambda }_1& 0\\ 0& -{\Lambda }_1^{-1}+{\Lambda }_2\end{array}\right]\left[\begin{array}{cc}I& {\Lambda }_1^{-1}\\ 0& I\end{array}\right]\left[\begin{array}{cc}S& 0\\ 0& S\end{array}\right]$ (7)

我们就得到了所需的对称正定块预处理器。

3.1.2. 针对对称不定块系统的预处理分析

通过分析，针对求解对称不定块系统一般以构造一个对称正定且容易求逆的预处理器为目标，加快MINRES方法的迭代求解效率。对于处理稠密的Toeplitz矩阵，通常选择使用基于离散正弦变换的τ近似来提升算法效率。

Huang等人[7]基于直接近似的思想，通过正弦变换对进行完τ近似的矩阵对角化并写成其分裂形式，使其正、负特征值的分布显露出来，在含有负特征值的对角矩阵前加负号，以实现矩阵的正定性。再利用舒尔补，将中间块矩阵进行精确分裂，所得块三角矩阵与块对角矩阵乘积的结构便于求逆。该预处理方法虽然通过数值实验证明了其有效性，但缺乏对于预处理矩阵的谱分析。其预处理方法还可以尝试使用循环矩阵近似代替τ矩阵近似，在某些情况下，特别是当矩阵近似平移不变时，循环近似的效果可能会更好，其可通过FFT实现快速对角化。

3.1.3. 针对单位矩阵与“Toeplitz-Like”结构矩阵相加系统的预处理形式及性质

关于方程组(4)，Huang等人[17]针对一个非对称的线性系统，分别基于合同变换和右乘逆变换提出了两种对称化方法。在得到等价系统后，采用基于离散正弦变换的τ近似，获得最终预处理器。

原非对称系统的系数矩阵为$M=I+\widehat{A}\widehat{B}$ ，其中

$\widehat{A}=-c_{0}A,c_0=\frac{\Delta t{c}_{\alpha }}{h^{\alpha }},$

A是通过离散Riesz分数阶导数所得到的对称负定的Toeplitz矩阵，

$\widehat{B}=\Delta tsI-\frac{{ϵ}^2}{2h^2}B$ ，$s>0$ 是稳定系数，

B是通过离散拉普拉斯算子得到的对称负定的三对角Toeplitz矩阵。$\widehat{A}$ 和$\widehat{B}$ 都是对称正定的Toeplitz矩阵。

第一种对称化方法为：由于$\widehat{B}$ 是一个自然的τ矩阵，通过在原始系统两边左乘矩阵${\widehat{B}}^{1/2}$ 可以得到一个对称Toeplitz系统，该系统的系数矩阵为

$\tilde{M}=I+{\widehat{B}}^{1/2}\widehat{A}{\widehat{B}}^{1/2}.$

对$\widehat{A}$ 进行τ近似得到一个对称正定且可逆的预处理器：

$\tilde{P}=I+{\widehat{B}}^{1/2}\tau \left(\widehat{A}\right){\widehat{B}}^{1/2}$ (8)

第二种对称化方法为：将原始系数矩阵右乘${\widehat{B}}^{-1}$ ，得到一个对称系统，其系数矩阵

$\overline{M}={\widehat{B}}^{-1}+\widehat{A}.$

同样地，对$\widehat{A}$ 进行τ近似得到一个对称正定且可逆的预处理器：

$\overline{P}={\widehat{B}}^{-1}+\tau \left(\widehat{A}\right)$ (9)

预处理矩阵${\tilde{P}}^{ }$ 和${\overline{P}}^{ }$ 具有以下谱性质。

定理1 ${\tilde{P}}^{-1}\tilde{M}$ 的特征值满足以下不等式

$1/2<\lambda \left({\tilde{P}}^{-1}\tilde{M}\right)<3/2$

定理2 ${\overline{P}}^{-1}\overline{M}$ 的特征值满足以下不等式

$1/2<\lambda \left({\overline{P}}^{-1}\overline{M}\right)<3/2$

3.1.4. 针对单位矩阵与“Toeplitz-Like”结构矩阵相加系统的预处理对比与分析

通过对比分析，Huang等人[17]的两种预处理方法都以将线性系统对称化以获得一个对称正定的系统为目标，再基于τ近似构造预处理器加速共轭梯度法(CG)的求解效率。从系统对称性的构造方式来看，第一种方法是利用${\widehat{B}}^{1/2}$ 的对称正定性，以及${\widehat{B}}^{1/2}\widehat{A}{\widehat{B}}^{1/2}$ 是$\widehat{A}$ 的合同变换，保证了$\tilde{M}$ 的对称正定性；第二种方法利用${\widehat{B}}^{-1}$ 和$\widehat{A}$ 各自的对称正定性，通过矩阵加法保证了$\overline{M}$ 的对称正定性。从计算复杂度来看，带有预处理器$\tilde{P}$ 的PCG方法需要反复处理${\widehat{B}}^{1/2}$ ，这导致在应用中调用多次DST，且初始化时还需计算$\widehat{B}$ 特征值的平方根。而有预处理器$\overline{P}$ 的PCG方法只涉及${\widehat{B}}^{-1}$ ，其结构允许预处理器应用和矩阵乘法均通过DST实现，且调用次数更少，完全避开了开平方运算。数值实验也表明第二种方法的CPU更短，因此第二种对称化方法构造的预处理器计算量更小，整体效果更好。对于第一种方法${\widehat{B}}^{1/2}$ 所带来的计算量，为避免显示构造平方根，我们建议可以在原系统两端左乘$L^{\text{T}}$ ，并利用$\widehat{B}$ 的Cholesky分解$\widehat{B}=L{L}^{\text{T}}$，其中$L$ 是一个稀疏的下三角双对角矩阵，记$y=L^{\text{T}}u$ ，所得对称正定系统为$\left(I+L^{\text{T}}\widehat{A}L\right)y=L^{\text{T}}b$，仍然基于τ近似构造预处理器$\widehat{P}=I+L^{\text{T}}\tau \left(\widehat{A}\right)L$。由于$L$ 稀疏且$\tau \left(\widehat{A}\right)$ 可快速对角化，整体仍能保持接近O (MlogM)的计算复杂度，与原方法相比，减少了DST的调用次数。

3.2. 针对高维情况的预处理

3.2.1. 针对对称不定块系统的预处理形式及性质

关于方程组(6)，Huang等人[6]基于线性近似与离散正弦变换构造了一个对称正定块预处理器。先将线性方程组进行等价变换，方便起见，我们将新系统简化为

$A{\tilde{u}}^{k+1}={\tilde{y}}^k$ (10)

其中

$A=\left[\begin{array}{cc}-cG& I\\ I& cG\end{array}\right],$

G是将多维空间分数阶导数离散化所产生的对称正定的块Toeplitz矩阵，记$\theta =\sqrt{c_1/c_2}$ ，$c=\sqrt{c_1{c}_2}$ ，$c_1={ϵ}^2$ ，$c_2=\Delta t$ 。为找到一个正定的预处理器，首先将系数矩阵$A$ 进行平方运算，得到

$A^2=\left[\begin{array}{cc}{c}^2{G}^2+I& 0\\ 0& c^2{G}^2+I\end{array}\right],$

其中0是与$G$ 大小相同的零矩阵。再对$A^2$ 进行开根操作，得到一个理论上与系数矩阵最接近的理想预处理器

$\overline{P}=\left[\begin{array}{cc}\sqrt{c^2{G}^2+I}& 0\\ 0& \sqrt{c^2{G}^2+I}\end{array}\right],$

理论分析表明${\overline{P}}^{-1}A$ 的特征值为±1。然而，$\overline{P}$ 中存在平方根运算难以处理，因此Huang利用线性化近似来代替平方根运算，得到中间近似预处理器

$\widehat{P}=\left[\begin{array}{cc}cG+I& 0\\ 0& cG+I\end{array}\right],$

通过理论分析可知${\widehat{P}}^{-1}\overline{P}$ 的所有特征值位于$\left[\sqrt{2}/2,1\right]$ ，这说明$\widehat{P}$ 与$\overline{P}$ 是谱等价的。但直接处理稠密对称Toeplitz矩阵$G$ 的运算(矩阵–向量乘法)的复杂度高达$O\left(M_i^2\right)$ ，于是Huang采用基于离散正弦变换的$\tau$ 矩阵来近似$G$ ，因此其矩阵–向量乘法可在$O\left(M_i\log M_i\right)$ 次运算中得到，最终预处理器为

$P=\left[\begin{array}{cc}c\tau \left(G\right)+I& 0\\ 0& c\tau \left(G\right)+I\end{array}\right]$ (11)

预处理器$P$ 可通过离散正弦变换求逆。记$S=S_{M_1}\otimes S_{M_2}\otimes \cdots \otimes S_{M_d}$ ，从而得到

$P=\left[\begin{array}{cc}S& 0\\ 0& S\end{array}\right]\left[\begin{array}{cc}{\displaystyle {\sum }_{i=1}^d\left(I_{M_i^{-}}\otimes \left(c{\Lambda }_{{\alpha }_i}+I\right)\otimes I_{M_i^{+}}\right)}& 0\\ 0& {\displaystyle {\sum }_{i=1}^d\left(I_{M_i^{-}}\otimes \left(c{\Lambda }_{{\alpha }_i}+I\right)\otimes I_{M_i^{+}}\right)}\end{array}\right]\left[\begin{array}{cc}S& 0\\ 0& S\end{array}\right]$

其中${\Lambda }_{{\alpha }_i}$ 为每个$\tau \left(G_{{\alpha }_i}\right)$ 的特征值。对于任何向量$v\in {ℝ}^{2M\times 1}$ ，将其拆分成${\left[v_1,v_2\right]}^{\text{T}}\in {ℝ}^{2M\times 1}$ ，其中$v_1,v_2\in {ℝ}^{M\times 1}$ 。计算$\tilde{y}=S{v}_1$ ，$\tilde{z}=S{v}_2$ ，该步骤可通过$d$ 维DST以$O\left(M\log M\right)$ 次运算实现。计算

$\overline{y}={\displaystyle {\sum }_{i=1}^d\left(I_{M_i^{-}}\otimes {\left(c{\Lambda }_{{\alpha }_i}+I\right)}^{-1}\otimes I_{M_i^{+}}\right)\tilde{y}}$ ，$\overline{z}={\displaystyle {\sum }_{i=1}^d\left(I_{M_i^{-}}\otimes {\left(c{\Lambda }_{{\alpha }_i}+I\right)}^{-1}\otimes I_{M_i^{+}}\right)\tilde{z}}$ ，需要$O\left(M\right)$ 次运算的线性复杂度。计算$y=S\overline{y}$ ，$z=S\overline{z}$ ，同第一步一样可通过$d$ 维DST以$O\left(M\log M\right)$ 次运算实现。最后输出${\left[y;z\right]}^{\text{T}}$ 。因此对于计算$P^{-1}v$ 需要$O\left(M\log M\right)$ 复杂度。

预处理矩阵$P$ 具有以下谱性质。

定理3 $P^{-1}A$ 的特征值满足

$\sum \left(P^{-1}A\right)\subset \left[-\frac{3}{2},-\frac{1}{2\sqrt{2}}\right]\cup \left[\frac{1}{2\sqrt{2},},\frac{3}{2}\right]$

3.2.2. 针对对称不定块系统的预处理分析

Huang等人[6]基于直接近似的思想，采用对系数矩阵取平方再开根的方式，将原系数矩阵的不定性转化为正定性。利用线性近似、τ近似以及矩阵的块对角结构降低了计算复杂度。然而对中间预处理器的块对角部分采用了线性近似的操作，虽然从结构上对根号部分进行了简化，但是却引入了线性误差，这种近似方式比较粗糙，再加上后续对Toeplitz矩阵进行τ近似，还会放大近似误差。对此，可以采用有理逼近(如Chebyshev-Padé逼近或Zolotarev最优有理逼近)来近似平方根部分，会比一阶线性近似更准确地拟合平方根运算。

对于Chebyshev-Padé逼近，设$\tau \left(G\right)$ 的特征值为${\lambda }_i\in \left[{\lambda }_{\min },{\lambda }_{\max }\right]$ ，且$\tau \left(G\right)=S\Lambda S^{\text{T}}$ ，其中$\Lambda =diag\left({\lambda }_i\right)$ ，S为正弦变换矩阵。构造有理函数$r_{m,n}\left(\lambda \right)\approx f\left(\lambda \right)=\sqrt{c^2{\lambda }^2+1}$ ，${\lambda }_{ }\in \left[{\lambda }_{\min },{\lambda }_{\max }\right]$ ，其中$r_{m,n}\left(\lambda \right)=P_m\left(\lambda \right)/Q_n\left(\lambda \right)$ 由Chebyshev级数截断并配凑前$m+n+1$ 项系数得到。则预处理器为$r_{m,n}\left(\tau \left(G\right)\right)$

$P=\left[\begin{array}{cc}{r}_{m,n}\left(\tau \left(G\right)\right)& 0\\ 0& r_{m,n}\left(\tau \left(G\right)\right)\end{array}\right].$

其中$r_{m,n}\left(\tau \left(G\right)\right)=S{r}_{m,n}\left(\Lambda \right)S^{\text{T}}$ 。若$f$ 解析，则误差阶

$\underset{{\lambda }_{ }\in \left[{\lambda }_{\min },{\lambda }_{\max }\right]}{\max }\left|f\left(\lambda \right)-r_{m,n}\left(\lambda \right)\right|=O\left({\left(\frac{{\lambda }_{\max }-{\lambda }_{\min }}{{\lambda }_{\max }+{\lambda }_{\min }}\right)}^{m+n+1}\right).$

构造时需计算${\lambda }_i$ ($O\left(M\log M\right)$ )和有理函数值($O\left(M\right)$ )。每步MINRES迭代应用$P_{CP}^{-1}$ 需两次DST ($O\left(M\log M\right)$ )和一次对角缩放($O\left(M\right)$ )，使得总复杂度为$O\left(M\log M\right)$ 。

对于Zolotarev最优有理逼近，令

$\kappa =\sqrt{\frac{c^2{\lambda }_{\max }^2+1}{c^2{\lambda }_{\min }^2+1}}>1,$

作变量替换

$x=\frac{c^2{\lambda }_{\max }^2+1}{c^2{\lambda }_{\min }^2+1}\in \left[1,{\kappa }^2\right],$

Zolotarev给出了次数为(n, n)的最优有理逼近$R_n\left(x\right)\approx \sqrt{x}$ 的显示形式：

$R_n\left(x\right)=\sqrt{\kappa }{\displaystyle \prod }_{j=1}^n\frac{x+c_j}{x+d_j},$

其中

$c_j=\kappa s{n}^2\left(\frac{2j-1}{2n}{K}^{\prime },k\right),$

$d_j=\kappa /c_j$ ，$s{n}^{ }$ 为Jacobi椭圆函数，$k=\sqrt{1-1/{\kappa }^2}$ ，$K^{\prime }=K\left(k^{\prime }\right)$ 。则$f\left(\lambda \right)=\sqrt{c^2{\lambda }^2+1}$ 的逼近为：

$f_{Zolo}\left(\lambda \right)=\sqrt{c^2{\lambda }_{\min }^2+1}{R}_n\left(\frac{c^2{\lambda }^2+1}{c^2{\lambda }_{\min }^2+1}\right),$

预处理器为

$P_z=\left[\begin{array}{cc}{f}_{Zolo}\left(\tau \left(G\right)\right)& 0\\ 0& f_{Zolo}\left(\tau \left(G\right)\right)\end{array}\right],f_{Zolo}\left(\tau \left(G\right)\right)=S{f}_{Zolo}\left(\Lambda \right)S^{\text{T}}.$

最大相对误差为

$\underset{\lambda \in \left[{\lambda }_{\min },{\lambda }_{\max }\right]}{\max }\frac{\left|f\left(\lambda \right)-f_{Zolo}\left(\lambda \right)\right|}{f\left(\lambda \right)}\le 2{\text{e}}^{-2n\log \left(\kappa +\sqrt{{\kappa }^2-1}\right)},$

即指数收敛。其计算复杂度与Chebyshev-Padé相同，也为$O\left(M\log M\right)$ 。

4. 总结

针对数值求解一维和高维空间分数阶Cahn-Hilliard方程，学者们围绕系数矩阵的结构和性质提出了一系列的预处理方法，有效提高了Krylov子空间法的迭代求解效率。本文整理了近几年出现的预处理方法，并从不同角度进行了对比分析。我们发现，一维空间分数阶Cahn-Hilliard方程在经过不同的离散化方法的作用下，会形成两类结构化方程组，一类为对称不定块系统，另一类为系数矩阵为单位矩阵与“Toeplitz-like”结构矩阵之和的系统；而对于高维空间分数阶Cahn-Hilliard方程目前的研究方向主要集中在对称不定块系统。对于一维和高维情况的对称不定块系统，我们最终都要构造一个对称正定且容易求逆的块预处理器以提高MINRES方法的求解效率，因此通常会基于其系数矩阵的块结构来构造谱等价的具有块对角或块三角结构的矩阵作为预处理器，在构造过程中典型的方法有基于正弦变换的对角化和舒尔补分裂。对于一维情况下系数矩阵为单位矩阵与“Toeplitz-like”结构矩阵之和的系统，若该系统是对称正定的，需要构造一个对称正定的预处理器以加速共轭梯度法(CG)的求解。若该系统非对称，则可尝试将其对称化，最常用的对称化方法是合同变换。结果表明，相比原有的方程组，以上方法预处理后的方程组的线性收敛性显著提升，Krylov子空间法的迭代求解效率也有所提高。

然而，目前已发表的预处理方法有限，且其他格式在高维情况下的推广及相关预处理技术仍有待探索，值得进一步的研究与改进。关于空间分数阶Cahn-Hilliard方程数值求解中结构化方程组的预处理方法研究，本文认为还可以从以下方面进一步探索：

一、构造计算复杂度更低的对称系统：采用合同变换对系统进行对称化容易引入矩阵平方根，增加计算复杂度，可以尝试使用Cholesky分解，避免以上问题，以降低计算复杂度。

二、构造更精确的对称Toeplitz近似矩阵：尝试采用循环近似代替基于离散正弦变换的τ近似，当矩阵近似平移不变时，循环近似的效果可能会更好，且还可通过FFT实现快速对角化。

基金项目

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

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

NOTES

^*通讯作者。