Jacobi 一万次也消不掉的误差,V-循环十次就解决 —— 多重网格的真正秘密

为什么 Jacobi 跑了一万次还结束不了#

工程师 B 在 1024×1024 的网格上用点 Jacobi 解压力 Poisson。两百秒过去，残差还在同一个数量级里抖。CFL 很小，网格正交，代码也没错。有些误差 Jacobi 永远消不掉。

同一个误差，换一层网格，频率就变了。V-循环是把这句话做成了步骤。

把 Pre/Post 在 1 和 5 之间拨，同一张 64 单元网格、同一个零右端项上，只有 V-循环的残差会更陡，多出来的松弛很快就不划算。

上面的橙色是 Jacobi。振荡变小，大鼓包还在。青色的对数残差几乎沿直线下降，橙色停在低频的平台上。

误差的频率 —— 谁先死,谁活到最后

设线性系统是 $A u = f$。当前猜测 $u^{(k)}$ 的误差是

一个向量。残差是 $r^{(k)} = f - A u^{(k)} = A e^{(k)}$。

加权 Jacobi 的误差传播里，傅里叶模 $\sin(j \pi x)$ 每扫一次变成

$\omega$ 是松弛系数，$j$ 是模的编号，$h$ 是网格间距。衰减系数的绝对值如下。

短波一次缩到三分之一。长波几乎不动。Jacobi 消不掉的，是铺满整张网格的缓坡。

在更粗的网格上,那段低频就是高频

多重网格在这里和普通松弛分开。$N$ 个单元上最长的波长是 $\lambda \sim N h$。搬到 $N/2$ 个单元上，同一段波只占一半的单元，看起来更短。说准确一点，细网格上最低的模，在粗网格上是中频，而那里的 Jacobi 擅长消中频。

再粗化一层，波长再短一截。下到 $\log_2 N$ 层，每个模总会在某一层变成高频。每层扫几次 Jacobi，理论上所有误差都会死掉。把低频留给迭代法去对照的记录是 从 Jacobi 到 BiCGSTAB 的收敛。

V-循环 —— 下去再上来的路#

把这个判断写成步骤，形状就是 V。

1. 预松弛 (pre-smoothing) — 细网格上 Jacobi $\nu_1$ 次。
2. 残差 — $r_h = f_h - A_h u_h$。
3. 限制 (restriction) — $r_{2h} = R\, r_h$。残差拿到粗网格。
4. 递归 — 在粗网格上用 V-循环解 $A_{2h} e_{2h} = r_{2h}$。网格够小就直接解。
5. 延拓 (prolongation) — $e_h = P\, e_{2h}$。修正抬回细网格。
6. 修正 — $u_h \leftarrow u_h + e_h$。
7. 后松弛 (post-smoothing) — 细网格上 Jacobi $\nu_2$ 次。

下去只带残差，上来只带修正。一个循环的花费是 $O(N)$：细网格扫一次，再加上几何级数变小的那些层。一个循环大约掉一位残差。

限制与延拓 —— 两层网格之间的汇率

两层之间的算子必须成对。一维顶点中心网格上最简单的一对是这样。

全加权限制: 粗点 $i$ 的残差是细点 $2i-1, 2i, 2i+1$ 的加权平均。

线性延拓: 对准粗点的位置原样拷贝，中间点取两个邻居的平均。

两者互为转置，差一个网格比常数($P = c R^T$)。V-循环的收敛论证里有这一对。不对称的一对可能发散。

用 Python 走一遍两层 V-循环#

用 V-循环解一维 Poisson $-u''(x) = f(x)$，$u(0)=u(1)=0$。从 64 个单元开始。

import numpy as np
 
class GridLevel:
    """一个网格层的状态:解、右端项、网格步长。"""
    def __init__(self, n):
        self.n = n              # 单元数
        self.h = 1.0 / n
        self.u = np.zeros(n + 1)
        self.f = np.zeros(n + 1)
 
def jacobi_relax(u, f, h, iters):
    """1D Poisson 的加权 Jacobi。两端 Dirichlet 零边界。"""
    n = len(u) - 1
    nxt = np.zeros_like(u)
    for _ in range(iters):
        nxt[1:n] = 0.5 * (u[0:n-1] + u[2:n+1] + h * h * f[1:n])
        u[1:n] = nxt[1:n]
    return u
 
def residual_norm(u, f, h):
    n = len(u) - 1
    r = np.zeros_like(u)
    inv = 1.0 / (h * h)
    r[1:n] = f[1:n] - (-u[0:n-1] + 2 * u[1:n] - u[2:n+1]) * inv
    return float(np.sqrt(np.mean(r[1:n] ** 2)))
 
def restrict_full_weighting(rh):
    n = len(rh) - 1
    n2 = n // 2
    r2 = np.zeros(n2 + 1)
    r2[1:n2] = 0.25 * rh[1:n-1:2] + 0.5 * rh[2:n:2] + 0.25 * rh[3:n+1:2]
    return r2
 
def prolong_linear(e2h):
    n2 = len(e2h) - 1
    n = 2 * n2
    eh = np.zeros(n + 1)
    eh[2:n:2] = e2h[1:n2]
    eh[1:n+1:2] = 0.5 * (e2h[0:n2] + e2h[1:n2+1])
    return eh
 
def v_cycle(u, f, h, n_pre=2, n_post=2):
    n = len(u) - 1
    if n <= 4:
        return jacobi_relax(u, f, h, 50)
    jacobi_relax(u, f, h, n_pre)
    inv = 1.0 / (h * h)
    r = np.zeros_like(u)
    r[1:n] = f[1:n] - (-u[0:n-1] + 2 * u[1:n] - u[2:n+1]) * inv
    r2 = restrict_full_weighting(r)
    e2 = np.zeros_like(r2)
    v_cycle(e2, r2, 2 * h, n_pre, n_post)
    e = prolong_linear(e2)
    u += e
    jacobi_relax(u, f, h, n_post)
    return u
 
# --- 运行 ---
n = 64
h = 1.0 / n
x = np.linspace(0, 1, n + 1)
f = (np.pi ** 2) * np.sin(np.pi * x)   # 精确解 u = sin(πx)
 
u_jac = np.zeros(n + 1)
u_mg = np.zeros(n + 1)
 
for k in range(20):
    jacobi_relax(u_jac, f, h, 4)
    v_cycle(u_mg, f, h, n_pre=2, n_post=2)
    print(f"step {k:3d} | jacobi r={residual_norm(u_jac, f, h):.3e} "
          f"v-cycle r={residual_norm(u_mg, f, h):.3e}")

一个 step 里 Jacobi 做 4 次，V-循环在细网格上也是 4 次(pre 2 + post 2)。花费接近，残差差很远。20 步之后，V-循环残差通常低于 $10^{-10}$，Jacobi 停在 $10^{-2}$ 附近。

表里的衰减在 $\omega = 2/3$ 上闭合。最快的模 $jh = 1$，$\cos(\pi) = -1$，系数是 $1 - 2\omega = -1/3$。绝对值就是表上的 $1/3$。只扫这个模十次，会掉到 $(1/3)^{10} \approx 1.7 \times 10^{-5}$。中频是表上的 $\approx 0.67$，把 0.67 乘 20 次得到 $0.67^{20} \approx 3.3 \times 10^{-4}$。最慢的因子是 $1 - O(h^2)$。$N = 64$，$h = 1/64$，于是 $h^2 = 1/4096 \approx 2.4 \times 10^{-4}$。一次能削掉的就是这个量级。

这和输出对得上。每个 step 的 Jacobi 是 4 次，20 个 step 就是 80 次。残差仍停在 $10^{-2}$ 附近。高频按 $1/3$ 衰减，在这之前已经死了。留下的 $10^{-2}$ 是表里的低频。同样 20 个 step 的 V-循环，细网格上 pre 2 次、post 2 次，残差通常在 $10^{-10}$ 以下。$10^{-2}$ 和 $10^{-10}$ 隔着八个数量级。

最粗的一层是 $n \le 4$ 时 Jacobi 50 次。64 个单元按 64、32、16、8、4 限制四次。全加权是 $1/4$、$1/2$、$1/4$。粗层花费之和是 $N/2 + N/4 + \cdots$，不超过 $N$。十次 V-循环是把这个和再乘十。一万次 Jacobi 在次数上就已经不是同一个量。

把同一式子放进最慢的模。$N = 64$、$jh = 1/N$、$\omega = 2/3$ 时，衰减大约是 $0.9992$。清单里的 80 次乘下去大约是 $0.94$。削掉一位大约要 2900 次。把同一式子放进 $N = 1024$，衰减大约是 $0.999997$，一位大约要 $7.3 \times 10^{5}$ 次。乘上一万次，这个模仍留在 $0.97$ 附近。1024×1024 上残差在 200 秒里停在同一个数量级，是同一个地方。这个式子是一维加权 Jacobi。工程师的网格是二维。一个循环掉一位的 V-循环，大约十次就能推进那些扫掠推不动的慢模。

$0.9992$ 不是测出来的残差，是把表里的式子放进最慢的模得到的。乘 80 次，大约只动百分之六。20 步之后 Jacobi 残差停在 $10^{-2}$ 附近，是衰减为 $1/3$ 的快模先死掉以后，只剩一个接近 1 的系数时的样子。$N = 1024$ 时同一个系数在一万次之后仍约 $0.97$，数量级不动。V-循环一轮大约掉一位，和这个接近 1 的系数不在同一个时间尺度上。若每轮掉一位，十轮大约是十位，和 20 步后残差低于 $10^{-10}$ 的印刷结果是同一种说法。粗层花费 $N/2+N/4+\cdots$ 不超过 $N$，十轮仍是细网格工作量的常数倍，不是一万次细网格扫描。$n \le 4$ 上的 50 次 Jacobi 只出现在 64、32、16、8、4 这条链的末端。全加权的 $1/4$、$1/2$、$1/4$ 不改变这个费用的几何级数。

它也不是总是快 —— 两个陷阱

V-循环不是魔法。实务里常断的地方有两处。

粗糙的系数。 各向异性或系数跳跃的 Poisson 上，Jacobi 不是合适的 smoother。在 $-\nabla \cdot (k(x) \nabla u) = f$ 里，$k$ 按数量级跳，标准 V-循环会把残差停在某一位。出路是 line-Jacobi、沿线排列的网格，或者代数多重网格(AMG)，层次从矩阵而不是从网格长出来。

不对齐的网格。 非结构网格或 AMR 上，细到粗不是把下标对半砍。聚团多重网格把细单元捆成粗单元，限制用体积加权，延拓用单元中心插值。实现很重，OpenFOAM 的 GAMG 走的就是这条路。

收敛率要和网格无关，smoother、算子、粗化比三样都得对上。缺一个，迭代次数就跟着 $N$ 慢慢涨。那时再看 W-循环或 F-循环。

标准 V-循环放到哪里为止#

$k$ 按数量级跳的时候，普通 Jacobi 不是 smoother，残差会粘在某一位上。换成 line-Jacobi 或 AMG。拉长边界层上的这种替换写在 线状预处理。

非结构网格和 AMR 不能靠把下标砍半来做粗层。OpenFOAM 的 GAMG 把单元捆起来。迭代次数开始跟着 $N$ 走，就试 W-循环或 F-循环。64 个单元、$n \le 4$ 上 Jacobi 50 次、20 步里 $10^{-10}$ 对 $10^{-2}$，只是在对齐的一维 Poisson 上读到的数。