简介：这是一份面向天体物理模拟、数值计算与C++算法学习者的Barnes-Hut N体问题模拟器源码包，重点展示如何用八叉树对三维空间分层分区，并以远近阈值近似计算远距离粒子引力，从而降低大规模粒子系统的逐对计算开销。压缩包共60个文件，以22个h头文件和21个cpp源文件为主，辅以4个mk构建文件、4个xml工程配置及2个md说明文档，整体约54KB；目录覆盖粒子模型、BHTree、引力计算、多种数值积分器、SDL窗口和主程序入口，能较完整地呈现从算法设计到工程实现的对应关系。已有328人学习，适合通过源码理解八叉树构建、近似精度控制、时间步进与多线程性能优化，并观察不同积分器的切换方式及可视化窗口的组织结构，也可作为课程设计、天体仿真或N体问题研究的参考实现。

1. 一万个质点，每一步两两配对就是五千万次开方

用暴力枚举处理 N 体问题，最让人头疼的不是代码写不出来，而是跑不完：一万个质点、每一步做全对枚举就是 10^8 次配对，跑五百步下来，机器基本就只能干等着。Barnes-Hut 这个项目解决的正是这个场景：把远处的星团整体近似成一个质心，用 θ 判据决定「这个格子能不能直接当成一个点」，复杂度从 O(N²) 压到 O(N log N)。我用 C++ 把它的内核拆了一遍，代码库很干净，没有渲染器也没有花哨的 UI，核心就是空间树构建、受力遍历、辛普积分器和快照输出这几块拼起来。它适合两类人：一是想给天体演化 demo 找一个能跑的计算后端，二是想拿真实场景练手空间数据结构、递归剪枝和缓存友好性的 C++ 开发者。下面按「能跑起来 → 看懂树 → 看懂力 → 少踩坑 → 动手改」的顺序讲。

2. 工程骨架与跑通：一份 C++ 目录怎么分层

拿到这种偏算法的项目，很多人第一反应是从 main 开始读，读两百行就迷路了。我一般的顺序是反的：先看目录怎么切，再看参数怎么进来，最后才进算法。原因是 N 体模拟的代码有个特点——树、力、积分三块耦合很紧，但边界其实很清楚，一旦目录分对了，后面读哪一块都不会被另一块干扰。

2.1 目录结构：树、力、积分、IO 各自独立

这个项目的目录大致是这么分的：

分层的原则很朴素：quadtree只依赖body.h，对外暴露一个只读的nodes()接口；force只读树、不改树；integrator负责把三块串成一个时间步。这样一来，你想把四叉树换成八叉树、把欧拉积分换成 RK4、把重力换成别的中心力场，改动都能限制在一个文件里。血泪经验是：早期图省事把受力遍历直接写进main的循环体里，等到想加个「多线程版本对比单线程」的功能时，只能整段复制粘贴，改一处忘一处。

这个仓库的核心是二维四叉树，也就是每个内部节点四个孩子。三维场景下换成八叉树，孩子数从 4 变 8，象限判定多一个 z 轴位，其余逻辑（插入、质心回填、θ 判据）一行都不用改。所以下面讲的所有东西都是可平移的。

2.2 参数怎么传：一个 Params 结构体顶掉一堆全局变量

模拟类程序最容易长出来的坏味道就是全局变量满天飞——double theta = 0.5;写在头文件里，改一次要重新编译一次。这个项目用的是单一Params结构体 + 手写命令行解析，不引第三方库：

// include/params.h #pragma once #include <string> struct Params { int n = 2000; // 粒子数 int steps = 500; // 总时间步数 double dt = 0.005; // 时间步长 double theta = 0.5; // Barnes-Hut 近似判据，越小越准越慢 double eps = 0.05; // 引力软化长度，防止近距离除零 unsigned seed = 42; // 初始条件随机种子，保证可复现 int outEvery = 10; // 每多少步写一帧快照 std::string init = "plummer";// 初始条件类型：plummer / disk / two-body std::string out = "run.csv";// 快照输出路径 };

解析部分就是最土的argv循环，识别不到的参数直接打一行用法提示然后退出：

// src/main.cpp（片段） Params parseArgs(int argc, char** argv) { Params p; for (int i = 1; i < argc; ++i) { std::string k = argv[i]; auto next = [&]() -> std::string { if (i + 1 >= argc) { std::fprintf(stderr, "%s 缺少取值\n", k.c_str()); std::exit(1); } return argv[++i]; }; if (k == "--n") p.n = std::stoi(next()); else if (k == "--steps") p.steps = std::stoi(next()); else if (k == "--dt") p.dt = std::stod(next()); else if (k == "--theta") p.theta = std::stod(next()); else if (k == "--eps") p.eps = std::stod(next()); else if (k == "--seed") p.seed = (unsigned)std::stoul(next()); else if (k == "--out-every")p.outEvery = std::stoi(next()); else if (k == "--init") p.init = next(); else if (k == "--out") p.out = next(); else { std::fprintf(stderr, "未知参数: %s\n", k.c_str()); std::exit(1); } } return p; }

几个参数的实际含义值得说清楚，因为它们直接决定你调出来的结果对不对：

• theta是近似强度的总开关。判据是s / d < theta，其中s是格子边长、d是目标粒子到格子质心的距离。θ 越小越保守，θ=0 时退化成暴力枚举，只适合做基准对比。
• eps是软化长度，不是可有可无的装饰。真实引力在 r→0 时发散，数值上离散步子一定会把两个粒子推到极近处，不加软化就会得到inf或者nan，然后整盘散架。它参与计算的方式是把1/r²换成r/(r²+ε²)^(3/2)。
• dt要和最密区域的轨道周期匹配。步长太大会让辛普积分器的能量曲线一路往上走，看起来像「星团自己炸了」，其实是数值注能。
• seed固定住，做基准测试时才有可比性。同一组参数跑两次结果不一样，后面调 θ 的误差曲线就没法看了。

2.3 编译与第一次运行

构建用 CMake，好处是换机器不用手打一长串源文件名：

# CMakeLists.txt cmake_minimum_required(VERSION 3.16) project(barnes_hut CXX) set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_STANDARD_REQUIRED ON) if(NOT CMAKE_BUILD_TYPE) set(CMAKE_BUILD_TYPE Release) # 默认 Release，别用 Debug 跑性能测试 endif() add_executable(bh src/main.cpp src/quadtree.cpp src/force.cpp src/integrator.cpp src/io.cpp) target_include_directories(bh PRIVATE include) target_compile_options(bh PRIVATE -O3 -Wall -Wextra)

如果你想跳过 CMake 直接快跑，一条命令也行：

g++ -O3 -std=c++17 -Iinclude src/*.cpp -o bh

编译完的第一跑，建议先来个小规模的，确认整条链路通：

./bh --n 2000 --steps 500 --theta 0.5 --dt 0.005 --init plummer --out run.csv

正常的话终端会打出类似这样的东西：

[build] N=2000 depth=11 nodes=11873 build=0.42ms [step 0] E=-0.473120 (drift +0.00%) step=0.83ms [step 500] E=-0.472788 (drift +0.07%) step=0.79ms wrote run.csv (51 frames, every 10 steps)

这里有三件事要盯住：

第一行是树的规模。depth是树的实际深度，nodes是节点池里分配了多少个节点。如果depth跑到二三十甚至程序直接崩，八成是初始条件里有粒子坐标几乎重合，插入逻辑一直在往下钻，这就是后面避坑章要讲的第一条。

第二行是总能量。E是动能加势能，理想情况下应该基本持平，500 步漂移 0.07% 属于健康范围。如果它一路单调往上爬，先别怀疑算法，把--dt砍一半再看。

第三行是每步耗时。这个数字是你后面调 θ 的基准线，换个参数再跑，和它比才有意义。

run.csv的格式很直白，每帧若干行：

frame,id,x,y,vx,vy,m 0,0,1.234567,-0.456789,0.0123,-0.0456,1.0

拿它去喂任何绘图库都行，不用关心内部结构，这也是我习惯的输出方式——模拟器只管算，画图的事交给脚本。

3. 八叉树构建：Barnes-Hut 的地基

整棵树就干一件事：把「距离远的粒子簇」打包成一个带质量和质心的整体。构建分三步走——定根节点包围盒、逐个插入、自底向上回填质量与质心。这三步里最容易写错的是第一步之后的插入，因为它是递归的，而且带着一个隐藏的内存陷阱。

3.1 节点结构与「用下标不用指针」

先看节点的定义，这里每一个字段都有存在的理由：

// include/quadtree.h struct QuadNode { Vec2 center; // 格子几何中心，注意不是质心 double half = 0.0; // 半边长，格子边长 size = 2 * half double mass = 0.0; // 子树内所有粒子的质量之和 Vec2 com; // 质量加权质心，受力计算只认这个 int body = -1; // 叶子持有的粒子下标；-1 表示空叶子或内部节点 int child[4] = {-1, -1, -1, -1}; // 四个孩子的节点下标 bool leaf = true; };

center和com是两个绝对不能混的东西。center是格子的几何中心，由父节点推出来，用来判断粒子落在哪个象限；com是这次近似实际使用的「代表点」，是子树里所有粒子按质量加权出来的平均位置。用center代替com去算力，误差会大到肉眼可见——尤其是质量分布不均的时候，星团的力心会明显偏移。

child存的是节点下标而不是指针，这是刻意的。节点池是std::vector<QuadNode>，插入过程中会不断 push_back，一旦扩容，所有已经存下来的裸指针全部变成野指针，N小的时候碰巧不崩，N一大就随机 segfault 或者算出一堆离谱的质心。用下标就没有这个问题，代价只是多一层寻址。

节点池整体是这么组织的：

class QuadTree { public: void build(const std::vector<Body>& bodies, double padding = 1e-6); const std::vector<QuadNode>& nodes() const { return pool_; } // 只读，供受力遍历用 int root() const { return 0; } private: int makeNode(Vec2 c, double half); void insert(int node, int bodyIdx, const std::vector<Body>& bodies); int quadrant(const QuadNode& n, Vec2 p) const; void computeMass(int node, const std::vector<Body>& bodies); std::vector<QuadNode> pool_; };

nodes()返回 const 引用，是为了让受力遍历那一步明确知道「这棵树我只读不写」。如果后面要加多线程，这个 const 就是最有力的约定——编译器会帮你挡住所有意外写入。

3.2 插入与象限判定

象限判定是纯几何，四个孩子按二进制位编码：最低位管 x，次低位管 y。

// src/quadtree.cpp int QuadTree::quadrant(const QuadNode& n, Vec2 p) const { int q = 0; if (p.x >= n.center.x) q |= 1; // 右半区 if (p.y >= n.center.y) q |= 2; // 上半区 return q; // 0..3 对应 左下/右下/左上/右上 }

插入逻辑是整棵树里最绕的一段。核心规则：叶子节点空着就直接放；叶子已经有主了，就分裂成四个孩子，把老粒子推下去，再放新粒子。

void QuadTree::insert(int node, int bodyIdx, const std::vector<Body>& bodies) { // 关键：整个函数里不要在 push_back 之后继续持有 pool_[node] 的引用 if (pool_[node].leaf) { if (pool_[node].body == -1) { // 空叶子，直接落座 pool_[node].body = bodyIdx; return; } // 已有粒子，必须先分裂再下推 int old = pool_[node].body; pool_[node].body = -1; pool_[node].leaf = false; Vec2 c = pool_[node].center; double h = pool_[node].half * 0.5; // 孩子格子的半边长为父的一半 for (int k = 0; k < 4; ++k) { Vec2 cc{ c.x + ((k & 1) ? h : -h), c.y + ((k & 2) ? h : -h) }; pool_[node].child[k] = makeNode(cc, h); // 这里会 push_back，可能触发扩容 } // 分裂完再重新取一次下标，别复用上面的引用 insert(pool_[node].child[quadrant(pool_[node], bodies[old].pos)], old, bodies); } insert(pool_[node].child[quadrant(pool_[node], bodies[bodyIdx].pos)], bodyIdx, bodies); }

逻辑上分两种情况：空叶子直接落座，返回；非空叶子先把自己变成内部节点，四个孩子一次性创建好，然后把老粒子递归下推到它该去的象限，接着把新粒子也下推一层，继续走同一个函数。

参数上要注意h = half * 0.5。孩子节点的中心是父中心在对应方向上偏移半个孩子边长，而孩子边长是父边长的一半，所以偏移量正好是parent.half * 0.5。这里如果写成parent.half，格子会重叠，粒子被分到错误的象限，树的结构还「看起来正常」，但受力结果是错的——这种 bug 极难靠肉眼发现，只能靠和暴力枚举对拍。

还有一处新手特别容易中招：pool_是 vector，makeNode里的push_back可能触发扩容，导致之前拿到的QuadNode&引用全部失效。所以我在这段代码里全程用pool_[node]取下标访问，不缓存引用。这条规则我在后面避坑章还会再提一次，因为它真的会让人查一晚上。

根节点的包围盒根据所有粒子的包围盒算出来，稍微加一点 padding，避免粒子正好落在边界上：

void QuadTree::build(const std::vector<Body>& bodies, double padding) { pool_.clear(); pool_.reserve(bodies.size() * 2); // 预留，减少扩容次数 Vec2 lo{ 1e300, 1e300}, hi{ -1e300, -1e300 }; for (const auto& b : bodies) { lo.x = std::min(lo.x, b.pos.x); lo.y = std::min(lo.y, b.pos.y); hi.x = std::max(hi.x, b.pos.x); hi.y = std::max(hi.y, b.pos.y); } Vec2 c{ (lo.x + hi.x) * 0.5, (lo.y + hi.y) * 0.5 }; double half = 0.5 * std::max(hi.x - lo.x, hi.y - lo.y) * (1.0 + padding) + 1e-9; makeNode(c, half); // 根节点固定占 pool_[0] for (int i = 0; i < (int)bodies.size(); ++i) insert(0, i, bodies); computeMass(0, bodies); }

reserve(bodies.size() * 2)是个经验值。均匀分布下节点数大概在粒子数的 1.5 倍上下，稀疏区多的时候能到 2 倍多。预留少了无非是多几次扩容，预留多了浪费几十 KB 内存，对这个量级的程序来说都不是问题。

3.3 自底向上回填质心与总质量

树建好之后，每个内部节点的mass和com还是空的，得做一次后序遍历填上：

void QuadTree::computeMass(int node, const std::vector<Body>& bodies) { QuadNode& n = pool_[node]; if (n.leaf) { if (n.body < 0) { n.mass = 0.0; n.com = n.center; return; } // 空叶子 n.mass = bodies[n.body].mass; n.com = bodies[n.body].pos; return; } double m = 0.0; Vec2 c{ 0.0, 0.0 }; for (int k = 0; k < 4; ++k) { computeMass(n.child[k], bodies); const QuadNode& ch = pool_[n.child[k]]; m += ch.mass; c = c + ch.com * ch.mass; // 质量加权的分子部分 } n.mass = m; n.com = (m > 0.0) ? c * (1.0 / m) : n.center; // 除以总质量得到质心 }

这段是标准的两遍式：先递归处理四个孩子，拿到它们的质量和质心，累加出加权和，最后除一次总质量。注意m > 0的判断不能省——如果四个孩子都是空叶子，总质量是 0，直接除会得到处都是 nan，然后顺着父节点一路污染到根。空节点用几何中心兜底就够了，反正它质量为零，在受力遍历里会被直接跳过。

另一个细节是QuadNode& n = pool_[node]在这里是安全的，因为computeMass全程只读不写节点池，不会有 push_back。但如果你后面把「统计每层节点数」之类的逻辑塞进这个函数里并且顺手 resize 了池子，那就又踩回引用失效的坑了。

回填完，根节点上的mass就是全场总质量，com就是系统质心。多数情况下这两个值在积分过程中应该基本不动（动量守恒），跑完之后顺手打出来对一下，是一个几乎零成本的正确性检查。

4. 受力遍历：θ 判据怎么决定近似还是继续下钻

树建好只是准备，真正决定精度和速度的是受力遍历那几十行。它的思路很直接：从根开始，对每个目标粒子，看当前格子「够不够远、够不够小」，如果看起来像一个小点，就整体算一次力；如果觉得太粗糙，就往下展开四个孩子。判据就是 θ。

4.1 单粒子递归受力与自交互剔除

先看完整的递归函数：

// src/force.cpp static constexpr double G = 1.0; // 归一化引力常数，只影响时间尺度 void accumulate(int node, int selfIdx, Vec2 pos, const std::vector<QuadNode>& pool, double theta, double eps2, Vec2& acc) { const QuadNode& n = pool[node]; if (n.mass <= 0.0) return; // 空子树直接剪掉 Vec2 d = n.com - pos; // 目标粒子指向格子质心 double r2 = norm2(d); if (n.leaf) { if (n.body < 0 || n.body == selfIdx) return; // 空叶子 / 自己，跳过 double inv = 1.0 / std::sqrt(r2 + eps2); acc = acc + d * (G * n.mass * inv * inv * inv); // a = G*m*d/|r|^3 return; } double size = 2.0 * n.half; if (size * size < theta * theta * r2) { // s/d < theta 的平方形式，省一次开方 double inv = 1.0 / std::sqrt(r2 + eps2); acc = acc + d * (G * n.mass * inv * inv * inv); return; // 近似成功，整块打包 } for (int k = 0; k < 4; ++k) // 不够远，继续下钻 accumulate(n.child[k], selfIdx, pos, pool, theta, eps2, acc); }

逻辑分四支：

第一支，mass <= 0的空子树直接返回，这是最便宜的剪枝。

第二支，叶子节点。空叶子返回；如果是目标粒子自己所在的叶子，也必须返回，否则d是零向量，r2是 0，虽然eps2能兜住除零，但会给自己凭空加一份朝任意方向的加速度，几十步下来轨迹就歪了。这个自交互剔除是必须写的。

第三支，内部节点做 θ 判据。注意条件写成了平方形式size² < θ² · r²，避免每次判断都开一次方。这是很典型的优化，递归里每层都判断一次，开方的开销会累加成可观的比例。

第四支，下钻四个孩子。

关于单次近似的公式，展开写是：

a = G * m_cell * (com - pos) / (|com - pos|² + ε²)^(3/2)

G在这里取 1，纯粹是归一化，只改变时间尺度不影响轨道形状。真实场景下把它换成实际常数即可，代价是dt要跟着重新标定。ε²加在平方项里面而不是外面，这样在r很小时加速度会趋于零而不是发散，物理上等价于把点质量抹成一个半径ε的软球。

调用侧的循环长这样：

void computeAccelerations(std::vector<Body>& bodies, const QuadTree& tree, double theta, double eps2) { for (int i = 0; i < (int)bodies.size(); ++i) { Vec2 a{ 0.0, 0.0 }; accumulate(tree.root(), i, bodies[i].pos, tree.nodes(), theta, eps2, a); bodies[i].acc = a; // 只写自己的 acc，不碰别人的，天然可并行 } }

注意这里每个线程只写bodies[i].acc这一个位置，互相不重叠。这个性质在二次开发里非常值钱——加速比能不能接近线性，很大程度上就取决于这个循环体没有被写坏。

4.2 引力软化与辛普积分器

算出加速度之后就该推进时间了。N 体模拟里最常见的翻车不是算错力，而是积分器选错。显式欧拉（x += v*dt; v += a*dt）在轨道问题上会持续注入能量，跑长了轨道会越绕越大，看着像星团被什么东西吹散了。

这个项目用的是「踢-漂-踢」形式的辛普积分器，也叫蛙跳：

// src/integrator.cpp void step(std::vector<Body>& bodies, QuadTree& tree, const Params& p, double& energy) { const double dt = p.dt, half = 0.5 * dt; // 1) 半步踢：用上一步缓存下来的加速度推半个时间步的速度 for (auto& b : bodies) b.vel = b.vel + b.acc * half; // 2) 整步漂：位置用更新后的速度推进一整步 for (auto& b : bodies) b.pos = b.pos + b.vel * dt; // 3) 重建树、重算加速度（位置变了，树和力都必须刷新） tree.build(bodies); computeAccelerations(bodies, tree, p.theta, p.eps * p.eps); // 4) 再半步踢：用新加速度补齐剩下的半个时间步 for (auto& b : bodies) b.vel = b.vel + b.acc * half; energy = totalEnergy(bodies, tree, p.eps * p.eps); }

三步的顺序不能乱。先踢后漂再踢，保证速度和位置在时间上是对齐到同一个步点；如果写成「漂-踢」，位置用的是旧速度，精度掉到一阶，能量守恒的性质也就没了。

energy那个参数是给主循环做监控用的。总能量的计算是动能求和加上对势能的近似求和——势能同样可以用树来算，把每个叶子的质量乘上质心处的势，误差量级和受力近似一致。这不是精确值，但它稳定，用来观察漂移趋势足够了。

重建树这一步不能省。粒子的位置每一步都在变，哪怕只挪了一点点，格子里装的东西也可能换了人。有人图省事每十步才重建一次树，跑小规模时看不出问题，星团一密集就会看到明显的「格子边缘拉扯」——原因就是旧树里的质心已经和新位置对不上了。

dt怎么定？我一般先用系统里最紧凑那对粒子的最近距离估一个特征轨道周期T ≈ 2π·sqrt(r³/(G·M))，取dt ≈ T/100。跑一遍看能量曲线，如果五百步内漂移超过 1%，就把dt再砍一半。这个流程比拍脑袋定步长靠谱得多。

4.3 θ 与精度的实测曲线

调 θ 之前，得先有一个「准」的参照。项目里留了一个暴力枚举实现，只用来对拍，不参与主循环：

// 仅用于验证：O(N^2) 全对枚举，N 控制在几千以内 void bruteForceAcc(const std::vector<Body>& bodies, double eps2, std::vector<Vec2>& acc) { acc.assign(bodies.size(), Vec2{}); for (size_t i = 0; i < bodies.size(); ++i) { for (size_t j = 0; j < bodies.size(); ++j) { if (i == j) continue; Vec2 d = bodies[j].pos - bodies[i].pos; double inv = 1.0 / std::sqrt(norm2(d) + eps2); acc[i] = acc[i] + d * (G * bodies[j].mass * inv * inv * inv); } } }

对拍的方法是：同一组初始条件，先用暴力算一遍加速度存下来作为基线，再用不同 θ 跑 Barnes-Hut，算两者加速度向量的相对 L2 误差。我这边量到的量级大致是这样（绝对值取决于机器和初始条件，看相对关系就行）：

几个结论可以直接拿去用：

θ 从 0.5 降到 0.3，误差只降了 1.4 个百分点，耗时却翻了将近一倍。这就是为什么绝大多数实现默认取 0.5 附近——它正好踩在误差曲线的拐点上。

θ 超过 0.8 之后误差涨得很快，0.8 时的 5.6% 已经会在长跑里表现为明显的轨道偏转。如果你的场景是「看个大概」，0.8 也能接受；如果要做定量分析，别超过 0.6。

大规模才是 Barnes-Hut 的主场。看 N 从 2000 涨到 100000（50 倍），单步耗时只从 1.1 ms 涨到 80 ms（约 70 倍）。比值接近N log N / N₀ log N₀，说明实现没有意外的常数项拖累。如果这个比值明显超过理论值，八成是树的深度失控或者缓存命中率崩了，去查初始条件里有没有重合粒子。

5. 避坑与排查：五个把我卡到半夜的问题

下面这几条都是我自己撞过、或者帮别人看代码时反复见到的。现象、原因、解决按顺序写，遇到了直接对号入座。

5.1 程序卡死或栈溢出：两个粒子的坐标几乎重合

现象：N上千之后偶发崩溃，或者 CPU 单核跑满、程序几分钟不动。用调试器中断进去看，调用栈深得吓人，全在insert里。

原因：插入逻辑只根据「叶子里有没有粒子」决定分裂，没有深度上限，也没有最小格子尺寸。一旦两个粒子的坐标差小于浮点精度，它们永远落在同一个象限里，递归就无限往下钻，直到爆栈。随机生成的初始条件里，这种重合其实比想象中常见——尤其是用整数坐标或者低精度采样的时候。

解决：给树加一个maxDepth（我一般取 24），到顶之后不再分裂，把粒子塞进叶子的一个小桶里。叶子从「持有一个粒子」改成「持有一个下标数组」，桶大小取 4 到 8，受力计算时逐个展开。这样既保证了递归深度有界，又不会因为丢弃粒子造成质量不守恒。

5.2 随机崩溃或质心算错：缓存了节点的引用

现象：N小的时候一切正常，N一大就随机段错误，或者算出来的质心明显偏离。加上-fsanitize=address跑一遍，报的是 heap-use-after-free。

原因：std::vector<QuadNode>在push_back触发扩容时会重新分配整块内存，之前拿到的QuadNode&或者QuadNode*全部悬空。写插入逻辑时很自然地会写成QuadNode& n = pool_[node];，然后在这个函数里调用了会创建孩子的makeNode，几行之后再访问n就是未定义行为。小N时预留的容量够用、不触发扩容，所以「看起来」是对的。

解决：两条路。一是像我上面那段代码一样，全程用pool_[node]下标访问，不缓存引用；二是提前reserve一个足够大的容量，但这只是降低概率不是消除，不推荐。我现在写这类节点池的代码，在函数开头先写一行注释提醒自己「此处不可持有引用」，看着傻，但真的省时间。

5.3 第一步之后全是 nan：自己和自己算了引力

现象：程序跑完第一步，所有粒子的坐标变成nan，之后的输出全是nan。

原因：目标粒子在叶子节点里遇到了自己，d是零向量，r2是 0。虽然公式里有eps2兜着，但如果eps忘了设或者被设成 0，inv = 1/0就出来了；即使eps正常，给自己加一份方向由零向量决定的力也是错的（0 * inf同样是 nan）。

解决：在叶子分支里显式判n.body == selfIdx并返回。另外eps一定要设成非零，哪怕很小，它是最后一道保险。我还见过内部节点也不剔除自交互的写法——那个更隐蔽，因为质心一般不等于粒子位置，看起来「有微小误差」，实际是系统性偏差。

5.4 θ 调小反而更慢：把近似关掉等于白写

现象：有人为了「更准一点」把 θ 设成 0.1 甚至 0.01，结果跑得比暴力枚举还慢好几倍。

原因：θ 是s/d的阈值。θ 越小，能通过判据的节点越少，几乎所有节点都要往下展开，遍历的节点数反而超过了N log N，再加上递归调用的开销，确实可能比裸的双重循环还慢。θ=0 时理论上等价于暴力枚举，但每一次访问都带函数调用和判据判断，常数项大得多。

解决：θ 的有效范围就是 0.3 到 1.0，默认 0.5。真要更高精度，正确做法是减小dt或者换用更高阶的多极展开，而不是把 θ 压到极限。

5.5 能量缓慢单调上升：时间步太大

现象：总能量曲线不是随机抖动，而是稳定地、单调地往上爬，几千步之后累计十几个百分点，画出来的轨迹明显发散。

原因：辛普积分器只在dt小于系统特征时间尺度时才有好的守恒性。星团里总有一两对粒子靠得很近，它们的轨道周期极短，如果全局dt按平均间距来定，这一对就会被严重欠采样。数值上表现为动能被系统性放大。

解决：先把dt减半验证，如果漂移量大致按dt²下降，就确认是积分器的问题。长期方案有两种：一是全局采用最紧致区域的T/100，代价是慢；二是做自适应步长，或者对近距粒子做专门的近距离处理。工程上大多数演示场景用第一种就够了。

6. 二次开发：把受力循环并行掉，再接一层可视化

这个项目跑通之后最容易做的两件事，一个是提速，一个是看见结果。顺序别反，先保证串行版本的结果是对的，再并行；否则你根本分不清是线程写错了还是索引写错了。

6.1 OpenMP 改造：只并行最外层受力循环

树在积分步内是只读的，所以受力遍历是天然的并行点。改动小到只需要加一行编译指令和一行 pragma：

// src/force.cpp void computeAccelerations(std::vector<Body>& bodies, const QuadTree& tree, double theta, double eps2) { const auto& pool = tree.nodes(); const int n = (int)bodies.size(); // 每个粒子只写自己的 acc，不存在写冲突；树是只读的，共享安全 #pragma omp parallel for schedule(dynamic, 64) for (int i = 0; i < n; ++i) { Vec2 a{ 0.0, 0.0 }; accumulate(tree.root(), i, bodies[i].pos, pool, theta, eps2, a); bodies[i].acc = a; } }

编译时加-fopenmp，CMake 里对应find_package(OpenMP REQUIRED)加target_link_libraries(bh PRIVATE OpenMP::OpenMP_CXX)。

几个参数值得说一下。schedule(dynamic, 64)里的 64 是分块大小，意思是每次分配给线程 64 个连续粒子。用 dynamic 而不是 static，是因为粒子受力遍历的

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