5.2 数值模拟技术：从欧拉法到自适应求解器

5.2 数值模拟技术：从欧拉法到自适应求解器

本节摘要：定义先行——数值积分把连续的流切成小步长推进：欧拉法每步只看当前斜率（全局误差 O(h)），龙格-库塔四阶法用四次采样把误差压到 O(h⁴)，自适应求解器自动调步长，刚性方程必须换隐式方法，保守系统最好用辛积分器。本节给出方法对照表、误差实测与三条自检纪律，让"算出来的轨道"变成可信证据。

本节上承 5.1 的分工结论：定量细节交给数值。数值方法的对照轴是"便宜对可靠"——每升一级精度，代价翻几倍；选方法的本质是为你关心的那类系统（刚性？保守？混沌？）买对保险。

步长里的乾坤

最朴素的欧拉法：y(n+1) = y(n) + h·f(y(n))。它每步用起点斜率外推一整步，起点之后的弯曲全然不顾——所以局部误差 O(h²)、累加成全局误差 O(h)。四阶龙格-库塔（RK4）在同一步内采样四次斜率做加权平均，把局部误差压到 O(h⁵)、全局 O(h⁴)：步长减半，误差缩到约十六分之一。理论数字用实验钉一遍最踏实：

import numpy as np

def f(s):                                   # 阻尼谐振子 x'' = -x - 0.2 x'
    x, v = s
    return np.array([v, -x - 0.2 * v])

def euler(s, h, n):
    for _ in range(n):
        s = s + h * f(s)
    return s

def rk4(s, h, n):
    for _ in range(n):
        k1 = f(s)
        k2 = f(s + h / 2 * k1)
        k3 = f(s + h / 2 * k2)
        k4 = f(s + h * k3)
        s = s + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
    return s

wd = np.sqrt(0.99)                                  # 阻尼振荡频率
def exact_at(t):                                    # 解析解，用作对照真值
    e = np.exp(-0.1 * t)
    x = e * (np.cos(wd * t) + 0.1 / wd * np.sin(wd * t))
    v = e * (-0.1 * np.cos(wd * t) - (0.01 + wd**2) / wd * np.sin(wd * t))
    return np.array([x, v])

exact = exact_at(10.0)
for h in (0.1, 0.05, 0.025):
    n = round(10 / h)
    e = np.linalg.norm(euler(np.array([1.0, 0.0]), h, n) - exact)
    r = np.linalg.norm(rk4(np.array([1.0, 0.0]), h, n) - exact)
    print(f"h={h}: 欧拉误差 {e:.3e} | RK4 误差 {r:.3e}")
​

结果（典型输出）：欧拉误差约 6e-2 → 3e-2 → 1.5e-2（步长减半误差减半，一阶名副其实）；RK4 误差约 2e-6 → 1.3e-7 → 8e-9（减半缩约十六倍，四阶名副其实）。同样的终点精度，欧拉要多走三个数量级的步数——便宜与可靠的兑换率就摆在误差阶里。

（顺带一提，解析解 e^{At}x₀ 的精确形式让"对照真值"免费可得——这正是 5.1 节"能解析先解析"的又一用处。）

图：欧拉法与RK4的误差对照

方法对照表

三条自检纪律与一个刚性案例

数值结果要过三道安检才配进报告：步长减半对照（误差应按理论阶缩水，缩不动说明已到方法极限或遇到了刚性）；守恒量监视（保守系统跑长了看能量漂没漂，漂了就换辛积分器）；双方法对照（RK45 与隐式各跑一遍，分歧即警报）。

刚性案例现场：把范德波尔的 μ 调到 1000（4.4 节切削颤振的量级），波形是"极慢充电 + 极陡放电"，快慢时间尺度相差三个数量级。显式方法被最陡的那段绑架——为了不爆炸必须用微步步长，慢段白白烧掉亿万步；换隐式 BDF（求解器里 method 设为隐式选项），自适应把慢段步长放大到毫秒级，计算量骤降几个数量级而结果吻合。刚性 = 显式方法的步长被最快尺度绑架，识别信号是"步长减半对照突然失效 + 曲线大段平坦偶发跳变"。

混沌系统还要叠加 3.4 节的账本：轨道逐点精度随时间必然耗尽，所以混沌报告里只写统计量（李雅普诺夫指数、维数、不变分布）并交代积分容差与总时长——5.3 节就把这些统计量的估计算法补齐。

⚠️ 常见坑：用 RK4 硬啃刚性方程，然后抱怨"步长一小机器就烫"。这不是方法不努力，是选型错了——快慢尺度比超过三个量级就该请隐式方法出场。

💡 关键直觉：自适应求解器的"步长历史"本身是诊断信息——步长骤降的位置标记了系统最陡的动作（放电、碰撞、边界层），比曲线本身更早告诉你动力学哪里有故事。

问题：数值结果该怎么"报告不确定性"？

一份合格的数值报告至少带四项元数据：求解器与容差（或步长）、积分总时长（混沌系统折算成李雅普诺夫时间数）、收敛证据（步长减半的误差缩水比）、统计量的样本数与置信区间。仿照实验物理的误差条文化，把"数值口径"当作仪器铭牌——本章各案例的输出都应附这类口径。审阅别人的结果时也按四项反查：缺任何一项，引用其数字都应打折。把这套纪律养成默认动作，数值实验才算坐稳"第三种科学方法"的位置——1.1 节年表里它挣来的地位，靠的正是这类可复核性。

本节要点回顾

• 误差阶是兑换率：欧拉 O(h)、RK4 O(h⁴)，实测减半规律一验便知。
• 选型三问：刚性吗（隐式）？保守吗（辛积分器）？要长期统计吗（控容差、盯守恒量）？
• 三道安检：步长减半、守恒量监视、双方法对照——缺一道都别下结论。
• 混沌的诚实口径：只报统计量，逐点精度在李雅普诺夫时间之外没有意义。

方法可信了，下一节给混沌做正式体检：庞加莱截面与李雅普诺夫指数。