Skip to content
Computational Lab.

Barnes-Hut Quadtree

Approximating the far field with a tree.

Advanced / tree / gravity / approximation

Barnes-Hut makes large N-body simulation affordable: build a quadtree over the bodies, and when a distant cell's size-to-distance ratio drops below theta, treat the whole cluster as one mass. The force-error metric compares the approximation against direct summation - the theta slider is the dial between speed and honesty.

Tagged: quadtree, n-body, theta criterion, spatial partitioning

Live simulation

The theta slider trades force accuracy for speed; the error metric audits it.

Speed

The theta slider trades force accuracy for speed; the error metric audits it.

Mathematics

Barnes-Hut approximates the gravitational effect of a distant cluster by its center of mass whenever the cluster subtends a small enough angle: s/d < theta, cell size over distance. The multipole acceptance criterion turns an O(n^2) pairwise sum into roughly O(n log n).

a_i = G sum over clusters near i of m_c (r_c - r_i) / (|r_c - r_i|^2 + eps^2)^(3/2)

accept cluster when s / d < theta; else subdivide

force error vs direct sum: err = mean |F_BH - F_direct| / |F_direct|

The theta slider is the honesty dial: theta near 0 forces the tree to open almost every node, converging to the direct sum (and the error metric agrees); theta near 1.4 opens almost nothing and the error shows it. Watching the measured error respond to theta is the calibration exercise every force-field code runs.

How it works

The quadtree is rebuilt every tick into flat typed arrays: four children per node, a body slot per leaf, and per-node mass and center of mass computed by one upward pass. Force queries walk from the root, opening nodes that fail the theta test.

Every 90 ticks the engine runs the direct O(n^2) sum on the same positions and reports the mean relative force error - the number that tells you the approximation is real.

Implementation

src/engine/experiments/barnes-hut.ts · 349 lines · the module running above

TypeScript · fixed-timestep engine

1import { BaseExperiment } from "../core/base-experiment";
2import type { ParameterDef, ThemeColors } from "../core/types";
3
4const CAPACITY = 320;
5
6/**
7 * Barnes-Hut: build a quadtree over the bodies each tick and approximate
8 * distant clusters by their center of mass when the node's size/distance
9 * ratio falls below theta. The accuracy metric compares the total force
10 * against a direct O(n^2) sum sampled periodically.
11 */
12export class BarnesHutExperiment extends BaseExperiment {
13 readonly id = "barnes-hut";
14
15 private count = 0;
16 private x = new Float32Array(CAPACITY);
17 private y = new Float32Array(CAPACITY);
18 private vx = new Float32Array(CAPACITY);
19 private vy = new Float32Array(CAPACITY);
20 private fx = new Float32Array(CAPACITY);
21 private fy = new Float32Array(CAPACITY);
22 // Quadtree in flat arrays: node -> children (4), center of mass, mass, bounds.
23 private nodeChild = new Int32Array(0);
24 private nodeMass = new Float64Array(0);
25 private nodeComX = new Float64Array(0);
26 private nodeComY = new Float64Array(0);
27 private nodeX = new Float64Array(0);
28 private nodeY = new Float64Array(0);
29 private nodeSize = new Float64Array(0);
30 private nodeBody = new Int32Array(0);
31 private nextBody = new Int32Array(CAPACITY);
32 private nodeCount = 0;
33 private directFx = new Float32Array(0);
34 private directFy = new Float32Array(0);
35 private forceError = 0;
36 private directTick = 0;
37
38 protected params(): ParameterDef[] {
39 return [
40 { key: "bodies", label: "Bodies", min: 60, max: 320, step: 20, defaultValue: 180 },
41 { key: "theta", label: "Opening angle theta", min: 0.2, max: 1.4, step: 0.1, defaultValue: 0.8 },
42 {
43 key: "showTree",
44 label: "Show quadtree",
45 options: [
46 { value: "off", label: "Off" },
47 { value: "on", label: "On" },
48 ],
49 defaultValue: "off",
50 },
51 ];
52 }
53
54 private readonly G = 9000;
55 private readonly softening = 4;
56
57 protected onReset(): void {
58 this.count = this.num("bodies");
59 for (let i = 0; i < this.count; i++) {
60 // Spiral galaxy initial condition.
61 const arm = i % 3;
62 const radius = 20 + (i / this.count) * Math.min(this.width, this.height) * 0.44;
63 const angle = radius * 0.02 + (arm * Math.PI * 2) / 3 + (this.rng() - 0.5) * 0.6;
64 this.x[i] = this.width / 2 + Math.cos(angle) * radius;
65 this.y[i] = this.height / 2 + Math.sin(angle) * radius;
66 const v = Math.sqrt((this.G * 60 * (i / this.count + 0.3)) / radius);
67 this.vx[i] = -Math.sin(angle) * v;
68 this.vy[i] = Math.cos(angle) * v;
69 }
70 this.directFx = new Float32Array(this.count);
71 this.directFy = new Float32Array(this.count);
72 }
73
74 protected onParameterChange(key: string): void {
75 if (key === "bodies") this.reset();
76 }
77
78 resize(width: number, height: number): void {
79 this.width = width;
80 this.height = height;
81 this.reset();
82 }
83
84 /**
85 * Insert one body into the quadtree. Leaves hold linked lists of bodies;
86 * a leaf with more than one body subdivides (depth-capped), so coincident
87 * positions never drop mass from the tree.
88 */
89 private insert(body: number): void {
90 let node = 0;
91 let depth = 0;
92 while (this.nodeBody[node] === -2) {
93 depth += 1;
94 const q = this.quadrant(node, this.x[body], this.y[body]);
95 if (this.nodeChild[node * 4 + q] === -1) {
96 this.nodeChild[node * 4 + q] = this.newNode(
97 this.childX(node, q),
98 this.childY(node, q),
99 this.nodeSize[node] / 2,
100 );
101 }
102 node = this.nodeChild[node * 4 + q];
103 }
104 this.placeInLeaf(node, body, depth);
105 }
106
107 private placeInLeaf(node: number, body: number, depth: number): void {
108 this.nextBody[body] = this.nodeBody[node];
109 this.nodeBody[node] = body;
110 if (depth >= 48 || this.listLength(node) <= 1) return;
111 // Subdivide: detach the list, mark internal, re-place every body.
112 let b = this.nodeBody[node];
113 const list: number[] = [];
114 while (b !== -1) {
115 list.push(b);
116 b = this.nextBody[b];
117 }
118 this.nodeBody[node] = -2;
119 for (const bb of list) {
120 const q = this.quadrant(node, this.x[bb], this.y[bb]);
121 if (this.nodeChild[node * 4 + q] === -1) {
122 this.nodeChild[node * 4 + q] = this.newNode(
123 this.childX(node, q),
124 this.childY(node, q),
125 this.nodeSize[node] / 2,
126 );
127 }
128 this.placeInLeaf(this.nodeChild[node * 4 + q], bb, depth + 1);
129 }
130 }
131
132 private listLength(node: number): number {
133 let n = 0;
134 for (let b = this.nodeBody[node]; b !== -1; b = this.nextBody[b]) n += 1;
135 return n;
136 }
137
138 private newNode(x: number, y: number, size: number): number {
139 if (this.nodeCount >= this.nodeMass.length) {
140 const grow = <T extends Float64Array | Int32Array>(arr: T, minLen: number): T => {
141 const next = new (arr.constructor as { new(length: number): T })(Math.max(minLen, arr.length * 2));
142 next.set(arr);
143 return next;
144 };
145 this.nodeChild = grow(this.nodeChild, 64).fill(-1, this.nodeCount * 4);
146 this.nodeMass = grow(this.nodeMass, 16);
147 this.nodeComX = grow(this.nodeComX, 16);
148 this.nodeComY = grow(this.nodeComY, 16);
149 this.nodeX = grow(this.nodeX, 16);
150 this.nodeY = grow(this.nodeY, 16);
151 this.nodeSize = grow(this.nodeSize, 16);
152 this.nodeBody = grow(this.nodeBody, 16);
153 }
154 const id = this.nodeCount++;
155 this.nodeChild[id * 4] = -1;
156 this.nodeChild[id * 4 + 1] = -1;
157 this.nodeChild[id * 4 + 2] = -1;
158 this.nodeChild[id * 4 + 3] = -1;
159 this.nodeMass[id] = 0;
160 this.nodeX[id] = x;
161 this.nodeY[id] = y;
162 this.nodeSize[id] = size;
163 this.nodeBody[id] = -1;
164 return id;
165 }
166
167 private quadrant(node: number, x: number, y: number): number {
168 const half = this.nodeSize[node] / 2;
169 const right = x > this.nodeX[node] + half;
170 const bottom = y > this.nodeY[node] + half;
171 return (right ? 1 : 0) + (bottom ? 2 : 0);
172 }
173
174 private childX(node: number, q: number): number {
175 return this.nodeX[node] + (q % 2 === 1 ? this.nodeSize[node] / 2 : 0);
176 }
177
178 private childY(node: number, q: number): number {
179 return this.nodeY[node] + (q >= 2 ? this.nodeSize[node] / 2 : 0);
180 }
181
182 private computeMass(node: number): void {
183 if (this.nodeBody[node] >= 0) {
184 let mass = 0;
185 let comX = 0;
186 let comY = 0;
187 for (let b = this.nodeBody[node]; b !== -1; b = this.nextBody[b]) {
188 mass += 1;
189 comX += this.x[b];
190 comY += this.y[b];
191 }
192 this.nodeMass[node] = mass;
193 this.nodeComX[node] = comX / mass;
194 this.nodeComY[node] = comY / mass;
195 return;
196 }
197 let mass = 0;
198 let comX = 0;
199 let comY = 0;
200 for (let q = 0; q < 4; q++) {
201 const child = this.nodeChild[node * 4 + q];
202 if (child === -1) continue;
203 this.computeMass(child);
204 mass += this.nodeMass[child];
205 comX += this.nodeMass[child] * this.nodeComX[child];
206 comY += this.nodeMass[child] * this.nodeComY[child];
207 }
208 this.nodeMass[node] = mass;
209 this.nodeComX[node] = mass > 0 ? comX / mass : this.nodeX[node];
210 this.nodeComY[node] = mass > 0 ? comY / mass : this.nodeY[node];
211 }
212
213 /** Accumulate force on body i by walking the tree with the theta criterion. */
214 private treeForce(i: number, node: number): void {
215 if (this.nodeMass[node] === 0) return;
216 // Leaves are always resolved exactly, body by body.
217 if (this.nodeBody[node] >= 0) {
218 for (let b = this.nodeBody[node]; b !== -1; b = this.nextBody[b]) {
219 if (b === i) continue;
220 const ddx = this.x[b] - this.x[i];
221 const ddy = this.y[b] - this.y[i];
222 const dd = Math.sqrt(ddx * ddx + ddy * ddy + this.softening * this.softening);
223 const force = (this.G * 1) / (dd * dd * dd);
224 this.fx[i] += ddx * force;
225 this.fy[i] += ddy * force;
226 }
227 return;
228 }
229 const dx = this.nodeComX[node] - this.x[i];
230 const dy = this.nodeComY[node] - this.y[i];
231 const d = Math.sqrt(dx * dx + dy * dy + this.softening * this.softening);
232 // Far internal node: one center-of-mass interaction.
233 if (this.nodeSize[node] / d < this.num("theta")) {
234 const force = (this.G * this.nodeMass[node]) / (d * d * d);
235 this.fx[i] += dx * force;
236 this.fy[i] += dy * force;
237 return;
238 }
239 for (let q = 0; q < 4; q++) {
240 const child = this.nodeChild[node * 4 + q];
241 if (child !== -1) this.treeForce(i, child);
242 }
243 }
244
245 private directForce(): void {
246 for (let i = 0; i < this.count; i++) {
247 let fxi = 0;
248 let fyi = 0;
249 for (let j = 0; j < this.count; j++) {
250 if (j === i) continue;
251 const dx = this.x[j] - this.x[i];
252 const dy = this.y[j] - this.y[i];
253 const d = Math.sqrt(dx * dx + dy * dy + this.softening * this.softening);
254 const s = (this.G * 1) / (d * d * d);
255 fxi += dx * s;
256 fyi += dy * s;
257 }
258 this.directFx[i] = fxi;
259 this.directFy[i] = fyi;
260 }
261 }
262
263 protected onUpdate(dt: number): void {
264 // Root bounds must enclose every body, or the tree's s/d criterion lies.
265 let minX = Infinity;
266 let minY = Infinity;
267 let maxX = -Infinity;
268 let maxY = -Infinity;
269 for (let i = 0; i < this.count; i++) {
270 if (this.x[i] < minX) minX = this.x[i];
271 if (this.y[i] < minY) minY = this.y[i];
272 if (this.x[i] > maxX) maxX = this.x[i];
273 if (this.y[i] > maxY) maxY = this.y[i];
274 }
275 const span = Math.max(maxX - minX, maxY - minY) * 1.05 + 1;
276 this.nodeChild = new Int32Array(this.count * 8 * 4).fill(-1);
277 this.nodeMass = new Float64Array(this.count * 8);
278 this.nodeComX = new Float64Array(this.count * 8);
279 this.nodeComY = new Float64Array(this.count * 8);
280 this.nodeX = new Float64Array(this.count * 8);
281 this.nodeY = new Float64Array(this.count * 8);
282 this.nodeSize = new Float64Array(this.count * 8);
283 this.nodeBody = new Int32Array(this.count * 8).fill(-1);
284 this.nodeCount = 0;
285 this.newNode(minX, minY, span);
286 for (let i = 0; i < this.count; i++) this.insert(i);
287 this.computeMass(0);
288
289 // Forces first, for every body, at this tick's positions...
290 for (let i = 0; i < this.count; i++) {
291 this.fx[i] = 0;
292 this.fy[i] = 0;
293 this.treeForce(i, 0);
294 }
295 // ...so the periodic audit compares direct and tree forces on the SAME
296 // positions - never a stale sample from before the bodies moved.
297 this.directTick += 1;
298 if (this.directTick % 90 === 0) {
299 this.directForce();
300 let err = 0;
301 let total = 0;
302 for (let i = 0; i < this.count; i++) {
303 err += Math.hypot(this.fx[i] - this.directFx[i], this.fy[i] - this.directFy[i]);
304 total += Math.hypot(this.directFx[i], this.directFy[i]);
305 }
306 this.forceError = total > 0 ? err / total : 0;
307 }
308 for (let i = 0; i < this.count; i++) {
309 this.vx[i] += this.fx[i] * dt;
310 this.vy[i] += this.fy[i] * dt;
311 this.x[i] += this.vx[i] * dt;
312 this.y[i] += this.vy[i] * dt;
313 }
314 }
315
316 render(ctx: CanvasRenderingContext2D, theme: ThemeColors): void {
317 ctx.fillStyle = theme.bg;
318 ctx.fillRect(0, 0, this.width, this.height);
319 if (this.str("showTree") === "on") {
320 ctx.strokeStyle = theme.line;
321 ctx.lineWidth = 1;
322 for (let node = 1; node < this.nodeCount; node++) {
323 ctx.strokeRect(this.nodeX[node], this.nodeY[node], this.nodeSize[node], this.nodeSize[node]);
324 }
325 }
326 for (let i = 0; i < this.count; i++) {
327 ctx.fillStyle = i % 24 === 0 ? theme.accent : theme.fg;
328 ctx.fillRect(this.x[i] - 1.2, this.y[i] - 1.2, 2.4, 2.4);
329 }
330 }
331
332 getMetrics() {
333 return {
334 bodies: this.count,
335 quadtreeNodes: this.nodeCount.toLocaleString("en-US"),
336 theta: this.num("theta").toFixed(1),
337 forceError: `${(this.forceError * 100).toFixed(2)} %`,
338 };
339 }
340
341 describe(): string {
342 return `${this.count} bodies under Barnes-Hut gravity: a ${this.nodeCount.toLocaleString("en-US")}-node quadtree approximates distant clusters within opening angle ${this.num("theta").toFixed(1)}, measured force error ${(this.forceError * 100).toFixed(2)}%.`;
343 }
344
345 entities(): number {
346 return this.count;
347 }
348}
349

Performance

Time per tick
O(n log n) per tick for the tree plus force walk
Space
O(n) tree nodes per tick
Notes
A spiral-galaxy initial condition makes clusters visible: turn on the quadtree overlay and watch the box boundaries follow the density.

Live numbers (FPS, frame time, steps per second, microseconds per step, entity count) are in the performance overlay on the simulation above: the pulse icon in its top-right corner.