简介：面向大规模粒子系统模拟的C++实现Barnes-Hut-Simulator，基于Barnes-Hut算法高效求解N体引力问题，适合星系演化、恒星运动等天体物理场景，也便于计算物理爱好者学习空间划分与数值积分。源码采用八叉树递归分区，以远近法则和质心近似降低引力计算复杂度，并集成Euler、RK、ADB等多种积分器，兼顾精度与效率。包内共60个文件，主体为22个头文件与21个C++源文件，涵盖粒子类、八叉树、N体模型、窗口显示与积分器模块，另有Makefile与CMake构建脚本、NetBeans工程配置、Markdown说明及打包脚本，压缩包仅54KB，结构清晰，附README与INSTALL说明，便于快速上手。已有327人学习/下载。读者可对照源码掌握八叉树构建、近似引力判定、并行性能优化思路及可视化实现，并通过切换Euler、RK、ADB等积分器理解数值模拟的时间推进方式；从NBody模型到SDL窗口渲染的完整代码路径，适合作为C++算法学习、课程设计或科研入门的参考工程。
1. Barnes-Hut-Simulator：N体问题为什么越跑越慢，以及怎么破
做N体模拟的同事十有八九遇到过这种尴尬：粒子数从几千涨到几万，暴力计算的模拟器直接变幻灯片，一步要算几分钟。N体问题的瓶颈从来不是公式难，而是每一对粒子之间都要算一次引力，复杂度是O(N²)，规模翻一倍，耗时就是四倍。Barnes-Hut-Algorithm的解法是把空间切成一棵四叉树，远处的一群粒子用它们的质心代替，复杂度降到O(N log N)。Barnes-Hut-Simulator就是冲着这个目标来的。这篇笔记会用一套可复现的2D实现讲清楚建树、受力计算、参数标定和常见的坑，适合正在写粒子模拟、星系演化或者分子动力学Demo的读者。
2. 算法原理：质心近似为什么能省出三个数量级
2.1 暴力算法的复杂度到底亏在哪
先看最基本的暴力做法。N个粒子两两计算引力，总共算N(N-1)/2次，渐近复杂度O(N²)。当N=1000时约50万对；N=10000时就是接近5000万对。每一步都要重算一遍，因为粒子在运动。5000万次重力计算在单线程Python里已经能跑到秒级以上，而模拟通常要跑几千步，完全不可接受。
Barnes-Hut的核心洞察是：引力具有“远处可近似”的特性。站在很远的地方看一个星系团，它引力的细节几乎不取决于内部每个星体的具体位置，只取决于总质量和质心位置。所以不需要对远处的每个粒子都精确计算，只要对远处的一团粒子算一个“合起来”的引力即可。这个近似造成的误差可以控制，而计算量从N²级别降到N log N级别，这就是它“有效解决N体问题”的根本原因。
2.2 θ判据：多远算“远”
“多远算远”不能拍脑袋，Barnes-Hut算法里用θ这一个阈值来控制。每次递归遍历树时，先看当前树节点往外看：设s是这个节点覆盖区域的边长，d是目标粒子到该节点质心的距离。如果s/d小于θ，说明这个节点在目标粒子看来“足够小、足够远”，可以直接用质心替代内部所有粒子；否则就继续往下展开，直到满足条件或者到达叶子节点。
if s / d < theta: # 远处一团粒子，用质量集中在质心上的等效粒子 a += G * node.mass * (node.cx - px) / (d^2 + eps^2) ^ 1.5 else: for child in node.children: compute(child)θ的取值直接影响精度和性能。θ越小近似越严格，θ=0时几乎等于完全展开，精度逼近暴力算法但速度也退回去；θ越大近似越激进，速度越快但轨道可能变形。实践中我一般先取0.5~0.7，追求视觉稳定用0.5，追求速度用0.7~0.9。后面会专门讲参数该怎么定。
2.3 四叉树：区域划分与递归路径
Barnes-Hut的树结构和普通的搜索树不同，它按空间划分而不是按数值划分。以2D为例，根节点覆盖整个仿真区域，插入粒子时先判断粒子落在当前节点的哪个象限，然后递归进入那个象限的子节点。当一个节点里同时出现两个位置不同的粒子，就把该节点分裂成四个子节点，然后把它们分别推下去。
树建好之后，每个节点都保存三样东西：子树总质量、子树质心坐标、四个子节点指针。注意这个质心是“质量的加权中心”，不是“几何中心”，两者在质量分布不均匀时有明显差别。后面代码中先更新质量再更新质心的顺序就是最容易踩坑的地方。
我一般把树的构建和受力计算分成两个阶段：每一步先根据当前粒子位置重建整棵树，再对每个粒子用递归遍历算受力。树和受力计算是同一个递归逻辑，但心思要分开。建树阶段保证每个节点信息正确，受力阶段保证近似条件正确，绝大多数bug都出在这两个阶段之间。
3. 构建四叉树：节点定义与质心上浮
3.1 节点数据结构与参数设计
四叉树的节点需要记录区域边界、质量和质心、子节点列表，以及叶子节点里存放的粒子索引。这里有一个容易被忽略的设计决策：叶子里存“多个粒子”还是“唯一粒子”。我选择让叶子节点持有一个粒子索引列表。当两个粒子位置完全重合时，直接放进同一个叶子的列表里，而不是继续分裂，这样从根源上避免递归深度失控。
class Node: """四叉树节点。x0,y0,x1,y1 是当前节点覆盖的方形区域。""" def __init__(self, x0, y0, x1, y1): self.x0, self.y0, self.x1, self.y1 = x0, y0, x1, y1 self.mass = 0.0 # 子树总质量 self.cx, self.cy = 0.0, 0.0 # 子树质心 self.children = [] # 内部节点：最多4个子节点 self.particles = [] # 叶子节点：存放粒子索引质量在Demo里统一取1.0，但代码结构完全支持每个粒子不同质量，只需要在插入时把new_mass从固定值换成粒子自己的质量数组。用类而不是字典的好处是属性访问快、代码直观，缺点是Python对象的开销偏高。真正高性能的做法是后面章节要说的“结构数组”方案，这里先把逻辑讲清楚。
3.2 插入粒子与递归分裂
插入的核心逻辑是：新粒子进入节点后，先更新节点的质量和质心，然后判断当前节点是叶子还是内部节点。叶子节点如果还没有粒子就直接存下来；已经有粒子且位置不同，就要分裂；如果位置重合，则保留在同一叶子中。分裂时要保证“每个粒子唯一归属一个子节点”，我统一用x < mid_x和y < mid_y判断西南象限，其余情况分别归入东南、西北、东北。
def insert(self, particle_index, xs, ys): # 先更新质量，再用旧质量加权更新质心，顺序不能反 new_mass = 1.0 self.cx = (self.cx * self.mass + xs[particle_index] * new_mass) / (self.mass + new_mass) self.cy = (self.cy * self.mass + ys[particle_index] * new_mass) / (self.mass + new_mass) self.mass += new_mass if self.is_leaf(): if not self.particles: self.particles.append(particle_index) return # 已有一个粒子：若坐标重合则合并到同一叶子，避免无限分裂 other = self.particles[0] if xs[particle_index] == xs[other] and ys[particle_index] == ys[other]: self.particles.append(particle_index) return old_list = self.particles self.particles = [] for idx in old_list + [particle_index]: self._push_down(idx, xs, ys) else: self._push_down(particle_index, xs, ys)这里的_push_down负责把粒子放到正确的子象限去。如果对应子节点还不存在就现场创建。需要注意的是，内部节点的particles列表永远是空的，因为一旦分裂后所有粒子都转移到了子节点中。
def _push_down(self, idx, xs, ys): mx = (self.x0 + self.x1) * 0.5 my = (self.y0 + self.y1) * 0.5 quad = self._quadrant(xs[idx], ys[idx], mx, my) child = self._get_child(quad, mx, my) child.insert(idx, xs, ys) def _quadrant(self, x, y, mx, my): # 等号严格归入右侧/上侧，保证同一个粒子只有一个归属 if x < mx and y < my: return 0 if x >= mx and y < my: return 1 if x < mx and y >= my: return 2 return 3分裂判断里有一个血泪经验：象限判断必须用“小于”和“大于等于”，不能一正一反都用等号。否则当粒子恰好落在边界线上时，同一个粒子可能被分到两个子节点，树结构立刻混乱。
3.3 批量构建与递归进深的边界
批量建树只需要循环插入即可。根节点覆盖的区域要提前框定，常见做法是先扫一遍所有粒子求x和y的min/max，再向外扩一点边距，防止粒子刚好落在边界上。
def build_tree(xs, ys): x_min, x_max = min(xs), max(xs) y_min, y_max = min(ys), max(ys) margin = 0.01 * max(x_max - x_min, y_max - y_min, 1e-9) root = Node(x_min - margin, y_min - margin, x_max + margin, y_max + margin) for i in range(len(xs)): root.insert(i, xs, ys) return root这里有一个深度隐患。虽然位置重合的粒子会被合并到同一叶子，但如果两个粒子距离极其接近但不重合，比如相差1e-12，树还是会一层层分裂下去，直到区域边长小于间距。最坏情况下递归深度能达到几百甚至几千，触发Python的递归上限。应对办法有两个：一是给insert加深度参数，超过一定深度就直接把粒子并进当前叶子不再分裂；二是在初始化时对粒子坐标做极小随机扰动。我建议两个都做，后面避坑章节还会细说。
4. 遍历树计算受力：核心代码与最小可运行模拟器
4.1 递归受力函数与θ使用细节
受力计算的递归思路和建树对称：从根节点出发，对每个节点先判断能否用质心近似，不能就进入子节点继续。每个粒子的受力计算单独跑一次遍历，所以每步要对N个粒子各自调用一次这个递归函数。叶子节点里可能存了多个粒子（重合或深度限制所致），遍历它保存的粒子列表时要注意跳过自己。
def compute_acceleration(self, idx, px, py, xs, ys, theta, eps, G=1.0): """对粒子idx求加速度。px,py是粒子位置，eps是软化参数。""" ax, ay = 0.0, 0.0 def walk(node): nonlocal ax, ay dx = node.cx - px dy = node.cy - py d2 = dx * dx + dy * dy if node.is_leaf(): for j in node.particles: if j == idx: continue # 自身不能对自己产生引力 fx = xs[j] - px fy = ys[j] - py r2 = fx * fx + fy * fy + eps * eps inv = 1.0 / (r2 * math.sqrt(r2)) ax += G * fx * inv ay += G * fy * inv else: s = max(node.x1 - node.x0, node.y1 - node.y0) if s * s < theta * theta * d2: # 远处节点，用质心近似；s/d < theta 等价于 s*s < theta*theta*d2 inv = 1.0 / ((d2 + eps * eps) * math.sqrt(d2 + eps * eps)) ax += G * node.mass * dx * inv ay += G * node.mass * dy * inv else: for child in node.children: walk(child) walk(self.root) return ax, ay这里有三个细节值得展开。第一，用s * s < theta * theta * d2代替s / d < theta，避免每次开根号的开销，虽然演示版看不出差别，但粒子数一上去就差很多。第二，软化参数eps必须加在距离平方上，而不是只在距离上取绝对值，否则两个粒子几乎重合时引力会趋向无穷大。第三，内部节点的node.cx和node.cy是子树质心，不是节点几何中心，算距离时千万别用错。
4.2 积分步进：Velocity-Verlet让能量不飘
受力算完后如何更新位置和速度，直接决定长期模拟是否稳定。最简单的欧拉法每步都会积累能量误差，跑几百步后系统明显膨胀。我用Velocity-Verlet方案，也就是kick-drift-kick两步式，它的优势是时间可逆，能量漂移比欧拉法低两个数量级。
def step(xs, ys, vxs, vys, tree, dt, theta, eps): """用kick-drift-kick更新一步。velocity verlet变体。""" n = len(xs) ax = [0.0] * n ay = [0.0] * n for i in range(n): ax[i], ay[i] = tree.compute_acceleration(i, xs[i], ys[i], xs, ys, theta, eps) for i in range(n): vxs[i] += ax[i] * dt * 0.5 vys[i] += ay[i] * dt * 0.5 xs[i] += vxs[i] * dt ys[i] += vys[i] * dt # 粒子位置变了，必须重建树，再用新位置算后半步加速度 tree = build_tree(xs, ys) for i in range(n): ax[i], ay[i] = tree.compute_acceleration(i, xs[i], ys[i], xs, ys, theta, eps) for i in range(n): vxs[i] += ax[i] * dt * 0.5 vys[i] += ay[i] * dt * 0.5 return tree第二步重新建树容易被新手忽略。粒子的位置已经更新了，旧树对应的空间划分是过时的，继续用旧树算受力就相当于用“上一帧的位置关系”算这一帧的力，轨道会越跑越歪。
4.3 最小可运行的星系盘Demo
把以上三个部分拼起来，就是一个能跑的2D引力模拟器。初始条件我用一个中心大质量粒子和一圈环绕的小粒子，速度近似圆轨道速度，加一点随机扰动。跑200步，每步打印动能、势能和总能量，方便直接观察能量是否漂移。
def main(): import math, random n = 500 xs = [0.0] * n ys = [0.0] * n vxs = [0.0] * n vys = [0.0] * n # 中央大质量粒子 xs[0], ys[0], vxs[0], vys[0] = 0.0, 0.0, 0.0, 0.0 big_mass = 500.0 # 其余粒子分布在半径0.5~2.5的圆盘上，速度约等于开普勒速度 for i in range(1, n): r = random.uniform(0.6, 2.5) phi = random.uniform(0, 2 * math.pi) xs[i] = r * math.cos(phi) ys[i] = r * math.sin(phi) v = math.sqrt(G * big_mass / r) * random.uniform(0.8, 1.2) vxs[i] = -v * math.sin(phi) vys[i] = v * math.cos(phi) tree = build_tree(xs, ys) dt, theta, eps = 0.01, 0.6, 0.01 for step in range(200): tree = step(xs, ys, vxs, vys, tree, dt, theta, eps) if step % 20 == 0: ke = 0.5 * sum(vxs[i]**2 + vys[i]**2 for i in range(n)) print(f"step={step:4d} KE={ke:.3f}")这段代码跑起来后，中心大质量及周围粒子会形成一个相对稳定的盘状结构。需要注意的是，这个简单Demo没有处理粒子之间的碰撞合并，也没有自适应时间步，所以跑几百步后边缘粒子可能逃逸，这属于正常现象。
5. Barnes-Hut实战避坑：5个常见问题与排查
5.1 粒子恰好落在象限边界上，树结构错乱
现象：模拟中个别粒子突然获得巨大加速度，甚至直接飞出场外，而且每次飞出的粒子都集中在坐标轴附近。
原因：象限判定时<=和>=混用，导致同一个粒子在分裂时被同时分进两个子节点。树结构因此出现“双亲”，受力计算时某些粒子被重复统计，引力瞬间翻倍。
解决：统一用“左小右大、下小上大”的规则。我在代码中写明x < mx归西南，x >= mx再按y继续分。修改后需要重置根节点重新建树，旧树不会自动修正。这个坑只靠肉眼看很难发现，建议专门构造一个(0.5, 0.5)边界上的粒子做单步调试。
5.2 质心计算顺序错误，星系被“撕裂”成两团
现象：粒子云整体在空间上分裂成两坨，并且两坨互相绕转，但质量分布看起来不对称，好像每个节点都偏向某一侧。
原因：更新质心时，先用新质量除以了新总质量，但分子里的旧质心还是乘以旧质量。也就是代码写成了cx = (cx * mass + x) / (mass + 1)，而正确的分子应该是cx * mass + x * 1。这个错误在小规模时只造成轻微偏差，粒子多时误差不断累积放大。
解决：把质心更新公式固定为“先算新的加权平均，再更新总质量”，或者反过来先更新质量再用旧质量分子。我在代码注释里特别标了顺序。自查方法：造两个质量相同的粒子放在(1,0)和(-1,0)，正确质心为(0,0)；如果结果偏向某一侧，就是这里写错了。
5.3 粒子重合造成递归深度爆炸
现象：模拟到中途突然报RecursionError: maximum recursion depth exceeded，或者程序直接闪退。
原因：两个粒子距离极近但不重合时，树会一直分裂，因为子节点边长不断减半，最终递归深度超过Python默认上限。更隐蔽的情况是粒子位置相等但不合并，insert永远在同一路径上分裂下去。
解决：我在insert里对重合粒子做了合并到同一叶子的处理。对于“极近但不重合”的情况，给树加最大深度限制，比如深度超过64层就直接把粒子追加到当前叶子的列表中，不再分裂。这样会损失一点局部精度，但保证了程序不崩。工程上更稳的做法是在读入数据时先对粒子坐标做去重或加最小间距扰动。
5.4 θ参数不当：精度和速度同时崩
现象：θ设为1.5以上时，粒子轨道明显变椭圆，能量快速增大；θ设为0.05时，N=2000粒子就卡得几乎不动，比暴力算法还慢。
原因：θ太大导致远处节点太容易被当作质点，引力方向误差无法控制。θ太小时，几乎所有节点都要递归展开，Barnes-Hut退化成了“用树遍历的暴力算法”，既没有精度优势也没有速度优势。
解决：θ选0.5~0.7作为默认区间。我一般用0.6跑大尺度演示，用0.4跑需要轨道细节的场景。调参时直接看两个指标：总能量漂移率和单步耗时。漂移率超过5%就把θ调小，单步耗时超过预期就把θ调大。这个平衡就是Barnes-Hut的算法玄学，没有公式能替代实测。
5.5 时间步长过大导致粒子逃逸
现象：前几步能量看起来还正常，几十步后粒子速度突然暴增，盘结构消失。
原因：dt过大，粒子在一个时间步内运动的距离已经超过它附近的树节点尺寸。它跨过大半个区域，受力计算时用的还是旧位置，近似失真严重。有人把这个问题误认为是θ太大，其实把θ调到接近0也救不回来。
解决：收缩dt。对半径为1~2的星系盘，dt取0.005到0.01比较稳。想省事可以加一个简单自适应：每步统计最大速度，如果max_speed * dt > 0.05 * cell_size_min，就把dt减半重跑本步。速度自适应能大幅减少人工调参次数，尤其适合粒子从密集到稀疏变化剧烈的模拟。
6. 进阶技巧：用能量守恒验证你的模拟器写对了
写完Barnes-Hut代码后，第一件事不是调θ，而是验证它“写对了”。我的习惯是先用小规模暴力法做基准：取N=200，分别跑暴力算法和Barnes-Hut（θ设为0.01），对比每一步的加速度，相对误差应低于1%。这能排除建树和受力逻辑里的隐蔽bug。之后再把θ调到0.6，对比此时Barnes-Hut与暴力结果，观察误差主要来自近似而非代码错误。前者证明代码正确，后者证明算法参数合理。
第二个验证指标是能量漂移率。取一个稳定的双星或圆盘初始条件，跑1000步，记录总能量。如果每步的相对能量漂移稳定在0.1%量级，说明积分器和树构建都没问题；如果出现能量单向单调增加或减小，优先检查softening是否加了、dt是否过大、质心更新顺序是否错误。我见过有人把速度更新写成了v += a * dt之后又额外加一次a * dt，计算量翻倍而能量照样飘，这种bug用能量曲线一眼就能看出来。
跑通后想提速，有几个立竿见影的方向。第一是改数据结构：把每个节点从Python对象换成扁平数组，用索引指向子节点，内存连续且GC开销小。第二是向量化：在受力计算层用numpy数组处理一批距离和受力，把Python递归拆成显式栈。第三是并行化：同一帧内不同粒子的受力计算互不依赖，把粒子分成几组丢给多线程即可；但需要注意所有线程共享同一棵只读树，不能一边遍历一边改树。我一般先做数据扁平化，再做numpy向量化，最后才上多线程。
我在调参上的教训是：用θ=0.6做视觉演示，用θ=0.3做科学分析，永远不要为了赶帧率把θ拉到1.0以上。宁可少放几百个粒子，也要保住轨道的定性正确。希望帮到你。
本文还有配套的精品资源，点击获取
