Barnes-Hut, from scratch

drag rotatectrl + scroll zoomWASDrotateR / F zoom

Introduction

I still remember the day I discovered numerical integration in my numerical physics course; it felt like a bright, almost magical moment. I suddenly realized that I finally could take my revenge as a first-year physics student trying to solve the three-body problem (not knowing it has no general solution)... I went straight to my computer and coded a small simulation ! The dance of those three spinning planets was printing stars on my eyes ! Three planets ! Then the Solar System ! Why not the whole galaxy !? I was incredibly enthusiastic but so naive... I broke my teeth on the brutal tyranny of $O(N^2)$.

A few years later, I discovered the Barnes–Hut algorithm, and finally, I defeated the tyranny... Here is how to do it :

---

The N-body Problem

The $N$-body problem is a classical mechanics system in which $N$ bodies interact through a given force. I'm going to have the audacity to remind you briefly how this works...

The dynamics of the system obey Newton's second law: the acceleration of each body depends on the force applied to it,

$$m_i \vec{a}_i(t) = \vec{F}_i$$

In the case of a gravitational system, the force acting on body $i$ depends on the positions of all the other bodies:

$$\vec{F}_i = \sum_{j \neq i}^N \frac{G m_i m_j}{||\vec{x}_i - \vec{x}_j||^2+\epsilon^2}\vec{u}_{ij}$$

Here, $\vec{u}_{ij}$ is the unit vector going from $\vec{x}_i$ to $\vec{x}_j$, $G$ is the gravitational constant, and $\epsilon$ is a softening (regularization) parameter introduced to avoid numerical divergences when bodies get very close to each other.

Once the forces are computed, we can integrate the system in time. An efficient (and underrated !) choice is the Verlet integrator:

$$\vec{x}(t+\Delta t) = 2 \vec{x}(t) - \vec{x}(t-\Delta t) + \vec{a}(t)\,\Delta t^2 + \mathcal{O}(\Delta t^4)$$

In principle, this allows us to compute the positions of all bodies at the next time step.

However, the computational cost of evaluating the forces quickly becomes prohibitive as $N$ increases. Indeed, the naive approach requires computing all pairwise interactions, leading to a complexity of $O(N^2)$, which becomes the major bottleneck for large systems.

---

Barnes–Hut Algorithm

This is where our friends Barnes and Hut come to the rescue (if they are indeed real people, to check...). The algorithm that took their name allows us, with a slightly more complicated implementation, to reduce the computational complexity to $O(N \log N)$.

The key idea is to exploit the fact that gravitational forces decrease with distance. Far-away groups of bodies can be approximated as a single "virtual" body, reducing the number of interactions that need to be computed.

The algorithm works by building an octree structure: space is recursively subdivided, and bodies are organized hierarchically according to their positions. Groups of distant particles are then approximated by virtual quantities (virtual mass and center of mass).

Everything is controlled by an hyperparameter $\theta$. At each time step, the algorithm proceeds as follows:

1. Build the Tree

algorithm · Building the treeloopifelsefunction↻calls itself

inputthe bodies

outputevery cell, with its virtual body

1function Build(cell)

2if the cell holds more than one body

3cut it into 4 sub-cells (8 in 3D)

4↻ do the same with each sub-cell

5the cell's virtual body ← total mass, at the centre of mass

6Build(the whole space)

2. Compute Forces

algorithm · Force on each bodyloopifelsefunction↻calls itself

inputthe tree, the angle $\theta$

outputthe force on each body

1function Force(cell, body)

2// $s$ = size of the cell, $d$ = its distance to the body

3if the cell is a single body, or $s/d < \theta$

4return the pull of the cell, as one body

5else

6for each sub-cell

7↻ add Force(sub-cell, body)

8for each body

9its force ← Force(the whole tree, body)

In Action

Once all forces are computed, we integrate the system in time using the Verlet integrator.

The $\log N$ factor in the complexity comes from the tree traversal and recursive subdivision of space. You can see this method as a 3D analog of a dichotomy algorithm.

2D Live Version

Here is a small galaxy simulation running live in 2D. You can click on a star to see it's interacting quadtree live !

Implementation and Visualization

This algorithm is particularly well-suited for large-scale simulations and high-performance computing. It is, in fact, a bit overengineered for small systems like $N = 300$.

For this reason, the implementation was done directly in C++, using object-oriented programming for clarity. Only the force computation was parallelized, as the tree construction itself is not trivially parallelizable.

For the largest simulations ($N = 100000$), we used an AWS EC2 instance with a ~100 CPUs.

Another challenge at this scale is visualization. I developed a small Processing-based pipeline for this purpose, but it remains relatively slow and limited. A more robust solution could be to explore visualization tools such as Blender or C++ renderer.

---

Hyperparameter, Error, and Complexity

The quality of the force calculation is entirely controlled by the hyperparameter $\theta$. Its impact on the "resolution" of interactions is easy to see in the interactive figure above: slide $\theta$.

plots

Impact of θ on the mean acceleration error (left) and on the computation time of one step (right).

As expected, decreasing $\theta$ improves accuracy, but also increases computational cost.

In practice, this trade-off is not a major issue. For large $N$, the system becomes highly chaotic, and obtaining precise quantitative agreement is often unrealistic. Since we are mainly interested in qualitative behavior, we can choose a value of $\theta$ that provides a good balance between accuracy and performance.

We can also compare the computational complexity of the naive approach, $O(N^2)$, with the Barnes–Hut method, $O(N \log N)$:

plot

Computational complexity comparison between naive and Barnes-Hut implementation.

---

Results

Finally, we can blow up our AWS budget and simulate the birth of a small galaxy; the whole point of all those efforts !

To initialize the system, we randomly place bodies within a flattened ellipsoid and assign them radial velocities according to the velocity distribution of a self-gravitating system.

Here is what we obtain:

Time evolution of a N=100000 gravitational system.

Outlooks

• Big Bang simulation
• BIGGEEEEEER
• Recode properly everything
• Geometric center or center of mass splitting ?

---