数值解法

欧拉法的局部截断误差是 O(h²)，但全局误差是 O(h)。为什么少了一阶？因为走到终点要走 T/h 步，每步 O(h²) 的误差累积 T/h 次，就变成了 O(h)。这个"局部减一阶等于全局"的规律对所有单步法都成立。右下那张 log–log 图直接把阶数量出来了：斜率就是 p，欧拉法应该接近 1。读数里还给了"步长减半误差变成几分之一"——一阶方法是 1/2，二阶是 1/4，四阶是 1/16。这是检验自己写的求解器有没有写对最快的办法：把步长减半，看误差是不是按理论倍数下降。如果不是，要么阶数写错了，要么已经撞到了机器精度的地板（把 h 拖到最左边你会看到误差不再下降反而上升，那就是舍入误差开始主导）。

RK4 每步要算四次右端函数，是欧拉法的四倍工作量——但它换回来的是四阶精度。看最后那个读数：要把误差降到 10⁻⁶，欧拉法需要的步数通常是 RK4 的几万倍。就算 RK4 每步贵四倍，总成本仍然差着几个数量级。"每步更贵但阶数更高"在几乎所有非刚性问题上都是划算的，这就是 RK4 成为默认选择的原因。不过请注意"阶数高"只在 h 足够小的时候才兑现——把 h 拖到最右边，四条线会挤在一起，那时候高阶方法并不比低阶的好多少。另外注意 Heun 和中点法都是二阶，但常数因子不同；阶数决定斜率，常数因子决定截距，两者都影响实际误差。

刚性问题的麻烦不在精度，而在稳定性。y′ = −λy 的精确解衰减得飞快，本身一点也不复杂；但显式欧拉只有在 h·λ < 2 时才不发散，λ 一大，h 就被逼得极小——不是为了算得准，而是为了不爆炸。把 λ 调到 200、h 调到 0.02，显式欧拉立刻振荡发散；换成隐式欧拉，同样的 h 稳稳当当。隐式欧拉的绝对稳定域包含整个左半平面（A-稳定），代价是每步要解一个方程——线性问题好办，非线性就要牛顿迭代。顺带一个反直觉的事实：RK4 的稳定界（≈2.78）只比欧拉（2）大一点点。也就是说，提高阶数几乎不能缓解刚性，真正的办法只有换成隐式方法。现实中的刚性问题（化学反应网络、电路、多尺度系统）几乎都要用 BDF、Radau 这类隐式求解器。

这是数值方法里最反直觉的一件事：一阶的辛欧拉在长时间积分上打败四阶的 RK4。默认参数（h = 0.1、120 个周期）下：显式欧拉的能量指数级上涨（漂移已是 10¹⁹ 量级，轨道直接螺旋飞出去）；RK4 只掉了 0.0094%，肉眼看不出来；辛欧拉在 0.476 到 0.526 之间来回振荡（±5%），但长期均值只动了 0.004%。关键差别不在幅度，而在有没有累积。把步长调到 0.3、周期数拉到 400 再比一次：RK4 的能量一路掉了 7.24%（误差随步数线性累积），而辛欧拉的振荡带虽然变宽到 ±15%，长期漂移仍然只有 0.14%——差了五十倍。原因是辛方法精确地保持了相空间的面积（辛结构），于是它虽然不精确求解原方程，却精确求解了一个非常接近的哈密顿系统。所以它的能量误差有界，而不是累积。这一点在天体力学里是决定性的：要算太阳系十亿年的演化，你需要的不是每一步都很准，而是误差不要累积。看"一圈之后的位置误差"这个读数：RK4 在单圈上确实远比辛欧拉准——单步精度和长期行为是两个不同的指标，不要用一个替代另一个。