参数化泊松方程的模型降阶预处理

Model order reduction preconditioning of parameterized Poisson equations

摘要

本文针对数值求解参数化泊松方程提出了一种基于模型降阶思想的预处理共轭梯度（Preconditioned Conjugate Gradient，PCG）算法。为加快所得离散线性方程组的求解速度，本文首先基于模型降阶思想构造了一般形式的预处理矩阵，并证明其对称正定性。在算法的off-line阶段，本文基于少量参数对应的全模型解数据，利用PCG算法并结合本征正交分解（Proper Orthogonal Decomposition，POD）方法生成了一组动态预处理矩阵。然后，在on-line阶段，本文利用动态预处理矩阵结合PCG算法建立所需算法。为了验证算法的性能，本文分别在单位矩形区域和L形区域上数值求解参数化泊松方程。结果显示，在相同精度条件下，算法的平均计算时间比标准共轭梯度算法快40倍以上。

Abstract

This paper aims at the fast numerical solution of parameterized Poisson equation.A preconditioned conjugate gradient （PCG） method based on the model order reduction method is proposed to speed up the solution of the obtained discrete linear system.First， the general formulation of preconditioning matrix based on model order reduction is designed， and the symmetry and positive definition of the matrix are proved.In the off-line stage of the method， by using very few solution data， the PCG algorithm combined with the proper orthogonal decomposition （POD） method are adopted to generate a set of dynamic preconditioning matrices.In the on-line stage， the MPCG algorithm is proposed by combining the PCG algorithm and these dynamic preconditioning matrices.To verify the performance of the method， parameterized Poisson equations on the unit rectangular and L-shaped domains are numerically solved.It is shown that， in comparison with the standard CG algorithm， the average computation time is speeded up by more than 43 times with the same calculation accuracy.

Graphical abstract

关键词

参数化泊松方程 / 模型降阶 / 预处理矩阵 / 共轭梯度

Key words

parameterized Poisson equation / model order reduction / precondition matrix / conjugate gradient（2020 MSC 65M60）

作为描述众多基本物理现象（如静电场、热传导、流体流动等）的关键数学模型，参数化泊松方程通过引入参数提升了模型的灵活性和适应性，能够高效处理具有复杂几何形状及存在多物理场耦合的问题^［1-2］。

目前，泊松方程的数值求解已有多种高精度方法，如有限元法、有限体积法及有限差分法等。这些方法的求解代价很高，通常涉及高达10⁶~10⁹个自由度的数值计算问题，需花费数小时甚至更长的CPU时间^［3-4］。特别地，对于参数化泊松方程而言，每当方程中的参数发生改变时就需要重新进行求解，每次都需花费巨大的计算成本。因此，如何使算法在满足精度的条件下实现快速求解是一个有意义同时也极具挑战性的问题。

模型降阶方法因其高效的在线计算效率日渐成为解决以上问题的关键方法^［5-6］。该方法通过构造低维问题来逼近高维全模型（full-model）问题来有效降低计算复杂度^［6］。为了保证计算精度，模型降阶方法往往需要在离线阶段（off-line）生成大量全模型解数据以构造低维逼近空间，代价非常大。针对此问题，近年来研究者尝试将模型降阶与机器学习相结合，降低参数化PDE的计算复杂度。值得注意的是，这样虽然能够避免降阶模型需要侵入式修改全模型的缺点，但无法减少离线阶段的计算成本^［7］，并可能导致在线计算精度降低。因此，如何将模型降阶思想与全模型解有效结合起来加速其中最耗时的离散方程组的求解便成为问题的关键。

对于大型线性方程组，主流的求解方法是预处理Krylov子空间法^［8-9］，其中的预处理矩阵的构造是加速算法收敛的关键^［10］。特别地，针对参数化PDE，离散得到的大型线性方程组需根据方程中参数的构造去设计相应的预处理矩阵，加速收敛。为解决模型降阶中的线性系统序列问题，Anzt等^［11］利用前一个分解作为初始猜测值进行更新，显著提高了算法的迭代求解效率。Carlberg等^［12］提出了结合Krylov子空间循环利用和本征正交分解（Proper Orthogonal Decomposition，POD）增强的共轭梯度（Conjugate Gradient，CG）方法，用于高效求解参数化PDE离散产生的系列方程组。Santo等^［13］提出了一种多空间降基方法，用于构建预处理器的粗空间部分，加速其收敛。

针对参数化泊松方程参数变化时需要多次快速求解的需求，本文构造了基于模型降阶的预处理矩阵的一般形式，在off-line阶段利用POD结合预处理共轭梯度算法（Preconditioned Conjugate Gradient，PCG）生成一组动态预处理矩阵，在on-line阶段将动态预处理矩阵与PCG算法结合起来，进而建立了基于模型降阶预处理的MPCG算法。由于在off-line阶段仅需少量的全模型解数据，本文方法相对标准模型降阶方法大幅减少了off-line阶段的计算量。

1 参数化泊松方程及其数值求解

1.1 方程及限元解

参数化泊松方程具有如下形式：

\left\{\begin{array}{l}-\mathrm{\Delta }\mathit{\phi }+\mu \mathit{\phi }=\mathit{f},\mathrm{i}\mathrm{n}\mathrm{\Omega },\\ \mathit{\phi }=\mathit{g},\mathrm{o}\mathrm{n}\partial \mathrm{\Omega }\end{array}\right.（1）

其中，\mathit{\phi }是待求未知函数，\mu \ge \mathrm{0}是参数，\mathit{f}是给定的源项，\mathit{g}是已知边界条件，区域\mathrm{\Omega }是单连通、有界的多边形开集。

有限元法是求解参数化泊松方程的常见方法^［3-4］，其具体求解步骤如下。

1） 将参数化泊松方程转换为变分形式，即求\mathit{\phi }\in U，满足：

a(\mathit{\phi },\mathit{v})=〈\mathit{f},\mathit{v}〉,\forall \mathit{v}\in V,

其中，

a(\mathit{\phi },\mathit{v})=(\nabla \mathit{\phi },\nabla \mathit{v})+\mu (\mathit{\phi },\mathit{v}),

U=H_g^{\mathrm{1}}(\mathrm{\Omega })=\{\mathit{\phi }\in H_{}^{\mathrm{1}}(\mathrm{\Omega }):\mathit{\phi }{|}_{\partial \mathrm{\Omega }}=\mathit{g}\},

V=H_{\mathrm{0}}^{\mathrm{1}}(\mathrm{\Omega })=\{\mathit{v}\in H^{\mathrm{1}}(\mathrm{\Omega }):\mathit{v}{|}_{\partial \mathrm{\Omega }}=\mathrm{0}\}。

2） 将求解区域划分为有限个单元，并构造有限元空间U_h和V_h。离散的变分形式为：求{\mathit{\phi }}_h\in U_h，使得：

a({\mathit{\phi }}_h,{\mathit{v}}_h)=〈{\mathit{f}}_h,{\mathit{v}}_h〉,\forall {\mathit{v}}_h\in V_h。

3） 组装刚度矩阵并施加边界条件，得到离散线性方程组：

{\mathit{A}}_h(\mu ){\mathit{\phi }}_h(\mu )={\mathit{f}}_h(\mu )（2）

其中，系数矩阵{\mathit{A}}_h(\mu )对称正定，方程组随参数\mu的改变而变化。

在参数化泊松方程求解过程中，对方程组（2）的求解最耗时。

1.2 PCG算法

CG是求解对称正定线性方程组的主流方法之一。它的基本思想是：构造一组共轭方向，然后沿着这些方向进行线性搜索，逐步逼近方程的解^［8］。PCG是CG算法的改进，通过引入预处理矩阵来改善系数矩阵的条件数，加速收敛。

对方程组（2），通过左预处理共轭梯度方法先将其转化为：

{\mathit{M}}^{-\mathrm{1}}{\mathit{A}}_h(\mu ){\mathit{\phi }}_h(\mu )={\mathit{M}}^{-\mathrm{1}}{\mathit{f}}_h(\mu )（3）

再对方程组（3）使用标准的CG算法求解，其中\mathit{M}是对称正定的预处理矩阵，即矩阵{\mathit{A}}_h(\mu )的近似。在每次迭代中，PCG算法需要计算{\mathit{M}}^{-\mathrm{1}}与残差向量的乘积。实际计算中通常将该步骤转化为求解如下方程组：

\mathit{M}\mathit{z}=\mathit{r}（4）

由于预处理矩阵\mathit{M}是{\mathbf{A}}_h(\mu )的近似，所以求解方程组（4）相当于近似求解：

{\mathit{A}}_h(\mu )\mathit{z}=\mathit{r}（5）

因此，预处理过程并不需要显式地给出预处理矩阵\mathit{M}。

当参数\mu发生改变时，方程组（2）中的系数矩阵{\mathit{A}}_h(\mu )相应改变。此时，传统预处理方法的预处理矩阵需要重新设计，难以实现高效求解。

1.3 POD模型降阶

基于POD的模型降阶方法是一种基于数据驱动的降阶技术，核心思想是通过捕捉系统的动态响应的主要特征来构造低维空间，逼近高维全模型的解、显著降低计算复杂度。该方法可以分为off-line阶段和on-line阶段两部分^［6］。

使用基于POD的模型降阶方法求解参数化泊松方程时，在off-line阶段首先通过有限元法获得方程（1）中n个参数对应的全模型解，形成快照矩阵\mathit{S}\in {\mathbb{R}}^{N_h\times n}，即：

\mathit{S}=\left(\begin{array}{cccc}{\mathit{\phi }}_h({\mu }_{\mathrm{1}}),& {\mathit{\phi }}_h({\mu }_{\mathrm{2}}),& \cdots ,& {\mathit{\phi }}_h({\mu }_n)\end{array}\right),

其中，{\mathit{\phi }}_h({\mu }_i)为参数{\mu }_i对应的全模型解。其次，对快照矩阵\mathit{S}进行奇异值分解，得到奇异值和一组N_h维正交的左奇异向量。然后，按照奇异值从大到小的顺序，选择前N个（N\le n）奇异值对应的左奇异向量，构成基矩阵，即POD基\mathit{Z}。本文将该过程记为\mathit{Z}=\mathrm{P}\mathrm{O}\mathrm{D}(\mathit{S},N)。最后，将\mathit{Z}作为低维空间的基向量构造降阶模型，以逼近方程（1）的全模型解。

在on-line阶段，对给定的新参数\mu，降阶模型只需要在基\mathit{Z}张成的低维逼近空间中进行求解，即：

{\mathit{\phi }}_h(\mu )\approx \mathit{Z}{\mathit{c}}_h(\mu )（6）

其中，{\mathit{c}}_h(\mu )为表示系数向量。将式（6）代入方程组（2），得：

{\mathit{A}}_h(\mu )\mathit{Z}{\mathit{c}}_h(\mu )\approx {\mathit{f}}_h(\mu )（7）

在方程组（7）中，{\mathit{A}}_h(\mu )和\mathbf{Z}的维数分别为N_h\times N_h和N_h\times N，且通常N_h\gg N，因而问题（7）为超定方程组。模型降阶方法一般在两侧同乘{\mathit{Z}}^{\mathrm{T}}，将问题（7）变成如下的适定方程组：

{\mathit{Z}}^{\mathrm{T}}{\mathit{A}}_h(\mu )\mathit{Z}\mathit{c}(\mu )\approx {\mathit{Z}}^{\mathrm{T}}{\mathit{f}}_h(\mu )（8）

求解方程组（8）得到表示系数向量，代入式（6），得：

{\mathit{\phi }}_h(\mu )\approx \mathit{B}(\mu ){\mathit{f}}_h(\mu )（9）

其中，\mathit{B}(\mu )=\mathit{Z}({\mathit{Z}}^{\mathrm{T}}{\mathit{A}}_h{(\mu )\mathit{Z})}^{-\mathrm{1}}{\mathit{Z}}^{\mathrm{T}}。

为保证降阶模型解的精度，需要基\mathit{Z}张成的空间能够足够逼近全模型解，这就要求在off-line阶段获得大量全模型解，并进行奇异值分解或设计复杂的自适应优选采样参数算法，两种方法都会导致off-line阶段的计算代价过大^［5-6］。

2 模型降阶预处理

2.1 模型降阶预处理矩阵

本文主要研究如何将模型降阶思想与PCG算法结合起来构造求解问题（1）的高效数值解法。为此，首先需要构造基于模型降阶的预处理矩阵。

方程组（2）的准确解为：

{\mathit{\phi }}_h(\mu )={\mathit{A}}_h^{-\mathrm{1}}(\mu ){\mathit{f}}_h(\mu )（10）

对比式（9）和式（10）不难发现，模型降阶其实是将\mathit{B}(\mu )作为{\mathit{A}}_h^{-\mathrm{1}}(\mu )的近似。我们并不能直接将{\mathit{B}}^{-\mathrm{1}}(\mu )作为PCG算法求解问题（2）时的预处理矩阵\mathit{M}，因为它要求预处理矩阵对称正定，而矩阵\mathit{B}(\mu )并不能保证这一点。为此，我们需要对\mathit{B}(\mu )进行修正。定义：

\mathit{M}(\mu )={\left(\mathit{B}(\mu )+\lambda \mathit{I}\right)}^{-\mathrm{1}}（11）

关于矩阵\mathit{M}(\mu )，我们有如下定理。

定理2.1 若矩阵{\mathit{A}}_h(\mu )对称正定，则对任意实数\lambda \mathrm{0}矩阵\mathit{M}(\mu )对称正定。

证明 因为{\mathit{A}}_h(\mu )对称正定，所以\mathit{Z}列满秩。对任意非零向量\mathit{y}，有

\mathit{y}{\mathit{Z}}^{\mathrm{T}}{\mathit{A}}_h(\mu )\mathit{Z}\mathit{y}={(\mathit{Z}\mathit{y})}^{\mathrm{T}}{\mathit{A}}_h(\mu )(\mathit{Z}\mathit{y})\mathrm{0}。

故{\mathit{Z}}^{\mathrm{T}}{\mathit{A}}_h(\mu )\mathit{Z}对称正定，其逆矩阵({\mathit{Z}}^{\mathrm{T}}{\mathit{A}}_h{(\mu )\mathit{Z})}^{-\mathrm{1}}也对称正定。因此，存在可逆矩阵\mathit{Q}，使得：

({\mathit{Z}}^{\mathrm{T}}{\mathit{A}}_h{(\mu )\mathit{Z})}^{-\mathrm{1}}={\mathit{Q}}^{\mathrm{T}}\mathit{Q}\mathrm{。}

则对任意非零向量\mathit{w}，成立：

\begin{array}{l}{\mathit{w}}^{\mathrm{T}}\mathit{B}(\mu )\mathit{w}={\mathit{w}}^{\mathrm{T}}\mathit{Z}({\mathit{Z}}^{\mathrm{T}}{\mathit{A}}_h{(\mu )\mathit{Z})}^{-\mathrm{1}}{\mathit{Z}}^{\mathrm{T}}\mathit{w}=\\ {\mathit{w}}^{\mathrm{T}}\mathit{Z}{\mathit{Q}}^{\mathrm{T}}\mathit{Q}{\mathit{Z}}^{\mathrm{T}}\mathit{w}=(\mathit{Q}{\mathit{Z}}^{\mathrm{T}}{\mathit{w})}^{\mathrm{T}}\mathit{Q}{\mathit{Z}}^{\mathrm{T}}\mathit{w}\ge \mathrm{0}\mathrm{。}\end{array}

故\mathit{B}(\mu )是对称半正定矩阵。由于\lambda \mathit{I}对称正定，所以\mathit{B}(\mu )+\lambda \mathit{I}也对称正定。从而矩阵\mathit{M}(\mu )对称正定。证毕。

根据定理2.1，恰当选择参数\lambda后，矩阵\mathit{M}(\mu )可以作为一个合适的预处理矩阵。需要注意的是，在\mathit{M}(\mu )中只有{\mathit{A}}_h(\mu )与模型参数\mu相关，其他均与\mu无关，从而可以在off-line阶段一次性准备好。在on-line阶段，只需计算({\mathit{Z}}^{\mathrm{T}}{\mathit{A}}_h{(\mu )\mathit{Z})}^{-\mathrm{1}}这个N阶方阵。通常N非常小，从而计算量是比较小的。

在PCG算法实现过程中，预处理矩阵只在求解方程组（4）时被使用。由于矩阵\mathit{B}(\mu )的维数为N_h\times N_h，在实现过程中并未直接形成\mathit{B}(\mu )，而是将\mathit{z}={\mathit{M}}^{-\mathrm{1}}(\mu )\mathit{r}做如下拆分：

\begin{array}{l}{\mathit{M}}^{-\mathrm{1}}(\mu )\mathit{r}=(\mathit{B}(\mu )+\lambda \mathit{I})\mathit{r}=\\ \mathit{B}(\mu )\mathit{r}+\lambda \mathit{r}=\\ \mathit{Z}[({\mathit{Z}}^{\mathrm{T}}{\mathit{A}}_h{(\mu )\mathit{Z})}^{-\mathrm{1}}({\mathit{Z}}^{\mathrm{T}}\mathit{r})]+\lambda \mathit{r}\end{array}（12）

按照（12）式进行计算，PCG算法运行过程中不需要显式形成\mathit{M}(\mu )，从而能够有效节约内存。算法的计算公式为：

(\mathit{\phi },\mathit{r},k)=\mathrm{P}\mathrm{C}\mathrm{G}(\mathit{A},\mathit{f},\mathit{Z},\lambda ,K,\epsilon ,{\mathit{\phi }}_{\mathrm{0}}),

具体实现方式如下。

算法1 PCG算法

输入： 系数矩阵\mathit{A}，右端向量\mathit{f}，POD基\mathit{Z}，预处理矩阵中参数\lambda，最大迭代步数K，容差\epsilon，初值{\mathit{\phi }}_{\mathrm{0}}；

输出： 解\mathit{\phi }，残差\mathit{r}，实际迭代步数k。

1） 初始化k=\mathrm{0}，残差{\mathit{r}}_{\mathrm{0}}=\mathit{f}-\mathit{A}{\mathit{\phi }}_{\mathrm{0}}，\mathit{C}={\mathit{Z}}^{\mathrm{T}}{(\mathbf{Z}\mathbf{A}\mathbf{Z})}^{-\mathrm{1}},{\mathit{z}}_{\mathrm{0}}=\mathit{C}({\mathit{Z}}^{\mathrm{T}}{\mathit{r}}_{\mathrm{0}}),{\mathit{p}}_{\mathrm{0}}={\mathit{z}}_{\mathrm{0}}；

2） for kK||{\mathit{r}}_k||/||\mathit{f}||\epsilon

3） a_k=({\mathit{r}}_k,{\mathit{z}}_k)/({\mathit{p}}_k,\mathit{A}{\mathit{p}}_k);

4） {\mathit{\phi }}_{k+\mathrm{1}}={\mathit{\phi }}_k+a_k{\mathit{p}}_k;

5） {\mathit{r}}_{k+\mathrm{1}}={\mathit{r}}_k-a_k\mathit{A}{\mathit{p}}_k;

6） {\mathit{z}}_{k+\mathrm{1}}=\mathit{C}({\mathit{Z}}^{\mathrm{T}}{\mathit{r}}_{k+\mathrm{1}})+\lambda {\mathit{r}}_{k+\mathrm{1}};

7） {\beta }_k=({\mathit{z}}_{k+\mathrm{1}},{\mathit{r}}_{k+\mathrm{1}})/({\mathit{z}}_k,{\mathit{r}}_k);

8） {\mathit{p}}_{k+\mathrm{1}}={\mathit{z}}_{k+\mathrm{1}}+{\beta }_k{\mathit{p}}_k;

9） k=k+\mathrm{1};

10） end。

2.2 动态预处理矩阵构造

给定模型参数\mu，若使用固定的\mathit{M}(\mu )作为PCG算法的预处理矩阵，则随迭代步数的增加，方程组（2）的相对残差将先快速下降，后缓慢下降，此时，若需满足精度要求，仍需花费较长时间。原因在于，PCG算法中每步迭代求解的方程组（4）应该是方程组（5）的近似求解，但由于右端项残差向量\mathit{r}随迭代不断发生改变，基于全模型解的数据构造的\mathit{Z}只能很好地逼近{\mathit{\phi }}_h(\mu )，而非方程组（5）的解。这样，当\mathit{f}与\mathit{r}差异变大时，就需要重新构造基矩阵\mathit{Z}，以保证预处理求解的方程组（4）始终是方程组（5）的近似求解。因此，在off-line阶段需要构造一组满足此需求的动态POD基，具体的构造方法见算法2，其中{\lambda }_i是第i个预处理矩阵需要的参数\lambda，p_i是第i个预处理矩阵使用的步数。需要注意的是，算法2在求解方程组{\mathit{A}}_h({\mu }_i){\mathit{\xi }}_i^{(k+\mathrm{1})}={\mathit{r}}_i^{(k+\mathrm{1})}时，并不直接求解，而是通过推导的方式得到，即：

\begin{array}{l}{\mathit{r}}_i^{(k+\mathrm{1})}={\mathit{A}}_h({\mu }_i){\mathit{\xi }}_i^{(k+\mathrm{1})}=\\ {\mathit{r}}_i^{(k)}-{\mathit{A}}_h({\mu }_i){\mathit{\phi }}_i^{(k+\mathrm{1})}=\\ {\mathit{A}}_h({\mu }_i){\mathit{\xi }}_i^{(k)}-{\mathit{A}}_h({\mu }_i){\mathit{\phi }}_i^{(k+\mathrm{1})}=\\ {\mathit{A}}_h({\mu }_i)({\mathit{\xi }}_i^{(k)}-{\mathit{\phi }}_i^{(k+\mathrm{1})}),\end{array}

因而，{\mathit{\xi }}_i^{(k)}-{\mathit{\phi }}_i^{(k+\mathrm{1})}是方程组的解。

算法2 动态POD的基构造算法

输入： n个参数对应的系数矩阵

\{{\mathit{A}}_h({\mu }_{\mathrm{1}}),{\mathit{A}}_h({\mu }_{\mathrm{2}}),\cdots ,{\mathit{A}}_h({\mu }_n)\}，

右端向量

\{{\mathit{f}}_h({\mu }_{\mathrm{1}}),{\mathit{f}}_h({\mu }_{\mathrm{2}}),\cdots ,{\mathit{f}}_h({\mu }_n)\}，

高精度数值解

\{{\mathit{\phi }}_h({\mu }_{\mathrm{1}}),{\mathit{\phi }}_h({\mu }_{\mathrm{2}}),\cdots ,{\mathit{\phi }}_h({\mu }_n)\}，

动态预处理矩阵数量m,POD基维数N,

参数\{{\lambda }_{\mathrm{0}},{\lambda }_{\mathrm{1}},\cdots ,{\lambda }_{m-\mathrm{1}}\},迭代步数\{p_{\mathrm{0}},p_{\mathrm{1}},\cdots ,p_{m-\mathrm{1}}\}；

输出：动态POD基\{{\mathit{Z}}_{\mathrm{0}},{\mathit{Z}}_{\mathrm{1}},\cdots ,{\mathit{Z}}_{m-\mathrm{1}}\}\mathrm{。}

1） 初始化k=\mathrm{0},初值{\mathit{\phi }}_{\mathrm{0}}=\mathrm{0},残差{\mathit{r}}_i^{(\mathrm{0})}={\mathit{f}}_h({\mu }_i),\mathrm{1}\le i\le n,{\mathit{\xi }}^{(\mathrm{0})}=\left(\begin{array}{cccc}{\mathit{\phi }}_h({\mu }_{\mathrm{1}}),& {\mathit{\phi }}_h({\mu }_{\mathrm{2}}),& \cdots ,& {\mathit{\phi }}_h({\mu }_n)\end{array}\right)；

2） for km-\mathrm{1}

3） 生成POD基{\mathit{Z}}_k=\mathrm{P}\mathrm{O}\mathrm{D}({\mathit{\xi }}^{(k)},N);

4） for i=\mathrm{1,2},\cdots ,n

5） 执行算法1，

\begin{array}{l}({\mathit{\phi }}_i^{(k+\mathrm{1})},{\mathit{r}}_i^{(k+\mathrm{1})},~)=\\ \begin{array}{cc}& \end{array}\mathrm{P}\mathrm{C}\mathrm{G}({\mathit{A}}_h({\mu }_i),{\mathit{r}}_i^{(k)},{\mathit{Z}}_k,{\lambda }_k,p_k,{\mathrm{e}}^{-\mathrm{20}},{\mathit{\phi }}_{\mathrm{0}});\end{array}

6） 生成方程组{\mathit{A}}_h({\mu }_i){\mathit{\xi }}_i^{(k+\mathrm{1})}={\mathit{r}}_i^{(k+\mathrm{1})}的解

{\mathit{\xi }}_i^{(k+\mathrm{1})}={\mathit{\xi }}_i^{(k)}-{\mathit{\phi }}_i^{(k+\mathrm{1})};

7） end

8） 令

{\mathit{\xi }}^{(k+\mathrm{1})}=\left(\begin{array}{cccc}{\mathit{\xi }}_{\mathrm{1}}^{(k+\mathrm{1})},& {\mathit{\xi }}_{\mathrm{2}}^{(k+\mathrm{1})},& \cdots ,& {\mathit{\xi }}_n^{(k+\mathrm{1})}\end{array}\right);

9） 更新k=k+\mathrm{1};

10） end。

2.3 MPCG算法

当输入新的模型参数\mu，再次对方程组（2）进行求解时，算法执行需要动态地选取POD基\{{\mathit{Z}}_{\mathrm{1}},{\mathit{Z}}_{\mathrm{2}},\cdots ,{\mathit{Z}}_{m-\mathrm{1}}\}。每次更新POD基后，由于预处理矩阵随之发生变化，需在更新预处理矩阵后重新调用PCG算法，以保证PCG算法的适用性，详见算法3，即基于模型降阶的动态预处理共轭梯度算法。

算法3 MPCG算法

输入： 系数矩阵{\mathit{A}}_h(\mu ),右端向量{\mathit{f}}_h(\mu ),最大迭代步数K,m个POD基\{{\mathit{Z}}_{\mathrm{0}},{\mathit{Z}}_{\mathrm{1}},\cdots ,{\mathit{Z}}_{m-\mathrm{1}}\}，每个POD基对应的迭代步数\{p_{\mathrm{0}},p_{\mathrm{1}},\cdots ,p_{m-\mathrm{1}}\},参数\{{\lambda }_{\mathrm{0}},{\lambda }_{\mathrm{1}},\cdots ,{\lambda }_m\}，容差\epsilon；

输出： 解{\mathit{\phi }}_h(\mu ),实际迭代步数k，残差{\mathit{r}}_h(\mu )\mathrm{。}

1） 初始化初值{\mathit{\phi }}_{\mathrm{0}}=\mathrm{0},迭代步数k=\mathrm{0},残差{\mathit{r}}_{\mathrm{0}}={\mathit{f}}_h(\mu )-{\mathit{A}}_h(\mu ){\mathit{\phi }}_{\mathrm{0}};

2） for j=\mathrm{1,2},\cdots ,m

3） 执行算法1，

\begin{array}{l}({\mathit{u}}_j,{\mathit{r}}_j,k_j)=\\ \begin{array}{cc}& \end{array}\mathrm{P}\mathrm{C}\mathrm{G}({\mathit{A}}_h(\mu ),{\mathit{f}}_h(\mu ),{\mathit{Z}}_{j-\mathrm{1}},{\lambda }_{j-\mathrm{1}},p_{j-\mathrm{1}},\epsilon ,{\mathit{u}}_{j-\mathrm{1}}),\end{array}

4） 更新迭代步数k=k+k_j;

5） 如果||{\mathit{r}}_j||/||{\mathit{f}}_h(\mu )||\epsilon或kK,则程序结束；

6） end

7） 如果||{\mathit{r}}_j||/||{\mathit{f}}_h(\mu )||\epsilon且k\le K

8） 执行算法1，

\begin{array}{l}({\mathit{u}}_{m+\mathrm{1}},{\mathit{r}}_{m+\mathrm{1}},k_{m+\mathrm{1}})=\\ \begin{array}{cc}& \mathrm{P}\mathrm{C}\mathrm{G}({\mathit{A}}_h(\mu ),{\mathit{f}}_h(\mu ),{\mathit{Z}}_{m-\mathrm{1}},{\lambda }_m,K-k,\epsilon ,{\mathit{u}}_m)\end{array};\end{array}

9） 更新迭代步数k=k+k_{m+\mathrm{1}}\mathrm{。}

在算法3中，前p_{\mathrm{0}}步使用POD基{\mathit{Z}}_{\mathrm{0}}和参数{\lambda }_{\mathrm{0}}构造预处理并进行迭代，紧接着在p_{\mathrm{1}}步使用POD基{\mathit{Z}}_{\mathrm{1}}和参数{\lambda }_{\mathrm{1}}构造预处理并进行迭代，以此类推。在更新POD基和参数时，本文采用重启思想以适应PCG算法中预处理矩阵的变化，即将上一步的解作为更新预处理矩阵后的初始猜测。若经过p_{\mathrm{0}}+p_{\mathrm{1}}+\cdots +p_{m-\mathrm{1}}步后计算精度仍未达到要求，则一直使用基于{\mathit{Z}}_{m-\mathrm{1}}和参数{\lambda }_m的预处理迭代，直到收敛为止。这样做的好处：此时残差虽然仍在变化，但已经降至比较小的程度，没有必要再重新构造新的POD基。

基于文献［14］的分析，PCG算法的预处理主要用于加速标准CG算法的初始迭代阶段，当残差已经降到比较小时，预处理加速效果变差，因而参数{\lambda }_m可以相对较大，此时的算法接近标准的CG方法。本文固定参数{\lambda }_m=\mathrm{1}。

3 数值算例

在方程（1）中取\mathit{f}=\mathrm{1}，边界条件\mathit{g}=\mathrm{0}，参数\mu \in [\mathrm{0,8}]，\mathrm{\Omega }分别考虑单位矩形区域和L形区域，如图1所示。对区域\mathrm{\Omega }，从图1出发经过一致加密建立三角形网格剖分，并选择标准线性有限元进行离散，建立全模型，然后用所得线性方程组的求解问题来验证本文算法的有效性。

3.1 算例1：矩形区域求解

在单位矩形区域上一致加密9次，离散方程组的未知量规模为261 121。在off-line阶段，从参数空间[\mathrm{0,8}]均匀采样9个参数，即\{\mathrm{0,1},\mathrm{2},\cdots ,\mathrm{8}\}，并将其作为方程（1）的参数输入得到全模型解。选取动态POD基个数m=\mathrm{4}，每个POD基的向量个数均为N=\mathrm{3}，迭代步数p_i均为5，动态预处理矩阵中的参数\{{\lambda }_{\mathrm{0}},{\lambda }_{\mathrm{1}},\cdots ,{\lambda }_{m-\mathrm{1}}\}均相同，并将其记为\lambda。

根据算法2，{\mathit{Z}}_{\mathrm{0}}的生成不受参数\lambda的影响。为了确定\lambda的值，首先生成POD基{\mathit{Z}}_{\mathrm{0}}，然后随机选择[\mathrm{0,8}]中的一个新参数（实际计算时取7.5），通过观察利用{\mathit{Z}}_{\mathrm{0}}和\lambda求解此参数对应的方程组（2）时\lambda的作用来选择合适的\lambda。图2给出了不同\lambda对应的相对残差随迭代步数的变化。可以看到，\lambda越小，初始迭代速度越快但伴随震荡，数值稳定性越差。因此，综合收敛速度和稳定性，最终取\lambda ={\mathrm{10}}^{-\mathrm{3}}。确定\lambda值后，根据算法2构造动态POD基即完成off-line阶段的准备工作。

在on-line阶段，从[\mathrm{0,8}]中随机选取100个数作为方程（1）的参数\mu，输入全模型生成对应方程组的系数矩阵{\mathit{A}}_h(\mu )和右端向量{\mathit{f}}_h(\mu )，然后用off-line阶段构造好的MPCG算法求解，将结果与标准CG算法进行对比。两种方法的停机准则均为\epsilon ={\mathrm{10}}^{-\mathrm{7}}或最大迭代步数k={\mathrm{10}}^{\mathrm{4}}。

图3为MPCG算法经过100次求解后满足精度所需的迭代步数、两种算法的迭代步数比值和求解时间比值的散点图。此外，表1给出了相应于MPCG算法的迭代步数的极值和均值等的统计结果。可以看到，100次随机求解后，MPCG算法达到精度需要的迭代步数范围约为7~13，时间约为0.095~0.214 s。在迭代步数方面，CG算法达到预定精度所需的步数平均约为MPCG算法的101倍。在计算时间方面，由于MPCG算法每一步迭代时需要处理预处理矩阵相关计算，所以每一步迭代比CG算法更耗时，但总求解时间平均仅为CG算法的\mathrm{1}/\mathrm{43}。因此，MPCG算法相较CG算法大幅提升了计算速度。

3.2 算例2：L形区域求解

在L形区域一致加密8次，离散方程组未知量规模为195 585。类似于算例1，在off-line阶段从参数空间[\mathrm{0,8}]均匀采样9个参数，将其作为方程（1）的参数输入，得到全模型解。选取动态POD基个数m=\mathrm{4}，每个POD基的向量个数均为N=\mathrm{3}，迭代步数p_i均为5，动态预处理矩阵中的参数\{{\lambda }_{\mathrm{0}},{\lambda }_{\mathrm{1}},\cdots ,{\lambda }_{m-\mathrm{1}}\}均相同，并取\lambda ={\mathrm{10}}^{-\mathrm{3}}。

在on-line阶段从[\mathrm{0,8}]中随机选取100个参数进行测试，计算结果如图4所示。此外，表2展示了相应的统计结果。虽然L形区域因解的正则性变差而导致求解更加困难，但经过100次随机求解后，MPCG算法达到预定精度需要的迭代步数范围仍然控制在7~13之间，求解时间控制在0.073~0.209 s之间，可见MPCG算法的收敛速度很快。此外，与CG算法相比，前者所需迭代步数平均仅为后者的\mathrm{1}/\mathrm{125}，求解时间平均仅为\mathrm{1}/\mathrm{59}，再次验证了MPCG算法在计算效率方面的提升。

4 结论

针对参数化泊松方程，本文结合模型降阶思想与PCG提出了MPCG算法，通过在off-line阶段仅使用少量全模型解数据构造高效预处理矩阵来降低计算成本、加快计算速度。单位矩形区域和L形区域上的数值算例表明，MPCG算法达到精度要求所需要的时间平均分别只有传统CG算法的\mathrm{1}/\mathrm{43}和\mathrm{1}/\mathrm{59}，验证了本文算法的高效性。

本文的算法也适用于一般对称正定参数化偏微分方程的快速求解。未来我们将对更多应用模型进行数值验证，并探索将其拓展应用于非对称情形。

参考文献

[1] Gao H， Sun L， Wang J X.PhyGeoNet： Physics-informed geometry-adaptive convolutional neural networks for solving parameterized steady-state PDEs on irregular domain ［J］.J Comput Phys， 2021， 428： 110079.

[2] Chen Y， Dong B， Xu J.Meta-MgNet： Meta multigrid networks for solving parameterized partial differential equations ［J］.J Comput Phys， 2022， 455： 110996.

[3] 王烈衡，许学军.有限元方法的数学基础［M］.北京： 科学出版社， 2004.

[4] Zienkiewicz O C， Taylor R L， Zhu J Z.The finite element method： Its basis and fundamentals ［M］.Oxford： Butterworth-Heinemann， 2013.

[5] Wu L， Azaïez M， Rebollo T C， et al.Certified reduced order method for the parametrized Allen-Cahn equation ［J］.Comput Math Appl， 2023， 134： 167-180.

[6] Quarteroni A， Manzoni A， Negri F.Reduced basis methods for partial differential equations ［M］.Cham： Springer， 2016.

[7] Zhang Q Y， Hu Q X， Wang H， et al.Fast solving of Darcy-Stokes equation by integrating model order reduction with machine learning ［J］.J Sichuan Univ （Nat Sci Ed）， 2025， 62（3）： 584-590.

[8] 张沁逸， 胡奇晓， 王皓， 等.融合模型降阶和机器学习的Darcy-Stokes方程快速求解［J］.四川大学学报（自然科学版）， 2025， 62（3）： 584-590.

[9] Hestenes M R， Stiefel E.Methods of conjugate gradients for solving linear systems ［J］.J Res Nat Bur Stand， 1952， 49（6）： 409-436.

[10] Saad Y， Schultz M H.GMRES： A generalized minimal residual algorithm for solving nonsymmetric linear systems ［J］.SIAM J Sci Stat Comput， 1986， 7（3）： 856-869.

[11] Saad Y.Iterative methods for sparse linear systems ［M］.Philadelphia： SIAM， 2003.

[12] Anzt H， Chow E， Saak J， et al.Updating incomplete factorization preconditioners for model order reduction ［J］.Numer Algorithms， 2016， 73（3）： 611-630.

[13] Carlberg K， Forstall V， Tuminaro R.Krylov-subspace recycling via the POD augmented conjugate-gradient method ［J］.SIAM J Matrix Anal Appl， 2016， 37（3）： 1304-1336.

[14] Santo N D， Deparis S， Manzoni A， et al.Multi-space reduced basis preconditioners for large-scale parametrized PDEs ［J］.SIAM J Sci Comput， 2018， 40（2）： A954-A983.

[15] Xu J， Zhu Y R.Uniformly convergent multigrid methods for elliptic problems with strongly discontinuous coefficients ［J］.Math Models Meth Appl Sci， 2008， 18（1）： 77-105.

基金资助

四川省自然科学基金(2023NSFSC0075)