Simulating Gravity at faster than N2

CS 482 Lecture, Dr. Lawlor

In a particle simulation, it's pretty easy to compute gravitational interactions directly from Newton's Law of Universal Gravitation: f = g m1 m2 / r2

g is the gravitational constant
m1 and m2 are the masses of the objects
r is the distance between the centers of the objects
f is the magnitude of the force

Since you want a vector force, you need to scale this magnitude f by a vector, which can be the normalized direction vector R between the object centers.  Or you can save one vector scale factor by using the non-normalize vector between centers R, and calculating F = R * g m1 m2 / r3, where the cube divides out the length of R during the rest of the scale factor.

The other problem in a simulation is particles will occasionally pass very close to one another, resulting in a distance r near zero, which gives a force near infinity.  To avoid this, the force is typically "softened" by dividing by 1/(a+r2), where a is a positive constant used to trim the maximum value.  This is equivalent to smoothing the point masses into small clouds.

Faster Neighbor Finding via Hashtable

One fairly simple way to get much better performance is to only consider nearby neighbors.  You can find your nearby neighbors by storing lists of particles by their integer coordinates in a hash table, and looking through your list to find your closest neighbors.  But there's a problem: the galaxy can't hold together anymore, even if we artificially inflate the gravitational constant, because the long-range forces that bind one side of the galaxy to the other are never evaluated.

Approximating Long-Range Interactions

The hashtable version gives the wrong answer because it ignores long-range interactions.  There are several neat ways to approximate these interactions.

A very simple approach is to sum all the particles into one "megaparticle", which will end up at the center of the galaxy.  We can then compute each particle's interaction with its local neighbors, and then with the megaparticle.  This doesn't even conserve momentum, because the influence of a particle on the megaparticle isn't symmetric, but at least the galaxy holds together!

A more accurate version of this is called Barnes-Hut: for neighboring particles, we compute full particle-particle interactions.  To interact with a group of distant particles, we lump the group into a megaparticle.  To accurately define which particles form a "group", we impose a spatial tree across the simulation, using the coarse levels of the tree to approximate distant regions, and increasingly fine levels for more nearby regions.  See my Barnes-Hut Gravity Tree code.

Here's some measured performance of gravity via Barnes-Hut:

n Particles Tree Build Lump Tree Gravity Export
20 0.3 0.015 0.31 3.2
200 1.2 0.03 5.4 23.5
2000 23.6 0.3 307.7 227.3
Scaling
 O(n) + big constant?
 O(n)
 O(n log n)
 O(n) with huge constant