承接上一篇。第 4 章把 Stokes 问题化为鞍点形式的变分问题，并落到代数系统 (4.10) 上；本篇继续第 5 章（完整 Navier–Stokes 方程的有限元处理）、第 6 章（鞍点系统的迭代解法）与附录 A（稀疏矩阵）。行文中穿插了我复习时产生的若干疑问，以及澄清的结果——其中有几处是我原先的理解出了偏差，一并记录下来。

5 Navier–Stokes 方程的有限元方法

第 4 章的变分理论是针对 Stokes 问题发展的，它是丢弃对流加速度之后所得的线性内核。然而真实流动携带惯性，而恰恰是对流项使流体力学既丰富又困难：它令方程非线性化，驱动向湍流的转捩，并且——一旦离散——摧毁了使 Stokes 系统如此易于处理的对称性。本章把有限元框架推广到完整的不可压缩 Navier–Stokes 方程：写出弱形式，记录关于存在唯一性的已知结论，讨论保持离散稳定的有限元空间，并且——对后文而言最为重要——说明非线性如何通过一列线性问题得到消解。这些线性问题中的每一个都是第 4 章那种鞍点系统，只是此时一般不再对称，其高效求解正是第 6 章的主题。

5.1 方程及其弱形式

在定常不可压缩情形下，把第 2 章的动量平衡中的物质导数展开，得

与 Stokes 问题 (4.2) 唯一的——但具有决定意义的——差别，是对流项 $\rho (\mathbf{u}\cdot \mathrm{\nabla })\mathbf{u}$ ，它关于未知速度是二次的。

像第 4 章那样作检验并加入对流贡献，弱问题除了 (4.3) 的粘性形式 $a$ 与约束形式 $b$ 之外，还要用到三线性形式

支配分析的是 $c$ 的两条性质。由 Hölder 不等式与 Sobolev 嵌入 $H^1\hookrightarrow L^4$ （对 $d\le 4$ 成立，故在 $d\in \{2,3\}$ 中可用），该形式有界，

常数 $N$ 只依赖于 $\mathrm{\Omega }$ 与 $\rho$ 。此外，对无散度的对流场，分部积分给出反对称性

恒等式 $c(\mathbf{w};\mathbf{u},\mathbf{u})=0$ 是一条乔装的离散能量平衡：它表明对流输运动量而不创造或销毁动能，这对连续问题的先验估计与离散的稳定性都不可或缺。由于有限元速度未必精确无散度，人们常把 $c$ 换成它的反对称部分 $\tilde{c}(\mathbf{w};\mathbf{u},\mathbf{v})=\frac{1}{2}(c(\mathbf{w};\mathbf{u},\mathbf{v})-c(\mathbf{w};\mathbf{v},\mathbf{u}))$ ，它恒满足 $\tilde{c}(\mathbf{w};\mathbf{u},\mathbf{u})=0$ 。

5.2 存在性与唯一性

非线性改变了适定性理论的性质。Stokes 问题对任意数据都唯一可解（定理 4.5），而 Navier–Stokes 问题只保证至少有一个解；唯一性则要求数据相对于粘性足够小。

存在性由 Galerkin 论证配合 Brouwer 不动点定理给出，唯一性由依赖反对称性 (5.5) 的压缩估计给出。引入 Reynolds 数 $\mathrm{Re}=UL/\nu$ （ $\nu =\mu /\rho$ 为运动粘度， $U,L$ 为特征速度与长度），条件 (5.6) 本质上就是关于 $\mathrm{Re}$ 的小性条件：低 Reynolds 数下流动唯一，且通常是定常层流；随着 $\mathrm{Re}$ 增大，解可能丧失唯一性，最终丧失定常性——这正映照着通向湍流的物理路径，也是高 Reynolds 数计算中数值困难的解析投影。

5.3 inf–sup 稳定单元

离散从第 4 章继承了这样一条要求：速度／压强对 $({\mathbf{V}}_h,Q_h)$ 必须关于 $h$ 一致地满足假设 4.8 的离散 inf–sup 条件。对流并不放松这一要求——它进一步约束速度空间，而压强稳定性仍由 inf–sup 支配。两族稳定单元贯穿全书：

• Taylor–Hood 族：速度用 $k$ 次、压强用 $k-1$ 次连续分片多项式（ $k\ge 2$ ），最低阶成员是三角形上的 ${\mathbb{P}}_2/{\mathbb{P}}_1$ ，或四边形上的 ${\mathbb{Q}}_2/{\mathbb{Q}}_1$ 。
• MINI 单元：连续 ${\mathbb{P}}_1$ 速度用一个内部泡函数加以富化，配以连续 ${\mathbb{P}}_1$ 压强；泡函数恰好提供了在最低阶满足 inf–sup 所需的那些速度自由度。

两种情形下速度空间都比压强空间"更富"，这正是离散散度得以满射的原因。相反，等阶插值（如 ${\mathbb{P}}_1/{\mathbb{P}}_1$ ）违反该条件并容许虚假压强模态——即第 6 章将从代数上分析的棋盘格不稳定性——除非加以稳定化（5.5 节）。

5.4 线性化：Picard 与 Newton 迭代

(5.3) 的离散对应物是一个非线性代数方程组，需迭代求解。两种标准线性化都把每一步化为 Oseen 型线性问题——即围绕一个已知速度场作对流增广的 Stokes 问题。

Picard（Oseen）迭代。 把对流冻结在当前迭代 ${\mathbf{u}}^k$ 上，求解 $({\mathbf{u}}^{k+1},p^{k+1})$ ：

${\mathbf{u}}^k$ 固定后问题是线性的。Picard 迭代全局收敛但只是线性收敛，且其压缩因子随 Reynolds 数增大而恶化；它稳健，是远离解时的首选。

Newton 迭代。 把映射 $\mathbf{u}\mapsto c(\mathbf{u};\mathbf{u},\mathbf{v})$ 在 ${\mathbf{u}}^k$ 处沿方向 $\delta \mathbf{u}$ 线性化，产生两项 $c(\delta \mathbf{u};{\mathbf{u}}^k,\mathbf{v})+c({\mathbf{u}}^k;\delta \mathbf{u},\mathbf{v})$ 。于是 Newton 步对增量 $(\delta \mathbf{u},\delta p)$ 求解

其中 $(r^k,s^k)$ 是当前迭代的残量。Newton 二次收敛但只是局部收敛，故实践中常先作几步 Picard 再切换到 Newton。

正是如此。 $F_{\nu }((v,p);(w,q))=0$ （对所有检验函数 $(w,q)$ 成立）是定常 Navier–Stokes 的弱非线性残量形式，而 Newton 法求解 $F_{\nu }=0$ 需要 $F_{\nu }$ 在当前迭代 $(v_k,p_k)$ 处的 Fréchet（Gâteaux）导数。按定义，它就是沿方向 $(u,r)$ 的方向导数：把当前迭代沿 $(u,r)$ 扰动一个小量 $\epsilon$ ，代入 $F_{\nu }$ ，对 $\epsilon$ 求导，再在 $\epsilon =0$ 处取值。

因此三组变量的角色是：

• $(v_k,p_k)$ ：当前的速度／压强迭代；
• $(u,r)\in \mathbf{V}\times Q$ ：每步待求的 Newton 增量，即上文 (5.8) 中的 $(\delta \mathbf{u},\delta p)$ ；
• $(w,q)$ ：一如既往的检验函数。

Newton 步即求解线性问题 $D{F}_{\nu }((v_k,p_k);(u,r),(w,q))=-F_{\nu }((v_k,p_k);(w,q))$ （对所有 $(w,q)$ ），再更新 $(v_{k+1},p_{k+1})=(v_k,p_k)+(u,r)$ 。之所以说这是"线性化"，正是因为该问题关于未知的 $(u,r)$ 是线性的。

推导中出现的两个线性化对流项 $((u\cdot \mathrm{\nabla })v_k,w)+((v_k\cdot \mathrm{\nabla })u,w)$ ，恰是三线性形式 $c$ 的线性化：由于 $\mathbf{v}=v_k+\epsilon u$ 同时出现在 $c$ 的前两个变量位置上，对 $\epsilon$ 求导时每一处各贡献一项。这也解释了 Newton 与 Picard 的分野：Picard 只冻结一个位置上的 $v_k$ （给出单独一项 $c(v_k;u,w)$ ），Newton 则对两处出现同时求导并保留两项——多出来的那一项换来的是局部二次收敛，而非 Picard 的线性收敛。

在两种情形下，每步所解线性问题的代数形式都是

其中 $A$ 是第 4 章那个对称的（矢量 Laplace）扩散矩阵， $N(\mathbf{w})$ 是围绕当前速度 $\mathbf{w}$ 由 $c(\mathbf{w};\cdot ,\cdot )$ 装配出的对流矩阵（对 Newton 而言， $F$ 中还多一个由 $c(\cdot ;\mathbf{w},\cdot )$ 来的类反应项）。由于对流是一阶、非自伴算子，Oseen 块 $F$ 不对称——与 (4.10) 中对称的 Stokes 块形成对照。于是 (5.9) 是一个非对称鞍点系统。第 6 章为对称 Stokes 系统发展的块预条件子可以推广到它，只是 Krylov 方法要由 MINRES 换成 GMRES。

5.5 对流主导流动的稳定化

标准 Galerkin 离散只在扩散控制对流时才准确。当网格 Péclet 数（单元 Reynolds 数） ${\mathrm{Pe}}_K=|\mathbf{w}|h_K/(2\nu )$ 超过 1 时——也就是在可负担的网格上作对流主导计算时——Galerkin 速度会发展出非物理振荡，与标量对流–扩散方程的情形完全一致。稳定化格式通过添加依赖网格的加权残量项来治愈它：这些项相容（对精确解为零），却能恢复对对流导数的控制。

流线迎风 Petrov–Galerkin（SUPG） 方法在动量方程上增添

即用流线导数 $\mathbf{w}\cdot \mathrm{\nabla }\mathbf{v}$ 去检验强残量，单元参数 ${\tau }_K\propto h_K/|\mathbf{w}|$ 。因为被加权的只是残量，精确解不受影响，最优精度得以保留，而添加项沿流线方向注入了恰好足够的数值扩散来压制振荡。相关的 Galerkin/最小二乘法与基于残量的变分多尺度法遵循同一原则。

另有两种稳定化常与 SUPG 联用：PSPG 用压强梯度去加权残量，从而彻底绕开 inf–sup 条件，容许方便的等阶插值——它正是第 6 章中代数稳定化 ${\mathbb{Q}}_1$ – ${\mathbb{Q}}_1$ 对的连续对应物；grad-div 稳定化添加 $\gamma {\int }_{\mathrm{\Omega }}(\mathrm{div}\mathbf{u})(\mathrm{div}\mathbf{v})\mathrm{d}x$ 以改善质量守恒与压强稳健性。对线性代数而言，稳定化把 (5.9) 中原本为零的压强–压强块填成一个对称半正定矩阵 $C$ ，从而给出下一章处理的一般形式 (6.1)，而并不改变所解系统的鞍点性质。

5.6 非定常流动的时间离散

对含时流动，动量方程保留局部加速度 $\rho {\partial }_t\mathbf{u}$ 。先作空间离散——线方法——把弱问题化为关于系数矢量的微分代数方程组

其中 $M$ 是速度质量矩阵。不可压缩约束不带时间导数，故 (5.11) 是指标为 2 的微分代数系统，压强扮演把速度约束在无散度流形上的代数乘子。

隐式时间积分器—— $\theta$ 格式（后向 Euler、Crank–Nicolson）或后向差分公式——在新时间层上处理刚性的粘性项与约束项。于是一步需要求解

其结构与 (5.9) 相同，只是速度块被 $M/\mathrm{\Delta }t$ 平移。这一平移使该块更接近对角占优、条件数更好，且 $\mathrm{\Delta }t$ 越小越好——正如第 6 章所示，这对迭代求解器及其预条件子是有利的。对流可以隐式处理（每步内作 Picard 或 Newton 求解），也可以在半隐式（IMEX）格式中显式处理，以对 $\mathrm{\Delta }t$ 的限制换取免去非线性求解。

求解完全耦合系统 (5.12) 之外的另一条路是投影（压强修正）方法：把每个时间步拆成速度更新，再跟一次恢复不可压缩性的压强 Poisson 求解。它以一列更简单的椭圆求解替换耦合的鞍点求解，在大规模非定常计算中被广泛使用。

5.7 小结

Navier–Stokes 方程的有限元处理建立在与 Stokes 问题相同的变分基础上，由三线性对流形式加以扩展，并受同一 inf–sup 条件约束。其特征性的困难是非线性，它在定常计算中由 Picard 或 Newton 迭代、在非定常计算中由隐式时间推进，消解为一列线性鞍点系统 (5.9)、(5.12)。这些系统一般非对称，且在细网格上规模极大；它们的高效迭代求解因此是任何实用 Navier–Stokes 求解器的关键所在。

6 鞍点系统的迭代解法

6.1 动机与问题设定

Stokes 方程——以及更一般地，任何以压强作为施行不可压缩约束的 Lagrange 乘子的流动模型——的有限元离散，导致具有特征性 $2\times 2$ 块结构的大型稀疏线性系统。稀疏性是有限元基函数局部支集的直接后果，使每个未知量只与少数几个邻居耦合。记 $u\in {\mathbb{R}}^n$ 为离散速度自由度矢量、 $p\in {\mathbb{R}}^m$ 为离散压强自由度矢量，这类系统取如下形式：

其中 $A\in {\mathbb{R}}^{n\times n}$ ， $B\in {\mathbb{R}}^{m\times n}$ ， $C\in {\mathbb{R}}^{m\times m}$ ，且通常 $m\le n$ 。

在假设 6.1 之下，矩阵 $\mathcal{A}$ 对称但不定：它同时具有正、负特征值（见下面的定理 6.3）。这一不定性是大部分数值困难的根源，并且排除了共轭梯度法的直接应用——后者要求系统是定的。

"鞍点系统"之名反映了 (6.1) 的变分来源。当 $C=0$ 时，解 $(u,p)$ 是 Lagrange 泛函

的唯一驻点，它对应于等式约束二次极小化问题

(6.3) 的一阶最优性（KKT）条件恰是 $C=0$ 时的 (6.1)。 $(u,p)$ 关于原变量 $u$ 是极小点、关于乘子 $p$ 是极大点，故 $\mathcal{L}$ 在此有一个鞍点，压强 $p$ 即施行离散不可压缩约束 $Bu=g$ 的 Lagrange 乘子。

6.2 代数结构与 Schur 补

由于 $A$ 是 SPD，特别地它可逆，于是速度未知量可以从 (6.1) 中消去。由第一块行解出

代入第二块行 $Bu-Cp=g$ ，得到只关于压强的约化系统

矩阵 $S$ 就是 $A$ 在 $\mathcal{A}$ 中的（负）Schur 补。它把速度与压强之间的全部耦合浓缩进一个 $m\times m$ 矩阵，在分析与预条件子构造中都居于核心地位。

证明。 对称性由 $A^{-1}$ 与 $C$ 的对称性立得。为证正定性，取 $q\in {\mathbb{R}}^m$ ， $q\ne 0$ ，则

因为 $A^{-1}$ 是 SPD、 $C$ 是 SPSD。取等号要求两项同时为零；特别地 $(B^{T}q{)}^T{A}^{-1}(B^{T}q)=0$ 迫使 $B^{T}q=0$ （因 $A^{-1}$ 为 SPD）。而 $B$ 满行秩，故 $B^T$ 单射， $B^{T}q=0$ 蕴含 $q=0$ ，矛盾。 $◻$

大方向对，但推导的顺序恰好相反，而且有两处需要更正。

第一， $S$ 并不是从假设的分解中"反解"出来的，而是直截了当消元的结果：由 (6.4) 解出 $u$ 、代回第二块行，剩下的关于 $p$ 的矩阵就是 $S$ 。分解 (6.6) 是事后对这一消元过程的识别； $L$ 的非对角块 $B{A}^{-1}$ 是把块乘法展开、与 $\mathcal{A}$ 的各块相比对而读出的（验证式 $B{A}^{-1}{B}^T-S=-C$ ），而非独立待解的未知量。

第二，我漏掉了使惯性论证得以成立的关键一环：引理 6.2（ $S$ 为 SPD）。没有它， $D$ 的惯性无从谈起。

第三，是 Sylvester 惯性律（Sylvester’s law of inertia），不是我写的"Silver"。它断言：两个实对称矩阵若合同（ $N=LM{L}^T$ ， $L$ 非奇异），则它们有相同的惯性 $\mathrm{In}(M)=(n_{+},n_{-},n_0)$ ——合同可以改变特征值本身，但不改变它们的符号。

于是完整链条是： $D=\mathrm{diag}(A,-S)$ 中 $A$ 为 SPD（贡献 $n$ 个正特征值）、 $-S$ 为负定（贡献 $m$ 个负特征值），故 $\mathrm{In}(D)=(n,m,0)$ ；而 $L$ 是单位下三角阵， $detL=1$ 非奇异，故 $\mathcal{A}=LD{L}^T$ 与 $D$ 合同，由 Sylvester 律 $\mathrm{In}(\mathcal{A})=(n,m,0)$ ： $\mathcal{A}$ 对称不定，无零特征值——这也再度确认了非奇异性，而且无需显式形成 $S$ （形成 $A^{-1}{B}^T$ 或 $S$ 一般会摧毁稀疏性）。

6.3 可解性与离散 inf–sup 条件

定理 6.3 保证了固定离散问题的可解性。但在有限元语境下，人们求解的是由网格尺度 $h$ 标记的一整族系统，真正的问题是解算子在 $h\to 0$ 时是否保持有界。这一稳健性由离散 inf–sup（LBB）条件支配。

这里的两个范数是相应有限元函数的自然能量范数： $A$ 与 $M_p$ 分别是 $a$ 与 $L^2$ 内积在两组基下的 Gram 矩阵，故 $\parallel v{\parallel }_A$ 以 $H^1$ 能量范数度量速度， $\parallel q{\parallel }_{M_p}$ 以 $L^2$ 范数度量压强。

不是，两者毫无关系。

Frobenius 范数是矩阵上的范数： $\parallel M{\parallel }_F=\sqrt{\sum _{ij}{M}_{ij}^2}=\sqrt{\mathrm{tr}(M^{T}M)}$ ，即把矩阵元当作一个长矢量取 Euclid 范数，纯属记账，与任何物理能量无关。

能量范数则是**矢量（或函数）**上的范数，由一个 SPD 双线性形式 $A$ （连续情形下为 $a(\cdot ,\cdot )$ ）诱导：

它之所以是真正的范数，正是因为 $A$ 是 SPD：正定性给出 $v\ne 0\Rightarrow \parallel v{\parallel }_A>0$ ，双线性与对称性经由 $A$ -内积 $⟨u,v{⟩}_A:=u^{T}Av$ 上的 Cauchy–Schwarz 给出三角不等式。

"能量"之名是字面意义上的物理陈述而非记号巧合：对本课程的 Stokes／弹性型问题， $a(\mathbf{u},\mathbf{u})={\int }_{\mathrm{\Omega }}2\mu \mathbf{D}(\mathbf{u}):\mathbf{D}(\mathbf{u})\mathrm{d}x$ 恰是（两倍的）粘性耗散率；在结构力学的历史源头处，线弹性问题的 $a(u,u)=\int \sigma :\epsilon \mathrm{d}x$ 就是变形 $u$ 中储存的应变能的两倍。离散地，由 Korn 不等式 $a({\mathbf{v}}_h,{\mathbf{v}}_h)=2\mu \parallel \mathbf{D}({\mathbf{v}}_h){\parallel }_{L^2}^2\approx \mu |{\mathbf{v}}_h{|}_{H^1}^2$ ，故 $\parallel v{\parallel }_A$ 与 $H^1$ 半范数等价。

这也解释了它在第 6 章反复出现的原因：CG 所极小化的正是能量范数（见 6.5 节），这与 Céa 引理在上一层的最佳逼近性质是同一个变分原理。

6.3.1 三种代表性单元对

离散 inf–sup 常数 $\beta$ 完全由速度与压强有限元空间的选取决定。在四边形网格上写 ${\mathbb{Q}}_k$ 表示每个坐标方向上 $k$ 次的连续分片（映射）多项式空间，三种定性不同的结局如下。

(i) 不稳定的等阶对 ${\mathbb{Q}}_1$ – ${\mathbb{Q}}_1$ 。 对速度各分量与压强使用同一连续双线性空间，是最自然、也最危险的选择。这一对直接违反 inf–sup 条件：存在非零压强矢量 $q$ ——即虚假压强模态——它对每一个离散速度都是"不可见的"，

在均匀网格上，典型的这类模态就是棋盘格模态：由于 ${\mathbb{Q}}_1$ 压强是节点型（连续）的，它对网格顶点 $(i,j)$ 赋以交替值 $(-1{)}^{i+j}$ 。

网格／基的策略是对的，但记号 $v{B}^{T}q$ 不成立——那不是一个合法的对象（ $v\in {\mathbb{R}}^n$ 是列矢量， $B^{T}q\in {\mathbb{R}}^n$ 也是，两者不能相乘）。正确的写法是转置：

最后一步用的是 (4.9) 中 $B$ 的定义 $B_{ki}=b({\phi }_i,{\psi }_k)$ 。于是

这才是要逐单元验证的对象。把 $q_h$ 展开为棋盘格模态、 $v_h$ 展开为 ${\mathbb{Q}}_1$ 速度基，在每个单元上算 ${\int }_K{q}_h\mathrm{div}{v}_h\mathrm{d}x$ ： $q_h$ 的 $\pm 1$ 节点值贡献在相邻单元之间成对相消——分片双线性的散度分辨不出波长为 $2h$ 的振荡，这就是"棋盘格模态与任何离散速度都不产生散度配对"的具体含义。

结论必须显式写出（这是我原先漏掉的）：既然对非零的 $q$ 有 $B^{T}q=0$ ，则 $B$ 不满行秩，假设 6.1(2) 被违反； $S=B{A}^{-1}{B}^T$ 奇异（因 $Sq=B{A}^{-1}(B^{T}q)=0$ ，即 $q$ 落在 $S$ 的核里），等价地离散 inf–sup 常数 $\beta =0$ 。这与引理 6.2 的证明严丝合缝：那里正是 $B^{T}q=0$ （ $q\ne 0$ ）破坏了 $S$ 的正定性。

(ii) 稳定的 Taylor–Hood 对 ${\mathbb{Q}}_2$ – ${\mathbb{Q}}_1$ 。 把速度次数提高一阶即可治愈不稳定性；在形状正则的网格族上，它满足 (6.7) 且常数 $\beta >0$ 与 $h$ 无关（ $d=2,3$ 皆然）。标准途径是 Fortin 判据，它把离散 inf–sup 条件的验证归结为构造单一的插值型算子。

(iii) 稳定化的 ${\mathbb{Q}}_1$ – ${\mathbb{Q}}_1$ 。 方便的等阶对可以由稳定化拯救：添加一个控制虚假模态、而对光滑解保持相容的压强相关项。经典例子是 Brezzi–Pitkäranta 稳定化，它在第二块行上增添双线性形式

在矩阵形式下，(6.13) 恰好贡献 (6.1) 中那个对称半正定块 $C$ ，使 $\mathcal{A}$ 的 $(2,2)$ 元变为 $-C\ne 0$ 。

从" $S^{\prime }=S+C$ 正定故非奇异"这一层说，事情就完结了；但值得再补两句，说明为什么恰好治得住棋盘格。虚假模态满足 $B^{T}q=0$ ，故 $Sq=0$ —— $q$ 正落在 $S$ 的核里，这才是 $S$ 单独奇异的方向；而稳定化项 $q^{T}Cq=\alpha \sum _K{h}_K^2\parallel \mathrm{\nabla }{q}_h{\parallel }_{L^2(K)}^2$ 对逐节点振荡的棋盘格模态取严格正值。也就是说， $C$ 恰在 $S$ 消失的那个方向上为正，这才把 $S+C$ 从奇异推成 SPD，而不是笼统地"加了个正则项"。

6.4 谱性质

迭代方法的收敛速率通常由所作用矩阵的谱——或预条件后算子的谱——支配。对未预条件的鞍点矩阵（ $C=0$ ），命题 6.8 给出：设 $0<{\lambda }_{min}(A)\le {\lambda }_{max}(A)$ 为 $A$ 的极端特征值、 ${\sigma }_{min}\le {\sigma }_{max}$ 为 $B$ 的极端奇异值（ ${\sigma }_{min}>0$ 因 $B$ 满行秩），则 $\mathcal{A}$ 的每个正特征值 ${\lambda }^{+}$ 与负特征值 ${\lambda }^{-}$ 满足

两个区间由原点分开：正特征值以 ${\lambda }_{min}(A)>0$ 为下界，负特征值以一个严格负的数为上界——之所以严格负，正是因为 ${\sigma }_{min}>0$ ，即 $B$ 满行秩（离散 inf–sup 性质）。对未预条件的 Stokes 矩阵，这两个区间在网格加密时会展开（ $\mathrm{cond}(A)=\mathcal{O}(h^{-2})$ ），谱分布愈发恶劣，迭代步数随网格加密而增长——这使 6.6 节的有效预条件不可或缺。

6.5 Krylov 子空间方法

在三维中，对 $\mathcal{A}$ 作稀疏 $LD{L}^T$ 直接分解会因填充而在内存与运算量上都变得不可行。大规模问题因此转向 Krylov 子空间方法，它们只通过矩阵–矢量乘积访问 $\mathcal{A}$ ，从而保持稀疏性（稀疏存储格式与该乘积内核见附录 A）。

6.5.1 投影方法与 Krylov 子空间

本节所有方法共享一个原则。求解 $\mathcal{A}x=b$ 的投影方法从 $m$ 维仿射搜索空间 $x_0+{\mathcal{K}}_m$ 中抽取近似 $x_m$ ，要求残量正交于 $m$ 维检验空间 ${\mathcal{L}}_m$ ：

当 ${\mathcal{L}}_m={\mathcal{K}}_m$ 时投影正交（Galerkin 条件），否则斜交（Petrov–Galerkin 条件）。

这些空间是嵌套的，并带有解释方法及其收敛性的多项式结构：每个 $x_m\in x_0+{\mathcal{K}}_m$ 形如 $x_m=x_0+q(\mathcal{A})r_0$ （ $\mathrm{deg}q\le m-1$ ），故残量为

于是 Krylov 方法隐式地构造了一个以 $p_m(0)=1$ 归一化的残量多项式，其收敛性取决于这样的多项式在 $\mathcal{A}$ 的谱上能被压到多小。

6.5.2 三种方法

• 共轭梯度（CG，用于 SPD 块）。 取 ${\mathcal{L}}_m={\mathcal{K}}_m$ 。Galerkin 条件说残量正交于 ${\mathcal{K}}_m$ ，等价地误差** $\mathcal{A}$ -正交**于搜索空间；对 SPD 的 $\mathcal{A}$ ，这就是能量范数下的最佳逼近 $\parallel x-x_m{\parallel }_{\mathcal{A}}=\underset{z\in x_0+{\mathcal{K}}_m}{min}\parallel x-z{\parallel }_{\mathcal{A}}$ ，与极小化二次能量 $\varphi (z)=\frac{1}{2}{z}^T\mathcal{A}z-b^{T}z$ 等价。优点：能量范数下最优、存储恒定、收敛只由 $\mathrm{cond}(\mathcal{A})$ 支配。缺点：只适用于 SPD 系统，在不定的 $\mathcal{A}$ 上失效。因此 CG 在本章只用于 SPD 子问题：速度块 $A$ 、Schur 补 $S$ ，以及 Uzawa 与块预条件迭代的内层求解。
• 极小残量（MINRES，用于对称不定系统）。 $\mathcal{A}$ 对称但不定，能量范数（以及 CG）不可用。改取 ${\mathcal{L}}_m=\mathcal{A}{\mathcal{K}}_m$ ，即极小化 $\parallel b-\mathcal{A}z{\parallel }_2$ 。对称性使 ${\overline{H}}_m$ 三对角，其 QR 因子经短递推更新，故每步代价恒定。预条件子 $\mathcal{P}$ 必须 SPD 才能保持对称性，此时 MINRES 在 ${\mathcal{P}}^{-1}$ 范数下极小化残量，收敛由 ${\mathcal{P}}^{-1}\mathcal{A}$ 的特征值支配：若它们落在 $[-a,-b]\cup [c,d]$ 中（ $a,b,c,d>0$ ），则 $2k$ 步后残量降低约

• 广义极小残量（GMRES，用于非对称预条件）。 非对称预条件子使 ${\mathcal{P}}^{-1}\mathcal{A}$ 非对称， ${\overline{H}}_m$ 退化为满上 Hessenberg 阵，Lanczos 缩短失效。同样的残量极小化定义了 GMRES。缺点是完整 Arnoldi 递推使工作量与存储随迭代指标线性增长，实践中须重启（GMRES( $k$ )）。

6.5.3 CG 的具体形式

定义。 对 SPD 矩阵 $A$ ，两个非零矢量 $u,v$ 称为 $A$ -共轭（或 $A$ -正交），若 $u^{T}Av=0$ ，即它们在 $A$ -内积 $⟨u,v{⟩}_A:=u^{T}Av$ 下正交——正是上文诱导能量范数的那个内积。一组方向 $\{p_0,p_1,\dots \}$ 互相 $A$ -共轭，若 $p_i^{T}A{p}_j=0$ （ $i\ne j$ ）。

为什么要这个性质。 CG 在 $x_0+{\mathcal{K}}_k$ 上极小化能量范数误差。若搜索方向 $p_0,\dots ,p_{k-1}$ 张成 ${\mathcal{K}}_k$ 且互相 $A$ -共轭，则在 ${\mathcal{K}}_k$ 上的极小化解耦为沿各个 $p_i$ 的 $k$ 个独立一维极小化——一旦沿某个方向走过，就再也不必回头重新优化。这一解耦正是 CG 能塌缩成短而廉价的递推、而不必像完整 GMRES／Arnoldi 那样作不断增长的正交化的全部原因。

导出 ${\beta }_k$ 。 把新方向取为残量被前一方向修正的形式 $p_{k+1}=r_{k+1}+{\beta }_k{p}_k$ ，并要求 $p_{k+1}$ 与 $p_k$ 共轭，即 $p_{k+1}^{T}A{p}_k=0$ ：

再用递推中已有的两条事实化简。由 $r_{k+1}=r_k-{\alpha }_{k}A{p}_k$ 得 $A{p}_k=(r_k-r_{k+1})/{\alpha }_k$ ，故

最后一步用了残量互相正交（这是 Galerkin 条件的推论，与共轭性一同归纳地证明）。又由步长公式 ${\alpha }_k=(r_k^T{r}_k)/(p_k^{T}A{p}_k)$ 得 $p_k^{T}A{p}_k=(r_k^T{r}_k)/{\alpha }_k$ ，代入即得

${\alpha }_k$ 恰好抵消，且不需要任何额外的矩阵–矢量乘积。

一个值得能答上来的细节： 为什么只让 $p_{k+1}$ 与 $p_k$ 共轭，就能保证整族方向互相共轭？这正是 Krylov 子空间结构发挥作用之处：归纳地 $p_i\in {\mathcal{K}}_{i+1}$ ，故 $A{p}_i\in {\mathcal{K}}_{i+2}$ ；而 $r_{k+1}\perp {\mathcal{K}}_k$ （残量对迄今整个搜索空间的 Galerkin 正交性），于是自动有 $r_{k+1}^{T}A{p}_i=0$ 对所有 $i<k$ 成立。也就是说，递推中只对 $p_k$ 施加共轭性，就白得了对所有先前方向的共轭性——这正是把一般的（稠密、 $O(m^2)$ 存储的）Lanczos 型 Galerkin 投影，变成 CG 的两项、 $O(1)$ 存储短递推的机制。

6.6 块预条件子

块分解 (6.6) 提示了由对角块 $A$ 与 $S$ 的（近似）构造预条件子。先分析理想版本，其中 $A$ 与 $S$ 被精确使用；它们给出可观的特征值聚集，是实用变体的理论蓝本。

证明的要点是： $q=0$ 的情形给出 $\lambda =1$ （对应 $\ker B$ 中的 $(n-m)$ 维解空间）； $q\ne 0$ 时消去 $v=(\lambda -1{)}^{-1}{A}^{-1}{B}^{T}q$ 并利用 $S=B{A}^{-1}{B}^T$ ，化为标量关系 $\lambda (\lambda -1)=1$ ，即 ${\lambda }^2-\lambda -1=0$ ，其根正是黄金分割值。极小多项式次数为 3，故 3 步终止。

保留非对角耦合则给出块上三角预条件子 ${\mathcal{P}}_T=\left(\begin{array}{cc}A& B^T\\ 0& -S\end{array}\right)$ 。它非对称，须与 GMRES 联用。此时 $\mathcal{A}{\mathcal{P}}_T^{-1}=\left(\begin{array}{cc}I& 0\\ B{A}^{-1}& I\end{array}\right)$ ，只有单一特征值 1 且 $(\mathcal{A}{\mathcal{P}}_T^{-1}-I{)}^2=0$ ，故 GMRES 至多 2 步终止（定理 6.13）。

实用版本。 理想预条件子不能直接使用：每次作用都要精确求解带 $A$ 与带稠密 $S$ 的系统，代价与解原问题相当。实用配方是把 $A^{-1}$ 与 $S^{-1}$ 换成谱等价且易于作用的近似 ${\hat{A}}^{-1}$ 与 ${\hat{S}}^{-1}$ ：

• 近似速度求解 ${\hat{A}}^{-1}$ ：由于 $A$ 是离散（矢量）Laplace 算子，单个多重网格 V-循环或代数多重网格（AMG）循环即可给出与 $A$ 谱等价、常数与网格无关的 $\hat{A}$ ，代价线性。
• 近似 Schur 补 ${\hat{S}}^{-1}$ ：由谱等价性，对 Stokes 问题，压强质量矩阵 $M_p$ （按粘度倒数缩放）是 $S$ 的网格无关近似； $M_p$ 良态，故几步 Chebyshev 或 Jacobi 扫描——甚至只取其对角——就足以廉价地作用 ${\hat{S}}^{-1}$ 。

当 $\hat{A},\hat{S}$ 与 $A,S$ 以网格无关的常数谱等价时， ${\hat{\mathcal{P}}}_D^{-1}\mathcal{A}$ 的特征值仍聚集在远离零的固定区间内，MINRES 界 (6.30) 便预言与 $h$ 无关的迭代步数。

6.7 Uzawa 方法及其变体

核心思想一句话：Uzawa 就是对 6.2 节的 Schur 补系统作 Richardson 迭代，只是包装得让你从不必显式形成 $S$ 。

第一步：回忆约化系统。 由 (6.5)，消去 $u$ 后得关于 $p$ 的 $m\times m$ SPD 系统 $Sp=c$ ，其中 $c:=B{A}^{-1}f-g$ （ $S$ 的 SPD 性由引理 6.2 保证）。

第二步：对该系统作 Richardson 迭代。 求解 $Sp=c$ 的一般定常（Richardson）迭代为

其中 $\omega >0$ 为松弛参数。由标准 Richardson 理论，迭代矩阵为 $I-\omega S$ ，谱半径小于 1 当且仅当 $0<\omega <2/s_{max}$ （ $s_{max}$ 为 $S$ 的最大特征值）。

第三步：不形成 $S$ ，改用一次速度求解替代。 关键的观察是

因此只要定义 $u^{k+1}:=A^{-1}(f-B^T{p}^k)$ ——字面意思就是：把当前压强迭代冻结为已知载荷，解一个 Stokes 型的速度问题——所需残量就只是 $B{u}^{k+1}-g$ ，全程不出现 $S$ 。这就给出

结构上，每个 Uzawa 步就是"给定当前压强猜测精确求解速度块，再沿减小散度约束违反量 $B{u}^{k+1}-g$ 的方向轻推压强"；但把 $u^{k+1}$ 重新消去后，它与 $S$ 上的纯 Richardson 在代数上完全等同。最优参数与收敛率随之由同一套 SPD 理论给出： ${\omega }^{\ast }=2/(s_{min}+s_{max})$ 使谱半径最小，此时渐近收敛因子为 $(\mathrm{cond}(S)-1)/(\mathrm{cond}(S)+1)$ 。

实用版本为何不同。 精确版本有两处不足：(i) 每个外迭代都要作一次完整的速度求解 $A^{-1}(\cdot )$ ，昂贵；(ii) 收敛随 $\mathrm{cond}(S)$ 增大而恶化（网格加密时它通常确实增大）。补救办法与前面"谱等价替身"的思路一致：把精确的 $A^{-1}$ 换成廉价近似求解 ${\hat{A}}^{-1}$ （例如一个 AMG V-循环），把压强更新用 ${\hat{S}}^{-1}$ （例如 $M_p$ 上的几步 Chebyshev 扫描）预条件。

采用与 6.6.3 节相同的谱等价构件 $\hat{A},\hat{S}$ ，非精确 Uzawa 以与网格无关的速率收敛。密切相关的 Arrow–Hurwicz 方法把内层速度求解换成单步梯度步，以更慢的收敛率换取更快的单步迭代。另一个稳健变体是增广 Lagrange 预条件：把 $(1,1)$ 块换成 $A_{\gamma }=A+\gamma B^T{W}^{-1}B$ ，其 grad-div 型项不改变 (6.1) 的解（因解处 $Bu=g$ ），却改善了 Schur 补的近似质量；代价是 $A_{\gamma }$ 更病态、更各向异性，需要专门设计的多重网格。

6.8 应用于 Stokes 系统

对 Stokes 问题， $C=0$ ，块 $A$ 是按粘度 $\nu$ 缩放的离散矢量 Laplace 算子， $B$ 是离散散度。决定性的性质是 Schur 补与压强质量矩阵之间的谱等价，对 Stokes 算子它取如下尖锐形式：

其中 $\beta$ 是定义 6.4 的离散 inf–sup 常数。上界常数为 1，下界为 $\beta$ ，对 inf–sup 稳定的单元对（如 Taylor–Hood $P_2$ – $P_1$ 族）两者都与 $h$ 无关。这就为选取

作为 Schur 补近似提供了依据。配合矢量 Laplace 算子的多重网格近似 ${\hat{A}}^{-1}$ ，块对角预条件子 ${\hat{\mathcal{P}}}_D$ 使 ${\hat{\mathcal{P}}}_D^{-1}\mathcal{A}$ 的特征值落在端点只依赖于 $\beta$ 与多重网格循环质量、而不依赖于 $h$ 的固定区间内。于是预条件 MINRES 迭代以与网格无关的步数收敛——这正是最优求解器的标志。同样的构件装配成块三角预条件子 ${\hat{\mathcal{P}}}_T$ 并配以 GMRES，常被观察到只需大约一半的迭代步数，与定理 6.12、6.13 的比较相符。