4.1 最小可用的三体积分器：从方程到毕达哥拉斯实验
4.1 最小可用的三体积分器：从方程到毕达哥拉斯实验
本节摘要：三体数值实验的第一步是造仪器：把牛顿方程翻译成"状态向量加导数函数加步进循环"的三件套。本节给出百余行的完整 RK4 积分器，讲清步长的两难与刻度方法，然后用它跑一遍著名的毕达哥拉斯三体问题——质量 3、4、5 的三个天体从直角三角形的三个顶点静止出发，看它们如何在近距交会中交换能量，直到最轻者被抛出、剩下两个结成双星。
先泼一盆冷水：动手写代码之前
动手之前有一个必须先立起来的观念：数值积分不是"让计算机算方程"，而是"造一台会产生系统性误差的机器，然后管住这台机器"。三体轨道对误差极度敏感（第 5 章定量），任何一步的侥幸都会被指数放大成完全不同的结局。所以本节的顺序是刻意的：先立接口与格式，再谈步长，先刻度后使用——每一件都是为了让后面的实验结论属于引力，而不属于算法。
好消息是门槛比想象中低。一台"最小但合格"的三体积分器只要三件东西：把十八个状态量排成一列的状态向量、按方程返回全部导数的导数函数、按选定格式往前推一步的步进循环。第 2 章的受力零件已完成导数函数最难的部分，剩下的是装配工作。
接口设计：三件套怎么分工
状态向量是约定：三个天体、平面运动，排列成"位置六连、速度六连"的十二元组。导数函数职责单一——吃进状态，吐出导数（位置导数是速度，速度导数是加速度，加速度来自第 2.1 节的受力零件）。步进循环是格式本体，本节用四阶龙格-库塔：每步在一步之内四次采样斜率、加权平均，步长折半时误差缩到约十六分之一。
接口分清的好处立刻显现：换格式（下一节的蛙跳）只动步进循环，换问题（限制性三体、四体）只动状态与导数函数——仪器是模块化的，误差来源也是模块化的。
动手实验：毕达哥拉斯三体的完整复盘
实验背景：1913 年布劳提出的问题至今仍是三体数值实验的"出厂自检"——三个天体质量取 3、4、5（恰是一组勾股数），从直角三角形顶点静止释放。它浓缩了三体数值实验的全部考点：近距交会、能量交换、结局判读，而且结局稳定可复现，刻度与验证两相宜。
操作：初值取质量 3 在上方、4 在左下、5 在右下（比例合适的直角三角形），速度全零；RK4 小步长推进到六十个时间单位，全程监测能量漂移与两两距离，识别"哪个天体最终独自远去"。
import math G = 1.0 # 状态：[位置6；速度6]，质量表单独存放 masses = [3.0, 4.0, 5.0] state = [1.0, 3.0, -2.0, -1.0, 1.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] def deriv(s): q, v = s[:6], s[6:] a = [0.0] * 6 for i in range(3): for j in range(3): if i == j: continue dx = q[2 * j] - q[2 * i] dy = q[2 * j + 1] - q[2 * i + 1] r = math.hypot(dx, dy) + 1e-9 f = G * masses[j] / r ** 3 a[2 * i] += f * dx a[2 * i + 1] += f * dy return list(v) + a def energy(s): q, v = s[:6], s[6:] T = 0.5 * sum(masses[k] * (v[2 * k] ** 2 + v[2 * k + 1] ** 2) for k in range(3)) U = 0.0 for i in range(3): for j in range(i + 1, 3): r = math.hypot(q[2 * j] - q[2 * i], q[2 * j + 1] - q[2 * i + 1]) U -= G * masses[i] * masses[j] / r return T + U def rk4_step(s, dt): k1 = deriv(s) k2 = deriv([s[i] + 0.5 * dt * k1[i] for i in range(12)]) k3 = deriv([s[i] + 0.5 * dt * k2[i] for i in range(12)]) k4 = deriv([s[i] + dt * k3[i] for i in range(12)]) return [s[i] + dt / 6 * (k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i]) for i in range(12)] dt, t = 1.0e-3, 0.0 E0 = energy(state) while t < 60.0: state = rk4_step(state, dt) t += dt E1 = energy(state) print("相对能量漂移:", abs((E1 - E0) / E0)) q = state[:6] for k in range(3): print("天体", k + 1, "（质量", masses[k], "）末位置约", [round(q[2 * k], 1), round(q[2 * k + 1], 1)], " 距质心:", round(math.hypot(q[2 * k], q[2 * k + 1]), 1))
结果：能量相对漂移在 10^{-6} 量级（dt 取 0.001、RK4 的正常水准）；末态的几何一目了然——质量 3 的天体孤零零地跑到几十个长度单位之外（持续远去），质量 4 与 5 相距一两个单位并保有相对速度，结成缠绕的双星。解读这个结局有三层。第一层，物理层：全程总能量、总动量、总角动量守恒，但能量在两两之间反复倒手，两次三体密近交会后，最轻者带走了正的"单飞能量"，剩下 4 与 5 落入束缚态——这正是第 5 章结局分类学里"分离结局"的标准剧本，布劳问题是最早的实物标本。第二层，数值层：三体交会时天体距离短到千分之一量级，固定步长的 RK4 能挺过来，靠的是 RK4 的高阶刚度；但把 dt 放大十倍重跑，末态会明显变形甚至换一个天体被抛出——误差被交会放大到足以改变结局，这就是"伪混沌"的具体形状。第三层，方法论层：判断"谁被抛出"不能只看末位置，要跟踪两两距离的时间序列，确认距离单调增长才算分离——单一时刻的快照会骗人。变式：把三个质量都改成 5 重跑，对称性让结局迟迟难产、交会次数翻倍，积分预算也得翻倍——对称初值反而更耗算力，这是扫描实验预算表的常见反直觉项。
图 4-2 毕达哥拉斯三体的开局与终局
步长两难与仪器刻度
步长的两难值得单独摊开：放大步长省机时，但近距交会处力变化陡峭，一步跨过去细节全丢；缩小步长保真，机时按比例上涨，而扫描实验动辄要跑百万条轨道。工程上的解法是自适应步长——每步估计局部误差，超阈值就折半、富余就加倍。RK4 家族的经典搭档是步长折半对比法：同一步用整步与两个半步各算一次，差值就是局部误差的估计，无需真解。
刻度是开工前的必做工序，标准动作有三步。第一步，拿二体椭圆当标尺：解析周期是现成的真值，积分出的周期误差随步长的缩小比例就是收敛阶（RK4 应当看到"步长折半、误差缩十六倍"）。第二步，验守恒量：能量漂移应随步长按同阶规律下降且无单向趋势。第三步，极限测试：用毕达哥拉斯这类含密近交会的案例压一遍，记录"结局不被误差翻转"的最大步长——它才是这台仪器的真实工作上限，比理论收敛阶实用得多。
⚠️ 最容易被略过又最致命的一条：积分器的收敛不等于轨道的可信。混沌系统里，误差随时间指数放大，"每步都收敛"只保证短期可信；宣称长程结论前，必须先估一估这条轨道的李雅普诺夫时间（下一章的道具），看看预测保质期够不够用。
本节要点回顾
- 三件套接口：状态向量、导数函数、步进循环各司其职，换格式与换问题互不牵连；
- RK4 特性：四阶收敛，步长折半误差缩十六倍；受力零件直接复用第 2 章的体检合格品；
- 布劳问题：质量 3、4、5 静止释放，中盘密近交会倒手能量，终局最轻者单飞、其余结双星——分离结局的标本兼仪器自检；
- 步长三步刻度：二体标尺量收敛阶、能量漂移验守恒、交会案例定工作上限；
- 冷水重申：收敛只保短期可信，长程结论前先问预测保质期——这道工序在下一站升级成主角。
仪器能跑了，但 RK4 的能量会"漏水"。下一站把这台机器的隐秘故障摆上台面，换上不漏水的辛格式，并做一次完整的误差对账。