5.3 数值分析：误差的艺术

5.3 数值分析：误差的艺术

本节摘要：计算机算数带来两类固有误差——截断误差（方法近似的偏差）与舍入误差（浮点表示的偏差），两者随步长反向变化，存在最优点。本节从浮点表示的现实讲起，区分"问题的病态性"与"方法的稳定性"，用数值微分实验展示两类误差的博弈，最后对比欧拉法与四阶龙格—库塔法的收敛阶，给出方法选择的工程准则。

0.1 加 0.2 不等于 0.3

每个工程师迟早撞上这堵墙：浮点数用二进制表示，0.1 在二进制下是无限循环小数，存进 64 位必须截断。误差微小，但灾难在于两个相近数相减——有效数字大量对消，相对误差暴涨。经典案例如下：

import math

print(0.1 + 0.2 == 0.3)          # False：二进制表示的天然缺陷
print(f"{0.1 + 0.2:.20f}")       # 0.30000000000000004441...

# 相近数相减的灾难：求 x^2 - 1 的根附近函数值的精度
x = 1.0000001
print(f"{x*x - 1:.16f}")         # 有效数字只剩约 7 位（本应有 16 位）
# 对策：代数改写 x^2 - 1 = (x-1)(x+1)，避开减法
print(f"{(x - 1) * (x + 1):.16f}")
# 同样的数学式，改写后精度找回大半——数值分析的第一课：公式不等价于算法
​

教科书公式 (1 - cos x) 除以 x 平方在 x 极小时也如此，改写成半个正弦平方再除（半角公式）即可痊愈。算式的形态决定精度，这是数值分析与数学分析的分界线。

病态与不稳定：两个不同方向的坏

数值坏消息有两个来源，必须区分。问题病态：问题本身对输入扰动敏感，输入变一点输出变很多——条件数衡量这一敏感度，任何算法都救不了病态问题，只能换问法（如正规方程改 QR 分解求解最小二乘）。方法不稳定：问题良好但算法放大误差——换方法就能救。用数值微分看清"博弈"结构：

def central_diff(f, x, h):
    return (f(x + h) - f(x - h)) / (2 * h)

import math
f, x0 = math.sin, 1.0
exact = math.cos(1.0)

print(f"{'h':>10} {'误差':>12}")
for h in [1e-1, 1e-3, 1e-5, 1e-7, 1e-9, 1e-11, 1e-13, 1e-15]:
    err = abs(central_diff(f, x0, h) - exact)
    print(f"{h:>10.0e} {err:>12.2e}")
# 输出形态：误差先随 h 减小而减小（截断误差主导，按 h^2 下降），
# 到 h 约 1e-5 处达到最小，随后反弹上升（舍入误差主导，按 1/h 放大）
# 两类误差的乘性组合存在最优点——中心差分的最优 h 约为机器精度的三分之一次方
​

截断与舍入误差的博弈曲线

收敛阶：方法的价格性能比

数值方法用"步长减半、误差降几倍"分级：一阶方法误差减半，二阶减四分之一，四阶减十六分之一。直觉上阶越高越好，但每步计算量也在涨，同等计算预算下谁划算要实际比。用微分方程数值解做对照实验（5.2 节的方程在此做精度试验台）：

def euler(f, y0, t0, T, n):
    # 一阶欧拉法：用左端点斜率的矩形近似（黎曼和思想，第 3 章）
    h = (T - t0) / n
    y, t = y0, t0
    for _ in range(n):
        y += f(t, y) * h
        t += h
    return y

def rk4(f, y0, t0, T, n):
    # 四阶龙格-库塔：四次斜率采样的加权平均
    h = (T - t0) / n
    y, t = y0, t0
    for _ in range(n):
        k1 = f(t, y)
        k2 = f(t + h/2, y + h*k1/2)
        k3 = f(t + h/2, y + h*k2/2)
        k4 = f(t + h, y + h*k3)
        y += h * (k1 + 2*k2 + 2*k3 + k4) / 6
        t += h
    return y

decay = lambda t, y: -y            # 测试方程：精确解 e^(-T)
T, exact = 1.0, math.exp(-1.0)

for n in [10, 20, 40]:
    e_err = abs(euler(decay, 1.0, 0, T, n) - exact)
    r_err = abs(rk4(decay, 1.0, 0, T, n) - exact)
    print(f"n={n:<4} 欧拉 {e_err:.2e}  RK4 {r_err:.2e}")
# 欧拉：n 翻倍误差减半（一阶）；RK4：n 翻倍误差降为十六分之一（四阶）
# 同样 n=10，RK4 已达 1e-6 量级精度——采样策略的差距是数量级的
​

稳定性：误差别滚雪球

精度之外还有稳定性：某些显式方法对刚性方程（快慢时间尺度并存，如化学动力学）会在某步长之后误差指数放大——不是不够准，是直接爆掉。显式欧拉法的稳定区间要求步长小于 2 除以特征值模长，刚性系统特征值差几个数量级，稳定步长被最快的尺度绑架，慢尺度被迫用极小步长硬磨。对策是隐式方法（每步解一个方程换取无条件稳定），化学、电路仿真清一色用隐式积分器。选型经验浓缩成一张表：

⚠️ 三条战场守则：其一，先查条件数再上算法，病态问题换问法；其二，别在代码里直接比较浮点相等，用容差区间；其三，生产环境用久经考验的求解器库，自己写的积分器只用于理解原理。

本节要点回顾

• 浮点不是实数：相近数相减毁灭有效数字，等价的数学式可以有天差地别的数值形态；
• 截断与舍入误差随步长反向变化，中心差分最优步长约为机器精度开立方；
• 病态与不稳定是两回事：前者是问题属性（条件数），后者是算法属性（换方法可救）；
• 收敛阶是方法的价格性能比：RK4 用四次采样换四阶精度，同等精度下计算量大幅占优；
• 刚性系统必须隐式方法，稳定性而非精度才是刚性场景的瓶颈。

误差已被驯服，最后一环是最优决策：优化理论——线性规划、凸优化与梯度法，下一节收官。