A region of space contains the number of N bodies. Each body has its own potential and a resulting force field. This could be charges causing an electrical field (Coulomb's law) or planets in space (Law of Gravity). The potential of the bodies can be summed up yielding a combined potential which depends on the location of each of the bodies and its physical properties (i.e. mass or charge). According the Newtons first law the bodies themself will experience an acceleration caused by the field. As a result there is a force field that depends on the locations of the n-nodies whilst the location of the body depends on the force field. Mathematically this problem can be described with a differential equation of second order.

When solving such problems one usually starts with a known initial distribution of the bodies and known velocities and accelerations. From the initial conditions the system can be advanced in time by using numerical computations and a suitable integration scheme. Typical problems of this kind are for instance the computation of galaxy evolution in astrophysics or the behaviour of molecules in biology. Usually the equations of motion are relatively simple. For instance for simulating a galaxy they can be formulated using the Law of gravity:

Die Bewegungsgleichungen sind in aller Regel recht einfach aufgestellt. Wenn man beispielsweise eine Galaxie simulieren möchte, dann lautet die Gleichung um die Kraft auf ein Teilchen der Gesamtmenge zu berechnen:

\begin{equation} m_i \ddot{\vec{r_i}} = G \sum_{j=1;j \ne i}^N \Big( m_i m_j \frac{ \vec{r_i} - \vec{r_j}} {|\vec{r_i} - \vec{r_j}|^3} \Big) \end{equation}with:

- \(m_i, m_j\) - Mass of the particles i and j
- \(\vec{r_i}, \vec{r_j}\) - Positional vectors of the particles
- \(G\) - Gravitational constant

Lets assume F_{ij} is the force acting between particles i and j. The total number of force calculations needed to compute the state of the system is N*(N-1) and according to Newtons third law every force has
an opposite forces equal to itself: F_{ij} = -F_{ji}. The total number of force calculations can then be reduced to:

The problem is of order O(N^{2}). If the number of particles double the number of calculations quadrupels. If the number of particles is increased by factor ten the number of calculations increases by a factor of 100 and so on... From this simple relation it is clear that computing the N-Body problem for large numbers of particles will quickly become very costly in terms of computational power. A more effective algorithm that scales better with increasing number of particles is needed.

The Barnes-Hut Algorithm describes an effective methood for solving n-body problems. It was originally published in 1986 by Josh Barnes and Piet Hut ^{[1]}. Instead of directly summing up all forces it is using a tree based approximation scheme which reduces the
computational complexity of the problem from *O(N ^{2})* to

The Barnes-Hut algorithm reduces the number of force calculations by grouping particles. The basic idea behind the algorithm is that the force which a particle
group excerts on a single particle can be approximated by the force of a pseudo particle located at the groups center of mass. For instance the force which
the andromeda galaxy excerts on the milky way can be approximated by a point mass located at the centre of the andromeda galaxy. There is no need
to integrate over all stars in the andromeda galaxy provided the distance between the two galaxies is large enough. This approximation is valid as long
as the distance from the point group to the particle is large and the radius of the group is small in relation to the distance between
the group and the particle. The ratio of group distance (d) and group radius (r) is called the *Multipole-Acceptance-Criterion (MAC)*:

If theta drops below a certain threshold the quality of the approximation starts to deteriorate resulting in larger errors. Barnes and hut use a value of θ=1. The general rule is: The smaller θ, the better the simulation results. The Barnes-Hut algorithm applies this basic principle by subdividing the space of the simulation domain recursively (see Image 2). This is described more indepth in the next chapter.

The following section will explain the algorithm. All explanations and examples are simplified and deal with a two dimensional domain. Extending the concept to a third dimension is pretty straightforward but not the scope of this article.

Lets start with a random distribution of points in a two dimensional domain. The first step of the Barnes-Hut algorithm is sorting the particles into a hierarchical tree structure. This structure is the so called quadtree (or octree in 3D). The algorithm starts by subdividing the domain into quadrants. Based on their location I will refer to these quadrants by the following names: NW, NO, SW, SO. In order to add a particle to the tree its quadrant is determined. If there is another particle in the same quadrant the quadrant is subdivided again. This is repeated until each particle is located in its own quadrant (however small this quadrant may be). Image 3 shows the result of this process for a random distribution of particles.

The creation of the quadtree is a recursive process. It can be described by the following pseudocode:

Function MainApp::BuildTree { ResetTree for all particles rootNode->InsertToNode(particle) end for } Function TreeNode::InsertToNode(newParticle) { if number of particles in this node > 1 { quad = GetQuadrant(newParticle); if subnode(quad) does not exist create subnode(quad) subnode(quad)->InsertToNode(newParticle) } else if number of particles in this node == 1 { quad = GetQuadrant(existingParticle); if subnode(quad) does not exist create subnode(quad) subnode(quad)->InsertToNode(existingParticle) quad = GetQuadrant(newParticle); if subnode(quad) does not exist create subnode(quad) subnode(quad)->InsertToNode(newParticle); } else if node is empty { store newParticle in node as existingParticle } Increase number of particles }

The tree has to be recreated for each step of the simulation. The computational cost of recreating the quadtree depends on the distribution of the particles (see also [2]). The cost of adding a particle to the tree is proportional to its distance from the root node. Distributions with many densly packed particles require more operations because the tree must be subdivided often to place all particles in their own subdivided quadrant.

Up until now the tree only contains the spacial distribution of the particles. The next step deals with calculating the mass distribution of the quadtree (or. octtree). It sonsists of computing the total mass of all particles contained in the child nodes of each tree node as well as their center of mass. This data must be computed for each tree node. The calculation can be implemented recursively by scanning all child nodes of each treenode. The following pseudo code should illustrate the algorithm:

Function TreeNode::ComputeMassDistribution() { if number of particles in this node equals 1 { CenterOfMass = Particle.Position Mass = Particle.Mass } else { // Compute the center of mass based on the masses of // all child quadrants and the center of mass as the // center of mass of the child quadrants weightes with // their mass for all child quadrants that have particles in them { Quadrant.ComputeMassDistribution Mass += Quadrant.Mass CenterOfMass += Quadrant.Mass * Quadrant.CenterOfMass } CenterOfMass /= Mass } }

This procedure is of order O(N log N) ^{[2]}.

You probably guessed it by now the last step of the algorithm is computing the forces excerted on the particles. This is what makes the Barnes-Hut algorithm faster than the direct summation of all mutual forces. Do you remember the images 1 and 2 ? This is the basic principle we are going to implement here. In order to computer the force excerted on a single particle the quadtree must again be iterated recursively. For each node θ is computed:

$$\theta = {d \over r }$$ where:- d - Size of the box
- r - Particle distance from the nodes center of mass

If θ lies below a certain threshold the force can be approximated by inserting the quadrant nodes mass and its center of mass in the law of gravity. The child nodes don't have to be summed up separately. A reasonable relation is θ < 1. If θ is larger than the threshold the quadrants effect can't be approximated and all of its child nodes have to be tested again. The iteration stops only after all nodes have been tested. The worst case scenario is having to test all nodes without finding any node that meets the MAC. In such a case the result is similar to summing up all mutual particle forces (θ=0). The iteration depth can be finetuned by adjusting θ. Animation 1 illustrates the influence of θ on the number of force computations.

Finally the pseude code for the force calculation:

Function MainApp::CalcForce for all particles force = RootNode.CalculateForceFromTree(particle) end for end Function force = TreeNode::CalculateForce(targetParticle) force = 0 if number of particle equals 1 force = Gravitational force between targetParticle and particle else r = distance from nodes center of mass to targetParticle d = height of the node if (d/r < θ) force = Gravitational force between targetParticle and node else for all child nodes n force += n.CalculateForce(particle) end for end if end end

By combining the Barnes-Hut-Algorithm with a suitable integration scheme (i.e. Adams-Bashforth) one can relatively easy compute a two dimensional galaxy collision:

You might also like: