企业官网建设流程全解析

我最早被这个题目绊住，是A同学拿一个半成品代码来问：为什么轨道一旦靠近其中一颗主星，几秒之后就彻底崩了。那个代码写的是平面限制性三体问题——两个大质量天体在同一个平面内绕公共质心匀速转动，第三个小质量天体只受它们的引力，反过来又不会影响前两者的运动。看起来是最简单的三体问题，但折腾了我将近三天。把它彻底跑通之后，我反而觉得这个模型是整个轨道力学里最值得做仿真的题目之一：五个引力平衡点的存在、转移轨道、借力飞行、还有混沌，都可以在一台普通笔记本上用几段代码复现出来，非常适合课程设计、航天科普，也适合第一次接触非线性动力学的人上手。

下面我不会讲太多推导，重点放在“为什么这么做”“代码怎么落地”“跑出来怎么判读”，把我踩过的坑一并写出来。

1. 这个问题到底在仿什么：限制性三体模型的物理边界

1.1 从N体问题到限制性三体问题的降维过程

真实的三体问题，三个天体质量差不多，彼此引力互相影响，那是一个真正的“混沌乐园”，数值上极其敏感，解析解更是想都别想。但工程和科研中大多数情况根本不是这种。比如一颗探测器绕地月系统运动，它对地球和月球的引力微乎其微，反过来却要被两者同时牵引；又比如某个小行星待在太阳和木星之间，也完全符合“小质量天体被两个大天体带着走”的场景。

限制性三体问题就是从这里来的一次漂亮降维：假设两个主天体质量分别为 m1 和 m2，它们的运动只受彼此引力影响，并且轨道是圆；第三个小天体质量 m3 可以忽略不计，只感受 m1、m2 的引力，不产生任何反作用。这样一来，主天体的运动可以先独立解出来，剩下的事情全部集中到小质量天体身上，自由度少了一大截。

更进一步的“平面”假设，是把小质量天体限制在 m1、m2 轨道所在的平面内。真实航天器当然有垂直分量，但很多关键现象，比如平动点稳定性、引力辅助转向、转移轨道设计，在平面模型里已经能展示出绝大部分物理本质。先在一个二维平面上把问题研究清楚，再往三维扩展，这是轨道力学里非常成熟的路径，不是偷懒。

1.2 旋转坐标系与无因次化的关键是“换视角”

我第一次自己实现这个模型时，第一版写的是惯性坐标系，就是太阳不动、行星动的那套视角。结果两条主星的轨迹本身就绕来绕去，小质量天体所受引力的方向也一直在变，画出来的图又乱又难以判断。后来才意识到，限制性三体问题里最重要的一个技巧，是跳到一个“跟主天体一起旋转”的坐标系里去看。

在这个旋转坐标系中，两颗主星是固定在 x 轴上的，一个在左侧，一个在右侧。这样引力的方向不再随时间变化，问题从时间依赖的“动问题”变成了近似“静电场”的结构。当然，天上不会掉馅饼，换到旋转系要付出代价——必须补上两个惯性力项：离心力和科里奥利力。别看到这两个词就跑，它们在方程里只是“多出来的两项”，但缺了它们，你算出来的轨道一定飞掉。

同时还有一招叫无因次化。把所有长度用两颗主星之间的距离去度量，把质量用总质量去度量，把时间用系统轨道角速度的倒数去度量。这样方程里不再出现像 3.84e8 米、6.67e-11 这些容易把数值搞炸的大数，只留下一个无量纲参数 μ，代表小主星质量占总质量的比值。比如一个接近真实地月系统的参数 μ=0.01215，整段代码只需要输入这一个常数，其余自动变成 0.1 到 1 之间的友好量级。

提醒：无因次化不是“可选的程序优化”，它决定了数值积分器的生死。你把 SI 单位直接带进 RK4，光量纲就差了十几个数量级，精度再高的算法也救不回来。

1.3 运动方程和守恒积分C

旋转坐标系下的无量纲运动方程，形式上非常紧凑：

• x'' - 2 y' = ∂Ω/∂x
• y'' + 2 x' = ∂Ω/∂y

其中伪势函数 Ω 是：

Ω = (x² + y²) / 2 + (1 - μ) / r1 + μ / r2

r1 是小天体到第一颗主星的距离，r2 是到第二颗主星的距离。真正给计算机写代码时，要把二阶方程拆成一阶方程组。设 vx = x'，vy = y'，于是有四个状态分量 x、vx、y、vy，右侧表达式如下：

• x' = vx
• y' = vy
• vx' = 2 vy + x - (1 - μ)(x + μ) / r1³ - μ(x - 1 + μ) / r2³
• vy' = -2 vx + y - (1 - μ)y / r1³ - μ y / r2³

这个系统里有一个非常重要的守恒量，通常记作 C：

C = x² + y² + 2(1 - μ) / r1 + 2 μ / r2 - (vx² + vy²)

它在任意时刻都应保持恒定。这个东西在后面的仿真调试里几乎是“照妖镜”：轨迹画得像模像样，但 C 漂得一塌糊涂，说明数值积分精度不够；C 保持得很稳，说明整条轨道大概率是可信的。所以每跑一个算例，我都会顺手把 C - C0 画出来看看，这个习惯帮我排掉了不少问题。

2. 数值积分选型：为什么我不能只做RK4

2.1 为什么不能只靠解析解

很多人第一次听说三体问题，都会问一句：为什么不能用公式直接算出轨道？因为限制性三体问题虽然已经是“限制性”的简化模型，但它依然没有通用解析解。能解析处理的情况极其有限，比如五个平动点本身的位置，以及平动点附近做线性化后的小振幅近似周期运动。真实轨道，无论是绕 L4 的大范围运动，还是从一颗主星附近飞向另一颗主星的转移，都只能靠数值积分一步一脚印地推进。

所以选对数值积分器，是这个仿真的第一道关卡。我在刚开始时直接用最朴素的欧拉法，那是为了验证方程推导有没有问题；但如果拿欧拉法去算周期轨道，哪怕步长已经取到 0.001，也会出现明显的“向外螺旋”假象。原因很简单：欧拉法每步都会引入系统性的能量误差，长时间累积后完全不可接受。

2.2 RK4固定步长：经典的起点

比欧拉法进一阶的自然是 RK4，四阶龙格-库塔方法。核心思想是在一个步长内取四个不同位置的斜率，然后加权平均，让局部截断误差达到 h 的五次方量级，全局误差是 h 的四次方量级。对大多数短时间演示来说，它已经足够用。

我自己写的一个固定步长 RK4 迭代函数是这样，非常简单：

def rk4_step(f, t, y, h, mu): k1 = f(t, y, mu) k2 = f(t + h/2, y + h/2 * k1, mu) k3 = f(t + h/2, y + h/2 * k2, mu) k4 = f(t + h, y + h * k3, mu) return y + h/6 * (k1 + 2*k2 + 2*k3 + k4)

步长怎么选？无量纲时间走一个主星轨道周期大约是 2π，约 6.28。如果你想看十来圈的轨道，dt=0.01 意味着每圈 628 步，总计算量只要几千步，毫秒级完成。但固定步长的隐患在于：小质量天体一旦飞到某颗主星附近，引力项的分母 r1 或 r2 变得很小，加速度瞬间暴涨。这时的局部误差会迅猛增加，而 RK4 不会自主察觉，只能靠人把 dt 整体调小。

实测下来，如果只是让小天体围绕 L4 点附近画几圈，dt=0.01 完全够；但如果要模拟飞掠小主星表面附近的超高速转向，同样的 dt 很容易在某一步跨过了“极剧烈引力变化区”，轨道直接弹飞。所以固定步长只适合“安全区”里的短时仿真。

2.3 RKF45变步长：实战主力

真正干活的时候，我更推荐自带自适应步长的 RKF45。它同时用四阶和五阶两种估计结果，比较它们之间的误差，从而决定当前步长是放大还是缩小。靠近主星时，误差估算变大，步长自动收紧；远离主星时，步长又自动放开。这比手动调固定步长聪明得多。

用 SciPy 的 solve_ivp 实现非常简单：

import numpy as np from scipy.integrate import solve_ivp mu = 0.01215 def rtbp(t, y, mu): x, vx, yy, vy = y r1 = np.sqrt((x + mu)**2 + yy**2) r2 = np.sqrt((x - 1.0 + mu)**2 + yy**2) ax = 2.0 * vy + x - (1.0 - mu) * (x + mu) / r1**3 - mu * (x - 1.0 + mu) / r2**3 ay = -2.0 * vx + yy - (1.0 - mu) * yy / r1**3 - mu * yy / r2**3 return [vx, ax, vy, ay] sol = solve_ivp(rtbp, [0, 20 * np.pi], [0.48785, 0.0, 0.86603, 0.0], args=(mu,), rtol=1e-10, atol=1e-12, method='RK45', dense_output=True)

这里的 rtol 和 atol 直接决定精度。我个人的经验是，rtol=1e-10、atol=1e-12 是一个“甜点区间”：既能保证十几圈轨道不飘，也不会慢到让人烦躁。如果再严到 1e-14，计算时间成倍增长，画面上的提升肉眼根本看不出来。如果你要算上百圈的长程演化，可以把 method 换成 DOP853，即高阶的 dormand-prince 8 阶方法，精度更高，耗时会多一些。

2.4 守恒量监控与辛积分的概念

就算用 RKF45，时间长了一样会有能量漂移。毕竟任何基于 RK 格式的方法都会在相空间中引入微小但持续的人工耗散。这里的“能量”我用守恒量 C 来代表。C 恒定的程度，比轨迹是否漂亮更能说明数值质量。

对绝大多数课程设计、科普演示来说，RKF45 + 严格容差已经非常够用。但如果你做的是长期轨道族研究，比如让小天体在 L4 附近绕几百上千圈，就要考虑换辛积分器。辛积分的核心概念是保持相空间体积不变，专门为哈密顿系统设计，能量误差不会持续线性增长，而会在某个小范围内来回震荡。

三种方法直观对比：

先说清楚，我不建议新手一上来就上辛积分。那需要理解生成函数、配置参数等一系列概念，对于第一次做平面限制性三体仿真的人来说太劝退了。先把 RKF45 用顺手，跑出稳定轨道，再去扩展精度，是更平滑的学习曲线。

3. 核心仿真代码实现：从零搭一套能直接跑的模型

3.1 参数、状态量与初始条件

写代码前先把约定铺好。两颗主星的质点圆心相对质心距离，在无量纲坐标下固定为：第一颗主星位于 (-μ, 0)，第二颗主星位于 (1-μ, 0)。当 μ=0.01215 时，前者在 (-0.01215, 0)，后者在 (0.98785, 0)。小天体的状态向量我统一用 [x, vx, y, vy] 这样的顺序，即 x 坐标、x 方向速度、y 坐标、y 方向速度。顺序其实没有标准，但一定要在代码里固定下来，不然后面索引记错就全是错乱。

初始条件的单位问题是最容易翻车的。所有量都是无量纲的，比如速度的单位是“主星间距除以时间单位”，不是 km/s。如果你从真实星历里拿到一个速度数值，直接代进去，结果一定不对。要先把真实量换算到无量纲量：长度除以两主星距离，时间乘以轨道角速度，速度按“单位长度 / 单位时间”缩放。

注意：方程里的速度是相对旋转坐标系的。假设你想让小质量天体在旋转坐标系里完全静止，那么初始速度 vx=0、vy=0 就对了。但这个“静止”对应惯性坐标系里的一个圆周运动速度，两者不能混着理解。

3.2 微分方程与积分器对接

完整可跑的代码如下，这部分是我实际在笔记本上调通的版本，加了注释方便对照：

import numpy as np from scipy.integrate import solve_ivp import matplotlib.pyplot as plt # 参数：接近真实地月系统的小主星质量占比 mu = 0.01215 def rtbp(t, y, mu): x, vx, yy, vy = y r1 = np.sqrt((x + mu)**2 + yy**2) r2 = np.sqrt((x - 1.0 + mu)**2 + yy**2) ax = 2.0*vy + x - (1.0 - mu)*(x + mu)/r1**3 - mu*(x - 1.0 + mu)/r2**3 ay = -2.0*vx + yy - (1.0 - mu)*yy/r1**3 - mu*yy/r2**3 return [vx, ax, vy, ay] # 初始状态：L4平动点附近 xL4 = 0.5 - mu yL4 = np.sqrt(3)/2.0 state0 = [xL4 + 1e-3, 0.0, yL4, 0.0] # 积分时长 20π，约为10个主星轨道周期 t_span = [0, 20 * np.pi] sol = solve_ivp(rtbp, t_span, state0, args=(mu,), rtol=1e-10, atol=1e-12, method='RK45', dense_output=True) # 旋转坐标系轨迹 plt.figure(figsize=(7, 7)) plt.plot(sol.y[0], sol.y[2], lw=0.8) plt.scatter([-mu], [0], s=80, marker='o', label='primary 1') plt.scatter([1-mu], [0], s=40, marker='o', label='primary 2') plt.xlabel('x') plt.ylabel('y') plt.axis('equal') plt.legend() plt.show()

dense_output=True 让解在任意时间点都能被插值，后面做惯性坐标系变换、画动画都很方便。注意 solve_ivp 的返回结果里，sol.y 是一个形状为 (4, N) 的数组，sol.y[0] 和 sol.y[2] 分别对应整条轨迹的 x 和 y 序列。索引顺序错位是高频 bug，我在这上面栽过两回。

3.3 可视化：旋转系、惯性系、能量监察

旋转系轨迹是最直接的展示方式，两条主星固定在画面里，小天体绕它们画出各种漂亮图案。但很多时候你想理解“实际看起来怎么样”，那就得把结果变换回惯性系。因为是圆轨道假设，无量纲角速度为 1，旋转矩阵很简单：

t = sol.t X = sol.y[0] * np.cos(t) - sol.y[2] * np.sin(t) Y = sol.y[0] * np.sin(t) + sol.y[2] * np.cos(t)

这样画出来的惯性系轨迹往往是复杂的编织曲线，看起来和旋转系里的结果完全不像同一个系统。我第一次画出来时都有点恍惚，这就是换参考系的威力：同样的物理过程，看问题的角度不同，直观感受天差地别。

能量监察代码也建议一开始就写好。计算守恒量 C 的序列，然后画它们相对初始值的偏差：

def energy_c(y, mu): x, vx, yy, vy = y r1 = np.sqrt((x + mu)**2 + yy**2) r2 = np.sqrt((x - 1.0 + mu)**2 + yy**2) return x**2 + yy**2 + 2*(1-mu)/r1 + 2*mu/r2 - (vx**2 + vy**2) E = np.array([energy_c(sol.y[:, i], mu) for i in range(sol.y.shape[1])]) plt.plot(sol.t, E - E[0]) plt.xlabel('time') plt.ylabel('C - C0')

这条偏差曲线应当在零附近做一个窄幅抖动。如果呈现明显的单调增长或下降趋势，说明必须在精度上做文章。

3.4 平动点求解与稳定域初探

五个平动点的位置是仿真里最常用来“开刀”的对象。四个中有两个是解析的：L4 在 (0.5 - μ, √3/2)，L5 在 (0.5 - μ, -√3/2)。它们和两颗主星构成等边三角形，这个结论在很多教材里都有，我自己证过一遍，其实就是让伪势的梯度为零后，发现 r1=r2=1 的等边构型是解。

剩下 L1、L2、L3 在 x 轴上，需要数值求根。用小段二分法就能解决。函数定义和求根如下：

def f_col(x, mu): r1 = abs(x + mu) r2 = abs(x - 1.0 + mu) # 避免 r1 或 r2 为零 if r1 < 1e-12 or r2 < 1e-12: return 1e10 return x - (1-mu)*(x + mu)/r1**3 - mu*(x - 1.0 + mu)/r2**3 def bisection(f, mu, a, b): fa = f(a, mu) for _ in range(80): c = 0.5*(a + b) fc = f(c, mu) if fa * fc < 0: b = c else: a = c fa = fc return 0.5*(a + b) # L1 位于两主星之间，区间取 (-mu*0.5, 1-mu*0.5) # L2 位于小主星外侧，区间取 (1-mu+1e-6, 10) # L3 位于大主星外侧，区间取 (-10, -mu-1e-6)

实际使用中，L1 求根区间最容易撞到 r2=0，也就是小主星质心位置本身。所以我在函数里加了保护逻辑，如果距离过小直接返回一个很大的数，避免二分法在奇点附近疯狂振荡。这个细节花了好久才发现。

求完平动点后，最简单的稳定性检验就是：把初值设为平动点，速度设为零，积分几个周期，看小天体是否一直“粘”在那里。五个点都不会漂移，因为它们本身是平衡点。接着把初值在某个方向上偏移 0.01，观察轨道的形态，就能直观感受到稳定域的形状。

4. 三个典型算例的物理复盘

4.1 算例1：L4平动点附近的周期轨道

第一个算例，把小天体放在 L4 附近，初始位置加一个 0.001 量级的偏移，速度为零。RKF45 积分 20π 后，旋转系里的轨迹会形成一个围绕 L4 的“花环状”闭合曲线。这个形状既不是圆也不是椭圆，而是两个本征频率叠加出来的类李萨如图形。

我做了几组不同偏移量的对比。偏移量越小，轨迹越紧贴 L4，接近线性理论预测的小椭圆；偏移量慢慢增大后，轨道变为更复杂的环状，但仍然困在 L4 附近的稳定域内。这背后和平面模型中 L4 点的稳定性判据有关：当质量比 μ 小于约 0.0385 时，L4/L5 在平面内是线性稳定的。真实地月系统参数 μ=0.01215 恰好落在这个区间。

这个算例最大的价值是给人信心：稳定性不是嘴上说说，而是真的能在数值模拟里看见“小天体绕着某个看不见的点转圈”。我在跑通这条轨道后，又顺手改了 μ=0.3 去测 L4，得到的轨迹很快就飞走了，正好反过来说明参数对稳定域的剧烈影响。

4.2 算例2：小主星附近的借力效应

第二个算例更有意思：把初始点设在大主星左侧，给一个朝向小主星区域的速度，让轨迹恰到好处地从小主星身边掠过。积分过程时，你会发现小质量天体经过小主星附近时，轨道完成一个非常锋利的转向，像台球撞库一样改变了运动方向。

从旋转坐标系看，小主星是固定在 x 轴右侧的。小天体飞过去时，引力会先把轨迹往小主星方向拉，然后由于速度太快又会“甩”出去，出射方向与入射方向完全不同。这就是引力辅助的微观机制：航天器并没有获得什么神秘能量，它只是在接近行星时借助行星的引力场改变了相对速度的方向；如果变换到惯性系，还能看到它相对大主星的轨道能量被显著改写。

复现这个算例的技巧是调整初始速度。我的做法是写一个小的循环，对 vx 从 0.3 扫描到 0.8，固定 vy，然后看哪些初值会让轨迹靠近小主星质心 0.05 以内的区域。挑中一个合适初值后，再把步长精度调高，让近距离散步的细节不被数值误差抹掉。

4.3 算例3：微小扰动导致的混沌轨道

第三个算例，我建议每一个做仿真的人都亲自跑一次。方法特别简单：取一个离 L4 较远、但又没远到立刻飞走的初值，生成两条初始条件相差 1e-8 的轨道。前几个时间单位内，两条轨迹几乎重合；到中途开始出现肉眼可见的分离，后面直接变成完全不同的缠绕图案。

可以这样检查分离速度：计算两条轨迹在相空间里的距离，每隔一段时间取一次值：

dist = np.sqrt((x1 - x2)**2 + (y1 - y2)**2)

如果距离呈指数级增长，这就是典型的混沌现象。而且要注意，这里加的 1e-8 差异，还没数值积分器本身的舍入误差大。换句话说，哪怕你不故意加扰动，只要积分时间足够长，数值误差本身就会让轨道偏离原方向。这就是为什么在长时间仿真里监控守恒量 C 变得极其重要——数值误差不是“小问题”，在混沌系统里，它会被系统本身成指数放大。

我在一个深夜里盯着两条轨迹发散的过程看了很久，有一种“书上说的敏感依赖终于变成了眼前画面”的直观冲击。如果你不信混沌，这个算例能说服你。

5. 常见问题排查与参数调试实录

5.1 一算就发散/NaN

这是我被问得最多的问题。轨道刚出发几步就飞掉，或者干脆出现 NaN，大概率不是方程写错，而是数值参数踩雷。我列一个排查顺序：

整套检查的次序也很重要：先确认初始条件没落在天体内部，再确认无因次化正确，最后才去怀疑积分器。很多时候问题出在最基础的单位上。

5.2 能量C漂移明显

守恒量 C 的理想表现是围绕零做微小振荡，而不是单调上涨或下跌。如果你画出的是明显的趋势线，第一反应不是怀疑物理，而是怀疑数值积分精度。

固定步长 RK4 最常见的症状就是 C 值随计算步数线性漂移。换成 RKF45 并设置严格容差后，漂移通常会被压到可接受范围。对于一个 10 个主星周期时长的仿真，C 的偏差最好控制在 1e-6 以下。如果做不到，先检查是不是经过了小主星的近距离区域，那里的引力梯度极大，是误差的主要来源。

如果想给观众演示一个“数字很漂亮”的能量曲线，我通常在画图时只显示 C-C0，并设置纵轴为 1e-8 量级，这样能看出微小的振荡细节。如果纵轴默认从零拉到很大，反而看不出趋势变化，容易被误导。

5.3 初始速度方向怪

这个坑特别隐蔽。很多初学者计算初始速度时，参考的是惯性系里的想法。比如想让探测器绕大主星画圆，就设一个切向速度，但因为整个方程组是在旋转坐标系里，科里奥利力项会带来额外的加速度，实际轨迹完全不是预期的样子。

旋转坐标速度和惯性坐标速度的换算，在平面模型中很简单：

• v_rot_x = v_inertial_x + y
• v_rot_y = v_inertial_y - x

注意这是无量纲角速度为 1 的情况。如果你的初始速度来自某个惯性系参考，先做这一步换算再放进 solve_ivp。实践中最稳妥的做法是：默认初始速度给 0，先看图里小天体怎么运动，再一点一点加速度，观察方向变化。这样能帮你建立起对科里奥利力的直觉。

5.4 时间尺度换算

很多人在最后一步发现自己画图的横轴含义理解错了。无量纲时间 1 个单位，对应的是系统轨道角速度的倒数，不等同于“一个轨道周期”。一个轨道周期是 2π 个无量纲时间单位，所以积分 20π 其实就是 10 个主星轨道周期。

如果你想把仿真结果对应到真实时间，需要一个缩放因子。假设双星系统的真实轨道周期是 T，那么 1 个无量纲时间单位等于 T/(2π) 天。比如一个周期约 27 天的系统，一个无量纲单位大约对应 4.3 天。这个换算在这些轨道可视化、与观测数据对比时特别有用。

6. 仿真结果背后：平面模型能用来干什么

6.1 教学、课程设计与科普实验

如果你的目标是课程设计或者科普演示，平面限制性三体仿真几乎是为这些场景量身定做的。代码量控制在两三百行以内，运行时间几秒到十几秒，一台普通笔记本就能轻松跑完。用它来做实验，我推荐几个直接能落地的方向：

• 验证 L1-L5 五个平动点的位置与稳定性，给每个点各跑一组偏置初值，对比轨迹形态
• 手动扫描初始速度，试图找出一条从大主星附近通往小主星附近的转移轨道，画出来会非常有成就感
• 展示混沌现象，把两条初始差异 1e-8 的轨道叠加在同一张图上
• 改变 μ，比较地月系统和近对称系统之间平动点稳定域的变化

这些实验不涉复杂代码，但每一步都能把轨道力学里最核心的概念“演”出来，比只看公式要直观得多。我在给A同学讲解时，就是用一张 L4 轨道图和一张混沌分离图让他快速建立了整体认知。

6.2 从平面模型向真实场景扩展的路线

平面限制性三体只是一个地基。真正做工程任务时，还需要把很多简化假设逐一拆掉：把 z 方向自由度打开，变成三维限制性三体问题；把主星圆轨道改成椭圆轨道，变成椭圆限制性三体模型；再把小质量天体对主星的反作用加回来，回到完整 N 体问题。每扩展一步，计算代价和调试难度都会上一个台阶。

我的建议是，先别急着一步到位。平面模型最大的优势在于“画得清、看得懂、跑得快”。你把平面版彻底吃透，理解旋转坐标系和守恒量 C 的意义后，再看三维模型和椭圆模型，会清楚每一步扩展到底买回了什么、留下了什么。

从我自己的经验看，空间任务轨道设计的前期物理分析，很多都在用这个级别的模型做初步概念验证和快速筛选。平面限制性三体仿真正是一个“肉眼可见”的物理实验室，让你在屏幕前亲手把天体力学那些耳熟能详的概念一个个测出来。

最后分享一个使用心得：每跑一个算例前，都把初始状态和参数存成一份记录，我习惯用一个字典写在代码头部，比如初始坐标、初始速度、积分时长、容差。后面回看调参过程时，能精确复现每一个版本。我在调 L4 轨道那几天反复改初始速度，多亏留了记录，否则早上算出的好轨道到下午就找不回来了。

如果你也想感受三体系统最迷人的一面，建议从 μ=0.01215 的参数开始，先复现 L4 稳定轨道，再复现混沌分离，两个算例跑通，你对这个模型的理解立刻就不一样了。