多重网格法

多重网格法（multigrid method）在不同尺度的网格之间传递误差信息，组合局部平滑与粗网格校正来求解离散方程。它针对一个常见现象：简单迭代能迅速消除网格上来回起伏的误差，却很难消除缓慢变化的误差；后者放到较粗网格上，会变成更容易处理的变化。

三个内部点上的粗网格校正

考虑 $- u^{″} = 2$ 、 $u (0) = u (1) = 0$ 。真解是 $u = x (1 - x)$ 。细网格间距 $h = 1 / 4$ ，三个内部点上的离散解也是 $(3 / 16, 1 / 4, 3 / 16)$ ，因为二次函数的中心二阶差分正好精确。

从零近似出发，残差为 $r = (2, 2, 2)$ 。用权重 $(1 / 4, 1 / 2, 1 / 4)$ 限制到粗网格唯一的中点，得到粗残差二。粗间距为 $H = 1 / 2$ ，误差方程是 $2 e_{H} / H^{2} = 2$ ，故 $e_{H} = 1 / 4$ 。线性插回细网格给出 $(1 / 8, 1 / 4, 1 / 8)$ 。校正后剩余误差为 $(1 / 16, 0, 1 / 16)$ ：粗层已消除平滑的主要部分，但仍需要细层平滑处理剩余起伏。

细网格精确离散解与粗层线性插值校正比较，中点完全校正而相邻点仍有剩余误差 — 从零近似出发的单次粗层校正。粗网格只有中点一个未知量，线性延拓不能完全表示细网格的二次曲线。

平滑器对不同误差频率的作用

对一维 Poisson 离散 $A_{h} = h^{- 2} tridiag (- 1, 2, - 1)$ ，加权 Jacobi 迭代为 $u^{n e w} = u + ω D^{- 1} (b - A_{h} u), D = 2 h^{- 2} I .$ 在无边界或局部 Fourier 分析中，误差模式 $e_{i} = e^{i θ i}$ 的放大因子是 $g (θ) = 1 - ω + ω \cos θ = 1 - 2 ω \sin^{2} (θ / 2) .$ 取 $ω = 2 / 3$ ，高频区间 $π / 2 \leq | θ | \leq π$ 上有 $| g | \leq 1 / 3$ ，一轮即可显著减小振荡；低频 $θ \to 0$ 时 $g \to 1$ ，衰减很慢。平滑器的任务因此是消除某类误差，不是单独把整个系统完全解好。

残差、限制与校正

设当前近似为 $u$ ，精确离散解为 $u_{*}$ ，误差 $e = u_{*} - u$ 。由 $A_{h} u_{*} = b$ 可得误差方程 $A_{h} e = b - A_{h} u = r .$ 虽然不知道误差，却知道右侧残差。将其限制到粗层，解 $A_{H} e_{H} = R r$ ，再作 $u \leftarrow u + P e_{H}$ 。 $R$ 是限制算子， $P$ 是延拓算子；加号来自这里误差定义为“精确减近似”。若定义反过来，校正符号也应一起改变。

对称正定问题中，取 $A_{H} = P^{T} A_{h} P$ 、 $R = P^{T}$ ，精确粗解使校正后误差满足 $P^{T} A_{h} (e - P e_{H}) = 0.$ 因此新误差在 $A_{h}$ 内积下与粗空间正交。将 $e = (e - P e_{H}) + P e_{H}$ 展开平方，交叉项为零，得到 $‖ e ‖_{A_{h}}^{2} = ‖ e - P e_{H} ‖_{A_{h}}^{2} + ‖ P e_{H} ‖_{A_{h}}^{2} .$ 粗层校正不会增加能量误差，并精确移除粗空间中能够表示的部分。实际常用限制与 $P^{T}$ 只差一个比例；同时按比例定义粗算子，得到同一校正。Netlib 的多重网格资料介绍了残差校正和层级循环。^[1]

V 循环与工作量

一次 V 循环在当前层作少量前平滑，将残差送到更粗一层递归求解，延拓校正后再作后平滑。粗到最小层时直接求解。若每层点数按固定比例减少、每层只递归一次，工作量满足 $T (M) = c M + T (M / 2)$ ，几何级数求和得到 $O (M)$ 。这是一轮循环的代价；总求解代价还取决于是否能获得与网格无关的收敛因子。

强各向异性、复杂系数或不合适的粗层可能使简单点 Jacobi 失效。线平滑、半粗化和代数多重网格分别针对不同结构。不能仅凭“有几个网格层”就保证线性复杂度或理想收敛。

一维 Poisson 的可运行 V 循环

代码使用 $2^{ℓ} - 1$ 个内部点、零端点、全加权限制、线性延拓及两次前后 Jacobi 平滑。传入右端和当前近似，返回新列表，不修改输入。一次循环需线性时间和线性总层级存储；这里重建列表以便展示，未作内存优化。

Python 3

import math


def vcycle_poisson(rhs, initial, h):
    b, u = list(map(float, rhs)), list(map(float, initial))
    n = len(b)
    if not n or len(u) != n or ((n+1) & n) != 0:
        raise ValueError("length must be 2**level - 1")
    h = float(h)
    if not math.isfinite(h) or h <= 0 or not all(math.isfinite(x) for x in b+u):
        raise ValueError("finite data and positive h required")

    def cycle(b, u, h):
        if len(b) == 1:
            return [b[0]*h*h/2]
        def smooth(v):
            padded = [0.0] + v + [0.0]
            return [(1/3)*v[i] + (1/3)*(h*h*b[i]+padded[i]+padded[i+2])
                    for i in range(len(v))]
        for _ in range(2):
            u = smooth(u)
        padded = [0.0]+u+[0.0]
        residual = [b[i]-(2*u[i]-padded[i]-padded[i+2])/(h*h)
                    for i in range(len(u))]
        coarse = [(residual[2*j]+2*residual[2*j+1]+residual[2*j+2])/4
                  for j in range((len(u)-1)//2)]
        correction = cycle(coarse, [0.0]*len(coarse), 2*h)
        extended = [0.0]+correction+[0.0]
        for i in range(len(u)):
            u[i] += (correction[(i-1)//2] if i % 2 else
                     (extended[i//2]+extended[i//2+1])/2)
        for _ in range(2):
            u = smooth(u)
        return u

    answer = cycle(b, u, h)
    if not all(math.isfinite(v) for v in answer):
        raise OverflowError("nonfinite multigrid state")
    return answer


n = 7
h = 1/(n+1)
u = [0.0]*n
for _ in range(8):
    u = vcycle_poisson([2.0]*n, u, h)
exact = [(i*h)*(1-i*h) for i in range(1, n+1)]
print(max(abs(a-b) for a, b in zip(u, exact)))  # small residual error

有解变式将右端改为零，并把任意初始误差设成全部可由粗层延拓表示的向量 $u = P w$ 。若不作平滑、只作一次精确 Galerkin 粗层校正，右端为 $- P^{T} A_{h} P w$ ，粗误差为 $- w$ ，校正后恰为零。它直接检验了符号、粗算子和延拓是否一致。

参考资料

↑ Richard Barrett et al. Templates for the Solution of Linear Systems，Multigrid Methods 与粗网格校正。

[netlib-1] Richard Barrett et al. Templates for the Solution of Linear Systems，Multigrid Methods 与粗网格校正。

[1]