8.1 运动微分方程的数值解法
8.1 运动微分方程的数值解法
本节摘要：4.2 节学会了欧拉与龙格-库塔的基本步法，本节补上工程级仿真剩下的三件事：怎么用收敛实验给精度定价，怎么用事件检测处理弹跳、停机这类"状态突变"，怎么用解析特例给数值代码当裁判。三件事凑齐，数值解才算能签字交付。
从步法到工艺
第四章演示过误差阶数的含义：步长砍半，龙格-库塔四阶法的误差缩到约十六分之一。本节把这句话变成可操作的工艺——收敛实验：同一个问题连跑几档步长，看误差按什么比例缩水。比例对了，代码没病；比例不对，多半是方程写错或步长越过了方法的舒适区。这比"结果看起来合理"硬得多，因为它检验的是方法的数学本性，不是眼睛的错觉。
第二件装备处理真实系统里的"突变"。落球触地、电路开关、振幅越限报警，状态变量在瞬间转向，普通的定步长积分要么糊弄过去、要么直接算飞。事件检测的思路是：积分器只管平滑段，每当状态逼近预设的"事件面"（比如高度归零）就精确截停、切换规则、重启积分。弹跳球是它的标准教学病例。
手术一：收敛实验给精度定价
选一台有解析解的病例当裁判：平方阻力落体其实积得动——速度等于收尾速度乘双曲正切。用精确解对照四档步长的龙格-库塔，误差比读数就是方法的体检报告。
# 手术一：收敛实验。平方阻力落体有解析解 v = vt*tanh(g*t/vt) import math m, g, c = 80.0, 9.8, 0.25 vt = math.sqrt(m * g / c) # 收尾速度 56.0 f = lambda t, v: g - (c / m) * v * abs(v) exact = vt * math.tanh(g * 2.0 / vt) # t=2s 的精确速度 def rk4(dt, T=2.0): t, v = 0.0, 0.0 for _ in range(int(round(T / dt))): k1 = f(t, v) k2 = f(t + dt/2, v + dt/2 * k1) k3 = f(t + dt/2, v + dt/2 * k2) k4 = f(t + dt, v + dt * k3) v += dt * (k1 + 2*k2 + 2*k3 + k4) / 6 t += dt return v prev = None for dt in (0.2, 0.1, 0.05, 0.025): err = abs(rk4(dt) - exact) ratio = f"{prev/err:5.2f}" if prev else " -" print(f"dt={dt:6.3f} 误差 = {err:.3e} 与上档误差比 = {ratio}") # dt= 0.200 误差 = 3.026e-07 与上档误差比 = - # dt= 0.100 误差 = 1.899e-08 与上档误差比 = 15.94 # dt= 0.050 误差 = 1.189e-09 与上档误差比 = 15.97 # dt= 0.025 误差 = 7.438e-11 与上档误差比 = 15.99 print("t=2s 精确速度 =", round(exact, 6), "m/s") # 18.837030
误差比一路咬在十六附近，四阶的本性验明正身；想再快，看表读成本：步长从零点一压到零点零二五，机器时间翻了四番，只换来三位小数——精度要够就好，别为第七位小数付四倍机时，这是仿真预算的第一课。
手术二：事件检测与弹跳停机
事件检测的示范台：三米高落球，恢复系数零点六，回弹高度不足一厘米视为"躺平"自动停机。每段飞行是平滑的自由落体，地面接触是事件，事件一触发就按恢复系数反转速度、重开一段。总时间与"无限次弹跳"的解析级数对照，停机阈值的物理含义立刻现形。
# 手术二：事件检测。地面穿越为事件，回弹高度不足 eps 停机 import math g, e, h0, eps = 9.8, 0.6, 3.0, 0.01 t_total = math.sqrt(2 * h0 / g) # 首段下落 h = h0 * e * e # 第一次回弹高度 bounces = 1 while h > eps: t_total += 2 * math.sqrt(2 * h / g) # 一段完整的升与落 bounces += 1 h *= e * e print("弹跳次数 =", bounces, " 总飞行时间 =", round(t_total, 4), "s") # 弹跳次数 = 6 总飞行时间 = 2.9473 s T_inf = math.sqrt(2 * h0 / g) * (1 + e) / (1 - e) print("无限次弹跳的理论总时间 =", round(T_inf, 4), "s") # 3.1298 print("停机阈值截掉的尾巴 =", round(T_inf - t_total, 4), "s") # 0.1825 # 第 n 次回弹高度预告：h0*e 的 2n 次幂 for n in (1, 3, 6): print(f"第 {n} 次回弹高度 =", round(h0 * e**(2*n), 5), "m") # 第 1 次回弹高度 = 1.08 m # 第 3 次回弹高度 = 0.13997 m # 第 6 次回弹高度 = 0.00653 m —— 已低于 1 厘米停机线
两组读数互为注脚：一厘米的停机线让仿真在六次弹跳后收工，比无限级数的理论总时间少零点一八秒——少的正是被阈值截掉的"微颤尾巴"。工程上要自己回答的问题由此浮出：停机阈值定多细，答案就贵到哪一级，阈值的选择从来不是纯技术决定。
💡 实验习惯：任何数值代码上线前，先给它找一道有解析解的特例当裁判（本节的阻力落体、上一章的单摆小角），收敛实验过关再上真实工况——没有裁判的仿真，错了都不知道该向谁申诉。
步长自适应：让代码自己定价
收敛实验靠人工档位，自适应步长把它交给代码：每步用两个精度的差值估计局部误差，超差就折半步长重走，富余就放大步长省机时。经济效果立竿见影——平滑段大步飞奔、转折点自动碎步，同样的精度指标，机时常能省下过半。下面用变步长演示这条逻辑的最简版本。
# 自适应步长最简演示：对阻尼落体按局部误差调步长 import math g, c, m = 9.8, 2.0, 1.0 f = lambda t, v: g - (c / m) * v def euler_step(t, v, dt): return v + f(t, v) * dt t, v, dt, tol, T = 0.0, 0.0, 0.4, 0.002, 3.0 big, small = 0, 0 while t < T: v_half = euler_step(t, v, dt / 2) v_two = euler_step(t, v_half, dt / 2) # 两半步 v_one = euler_step(t, v, dt) # 一整步 err = abs(v_two - v_one) if err > tol and dt > 0.01: dt /= 2; small += 1 # 超差折半 else: v = v_two; t += dt; big += 1 if err < tol / 10 and dt < 0.4: dt *= 2 # 富余翻倍 print("整步推进次数 =", big, " 步长折半次数 =", small) # 整步推进次数 = 106 步长折半次数 = 5 print("t=3s 时 v =", round(v, 4), "m/s 解析值 =", round(g / (c/m) * (1 - math.exp(-(c/m)*3)), 4)) # t=3s 时 v = 4.8934 m/s 解析值 = 4.8879 m/s（阻尼衰减时间常数 0.5 s）
折半集中在起始的速度猛增段，之后步长一路放大——自适应的价值不是更准，而是"把机时花在刀刃上"；末速与解析值差在容差量级，正是 tol 定价的精度。真实求解器还带误差控制器与事件检测的联动，本节给的是它的骨架。
要点复盘：
- 收敛实验是体检报告：误差比咬住方法的阶数，代码有病当场现形；
- 精度按需购买：步长压四倍只换三位小数，仿真预算花在刀刃上；
- 事件检测管突变：平滑段积分、事件面截停、规则切换重启，弹跳与停机不再糊账；
- 阈值是物理决定：停机线截掉的尾巴要用理论级数估价，精度与成本一起谈判。
方程解得动了，交付的规矩还没立。下一节把受力分析图与建模流程标准化，让计算书经得起陌生人的复核。