前阵子帮一个做物理仿真的小团队看代码他们折腾了将近三个月把一套三体数值积分从单机跑到了多卡集群结果轨迹还是对不上参考结果。我打开仓库做的第一件事不是看并行部分而是打印初始总动量——三个天体的动量加起来不是零整个系统在以一个恒定的速度朝某个方向平移。相对构型全对绝对坐标全错。他们追了三个月的算力不够根子上是一个谁都没检查的初值问题。这个故事我后面还会提到。它几乎把三体与三算这件事讲透了三体问题没有解析解只能靠数值近似往前推而数值近似的质量由三个互相牵制的东西决定——算据你喂进去的初始条件和数值精度、算法你怎么推进时间、算力你有多大的机器、怎么把它用起来。我把这三样合称三算。任何一个环节短了另外两个环节的投入都会按比例打水漂和木桶装水一个道理。这篇写给两类人。一类是被三体问题吸引、想自己动手写一套仿真但不知道坑在哪的另一类不一定关心天体力学但手上经常处理需要长期演化的复杂系统——分子动力学、流体、控制回路、种群模型都算。三体只是一个把问题放大到极致的沙盘三算这笔账在任何长时程仿真里都躲不掉。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)代入两组数字看差别初值相对误差 e0e_target 0.1 时的预测窗口1e-16双精度满精度约 34.5 个时间单位1e-12传感器中等精度约 25.3 个时间单位1e-6六位有效数字实测值约 11.5 个时间单位1e-3粗测或估算值约 4.6 个时间单位这张表最值得看的地方不是具体数值而是它的对数关系。精度从 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, G1.0, eps0.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, G1.0, eps0.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, G1.0, eps0.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 这类方法就是这么干的效果对短时仿真很好。但这里有个硬约束变步长会破坏辛结构。辛积分器的长期能量守恒性质建立在固定步长准确说是固定时间映射的前提上。步长一变每一步的映射都不一样辛性就丢了能量又开始漂移。你想同时要自适应应对近距遭遇和长期能量守恒数学上做不到。工程上的取舍通常是三选一方案适用场景代价固定步长辛积分 足够小的 dt长期定性研究、系综统计全程都在为偶发的近距离遭遇付代价慢自适应 RK短时高精度、单次遭遇事件的精确刻画长期能量漂移跑不长时间变换辛积分既要长跑又要处理大动态范围实现复杂需要构造合适的时间变换函数大多数人的实际需求落在第一行。选一个能覆盖最短特征时间的固定步长接受平时浪费一点的事实。经验做法是先跑一小段记录最小天体间距估算该间距下的轨道周期取周期的 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 加速只在几万原子以上才有明显收益而三体这种 N3 的问题多线程都没必要开。4.2 batched 思路几千个三体打包成张量一起推真正需要算力的场景不是跑一条三体轨迹而是跑一万条不同初值的三体轨迹看统计分布或者扫参数空间找周期解。这时候并行维度不该在系统内部而应该在系统之间。把 B 个独立的三体系统打包成形状(B, 3, 3)的张量一次推进所有系统def accel_batched(pos, mass, G1.0, eps0.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和velB10000、N3、双精度的话两个张量一共 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 四类症状的归因对照跑出问题的时候哪里错了和看起来哪里错了经常不是一回事。这张表是我在实际项目里反复用到的一张对照按症状反查根因能省掉大量瞎猜症状首要怀疑对象验证方法能量持续单调下降/上升算法非辛或实现错误换 velocity Verlet 重跑看是否变成有界振荡能量守恒但整体轨迹平移漂移算据初始总动量不为零打印sum(mass[:,None]*vel, axis0)应接近 0近距离遭遇后突然飞出去算法步长过大或软化不足记录每步最小间距和步长对比同参数两次运行结果不一致算力并行归约顺序不确定固定单线程、固定种子重跑对照长时间跑后能量缓慢单向漂移但幅度很小软化前后势能定义不一致检查力函数和能量函数的 eps 形式是否一致第一行和第五行最容易混。前者是算法结构性问题漂移速度基本恒定跑多久漂多久后者是定义不一致漂移量通常很小但一直在换掉软化长度后漂移量会变。区分方法是跑两段不同长度的仿真看漂移量是不是和时长成正比。第二行值得单独说。总动量不为零会导致整个系统以一个恒定速度平移。这个平移不影响相对构型所以如果只做可视化或者只看相对距离你完全发现不了。但绝对坐标全错一旦你的判据里包含绝对位置比如判断某个天体是否进入某个区域结论就是错的。检查方法一行代码总动量应该是零或者至少相对于典型单体动量小几个数量级。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(axis0) / M vcm (mass[:, None] * vel).sum(axis0) / M return pos - com, vel - vcm改完之后跑 2000 步轨迹和参考结果的分叉时刻推迟到了 800 步左右剩下的差异才是真正的数值精度差异。整条链路花了大概两天。前面三个月的并行优化对精度没有一丁点帮助——这就是三算失衡的典型代价算力这一项被拉满算据和算法这两项是空的整体产出为零。6. 自己动手三个能立刻跑起来的验证实验6.1 用 8 字形周期解当基准尺三体问题虽然一般不可解但存在一些特殊的周期解其中最著名的是三个等质量天体在一条 8 字形轨道上追逐的解。它被广泛引用初值也公开非常适合当基准测试。归一化取G 1三个质量都是 1初值如下位置和速度都在平面内z 分量为 0pos0 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 节的周期性回归和能量振荡幅度。三个实验跑完大概一下午。但这三个小时建立起来的直觉比读十篇讲数值方法的文章都管用。三体这个题目最迷人的地方也在这儿——它的复杂度不是藏在高深的数学里而是藏在最朴素的浮点运算和一步一脚印的时间推进里你随时可以自己动手验证每一个结论。我在实际项目里带新人的时候一般都会让他们先把这三个实验跑一遍跑完再碰真实的仿真代码能省掉后面一大半的排查时间。
列了顺序)
