三体仿真中的三算平衡：算据、算法与算力如何影响数值积分稳定性

前阵子帮一个做物理仿真的小团队看代码，他们折腾了将近三个月，把一套三体数值积分从单机跑到了多卡集群，结果轨迹还是对不上参考结果。我打开仓库做的第一件事不是看并行部分，而是打印初始总动量——三个天体的动量加起来不是零，整个系统在以一个恒定的速度朝某个方向平移。相对构型全对，绝对坐标全错。他们追了三个月的"算力不够"，根子上是一个谁都没检查的初值问题。

这个故事我后面还会提到。它几乎把"三体与三算"这件事讲透了：三体问题没有解析解，只能靠数值近似往前推，而数值近似的质量由三个互相牵制的东西决定——算据（你喂进去的初始条件和数值精度）、算法（你怎么推进时间）、算力（你有多大的机器、怎么把它用起来）。我把这三样合称"三算"。任何一个环节短了，另外两个环节的投入都会按比例打水漂，和木桶装水一个道理。

这篇写给两类人。一类是被三体问题吸引、想自己动手写一套仿真但不知道坑在哪的；另一类不一定关心天体力学，但手上经常处理"需要长期演化的复杂系统"——分子动力学、流体、控制回路、种群模型都算。三体只是一个把问题放大到极致的沙盘，三算这笔账在任何长时程仿真里都躲不掉。

1. 三体难的不是"三个"，是"长跑"

1.1 解析解消失之后，剩下的全是数值手艺

两体问题有漂亮的解析解，开普勒轨道，椭圆、双曲线写得清清楚楚，代入时间就能算出位置。加进第三个天体，这套东西立刻崩掉。牛顿本人在这上面耗了大量精力，后来庞加莱在十九世纪末给出了一个让人不太舒服的结论：三体问题一般不存在用初等函数表达的守恒量，轨道对初始条件敏感依赖。这不是"暂时还没算出来"，而是问题本身的性质。

这个结论对工程实践的意义非常直接：你不可能通过"更聪明地推导公式"绕开数值积分。想要轨迹，就得一步一步往前推，每一步都带着浮点误差，误差还会被系统本身的动力学放大。所以三体仿真的质量，本质上不取决于你多懂天体力学，而取决于你把数值这套手艺做到什么程度。

我用一个类比来解释敏感依赖。把三个天体想象成三颗在山脊上滚的球，山脊两边都是陡坡。初始位置差一根头发丝，滚下去几十米后，一个落到左边的谷底，一个落到右边的谷底，轨迹完全分岔。你没法说哪个是"错"的，因为输入本来就不同。数值仿真面对的就是这种局面——双精度的初始误差大约在 1e-16 这个量级，它迟早会被放大到肉眼可见。

1.2 三算的分工：谁负责什么，谁拖累谁

算据、算法、算力这三者不是并列关系，是有明确的职责边界和失效模式的。

• 算据管"输入有多真"。它包括初始位置速度的有效位数、单位制是否归一化、坐标系是否选对（质心系还是实验室系）、存储用的浮点格式。算据是天花板——输入精度是 1e-6，后面用什么算法都救不回来。
• 算法管"怎么往前推"。是欧拉、龙格库塔还是辛积分，步长固定还是自适应，要不要软化。算法决定了误差是"有界振荡"还是"单调漂移"，这两种失效模式的后果差别巨大。
• 算力管"单位时间能推多少步、扫多少组初值"。它对单条轨迹的精度几乎没有帮助，但对统计性质的研究、对参数扫描、对寻找周期解是决定性的。

它们之间的拖累关系很典型。算力不足的时候，人会本能地放大步长来换取速度，这一步直接把算法精度毁掉；算法选错导致能量漂移，人会以为是时间不够长、机器不够快，于是加机器，加了还是漂；算据不准导致结果发散，人会去怀疑并行归约顺序、怀疑浮点不确定性，排查方向整个歪掉。我在实际项目里见过太多次这种循环：问题出在 A，人在 B 上花了两周。

一个实用的判断顺序：先查守恒量（能量、动量、角动量），再查精度格式，最后才动并行。守恒量不对的话，后面两项查了也白查。

2. 算据：初值里那一位小数，滚到最后就是另一条轨迹

2.1 用 Lyapunov 时间算一算你的预测窗口有多长

误差怎么增长，可以用一个粗略但足够用的公式估：e(t) ≈ e0 · exp(λt)，λ 是系统的最大 Lyapunov 指数。等质量三体系统在归一化单位下，λ 的量级大概在 0.1 到 1 之间，具体取决于构型。取 λ = 1 做保守估计。

那么误差从 e0 放大到 e_target 需要的时间是：

t_max = (1/λ) · ln(e_target / e0)

代入两组数字看差别：

这张表最值得看的地方不是具体数值，而是它的对数关系。精度从 1e-6 提到 1e-16，整整十个数量级，预测窗口只从 11.5 延长到 34.5，三倍。这就是为什么"提高初值精度"这条路收益递减得厉害。大气模型也是同一笔账，观测精度再涨几个数量级，预报窗口也就延长一倍出头，这不是算力问题。

所以我给团队的建议通常是：先明确你要的预测窗口有多长。如果只需要跑 10 个时间单位看构型变化，初值精度 1e-10 绰绰有余，把精力花在算法和可视化上。如果要做长时间的统计性质（比如碰撞率、逃逸率），那就别指望单条轨迹有意义，直接做系综平均，几万条轨迹跑起来看分布，这时候算力才真正开始值钱。

2.2 归一化没做对，双精度也只跑出单精度的效果

这是我见过最隐蔽的算据问题。很多人写第一版三体代码时，直接拿真实天文单位往里塞：距离用米（1e11 量级），质量用千克（1e24 到 1e30），引力常数 6.67e-11。三个数量级跨度巨大的量乘在一起，浮点有效位立刻被吃掉一大截。

浮点数不是均匀分布的。双精度有 52 位尾数，能保证大约 15 到 16 位十进制有效数字，但这个"有效"是相对量级而言的。你在 1e24 量级的数上加一个 1e10 量级的修正，那个修正直接掉进尾数下面被截断。位置更新的增量很小而位置本身很大时，同样的事情会发生。

标准做法是归一化到自然单位：

• 令引力常数G = 1
• 让天体质量都在 O(1) 附近（比如用总质量归一）
• 让典型轨道半径在 O(1) 附近（比如用初始平均间距归一）
• 时间单位随之确定为sqrt(L³/(G·M))量级

归一化之后，所有参与运算的量都在 1 附近晃，浮点的相对精度用满，也不会出现大数吃小数。这一步不做，后面所有关于精度的讨论都是空谈——你用双精度存的东西，实际精度可能只有单精度。

顺手贴一段归一化后的基础实现，后面几节都会用到：

import numpy as np def accelerations(pos, mass, G=1.0, eps=0.0): """pos: (N, 3), mass: (N,), 返回 (N, 3) 加速度""" n = pos.shape[0] acc = np.zeros_like(pos) for i in range(n): for j in range(i + 1, n): d = pos[j] - pos[i] r2 = np.dot(d, d) + eps * eps inv_r3 = 1.0 / (r2 * np.sqrt(r2)) acc[i] += G * mass[j] * d * inv_r3 acc[j] -= G * mass[i] * d * inv_r3 return acc def total_energy(pos, vel, mass, G=1.0, eps=0.0): ke = 0.5 * np.sum(mass[:, None] * vel * vel) pe = 0.0 n = len(mass) for i in range(n): for j in range(i + 1, n): r = np.sqrt(np.dot(pos[j] - pos[i], pos[j] - pos[i]) + eps * eps) pe -= G * mass[i] * mass[j] / r return ke + pe

注意eps在能量函数里也要用同样的形式，不然你算出来的能量漂移里混进了"软化前后势能不一致"的假信号，会误判算法有问题。这是个非常常见的自坑。

3. 算法：长时程要的是结构正确，不是单步精确

3.1 辛积分器为什么能在长跑里赢过高阶龙格库塔

刚接触数值积分的人几乎都有同一个直觉：阶数越高越准。RK4 是四阶，欧拉是一阶，那当然用 RK4，甚至有人上 RK45 自适应。单步看这个判断没错，跑一万步就完全反过来了。

原因在于哈密顿系统的几何结构。三体问题在没有耗散的情况下是一个哈密顿系统，它的相空间演化保持辛结构，等价地说，能量在真实轨道上是严格守恒的。RK4 不是辛方法，它在每一步引入的误差对能量有一个系统性的偏向，不是随机的正负抵消，而是每次都朝同一个方向偏一点。跑得越久，能量单调漂移越明显，通常是往下掉（数值耗散），表现为轨道慢慢向内收缩，天体一个个掉进中心。

辛积分器（最常用的是 velocity Verlet，也叫 leapfrog）做的事情不一样。它是二阶的，单步精度不如 RK4，但它精确保持辛结构，能量误差不累积漂移，而是在真值附近做有界振荡，振幅正比于O(dt²)。这个性质在长时程仿真里是决定性的。

def velocity_verlet(pos, vel, mass, dt, steps, G=1.0, eps=0.0): pos = pos.copy() vel = vel.copy() acc = accelerations(pos, mass, G, eps) energies = np.empty(steps) for k in range(steps): vel += 0.5 * dt * acc # 半步速度 pos += dt * vel # 整步位置 acc = accelerations(pos, mass, G, eps) vel += 0.5 * dt * acc # 再半步速度 energies[k] = total_energy(pos, vel, mass, G, eps) return pos, vel, energies

这段代码有两个地方值得盯。第一，加速度每步只算一次（重用上一步末尾算出的acc），这是 velocity Verlet 效率高的关键，别写成算两次。第二，半步速度的顺序不能颠倒，这是算法定义的一部分，写反了时间可逆性就没了，后面 6.3 节的自检方法会直接把它抓出来。

实测对比很直观。同样跑 10000 步，RK4 的能量相对误差大概会漂到 1e-4 甚至更差，而且是单调的；velocity Verlet 的能量相对误差在一个固定带宽里来回晃，量级取决于步长。步长减半，振荡幅度降到四分之一。这个"有界"两个字，就是辛积分器在长跑里的全部价值。

3.2 自适应步长和辛结构，只能选一个

三体系统有个讨厌的特性：大部分时间三个天体离得挺远，力变化平缓，大步长完全够用；偶尔发生近距离遭遇，两个天体擦身而过，加速度按1/r²暴涨，固定步长直接跨过去，一步就把两个天体甩到不知道哪里。

绝大多数人的第一反应是上自适应步长——误差大的时候自动缩小步长，误差小的时候放大。RK45 这类方法就是这么干的，效果对短时仿真很好。但这里有个硬约束：变步长会破坏辛结构。辛积分器的长期能量守恒性质，建立在固定步长（准确说是固定时间映射）的前提上。步长一变，每一步的映射都不一样，辛性就丢了，能量又开始漂移。你想同时要"自适应应对近距遭遇"和"长期能量守恒"，数学上做不到。

工程上的取舍通常是三选一：

大多数人的实际需求落在第一行。选一个能覆盖最短特征时间的固定步长，接受"平时浪费一点"的事实。经验做法是：先跑一小段，记录最小天体间距，估算该间距下的轨道周期，取周期的 1/500 到 1/1000 作为步长。别照抄别人的 dt 数值，归一化方式不同，数值没有可比性。

3.3 软化长度是个躲不掉的误差源

为了不让近距离遭遇把积分炸掉，常见做法是引入软化长度eps，把1/r²换成1/(r² + eps²)。这个操作在星系模拟、暗物质模拟里几乎是标配，因为那里天体被当作质量团块，本来就没有点粒子意义。但对三体这种"点质量"问题，软化是实打实地改了物理模型。

改了多少，可以估。当r >> eps时，1/(r²+eps²) ≈ 1/r² · (1 - eps²/r²)，相对偏差是(eps/r)²量级，影响很小。当r ~ eps时，力的误差是百分之几十，轨迹就完全不是原来那条了。所以软化长度的选取要和你的研究目标绑定：

• 只关心系统整体的定性行为（会不会散开、会不会形成束缚对），软化取到初始平均间距的 1% 到 5%，够用。
• 关心精确的碰撞时刻、逃逸速度，软化必须小于最小间距，宁可用小步长硬扛。
• 关心系综统计，软化会系统性改变碰撞截面，这时候软化长度必须当作模型参数写进论文里，不能当成数值技巧藏着。

还有一个实操细节前面提过：加软化之后，能量函数必须用同一个软化形式计算。很多人只改了力函数，没改势能，跑出来能量一直在漂，然后开始怀疑辛积分器不灵。这种自坑我见过不止一次，排查起来还挺费时间，因为代码看起来"没什么问题"。

4. 算力：三体本身不需要并行，需要并行的是别的东西

4.1 小系统上 GPU 为什么反而慢一截

N = 3 的三体系统，每步的力计算是 3 对相互作用，总共 9 次（考虑对称性其实算 3 次）。这个计算量在 CPU 上不到一微秒。而一次 GPU kernel launch 的固定开销，实测在 5 到 10 微秒之间。也就是说你为了算 3 次浮点运算，付了 10 微秒的调度成本，反而比 CPU 慢一个数量级。

这个反直觉的结论我特别想强调，因为"上 GPU 加速"几乎是所有人遇到性能问题时的第一反应。GPU 的算力优势建立在"海量同构任务并行"上，几万个线程同时干活才能填满那些流处理器。单个三体系统是完全串行的依赖链——算完力才能更新速度，更新完速度才能更新位置，下一步的力又依赖新位置。这条链上的并行度只有 3 个天体，GPU 大部分单元在闲着。

经验阈值：单个系统的粒子数低于几百，上 GPU 基本都是亏的，除非你能把多个系统拼成一批。这也解释了为什么分子动力学软件里 GPU 加速只在几万原子以上才有明显收益，而三体这种 N=3 的问题，多线程都没必要开。

4.2 batched 思路：几千个三体打包成张量一起推

真正需要算力的场景不是"跑一条三体轨迹"，而是"跑一万条不同初值的三体轨迹，看统计分布"或者"扫参数空间找周期解"。这时候并行维度不该在系统内部，而应该在系统之间。把 B 个独立的三体系统打包成形状(B, 3, 3)的张量，一次推进所有系统：

def accel_batched(pos, mass, G=1.0, eps=0.0): """pos: (B, N, 3), mass: (N,), 返回 (B, N, 3)""" diff = pos[:, None, :, :] - pos[:, :, None, :] # (B, N, N, 3) r2 = np.sum(diff * diff, axis=-1) + eps * eps # (B, N, N) idx = np.arange(pos.shape[1]) r2[:, idx, idx] = np.inf # 排除自相互作用 inv_r3 = 1.0 / (r2 * np.sqrt(r2)) inv_r3[:, idx, idx] = 0.0 return G * np.einsum('bij,bijd->bid', mass[None, :] * inv_r3, diff)

B 取一万、N 取 3 的时候，diff这个中间张量是 10000 × 3 × 3 × 3 = 27 万个元素，几个 MB，塞得进缓存。整个计算完全向量化，没有 Python 层的循环，在 GPU 上能跑满。推进逻辑照样用 velocity Verlet，只是所有操作对 batch 维度广播。

这里有几个实测下来很有用的调整：

• 控制中间张量的生命周期。diff和r2都是临时量，用einsum一步算完可以让框架少落一次显存。B 再大一点的时候，这两个中间张量的显存占用会变成瓶颈。
• batch 内系统要同构。所有系统都是三体、质量相同、步长相同，这样才能完全向量化。如果混着两体和四体，就得 padding 或者分组跑，效率立刻掉下来。
• 输出只存关键量。一万个系统 × 十万步 × 完整轨迹 = 存不下。实践里只存采样时刻的快照，或者只存统计量（能量漂移、最终逃逸数、最大间距）。这个决定要在写代码之前定好，事后改起来很痛苦。

4.3 混合精度能省一半带宽，但要看你在算什么

大批量扫描的时候，真正的瓶颈往往不是算力，是显存带宽。每步都要读写pos和vel，B=10000、N=3、双精度的话，两个张量一共 10000×3×3×8×2 ≈ 1.4 MB，跑起来带宽吃满。

这时候把位置速度降到 float32，带宽直接减半，理论上能快接近一倍。但要用对地方：

• 做参数扫描找周期解，只关心"这组初值会不会回到初始构型"这种定性判断，float32 完全够，因为判断标准是 O(1) 量级的阈值，1e-7 的精度绰绰有余。
• 做单条轨迹的长期演化，别用 float32。单精度尾数 23 位，相对精度 1e-7，按 2.1 节的公式算，预测窗口只有 16 个时间单位左右，比双精度差了一半还多。
• 混合精度是个折中：力计算用 float32（这是计算密集的部分，省时间），位置累积用 float64（这是精度敏感的部分，保窗口）。实现上就是每次算完力之后 cast 回 float64 再累加，代价是一次类型转换，通常在带宽上还是赚的。

还有一点容易忽略：float32 在归一化没做好的情况下会灾难性失效。如果位置量级在 1e8，单精度的绝对精度只有 10 左右，整个仿真直接失去意义。所以 2.2 节那一步永远排在最前面。

5. 三算失衡的排查手册

5.1 四类症状的归因对照

跑出问题的时候，"哪里错了"和"看起来哪里错了"经常不是一回事。这张表是我在实际项目里反复用到的一张对照，按症状反查根因，能省掉大量瞎猜：

第一行和第五行最容易混。前者是算法结构性问题，漂移速度基本恒定，跑多久漂多久；后者是定义不一致，漂移量通常很小但一直在，换掉软化长度后漂移量会变。区分方法是跑两段不同长度的仿真，看漂移量是不是和时长成正比。

第二行值得单独说。总动量不为零会导致整个系统以一个恒定速度平移。这个平移不影响相对构型，所以如果只做可视化或者只看相对距离，你完全发现不了。但绝对坐标全错，一旦你的判据里包含绝对位置（比如"判断某个天体是否进入某个区域"），结论就是错的。检查方法一行代码：总动量应该是零，或者至少相对于典型单体动量小几个数量级。

5.2 一次完整排查链路的现场复现

回到开头那个团队。他们的症状是"跑大约 500 步后轨迹开始和参考结果分叉，2000 步后完全对不上"。完整的排查过程是这样走的：

第一步，确认参照物本身。他们拿来对比的"参考结果"是另一套代码跑出来的。先让两套代码都用极小步长（dt 缩小 100 倍）跑 200 步，看是否一致。结果一致，说明公式和符号没写错，问题出在数值层面，不是物理层面。这一步很关键，很多人跳过它，直接在"我的代码是不是错了"和"参考是不是错了"之间反复横跳。

第二步，查守恒量。打印每步的总能量、总角动量、总动量。能量开始漂移的时刻和轨迹分叉的时刻对上了，都在第 500 步附近。总角动量完全守恒（因为是有心力，理论上应该守恒，这是个独立的一致性检查）。总动量不守恒——第一个信号出现了。

第三步，追总动量。总动量按理应该是个常量，不守恒说明每步都在被注入。逐行追下去，发现问题在速度更新：他们用了半隐式欧拉，但位置更新用的是旧速度而不是新速度。这个写法本身有动量守恒的版本，他们混用了两种。修正之后总动量稳定守恒，但轨迹还是在 500 步附近分叉。

第四步，换算法。既然能量在漂，先换成 velocity Verlet。能量从单调漂移变成有界振荡，振荡幅度符合O(dt²)预期。但轨迹仍然在 500 步左右分叉。

第五步，回头查算据。到这里只剩下一种可能了。打印初始总动量——果然不为零，量级还不小。他们从外部文件读初值，那个文件里的速度只给了六位有效数字，加起来不等于零，"不等于零"的那部分就是平移速度。修正方法很简单：在初始化之后做一次质心系变换，把总动量和质心位置都归零。

def normalize_to_com(pos, vel, mass): """平移到质心系并把总动量清零""" M = mass.sum() com = (mass[:, None] * pos).sum(axis=0) / M vcm = (mass[:, None] * vel).sum(axis=0) / M return pos - com, vel - vcm

改完之后，跑 2000 步轨迹和参考结果的分叉时刻推迟到了 800 步左右，剩下的差异才是真正的数值精度差异。

整条链路花了大概两天。前面三个月的并行优化，对精度没有一丁点帮助——这就是三算失衡的典型代价：算力这一项被拉满，算据和算法这两项是空的，整体产出为零。

6. 自己动手：三个能立刻跑起来的验证实验

6.1 用 8 字形周期解当基准尺

三体问题虽然一般不可解，但存在一些特殊的周期解，其中最著名的是三个等质量天体在一条 8 字形轨道上追逐的解。它被广泛引用，初值也公开，非常适合当基准测试。

归一化取G = 1，三个质量都是 1，初值如下（位置和速度都在平面内，z 分量为 0）：

pos0 = np.array([ [ 0.97000436, -0.24308753, 0.0], [-0.97000436, 0.24308753, 0.0], [ 0.0, 0.0, 0.0], ]) v3 = np.array([-0.93240737, -0.86473146, 0.0]) vel0 = np.array([-v3 / 2, -v3 / 2, v3]) # 总动量恰好为零 mass = np.array([1.0, 1.0, 1.0])

这组初值值得留意两件事。第一，总动量m·v加起来正好是零，这是解本身的性质，你的实现如果算出来不是零，说明初始化或者单位制有问题。第二，周期大约是 6.3259 个归一化时间单位，跑完一个周期后三个天体应该回到出发位置。这个"回到原位"的判据非常好用——不需要外部参考数据，自己就能验证。

我一般这么做：用 velocity Verlet，步长从 1e-2 开始，每个周期采样一次位置，跑十个周期。看|pos(T) - pos(0)|随周期数的变化。好的实现里这个偏差应该线性增长（浮点误差累积），不会指数增长。如果看到指数发散，说明要么算法不是辛的，要么步长太大已经跳出了稳定区。

6.2 单双精度对照跑一遍

这个实验的目的是把 2.1 节那张表变成你自己的实测数字。用同一组初值，一次用 float64，一次用 float32，打印每一步的总能量，画能量漂移曲线。

预期结果：float64 的能量振荡幅度在 1e-12 量级附近，float32 的在 1e-5 量级附近，比值大致对应两种格式的机器精度之比（约 1e-16 对 1e-7）。如果实测的比值远小于这个数，说明你的算法本身的误差（截断误差）已经盖过了浮点误差，还需要减小步长；如果比例符合，说明浮点误差是主要来源，再加步长收益不大。

这个实验还有一个副产品：你会直观地看到"预测窗口"这个概念。float32 跑出来的轨迹，在前一段和 float64 几乎重合，到了某个时刻突然分开，然后再也不重合。分开的那个时刻，就是你的精度能支撑的预测窗口。把它和 2.1 节公式算出来的理论值对一下，通常能对上量级，差个两三倍很正常，因为 Lyapunov 指数是估的。

提醒一句：做这个对照实验时，两次运行必须用完全相同的步长、完全相同的算法，只改存储精度。改两处你就分不清是谁的贡献了。

6.3 时间可逆性：一个能抓出实现 bug 的自检

velocity Verlet 有一个数学上的好性质：时间可逆。把速度取反，用同样的步长往回跑相同的步数，理论上应该精确回到起点。浮点实现做不到"精确"，但偏差应该保持在一个很小的量级，并且随步数多项式增长，不是指数增长。

具体操作：

pos, vel, _ = velocity_verlet(pos0, vel0, mass, dt, N) pos_back, vel_back, _ = velocity_verlet(pos, -vel, mass, dt, N) print(np.abs(pos_back - pos0).max()) # 应该很小

这个自检的威力在于，它能抓住很多"看起来没问题"的实现错误。比如半步速度的更新顺序写反、位置更新用了旧的加速度、力的符号搞错——这些错误在短时间的轨迹上看不出来，但会破坏可逆性，跑一遍反向立刻暴露。我自己写三体代码时，这个检查是每次改完力函数必跑的。

需要说清楚它的边界：时间可逆性检验的是实现的正确性，不是精度的好坏。一个数学上正确的 velocity Verlet 实现，不管步长多大，反向后都能回到原点附近（偏差量级取决于步长）。所以它能抓 bug，但不能告诉你"我的步长够不够"。判断步长够不够，还得看 6.1 节的周期性回归和能量振荡幅度。

三个实验跑完，大概一下午。但这三个小时建立起来的直觉，比读十篇讲数值方法的文章都管用。三体这个题目最迷人的地方也在这儿——它的复杂度不是藏在高深的数学里，而是藏在最朴素的浮点运算和一步一脚印的时间推进里，你随时可以自己动手验证每一个结论。我在实际项目里带新人的时候，一般都会让他们先把这三个实验跑一遍，跑完再碰真实的仿真代码，能省掉后面一大半的排查时间。