4.2 Runge-Kutta 家族

4.2 Runge-Kutta 家族

本节摘要：显式欧拉法一阶精度且在振荡系统上能量单调膨胀；RK4 用四次函数评估换取四阶精度，成为通用默认；嵌入对方法（RK45）让步长随解的形态自适应。本节在谐振子试金石上完整复现能量漂移现象，并剖析自适应步长的控制回路。

一、案发现场：欧拉法给弹簧注入能量

谐振子 y'' = -y 化为一阶系统（位置 y、速度 v）后，显式欧拉一步只是"沿当前切线走"：

import numpy as np

def euler_step(f, t, y, h):
    return y + h * f(t, y)

osc = lambda t, y: np.array([y[1], -y[0]])   # y[0]=位置, y[1]=速度

h, T = 0.05, 50.0
n = int(T / h)
y = np.array([1.0, 0.0])
for k in range(n):
    y = euler_step(osc, k*h, y, h)

energy = 0.5 * (y[0]**2 + y[1]**2)           # 初能量 = 0.5
print(f"t={T} 时能量: {energy:.4f}")          # 典型 1.6+，能量翻三倍！
​

取证：每一步欧拉都把圆轨道上的点移到切线上——切线永远落在圆外。单步误差 O(h²) 看似无害，但误差方向系统性向外，数十步内振幅肉眼可见地膨胀。这不是 bug，是一阶方法对保守系统的结构性缺陷：它把哈密顿系统当成了能量源。

二、RK4：四次采样换四阶精度

RK4 在每步内采样四个斜率（起点、两个中点试探、终点），加权平均出等效的"高阶切线"：

import numpy as np

def rk4_step(f, t, y, h):
    k1 = f(t, y)
    k2 = f(t + h/2, y + h/2 * k1)
    k3 = f(t + h/2, y + h/2 * k2)
    k4 = f(t + h, y + h * k3)
    return y + h/6 * (k1 + 2*k2 + 2*k3 + k4)

osc = lambda t, y: np.array([y[1], -y[0]])

h, T = 0.05, 50.0
y = np.array([1.0, 0.0])
for k in range(int(T/h)):
    y = rk4_step(osc, k*h, y, h)
print(f"RK4 t=50 时能量: {0.5*(y[0]**2+y[1]**2):.6f}")   # 约 0.5000x，漂移极小
​

同样的步长，RK4 的能量漂移比欧拉小五个数量级以上。代价是每步四次函数评估（欧拉一次）——换算成"每次评估买的精度"，RK4 在步长不太小时几乎是纯赚。这就是它从手算时代服役到今天的原因：在精度阶敏感、问题不刚性、函数评估不太贵的交集里，它仍是默认答案。

慢镜头看一个隐藏问题：RK4 的能量漂移小但不为零，且随时间线性累积。把 T 拉到几千、几万，相位误差逐渐显形——数值解的周期比真值略长，长时间后峰值出现在错误时刻。天体轨道、分子动力学的长期能量追踪对这类漂移零容忍，对症的药是辛积分器（保持哈密顿结构的保结构方法），它不追求单步高精度，而保证能量误差有界振荡而非漂移。

三、自适应步长：让积分器自己踩油门刹车

RK45（Dormand-Prince）同时算四阶和五阶两个版本，两者之差就是现成的局部误差估计——不需要真解就能知道"这步走得多准"：

from scipy.integrate import solve_ivp
import numpy as np

osc = lambda t, y: [y[1], -y[0]]

sol = solve_ivp(osc, [0, 50], [1.0, 0.0], method='RK45',
                rtol=1e-10, atol=1e-12, dense_output=True)
print(f"RK45 步数: {sol.t.size-1}")
print(f"t=50 能量误差: {abs(0.5*(sol.y[0,-1]**2 + sol.y[1,-1]**2) - 0.5):.2e}")

# 对比刚性逼问下显式方法的步长（下一节主题的预告）：
sol2 = solve_ivp(osc, [0, 50], [1.0, 0.0], method='RK45', rtol=1e-6)
print(f"宽松容差下步数: {sol2.t.size-1}")
​

控制回路拆解：局部误差估计与容差比较——超了就把步长乘一个小于 1 的因子重算（刹车），远小于容差就放大步长（油门）。增益因子按误差阶数设定，保证调节平稳不振荡。两个容差的分工要读准：rtol 管相对精度（跟随解的量级），atol 兜住解过零时的绝对精度，只设 rtol 会让过零段精度失控。

图 4.2-1 自适应步长控制回路

四、工程清单：solve_ivp 使用要点

from scipy.integrate import solve_ivp
import numpy as np

def rhs(t, y):
    return [y[1], -0.1*y[1] - np.sin(y[0])]   # 阻尼摆

sol = solve_ivp(rhs, [0, 100], [2.0, 0.0], method='RK45',
                rtol=1e-8, atol=1e-10, dense_output=True,
                max_step=0.5)
print(sol.success, sol.message)
​

清单四条：①rhs 的签名是 (t, y)，与手写循环的习惯 (y, t) 相反，顺序写反是高频事故；②一定检查 sol.success，积分失败时结果数组不可用；③dense_output=True 提供连续插值，避免事后自己插值引入误差；④max_step 防止积分器跨过物理上重要的事件（碰撞、开关切换）。

💡 关键直觉：rtol 从 1e-6 收紧到 1e-10，步数大约变为原来的 1.6 倍（四阶方法的四分之一次方关系）——精度便宜，但先把量纲和 atol 设对，比一味收紧 rtol 有效得多。

本节要点回顾

• 欧拉的结构性原罪：切线永远落在圆外，保守系统能量单调膨胀；教学价值高，生产用途低
• RK4 的汇率：四倍评估换四阶精度，"每次评估的精度"在非刚性场景下几乎纯赚
• 长期能量体检：跟踪守恒量是 ODE 结果验证的第一手段，漂移线性累积是常态，需要保结构方法时别硬扛
• 嵌入对自适应：两个精度阶之差 = 免费的误差估计器；rtol/atol 分工必须理解
• solve_ivp 纪律：签名顺序、success 检查、dense_output、max_step 四件事做完才算交付

下一节处理显式方法的终极噩梦：刚性系统——步长被稳定性而不是精度绑架的疑难案件。