4. 线性方程组求解异常

在计算物理中，大量问题最终都归结为求解形如 A·x = b 的线性方程组，例如有限差分离散后的泊松方程、隐式时间推进中的线性化系统、最小二乘拟合的正规方程等。求解器返回的结果有时看起来"算完了"，但数值完全不可信——这通常不是代码写错了，而是矩阵本身出了问题：要么接近奇异，要么严重病态。本节要解决的核心问题是：如何识别这两类异常，以及如何用正则化手段把不可解或不可信的问题转化为可稳定求解的问题。

一、先分清"奇异"与"病态"

奇异矩阵指行列式严格为零、存在非零向量 v 使 A·v = 0 的矩阵，此时解要么不存在，要么有无穷多个。病态矩阵则行列式不为零、理论上解唯一，但条件数 κ(A) = ‖A‖·‖A⁻¹‖ 极大，输入中极小的扰动（舍入误差、测量噪声）会被放大成解的巨大偏差。二者的关系是：奇异是病态的极端情形，而病态是"几乎奇异"。工程上更常见、也更隐蔽的是病态——因为求解器不会报错，只会安静地给出错误答案。

一个直观的判据是：若 κ(A) ≈ 10^k，而双精度浮点约有 16 位有效数字，则解大约只剩 16 − k 位可信数字。当 k 超过 12～14 时，结果基本失去意义。因此排障的第一步不是急着换求解器，而是先量化矩阵的"坏程度"。

二、检测：三个层次的诊断手段

第一层是求解器自身的信号。直接法（LU、Cholesky）在遇到零主元时会失败或给出警告；迭代法（CG、GMRES）则表现为残差停滞不降、达到最大迭代数仍未收敛。注意：残差小不等于解准确。对病态系统，‖b − A·x‖ 可以很小而 ‖x − x_true‖ 很大，这是最容易被误判的陷阱。

第二层是条件数估计。不必显式求逆，用 numpy.linalg.cond(A) 或 LAPACK 的 dgecon 即可在 O(n³) 内给出 1-范数条件数的估计。对大规模稀疏矩阵，可用 scipy.sparse.linalg.onenormest 配合近似逆来低成本估计。

第三层是残差回代与解扰动检验。计算相对残差 ‖b − A·x‖/‖b‖，再人为给 b 加一个 1e-10 量级的相对扰动重解一次，比较两次解的相对差异。若解的相对变化远大于扰动量（例如扰动 1e-10 导致解变化 1e-2），则病态确认无疑。

三、正则化的实际做法

对病态系统，最常用的手段是 Tikhonov 正则化：把 A·x = b 替换为最小化 ‖A·x − b‖² + λ‖x‖²，等价于求解

(AᵀA + λI)·x = Aᵀb

其中 λ > 0 是正则化参数。它的作用是把 AᵀA 的最小特征值从接近 0 抬到至少 λ，从而把条件数压到可控范围。代价是引入偏差：λ 越大解越平滑、越稳定，但偏离真实解越多。因此 λ 的选择是核心，常用 L 曲线法（在 log‖A·x−b‖ 对 log‖x‖ 的曲线上取拐点）或广义交叉验证（GCV）。

另一种手段是截断奇异值分解（TSVD）：对 A = UΣVᵀ，只保留前 k 个较大的奇异值，把 σ_i < ε 的项直接丢弃。这相当于把解投影到数值稳定的子空间上，适合奇异值谱有明显"断崖"的问题。

下面给出一个可直接运行的检测与正则化示例，用 Hilbert 矩阵（经典病态矩阵）构造问题：

import numpy as np

n = 10
A = np.array([[1.0/(i+j+1) for j in range(n)] for i in range(n)])
x_true = np.ones(n)
b = A @ x_true

# 1) 诊断
kappa = np.linalg.cond(A)
print("cond(A) =", kappa)          # 约 1.6e13，严重病态

x_direct = np.linalg.solve(A, b)
res = np.linalg.norm(b - A @ x_direct) / np.linalg.norm(b)
err = np.linalg.norm(x_direct - x_true) / np.linalg.norm(x_true)
print("直接解: 相对残差=%.2e, 相对误差=%.2e" % (res, err))
# 典型输出：残差 ~1e-16，误差 ~1e0 —— 残差极小但解完全错误

# 2) Tikhonov 正则化
lam = 1e-6
x_reg = np.linalg.solve(A.T @ A + lam*np.eye(n), A.T @ b)
err_reg = np.linalg.norm(x_reg - x_true) / np.linalg.norm(x_true)
print("正则化解: 相对误差=%.2e" % err_reg)   # 误差显著下降

# 3) TSVD
U, s, Vt = np.linalg.svd(A)
k = 6
x_tsvd = Vt[:k].T @ ((U[:, :k].T @ b) / s[:k])
print("TSVD 解: 相对误差=%.2e" % (np.linalg.norm(x_tsvd - x_true)/np.linalg.norm(x_true)))

这段代码刻意展示了那个最危险的陷阱：直接解的相对残差只有 1e-16 量级，看起来"完美收敛"，但相对误差却是 O(1)。如果只看残差就下结论，排障必然失败。加入 λ = 1e-6 的 Tikhonov 项后，误差通常能降到 1e-2 量级；TSVD 取 k = 6 也能得到相近效果。具体数值随 n 和 λ 变化，这里给出的是量级判断，实际使用时应用 L 曲线扫描 λ 从 1e-12 到 1e-2 的对数网格，观察误差与残差的权衡拐点。

四、工程上的补充手段

除了正则化，还有两类常用策略。其一是预处理：对迭代法使用不完全 LU 分解（ILU）或对角缩放作为预条件子，降低有效条件数，让 CG/GMRES 真正收敛。其二是变量缩放与方程平衡：当矩阵各行量纲差异巨大（如同时含 1e-9 和 1e6 量级的系数）时，先做行、列均衡（如 scipy.sparse.linalg.spsolve 内部的 equilibration），往往能把条件数降低几个数量级，且不引入任何偏差。

最后要建立一条固定的验证习惯：任何线性求解完成后，至少做三件事——计算相对残差、估计条件数、对右端项施加 1e-10 相对扰动重解并比较解的相对变化。三者结合，才能把"算完了"和"算对了"区分开。对病态问题，还应把正则化参数 λ 与最终误差一并记录，便于复现和后续调参。