4. 多体引力运动的数学表述
要把“多个天体在相互引力下如何运动”变成计算机能处理的问题，第一步是写出严格的数学表述。本节围绕四个相互关联的要点展开：N体问题的运动方程、质心坐标系与相对坐标系的选择、引力常数与单位制的统一，以及数值求解的基本思路。它们共同构成从物理定律到可执行程序的桥梁。
一、N体问题的运动方程
设系统中有 N 个质点，第 i 个质点的质量为 mi，位置矢量为 ri（相对某一惯性参考系），速度为 vi。根据牛顿万有引力定律与第二定律，第 i 个质点受到其余所有质点的引力合力为：
mi d²ri/dt² = G Σj≠i mi mj (rj − ri) / |rj − ri|³
两边约去 mi，得到每个质点的加速度：
a_i = G * Σ_{j≠i} m_j * (r_j - r_i) / |r_j - r_i|^3这个形式非常重要：加速度只依赖其他天体的质量与相对位置，与自身质量无关（引力质量与惯性质量相等）。这也意味着，若把方程写成状态空间形式，每个天体需要 6 个状态量（三维位置 + 三维速度），N 个天体共 6N 个一阶微分方程：
dr_i/dt = v_i
dv_i/dt = G * Σ_{j≠i} m_j * (r_j - r_i) / |r_j - r_i|^3需要注意分母中的 |rj − ri|³：分子提供方向，分母三次方中的一次用于抵消距离平方反比，另外两次来自单位矢量归一化。实际编程时通常先算距离平方 d²，再算 invd3 = 1/(d²·√d²)，可减少一次开方并避免重复计算。
二、质心坐标系与相对坐标系
N 体系统整体存在平动，若直接在惯性系中积分，质心会以恒定速度漂移，数值上容易积累大数相减的误差。因此常把参考系原点移到系统质心：
Rc = (Σ mi ri) / (Σ mi)
在质心系中，Σ mi ri′ = 0，质心静止于原点。这样位置量级由系统尺度决定，而不是由绝对坐标决定，数值稳定性更好。质心系适合研究整体角动量、总能量守恒以及孤立系统的长期演化。
相对坐标系则把某个天体（如中心恒星）固定在原点，其他天体用相对它的位置表示。它适合行星系统、卫星系统等“一个大质量中心 + 多个小质量伴随体”的场景，物理图像直观，也便于与观测数据（如相对中心天体的轨道根数）对应。但要注意：中心天体并非严格惯性系，若伴随体总质量不可忽略，需引入间接项或改用质心系。
三、引力常数与单位制统一
引力常数 G ≈ 6.674×10⁻¹¹ m³·kg⁻¹·s⁻²。在太阳系尺度下，用米和千克会让数值跨度极大（10¹¹ 米量级的位置、10³⁰ 千克量级的质量），浮点运算的有效位数容易被浪费。常用做法是采用天文单位制：长度用 AU，质量用太阳质量 M☉，时间用天或年。此时 G 的数值随之改变。
若取 AU、M☉、年为单位，则 G ≈ 4π² ≈ 39.478（因为地球公转周期约 1 年、轨道半径约 1 AU，由开普勒第三定律 GM = 4π²a³/T² 可得）。若时间用天，则 G ≈ 0.01720209895² × ... 实际常用值约为 2.959×10⁻⁴ AU³·M☉⁻¹·d⁻²。关键是：位置、速度、质量、G 必须属于同一单位制，否则加速度量纲错误，轨道会立刻发散或塌缩。
| 单位制 | 长度 | 质量 | 时间 | G 近似值 |
|---|---|---|---|---|
| 国际单位制 | m | kg | s | 6.674×10⁻¹¹ |
| 天文单位制（年） | AU | M☉ | 年 | 4π² ≈ 39.478 |
| 天文单位制（天） | AU | M☉ | 天 | ≈ 2.959×10⁻⁴ |
实操中建议在程序入口统一声明单位制，所有输入参数先换算到该单位制再进入积分器，输出时再换回习惯单位。速度单位也要一致：AU/年 与 AU/天 相差约 365 倍，混用会导致轨道周期错误。
四、数值求解的基本思路
N 体方程没有一般解析解，只能数值积分。基本思路是：把二阶方程降为一阶状态向量 y = [r₁, v₁, r₂, v₂, …, rN, vN]，写成一阶常微分方程 dy/dt = f(t, y)，再用积分器逐步推进。每一步的核心计算是“由当前所有位置算加速度”，复杂度为 O(N²)；若 N 很大，可用树算法或快速多极子降到 O(N log N)，但小规模教学示例用直接求和即可。
积分器选择上，短时演示可用四阶龙格–库塔（RK4），每步需四次加速度求值，精度较高；长期轨道演化更推荐辛积分器（如 leapfrog / 速度 Verlet），它保持相空间体积、能量误差有界振荡而不长期漂移。步长 h 需远小于系统最短特征时间（如最内层轨道周期），可用自适应步长在近距交会时自动缩小。
上图把四个环节串成一条可执行的链路：物理模型给出加速度公式，参考系决定坐标如何取，单位制保证量纲一致，状态向量与积分循环完成推进，最后用能量和角动量守恒做结果校验。实际编码时，建议先写一个两体圆轨道测试用例：取 G=4π²、中心质量 1、轨道半径 1，初速度 2π，积分一年后位置应回到起点附近，以此确认单位制、参考系和积分器都没有接错。
总结来说，N 体问题的数学表述并不复杂，难点在于把方程、参考系、单位制和积分器作为一个整体来设计。方程决定物理正确性，参考系决定数值稳定性，单位制决定量纲一致性，积分器决定长期精度。四者对齐之后，后续的轨道计算、误差分析和可视化才有可靠基础。