4.2 辛积分器与误差对账：别让算法替引力做主

4.2 辛积分器与误差对账：别让算法替引力做主

本节摘要：RK4 积分器有个隐疾：长程运行时系统能量单向漂移，轨道被算法慢慢"加热"，最终呈现的可能是伪混沌。辛积分器（蛙跳、维雷特）凭几何性质把能量误差锁成有界振荡，是长程三体实验的标配。本节做一次完整的对账实验，量出两种格式的差异，并给出格式选择的工程准则。

一条会自己变胖的轨道

把上一节的二体验证实验换个跑法：不跑三十圈，跑三千圈，然后盯着椭圆的半长轴看。用 RK4、固定步长，你会看到一件怪事——椭圆在缓慢地"变胖"，能量沿着一个方向单调爬升，像被一只看不见的手持续加热。这条轨道没有任何物理理由扩散，牛顿二体解封在圆锥曲线里一动不动。变胖的是算法：RK4 每一步都产生极小的误差，这些误差在时间上不平均、正负不抵消，累积成单向漂移。

对短程实验（上一节的六十个时间单位）这点漂移无伤大雅；对长程实验（星团演化要模拟百万年时间单位、稳定性研究要跑十万圈周期轨道）它是致命的——漂移积少成多，会把规则轨道推过稳定边界，让你把"算法造成的解体"误读成"引力的判决"。三体数值史上不乏这种乌龙：早期一些"三体系统必然解体"的结论，后来被证明部分是积分器加热的产物。

辛格式的原理：不更准，但不使坏

十九世纪中叶打下的哈密顿框架（第 2 章）在这里等了半个世纪的工程回报。哈密顿系统的相空间流动有一个几何性质——保辛：它保持相空间体积与一种叫辛结构的角度关系。RK4 这类通用格式破坏这个性质，误差在能量方向上随机游走再叠加系统性漂移；辛积分器每一步都精确保持这个几何，因此它实际上精确求解的是"另一个与真系统只差微小项的哈密顿系统"（所谓影子哈密顿）——能量误差不再漂移，只在真值附近做有界振荡。

工程上最常用的辛格式是蛙跳/维雷特格式，它把一步拆成三拍：先按引力推半步速度（kick），再按新速度推整步位置（drift），最后再用新位置的引力补半步速度（kick）。代码比 RK4 还短，每步只需一次受力评估（RK4 要四次），又快又保结构，代价是阶数只有二阶、且固定步长下处理高偏心率轨道时近距点误差偏大。

图 4-3 两种格式的能量误差对账单

动手实验：给两种格式做一次对账

实验背景：量化的结论才叫对账。同一二体椭圆、同一步长，分别用 RK4 与蛙跳跑长程，记录能量误差随圈数的曲线，检验"单向漂移"与"有界振荡"的定性差异。

操作：取偏心率 0.3 的椭圆，两种格式各跑两千圈，每十圈记录一次相对能量误差；步长取周期的百分之一（两种格式同价对比）。代码如下：

import math

mu = 4 * math.pi ** 2
a, e = 1.0, 0.3
T = 2 * math.pi * math.sqrt(a ** 3 / mu)
dt = 0.01 * T

def acc(p):
    r = math.hypot(p[0], p[1])
    f = -mu / r ** 3
    return [f * p[0], f * p[1]]

def energy(p, v):
    return 0.5 * (v[0] ** 2 + v[1] ** 2) - mu / math.hypot(p[0], p[1])

# 初值：远日点出发，速度纯切向
p0 = [a * (1 + e), 0.0]
v0 = [0.0, math.sqrt(mu * (1 - e) / (a * (1 + e)))]
E0 = energy(p0, v0)

def run_rk4(n_cycles):
    p, v = list(p0), list(v0)
    worst = 0.0
    per = int(round(T / dt))
    for _ in range(per * n_cycles):
        a1 = acc(p)
        v2 = [v[i] + 0.5 * dt * a1[i] for i in range(2)]
        p2 = [p[i] + 0.5 * dt * v2[i] for i in range(2)]
        a2 = acc(p2)
        v3 = [v[i] + 0.5 * dt * a2[i] for i in range(2)]
        p3 = [p[i] + 0.5 * dt * v3[i] for i in range(2)]
        a3 = acc(p3)
        v4 = [v[i] + dt * a3[i] for i in range(2)]
        p4 = [p[i] + dt * v4[i] for i in range(2)]
        a4 = acc(p4)
        v = [v[i] + dt / 6 * (a1[i] + 2 * a2[i] + 2 * a3[i] + a4[i]) for i in range(2)]
        p = [p[i] + dt / 6 * (v2[i] + 2 * v3[i] + 2 * v4[i] + v[i]) for i in range(2)]
    return abs((energy(p, v) - E0) / E0)

def run_leapfrog(n_cycles):
    p, v = list(p0), list(v0)
    a_now = acc(p)
    per = int(round(T / dt))
    for _ in range(per * n_cycles):
        v = [v[i] + 0.5 * dt * a_now[i] for i in range(2)]
        p = [p[i] + dt * v[i] for i in range(2)]
        a_now = acc(p)
        v = [v[i] + 0.5 * dt * a_now[i] for i in range(2)]
    return abs((energy(p, v) - E0) / E0)

for n in (10, 100, 1000, 2000):
    print("圈数", n,
          " RK4 相对误差:", f"{run_rk4(n):.2e}",
          " 蛙跳相对误差:", f"{run_leapfrog(n):.2e}")
​

结果：RK4 的误差随圈数单调上升，从十圈的 10^{-8} 量级爬到两千圈的 10^{-6} 到 10^{-5} 量级，方向恒定；蛙跳的误差起点稍高（二阶格式，同为 10^{-7} 量级），但从十圈到两千圈始终在同一水平附近往复，没有方向、没有增宽。解读：这张对账单就是辛性质的实物形态——RK4 每圈欠一点账，两千圈下来账面可见；蛙跳每步也在欠账，但它欠的账永远绕着零点转圈，永不结转。注意"起点稍高"这一细节：辛格式不是全面碾压，短程高精度需求（比如精确定位某次交会的时刻）下 RK4 依然是好选择，它的优势在长程守恒而非单步精度。变式：把偏心率加到 0.9 重跑——蛙跳的误差包络明显增宽（近距点力太陡，固定步长吃力），此时要么缩小步长，要么换自适应格式：辛性质在高偏心率下不是免死金牌，这是固定步长辛格式最著名的软肋。

格式选择的工程准则

把选择逻辑压成一张表：

最后一个格子多说一句：如今的实践者很少从零写积分器，天体力学社区有几款成熟的开源 N 体库（如天文学界常用的 rebound 一类工具），内置多种辛格式、自适应方案与碰撞正则化，本册教手写积分器不是为了替代它们，而是为了让你在使用它们时能看懂配置项背后的取舍——知道什么时候该怀疑一份数据，比会用软件重要得多。

本节要点回顾

• RK4 的隐疾：非辛格式长程运行时能量单向漂移，规则轨道会被算法"加热"越过稳定边界，产出伪混沌；
• 辛格式的底牌：精确保持哈密顿几何，实际求解影子哈密顿，能量误差有界振荡、永不结转；
• 对账实验：同价对比两千圈，RK4 漂移、蛙跳振荡——误差有没有"方向感"是最快的第一判据；
• 软肋明确：固定步长辛格式在高偏心率下包络增宽，近距交会需自适应或正则化补位；
• 选择准则：短程精度用通用高阶，长程守恒用辛，交会用自适应，规模化用成熟 N 体库——按场景点菜，别背通用答案。

仪器验讫，可以出片了。下一站把连续轨迹压成一张可判读的照片：庞加莱截面，混沌的第一张证件照。