简介：这份C++工程面向天体物理模拟与数值算法学习者，围绕Barnes-Hut算法高效求解N体问题。它用八叉树对三维空间分层划分，并按远近阈值以节点质心和总质量近似远处粒子团，从而把传统逐对引力计算的复杂度显著降低；从粒子模型、八叉树构建到引力计算与时间推进，完整呈现了近似算法的实现链路，代码中还能看到欧拉、龙格-库塔、Adams-Bashforth等积分器与SDL可视化窗口的组织方式，适合研究多体模拟、空间数据结构与性能优化。压缩包共60个文件，约54KB，以21个cpp源文件和22个h头文件为主体，另有4个mk构建脚本、工程配置与说明文档，覆盖粒子模型、八叉树、引力计算、积分推进和窗口渲染等模块。已有329人学习，可作为课程设计、算法复现和C++工程实践的可编译底稿，便于对照理解近似精度与计算效率的权衡。

1. 一万个天体、每一步五千万对引力：Barnes-Hut 省掉的是哪部分

一万个天体放进二维方框做引力演化，每一步要算多少对引力？N(N-1)/2 ≈ 5×10⁷。单步还能忍，一套轨道跑十万步就是 5×10¹² 次开方和除法，一台工作站从天亮跑到天黑。Barnes-Hut 算法不改变物理，它把「远处的天体团」折叠成一个质点，让每个粒子的受力计算从一万次配对降到几百次，整体复杂度从 O(N²) 压到 O(N log N)。

这套 Barnes-Hut Simulator 的内核就三件事：四叉树（三维是八叉树）、开角判据、蛙跳积分器。它拿的是「可接受的近似误差换十倍以上速度」这笔账。适合写星团演化、暗物质晕塌缩、分子动力学粗粒化、粒子群可视化的人；不适合需要长时间精确轨道的双星系统，也不适合 N<1000 的场景——那种规模直接求和反而更快，还少一层玄学。

下面按「树怎么建 → 力怎么算 → 时间怎么走 → 哪里会崩」的顺序拆开，代码用 Python + NumPy 写，重点在结构而不在语言，换成 C++ 或 Rust 时一一对应即可。

2. 从粒子表到四叉树：Barnes-Hut 的节点里到底存什么

2.1 一个节点必须携带的五样东西

四叉树节点的最小字段集是：子树总质量 mass、子树质心 (mx, my)、节点半边长 half、四个子节点指针 kids、以及叶节点持有的粒子下标数组 body。五个字段各自对应一件事，少一个都要在后面付代价。

mass 和 mx/my 是远场近似的全部依据。当判据认定某个节点「足够远」，我们用它整体的质心和一个等效质量代替它内部所有粒子，这一步等价于多极展开只保留零阶（单极）项，误差量级由下一节的开角判据控制。half 提供空间尺度：判据里的 s 就是 2×half，也就是节点边长。kids 决定还能不能继续细分，body 决定叶子上要老老实实做几次精确配对。

这里最容易踩的坑是把「节点几何中心 (cx, cy)」当成质心用。几何中心是建树时切分出来的方框中心，质心是质量加权平均，两者在节点内质量分布不对称时可以差出半个边长。用错的那个，力会有系统性偏差，而且不会报错，只会让轨道慢慢漂——这种血泪经验后文还会专门讲。

import numpy as np class QuadNode: __slots__ = ("cx", "cy", "half", "mass", "mx", "my", "kids", "body") def __init__(self, cx, cy, half): self.cx, self.cy, self.half = cx, cy, half self.mass = 0.0 # 子树总质量 self.mx, self.my = 0.0, 0.0 # 子树质心（质量加权），不是几何中心 self.kids = None # None 表示叶节点 self.body = None # 叶节点持有的粒子下标数组

__slots__在这里不是炫技：N=10⁵ 时节点数大约是粒子数的两到三倍，每个节点省掉一个__dict__，内存占用和缓存命中率都会明显不同。这一点在 Python 里只是小优化，在 C++ 里换成紧凑数组就是数量级差别。

2.2 自顶向下按象限分治，别逐个插入

建树有两种写法。一种是逐个粒子插入，每遇到已占用的叶子就分裂；另一种是先算全体包围盒，再按象限把粒子下标数组递归切分。前者写起来顺手，后者快得多，而且在 Python 这种解释执行的环境里差距会被放大。

原因在于下标数组的切分是一次向量化的布尔索引，一次处理一整批粒子；逐个插入则是 N 次指针跳跃，每次都要走一遍从根到叶的路径。用分治写法，建树整体是 O(N log N)，递归深度约 log₄N。

def build_tree(pos, mass, leaf_cap=1, max_depth=24): lo, hi = pos.min(axis=0), pos.max(axis=0) center = 0.5 * (lo + hi) half = 0.5 * float(np.max(hi - lo)) + 1e-9 # 取最大跨度，保证是正方形 root = QuadNode(center[0], center[1], half) _build(root, np.arange(len(pos)), pos, mass, leaf_cap, max_depth) return root def _build(node, idx, pos, mass, leaf_cap, depth): w = mass[idx] node.mass = float(w.sum()) node.mx = float((w * pos[idx, 0]).sum()) / node.mass node.my = float((w * pos[idx, 1]).sum()) / node.mass # 两个终止条件：粒子数不超过叶容量，或已到深度上限（重合点保护） if len(idx) <= leaf_cap or depth == 0: node.body = idx return rel = pos[idx] - np.array([node.cx, node.cy]) # 象限编号：x 正负 + 2 * y 正负 → 0..3，对应左下/右下/左上/右上 quad = (rel[:, 0] > 0).astype(np.int8) + 2 * (rel[:, 1] > 0).astype(np.int8) h = node.half * 0.5 offsets = ((-h, -h), (h, -h), (-h, h), (h, h)) node.kids = [] for q in range(4): sub = idx[quad == q] ox, oy = offsets[q] child = QuadNode(node.cx + ox, node.cy + oy, h) if len(sub): _build(child, sub, pos, mass, leaf_cap, depth - 1) node.kids.append(child)

代码逻辑顺着读：先把这个节点负责的那批粒子的总质量和质心算出来（这是后面前序遍历力的前提），再判断要不要继续分。要分就把粒子按所在象限打标签，四个子节点各自拿走一批，递归下去。

参数上有三个可以调。leaf_cap是叶容量，取 1 时树最深、精度最高、节点数最多；取 4~8 可以显著减少节点总数和递归开销，代价是叶子内部仍然要做几次精确配对——因为叶子上的粒子没法再近似了。实际项目里我一般先用leaf_cap=1跑通，确认精度对得上再往上调。max_depth=24是重合点的后悔药：如果两个粒子坐标完全一致，象限切分永远分不开它们，没有深度上限就会一路递归到栈溢出。取 24 意味着最小格子边长是初始边长的 2⁻²⁴，早就低于双精度在这个尺度上的分辨能力了，到此为止是合理的。

还有一处容易忽略的细节：根节点的half用的是max(hi - lo) / 2而不是分别取 x、y 跨度。因为四叉树的子节点是正方形，根节点不是正方形的话，分割出来的子节点边长会不一致，开角判据里的 s 就失去了统一含义。

2.3 质心必须由子节点回算，不能靠插入时累加

上面这份实现里，质心是在每次_build入口用mass[idx]直接对这批粒子求加权平均的，所以严格来说不存在「回填」这一步。如果你采用逐个插入的写法，就必须在每次分裂时把父节点原有的质量重新分配到子节点，并在整棵树建完后做一次后序遍历自底向上重算。

两种写法的结果应当一致，但逐个插入的写法在分裂时特别容易漏掉「把已有那个粒子的贡献搬到新子节点」这一步，症状是父节点的 mass 对、质心偏，越靠近树顶偏得越离谱，而粒子少的地方又看不出来。所以我更倾向分治写法：质心由定义直接算出来，不依赖任何中间状态的正确性。

代价是每个节点都要对它那一批粒子做一次加权求和，整体是 O(N log N) 而不是 O(N)。这个代价在 N=10⁵ 时大约是建树总时间的 30%~40%，可以接受；真要抠，就在递归时把子节点已经算好的 mass 和 mx/my 传上来合并，回到 O(N)。

2.4 三条建完就该跑的自检断言

树建错的典型表现不是崩溃，是「跑得动但结果不对」。所以在进入力的计算之前，加三条断言把它摁住：

def check_tree(root, pos, mass, tol=1e-10): total_m = root.mass assert abs(total_m - mass.sum()) < tol * mass.sum(), "总质量不守恒" cx = float((mass * pos[:, 0]).sum()) / mass.sum() cy = float((mass * pos[:, 1]).sum()) / mass.sum() assert abs(root.mx - cx) < tol and abs(root.my - cy) < tol, "根节点质心不对" # 总质量为零的空节点不该被当成有效源 assert total_m > 0.0, "全质量为零"

第一条查的是质量有没有在分裂过程中丢失或重复计入。第二条查的是「质心 vs 几何中心」这个最容易写错的地方——如果mx恰好等于root.cx，八成是把几何中心赋过去了。第三条是为后面力计算的分母兜底：mass 为零的节点在近似公式里会给出 0 加速度，看起来无害，但如果它同时又是某个粒子的唯一父节点，说明这颗粒子根本没进树。

3. 开角判据与受力遍历：一次递归同时买到精度和速度

3.1 s/d < θ 到底近似掉了什么

对一个距离为 d、边长为 s 的节点，Barnes-Hut 用它的质心代替内部所有质量。这个近似展开到单极项就截止，丢掉的是偶极项（质心定义保证它为零）和四极项及更高阶。四极项的相对大小大约是 (s/d)² 量级，所以判据写成 s/d < θ 时，力的相对误差大致按 θ² 走。

这就解释了两件事。第一，θ 从 0.5 放宽到 1.0，误差不是翻倍而是涨到四倍左右，实测里经常看到「θ 调到 0.8 就明显不对了」。第二，θ 趋近 0 时这个近似退化成精确配对，算法复杂度回到 O(N²)——所以 θ 不是「越小越好」的精度旋钮，它是精度和速度之间唯一的那根杠杆。

3.2 单粒子受力遍历的完整实现

遍历从根开始。对当前节点：如果它是叶子，就对叶内每个粒子做精确配对；如果它是内部节点且满足开角判据，就用它的 mass 和质心算一次近似力然后返回；否则递归进四个子节点。整个函数是尾递归形式，深度受max_depth约束。

def accelerations(root, pos, mass, G=1.0, theta=0.5, eps=0.01): acc = np.zeros_like(pos) for i in range(len(pos)): _walk(root, i, pos, mass, G, theta, eps, acc) return acc def _walk(node, i, pos, mass, G, theta, eps, acc): if node.mass <= 0.0: return dx = node.mx - pos[i, 0] dy = node.my - pos[i, 1] # 叶子节点：逐个粒子精确计算，不做近似 if node.kids is None: for j in node.body: if j == i: continue ex = pos[j, 0] - pos[i, 0] ey = pos[j, 1] - pos[i, 1] s2 = ex * ex + ey * ey + eps * eps # 软化因子在这里进分母 inv = G * mass[j] / (s2 * np.sqrt(s2)) acc[i, 0] += ex * inv acc[i, 1] += ey * inv return d2 = dx * dx + dy * dy d = np.sqrt(d2) # 内部节点：s = 2 * half 是节点边长，s/d < theta 时整体近似 if d > 0.0 and (2.0 * node.half) / d < theta: s2 = d2 + eps * eps inv = G * node.mass / (s2 * np.sqrt(s2)) acc[i, 0] += dx * inv acc[i, 1] += dy * inv return for child in node.kids: _walk(child, i, pos, mass, G, theta, eps, acc)

逐段看逻辑。开头的node.mass <= 0.0是空节点早退，用<=而不是==是为了顺手挡掉负质量这种不可能但一旦出现就会静默污染全局的输入。叶子分支里j == i的跳过不能省：叶节点持有的是下标数组，自身一定在里面，不跳过就是自相互作用，会得到一个巨大的虚假加速度。近似分支里的d > 0.0是防除零，d恰好为零意味着粒子和节点质心重合，这在实际演化里出现概率极低但一旦出现就是 NaN 扩散。

近似公式acc += G * M * Δr / (d² + ε²)^{3/2}是标准的 Plummer 软化形式。注意 ε 是加在 d² 上而不是加在 d 上，也不是加在最终的力上——加的位置不同，近距离行为完全不同，这一点在第 5 章会展开。

3.3 θ 的取值：一张可以直接抄的对照表

表里的误差是量级估计，不是可以直接引用的数字——它和粒子分布、软化因子、节点形状都有关系。稳妥做法是：拿自己的初始条件，用直接求和（下一节的脚本）算一遍精确加速度，再扫一遍 θ，画一条「θ — 相对误差」曲线，选误差落在你可接受范围内、且访问节点数已经明显下降的那个拐点。我一般从 0.5 起步，因为这个值在多数分布下误差已经压到个位数百分比，而节点访问量相比精确配对少了两个数量级。

还有一点：θ 是全局常数还是逐节点变化，对结果影响不小。有些实现会根据节点里粒子数动态放宽 θ（粒子多、距离远就更敢近似），这在密度反差极大的分布里有效，但它会让误差不再有统一量纲，做精度标定时很难说清是哪个 θ 在起作用。前期建议先固定全局 θ，把其他环节调稳再考虑动态策略。

3.4 软化因子 ε 不是精度参数，是物理参数

ε 的作用是把 1/r² 的奇点抹平。在真实的星团里，两颗星靠得足够近时会有非引力相互作用或者直接碰撞，软化就是在粗粒化掉这个过程。所以 ε 的选择应该跟你的物理问题挂钩：模拟的是稀疏星团，ε 取平均间距的 1/10 到 1/100 就够；模拟的是致密系统，ε 太小会让近距离交会主导整个演化，时间步被迫压得很小。

经验做法是先估一下初始分布的平均最近邻距离，把 ε 设在它的 1/10 附近，然后跑短轨道看能量曲线的斜率。如果能量快速漂移，先怀疑时间步，再怀疑 ε；如果双星轨道被异常散射，那基本是 ε 太小。

4. 时间积分：把加速度接成能跑十万步的轨道

4.1 为什么显式 Euler 在 N 体里必然翻车

力的计算再准，积分器选错也白搭。显式 Euler 每一步用当前加速度同时更新速度和位置，局部误差 O(dt²)，而且是系统性地往系统里注入能量——在引力这种长程、近距离交会频繁的势场里，几万步之后能量偏移足以让整个系统「热膨胀」到面目全非。这不是精度不够的问题，是方法本身没有守恒结构。

蛙跳积分（leapfrog，KDK 形式）属于辛积分器，它精确保持相空间体积，能量误差不在时间上单调增长，而是在一个有界区间里振荡。这个性质对 N 体模拟是决定性的：你不需要很小的时间步就能让总能量在十万步后仍然待在同一条带子里。代价是它只有二阶精度，补偿手段是「半步 kick、整步 drift、再半步 kick」。

4.2 KDK 三步和树的重建时机

def leapfrog_step(pos, vel, mass, dt, acc, cfg): vel += 0.5 * dt * acc # 半步 kick pos += dt * vel # 整步 drift root = build_tree(pos, mass, cfg["leaf_cap"], cfg["max_depth"]) acc = accelerations(root, pos, mass, G=cfg["G"], theta=cfg["theta"], eps=cfg["eps"]) vel += 0.5 * dt * acc # 再半步 kick return pos, vel, acc

顺序不能换。第一步用旧位置的加速度把速度推半步，第二步用这个中间速度把位置推一整步，第三步在新位置重算加速度再补半步速度。这样速度和位置在时间上错开半步，正好构成一个辛映射。

注意build_tree和accelerations都必须在pos更新之后调用——树是位置的函数，位置变了树就作废了。这一步是整段代码里最容易被「优化」错的地方：有人为了省时间在多个步长内复用同一棵树，结果粒子已经移出原来的节点范围，近似判据算出来的 s/d 完全是错的，误差不会报错，只会让轨道莫名其妙地发散。

acc在步与步之间传递，省掉一次重复计算。第一次调用前需要单独算一次acc，从那之后每一次的返回值都是下一步要用的旧加速度。

4.3 能量漂移是唯一可信的黑匣子读数

N 体模拟出问题时，画面往往还「看着像那么回事」，只有能量曲线会说实话。总能量 E = KE + PE，KE 是 Σ½mv²，PE 是所有粒子对的 -Gmᵢmⱼ/rᵢⱼ（带同样的软化）。

def energy(pos, vel, mass, G=1.0, eps=0.01): ke = 0.5 * np.sum(mass * np.sum(vel * vel, axis=1)) # 势能用直接求和，O(N²)，只在 N 较小或做基准时用 dx = pos[:, None, 0] - pos[None, :, 0] dy = pos[:, None, 1] - pos[None, :, 1] r2 = dx * dx + dy * dy + eps * eps np.fill_diagonal(r2, np.inf) # 排除自相互作用 pe = -0.5 * G * np.sum(mass[:, None] * mass[None, :] / np.sqrt(r2)) return ke + pe

逻辑上有两个点要抠。np.fill_diagonal(r2, np.inf)把自配对的距离置成无穷大，对应的贡献是 0，同时避免了对角线上的除零。前面的 0.5 是因为mass[:,None] * mass[None,:]这张 N×N 矩阵把每一对算了两遍。

软化因子必须和力的计算里用的完全一致。如果力里用(d²+ε²)^{3/2}、势能里用1/d，能量曲线会呈现出「看起来在漂移」的假象，实际上是两套势场在打架，这种不一致比数值误差更难查，因为它不随步长减小而收敛。

判断标准很简单：把 dt 减半再跑一遍，如果能量漂移的形状和幅度都明显改变，说明是时间步问题；如果几乎不变，说明误差来自树的近似（θ）或者积分器的守恒性质，跟 dt 无关。这一步判断能省掉大量盲目调参。

4.4 自适应时间步：让近距离交会自己决定节奏

固定 dt 的问题在于它必须迁就最极端的那个事件。一套星团里大部分粒子运动平缓，个别粒子偶尔会有近距离交会，用统一的小步长跑全程是巨大的浪费。

常见做法是 Aarseth 判据的变体：每个粒子按自己的加速度算一个特征时间 τᵢ = η·sqrt(ε/|aᵢ|)，取全局最小值作为这一步的 dt。η 一般取 0.01~0.05，它控制的是「一步之内速度变化占特征速度的比例」。这么做能自动在系统剧烈变化时收缩步长，平静时放大。

代价是 dt 变成时变的，严格的辛性质会打折扣——不过只要 dt 变化缓慢，能量误差仍然能保持在有界振荡而不是单调漂移，这一点在实践中比理论上的严格性更重要。如果观察到能量曲线变成阶梯状上跳，多半是 dt 变化太剧烈，把 η 调小或者对 dt 做平滑（比如新 dt 不超过旧 dt 的 1.2 倍）就能压住。

5. 排查手册：θ、ε 和退化树导致的五类典型故障

5.1 力算错的三类症状，以及该先看哪里

Barnes-Hut 的 bug 有一个共同特点：不崩溃。程序照跑，图照出，只有和基准对不上。所以排查顺序要固定下来，别东一榔头西一棒子。

第一顺位是能量漂移的形状。单调上升或下降，查积分器和 dt；有界振荡但幅值偏大，查 θ 和 ε；出现台阶状突变，查自适应步长。

第二顺位是加速度的逐粒子对比。拿 N=200 的小规模，用直接求和算出精确加速度，和树方法逐分量相减，看误差是均匀分布还是集中在某些粒子上。均匀分布是 θ 的问题，集中在少数粒子（通常是靠近节点边界的）是判据实现的问题，集中在质量最大的那几个粒子上是质心算错了。

第三顺位才是看轨迹图。因为轨迹是最后才被误差污染的东西，从它反推原因效率很低。

5.2 五条翻车记录

一、跑几十步后位置全是 NaN

现象：pos里出现nan，能量变成-inf，日志里没有任何异常。

原因：初始条件里有两个粒子坐标完全相同（或者近到双精度在这个尺度上分辨不出）。象限判据用的是> 0，两个点落在同一侧，永远分不开，递归一路到底直到栈溢出；即使有深度上限兜住，到了力的计算阶段 d≈0 又会除零。

解决：建树时加max_depth兜底（本文实现里已有），力计算里所有分母都过一遍软化因子 ε，初始化生成粒子时对重复坐标加 1e-10 量级的随机扰动。三条里最不该省的是第一条，因为它是唯一能防住栈溢出的。

二、同一个 θ，别人的代码误差 1%，我的误差 8%

现象：对照直接求和，相对误差差了好几倍，但代码逻辑逐行比对看不出区别。

原因：θ 的比较方式没有统一约定。常见写法有三种：s/d < θ、s²/d² < θ²（等价）、以及用节点对角线代替边长（差一个 √2）。更隐蔽的是 s 取边长还是半边长——差一个 2 的因子，等效于 θ 差了 2 倍，误差按 θ² 走就是 4 倍差距。

解决：在代码里把约定写死并注释清楚（本文用的是「s = 2 × half 即边长，判据 s/d < θ」），然后用同一份初始条件、同一套直接求和基准，扫一遍 θ 画出自己的误差曲线。别拿别人的 θ 数字直接用。

三、总质量守恒，但根节点的质心慢慢偏离

现象：root.mass一直对，root.mx和全体粒子的加权平均对不上，偏差随时间增大。

原因：用了节点的几何中心 (cx, cy) 代替质心做远场近似。几何中心在节点内质量分布对称时恰好等于质心，所以小规模测试很容易蒙混过关；一旦分布不对称，偏差立刻显现，而且它表现为力的系统性偏移，不会报错。

解决：每个节点的 mx/my 必须由(mass * pos).sum() / mass.sum()算出来，并且加一条断言：根节点的 mx 与全体粒子的加权平均之差应当在 1e-10 量级。如果这两个数恰好相等而root.cx也相等，反而要警惕——那说明初始分布太对称，没测出问题。

四、双星靠得很近时被弹飞，能量突然跳一截

现象：两条轨道接近到最小距离时，其中一个粒子的速度异常增大，之后整团天体的能量曲线抬升一个台阶。

原因：软化因子 ε 只加在了力的大小上，势能里用的还是1/r；或者力和势能虽然都加了 ε，但形式不同（一个用(d²+ε²)^{3/2}，一个用d²+ε²）。两套势场不一致，在近距离时表现出的等效能量不守恒。

解决：力和势能统一使用 Plummer 形式——力GmΔr/(d²+ε²)^{3/2}，势能-Gm/sqrt(d²+ε²)。改完之后能量曲线的台阶应当消失，剩下的平滑漂移才是数值误差。

五、总能量单调上升，系统整体「变热」

现象：几万步之后，原本的团块逐渐散开，速度弥散变大，能量曲线一路向上。

原因：多半是 dt 太大遇上近距离交会。也有可能是 θ 太松，导致树近似在每步引入的误差方向不是随机的，而是有偏的（近似总是低估远场的力，长期累积就成了净能量注入）。

解决：先把 θ 从 0.5 收到 0.3 重跑，如果漂移速度明显下降，问题在树近似；如果几乎没变，把 dt 减半再跑，漂移幅度应随之减半（二阶方法）。两个都试过还漂，检查自适应步长是否变化太剧烈。

5.3 一个三十行的交叉验证脚本

上面所有排查都指向同一件事：你得有一个「随时能拿来对答案」的精确解。直接求和写起来很短，值得常备。

def direct_accelerations(pos, mass, G=1.0, eps=0.01): dx = pos[:, None, 0] - pos[None, :, 0] # dx[i, j] = x_i - x_j dy = pos[:, None, 1] - pos[None, :, 1] r2 = dx * dx + dy * dy + eps * eps np.fill_diagonal(r2, np.inf) # 排除自配对 inv = G / (r2 * np.sqrt(r2)) # 1/(r²+ε²)^{3/2} ax = -np.sum(mass[None, :] * dx * inv, axis=1) ay = -np.sum(mass[None, :] * dy * inv, axis=1) return np.stack([ax, ay], axis=1) def compare(pos, mass, root, theta, eps): a_tree = accelerations(root, pos, mass, theta=theta, eps=eps) a_ref = direct_accelerations(pos, mass, eps=eps) rel = np.linalg.norm(a_tree - a_ref) / np.linalg.norm(a_ref) worst = np.max(np.linalg.norm(a_tree - a_ref, axis=1) / np.linalg.norm(a_ref, axis=1)) return rel, worst

符号要留意：dx[i, j] = x_i - x_j，而作用在 i 上的加速度指向 j，所以整体带一个负号。这个负号和前面树遍历里的写法是等价的（那边写的是pos[j] - pos[i]，不带负号），两种约定都对，但混用会让误差看起来像 200% 而不是 2%。

worst这个指标比全局相对误差更有用。全局误差被大量平缓的远场粒子摊薄了，看起来总是很漂亮；worst 会告诉你最差的那个粒子错了多少，而轨道演化恰恰是被最差的那几个粒子带偏的。我一般要求 θ=0.5 时 worst 控制在 5% 以内。

6. 进阶：给质心项补上四极矩，以及我踩过的那个坑

单极近似丢掉的是四极项。如果你想在同样的 θ 下把误差压下来，最划算的改法不是调小 θ（那是指数级的代价），而是给每个节点额外存两个二阶矩分量，力公式里补一项：

四极修正的形状是：把G·M·Δr/(d²+ε²)^{3/2}换成G·M·Δr/(d²+ε²)^{3/2} + 修正项(Δr, Q)。修正项里包含Δr的三次方组合和 Q 的分量，具体形式按二维展开写开有六项，实践中多数人直接调用现成的多极展开库或者查表实现，没必要手推。

我更想讲的是另一样东西：在 Python 里用对象指针建树，N 超过 5×10⁴ 之后，耗时的大头不是浮动运算而是属性访问和缓存未命中。那次翻车是把一份能跑 N=10⁴ 的代码直接拿去跑 10⁵，以为复杂度是 O(N log N) 所以只慢十倍，结果慢了将近四十倍，而且 profiler 显示 70% 的时间花在_walk的属性读取上。

改法是把树压成几个并行的 NumPy 数组：mass[]、mx[]、my[]、half[]、first_child[]、child_count[]

本文还有配套的精品资源，点击获取