3.3 实战：大规模稀疏方程组求解

3.3 实战：大规模稀疏方程组求解

本节摘要：把二维泊松方程用五点差分离散成线性方程组，得到五对角稀疏矩阵。同一台机器上，稠密路线用 numpy.linalg.solve，稀疏路线用 spsolve 与 cg 加预条件，网格从 50 乘 50 推到 1000 乘 1000，实测内存与耗时。你会看到稠密路线在哪个规模点崩溃，稀疏路线如何一路跑到百万未知数，并验证解的误差按步长平方收敛。本节还给出可直接复用的误差检验函数与求解路线选型对照表，是偏微分方程数值求解的起步模板，稍加改动即可迁移到热传导等方程。

本节目标

阅读完本节，你应当能够：

1. 把二维泊松方程用五点差分离散成五对角线性方程组；
2. 用 diags 一次构造五对角稀疏矩阵，并说明每条对角线的含义；
3. 用内存与耗时两个维度对比稠密与稀疏求解路线，解释差距的来源；
4. 配置 cg 加不完全 LU 预条件，并观察收敛步数的变化；
5. 用已知解析解检验数值解，验证二阶收敛。

一、问题与直觉

前面两节的工具都认识了，但散着用和合起来用是两回事。这一节我们做一件具体的事：解一个二维泊松方程。它是热传导、流体、电磁场模拟里最常见的方程之一，选它当案例有三个理由：离散后矩阵是典型的稀疏结构；问题本身有解析解，误差可以精确检验；构造过程只依赖本章的知识，不需要额外背景。

问题设定如下。单位正方形区域内，未知函数 u 的拉普拉斯算子取负，等于已知函数 f，边界上 u 恒为零。我们取 f 为 2 派平方乘 sin 派 x 乘 sin 派 y，那么真解就是 sin 派 x 乘 sin 派 y。有了真解，数值解准不准一目了然——这是检验数值方法的黄金手段。

边界取零是刻意为之：齐次边界让离散矩阵的结构干净，对角占优也有保证，我们能把注意力全部放在求解本身。等你会跑了，把右端项换成自己的数据、把边界条件改成非零，代码框架一个字都不用动。

先跑通一条最小链路，把整个流程走一遍，再上规模。

import numpy as np
from scipy.sparse import diags
from scipy.sparse.linalg import spsolve

N = 50                       # 网格内部点数，未知数共 N 的平方
h = 1.0 / (N + 1)

# 五点差分矩阵：主对角 4，四条次对角负一
offsets = [0, 1, -1, N, -N]
values = [4, -1, -1, -1, -1]
A = diags(values, offsets, shape=(N*N, N*N), format='csr') / (h * h)

# 右端项：f 在内部网格点上的取值
x = np.linspace(0, 1, N + 2)[1:-1]
X, Y = np.meshgrid(x, x)
F = 2 * np.pi**2 * np.sin(np.pi * X) * np.sin(np.pi * Y)

u = spsolve(A, F.ravel()).reshape(N, N)
print(u[N // 2, N // 2])     # 中心值，理论上是 1.0
​

跑完打印中心点的值，和理论值 1.0 对照，量级对得上就说明链路是通的。

二、核心原理

2.1 从连续方程到线性方程组

离散的思路是把连续区域切成网格。内部每个点上的二阶导数用相邻点的差分近似：二阶导约等于两边值减两倍中间值，再除以步长平方。把网格按行拉直成一维索引，每个未知数只出现在自己、上下、左右五个位置上的方程里——这就是五对角矩阵的来源。

写清楚一点：在第 i 行第 j 列的网格点上，方程左边是四倍的自身值减去上邻居、下邻居、左邻居、右邻居的值，再除以步长平方，右边是 f 在该点的值。边界点上的未知数因为边界为零而直接消失，所以内部点之间只有五条连接，矩阵的带宽也只有五条对角线。

主对角元素是四除以步长平方，四条次对角是负一除以步长平方。步长越小，对角占优越强，矩阵性质越好，但未知数按平方级增长，规模压力随之而来。

2.2 五对角矩阵的构造

diags 的 offsets 参数一次搞定：0 是主对角，正负一是上下邻居，正负 N 是隔行邻居。构造过程零循环，矩阵本身不占额外内存。非零个数约五乘 N 平方，而稠密版本要 N 的四次方个元素——N 等于 1000 时，稀疏约五百万非零，稠密是十万亿个浮点数，差出五个数量级。上一节讲过 diags 的用法，这里直接看出它的价值：一条语句，一个五对角矩阵。

矩阵搭好后值得花一行代码确认它的性质：对称性用矩阵减转置的最大值检查，正定性在小规模上算几个特征值验证。泊松矩阵天然对称正定，cg 因此可用；确认这一步，后面选求解器才不会选错。顺带的好处是，这种规则的带状结构让稀疏直接法的分解填充非常有限，所以 N 在一两百时 spsolve 依旧又快又稳，不必急着上迭代法。

图：稀疏与稠密存储内存对比

2.3 两条路线的实测对比

用下面这段代码分别在 N 等于 50、100、200 时跑两条路线：稠密路线先 toarray 再用 numpy.linalg.solve，稀疏路线直接 spsolve。

import time

def dense_route(A, b):
    Ad = A.toarray()
    t0 = time.perf_counter()
    x = np.linalg.solve(Ad, b)
    return x, time.perf_counter() - t0, Ad.nbytes

def sparse_route(A, b):
    t0 = time.perf_counter()
    x = spsolve(A, b)
    return x, time.perf_counter() - t0, A.data.nbytes + A.indices.nbytes + A.indptr.nbytes
​

我按常见笔记本配置把趋势列在下面，具体数字会因机器而异，但相对关系稳定：

这张表不是让你背数字，而是建立三个直觉：稠密内存按未知数平方涨，稀疏按非零个数线性涨；稠密求解耗时按未知数三次方涨，稀疏按迭代次数乘非零个数涨；两条路线的交叉点大概在万阶附近，过了这个点再犹豫就是白烧内存。

注意规律：网格每放大一倍，稠密路线的内存涨四倍，求解耗时涨约六十四倍——因为未知数是 N 的平方，稠密求解复杂度是未知数的三次方；而稀疏路线内存只按非零个数线性涨，迭代法每步成本与非零个数成正比。这就是大规模数值计算必须稀疏化的根本原因。

把复杂度写全：稠密路线是 N 的六次方量级，稀疏路线是 N 的平方乘以迭代次数。N 从 100 涨到 1000，稠密路线涨一千万倍，稀疏路线涨一百倍出头——前提是迭代次数不随规模恶化，这正是预条件要守住的东西。

2.4 迭代求解与预条件

把 N 推到 1000 时，spsolve 的直接分解开始吃力，我们换 cg 迭代。先不做预条件，再挂上 spilu，对比收敛步数：

from scipy.sparse.linalg import cg, spilu

b = F.ravel()
x0, info0 = cg(A, b, tol=1e-8, maxiter=5000)          # 无预条件

M = spilu(A)
Mx = lambda v: M.solve(v)
x1, info1 = cg(A, b, M=Mx, tol=1e-8, maxiter=5000)    # 带预条件

cg 的返回值里，info 等于 0 表示收敛，非零表示撞上最大迭代次数。这个返回值值得每次检查——迭代法失败是静默的，解照样给你，只是可能差得离谱。
​

经验上，同样的容差，无预条件可能要数千步，带预条件几十步就收敛。原因在于预条件把系数矩阵的条件数大幅压低，而 cg 的收敛速度直接受条件数支配。迭代法省内存的代价是收敛不确定性，预条件正是用来买回这份确定性的。

容差的选择也有讲究：tol 设到 1e-6 对多数工程问题够用，设到 1e-12 每步收敛都会变慢，最后几位数字的精度未必值得等。gmres 还有 restart 参数，内部每迭代若干步重启一次基向量，防止存储无限增长；默认值在多数问题里不用动，真遇到不收敛再试着调小。

2.5 误差检验：二阶收敛

最后检查解的质量。真解已知，数值解与真解的最大绝对误差随网格加密的变化规律是：步长减半，误差缩到约四分之一，这叫二阶收敛。验证方法：分别算 N 等于 50 和 100 的最大误差，比值应该在 4 附近。

def max_error(N):
    h = 1.0 / (N + 1)
    A = diags([4, -1, -1, -1, -1], [0, 1, -1, N, -N],
              shape=(N*N, N*N), format='csr') / h**2
    x = np.linspace(0, 1, N + 2)[1:-1]
    X, Y = np.meshgrid(x, x)
    F = 2 * np.pi**2 * np.sin(np.pi * X) * np.sin(np.pi * Y)
    u = spsolve(A, F.ravel()).reshape(N, N)
    u_exact = np.sin(np.pi * X) * np.sin(np.pi * Y)
    return np.max(np.abs(u - u_exact))

e50 = max_error(50)
e100 = max_error(100)
print(e50, e100, e50 / e100)   # 比值应接近 4
​

如果比值明显偏离 4，说明程序里有 bug——网格没对齐、边界处理错了、或者右端项采样位置不对。比值落在 3 到 5 之间都算正常，偏离一个数量级以上就必须查。误差检验是数值实验里最划算的体检，比盯着残差可靠得多——残差小只能说明方程被满足了，误差小才说明方程本身建对了。

三、工程实践要点

3.1 路线选型对照

这套流程稍加改动就能迁移到别的方程：热传导方程做隐式欧拉，每个时间步解的矩阵是单位阵减时间步长乘拉普拉斯矩阵，五对角结构不变；换成诺伊曼边界，只需把边界行的主对角与邻居系数改一改，骨架依然成立。学会一个泊松问题，等于学会半类椭圆与抛物方程的离散求解套路。

3.2 常见问题

问：为什么 N 等于 200 时稠密路线直接被杀？
答：4 万未知数的稠密矩阵要 12.8 GB，超过多数机器的物理内存后操作系统开始换页，程序要么极慢，要么被内存保护机制终止。这不是代码问题，是存储方案的必然结果。

问：迭代法每次算出来的解不一样？
答：cg 对对称正定矩阵在浮点误差范围内是确定的；gmres 的收敛轨迹依赖初始猜测。给迭代器固定初始向量与容差，结果可复现。

问：为什么先跑通再上规模的顺序这么重要？
答：小网格上矩阵小、求解快、误差可检验，任何 bug 都会在误差比值上现形；直接上大网格，内存和时间成本把排查次数限制到可怜的一两次。先小后大，是数值实验里性价比最高的习惯。

⚠️ 常见坑：对稀疏矩阵调用 numpy 的 solve 或 inv，numpy 会尝试把输入转成稠密数组，小矩阵看不出问题，大矩阵直接内存耗尽。判断一个操作是否安全，就看它是否要求稠密输入。
💡 关键直觉：网格加密一倍，稠密路线的成本涨约六十四倍，稀疏路线只涨约四倍。这个数量级差异，决定了"先离散成稀疏结构，再用迭代法求解"是大型数值模拟唯一现实的路线。

核心回顾

• 要点一：五点差分把泊松方程离散成五对角方程组，主对角与四条次对角分别对应自身和上下左右邻居。
• 要点二：diags 用偏移量一次构造，零循环，五对角结构一目了然。
• 要点三：稠密存储的规模墙出现在万阶附近，内存按未知数平方增长。
• 要点四：稀疏路线把内存压到非零个数量级，百万未知数约 80 MB。
• 要点五：网格放大一倍，稠密求解耗时涨约六十四倍，稀疏路线只按比例增长。
• 要点六：十万阶以上优先迭代法，spilu 预条件把收敛步数从上千压到几十。
• 要点七：用解析解做误差检验，步长减半误差缩到四分之一，验证二阶收敛。

线性代数与稀疏计算到这里就齐了。第 4 章我们进入信号与图像处理，傅里叶变换与滤波器会用到本章的矩阵直觉，但视角从解方程转向频率。