I'm trying to solve a 5-body problem (for the sake of the argument it can be any given $n$) numerically, yet every initial condition I choose seems to eventually drive the bodies apart to infinity or make them collide altogether.

I want them to move somewhat "periodically", meaning to neither vanish to infinity nor to collide.

I understand that the total energy has to be negative, but what else can I try to make it happen, other than making numerous attempts, each time trying to change the initial conditions a little bit?


I would choose highly symmetric initial conditions which are easily shown to be periodic. Then I would perturb the system by small steps.

For example, four identical particles (of mass $m$) moving on a circle and interacting through Newton gravitation alone is a solution of the equations of motion.

$$ \vec{x}_i = R \left(\begin{array}{c}\cos\left(\omega t +\delta_i\right) \\ \sin\left(\omega t + \delta_i\right) \end{array}\right) \, , \, \delta_i = (i-1)\frac{\pi}{2} \, ,$$

($i = 1,2,3,4$) solves Newton equations if

$$ \omega^2 = \frac{\sqrt{2}+4}{4 \sqrt{2}} \frac{G m}{R^3} \, . $$

You can use this as initial conditions and start by checking your numerical code. Then you can slightly change $\omega$, $m$, $\delta_i$, etc.

I did this for $4$ particles because it's easier. I leave the $5$ particles case to you.

  • $\begingroup$ Can you please explain how you found the frequency? $\endgroup$ – Talco Aug 18 '15 at 16:10
  • 1
    $\begingroup$ $\omega$ is found by writing down Newton's equations $\vec{F}_i = m \vec{a}_i$. Since each particle is moving in a circle, the right-hand side is directed towards the centre of the circle and is $m v^2 /R = m R \omega^2$. The left-hand side is found by adding up the forces exerted by the three other particles. Because of the particular choice of $\delta_i$, these are distributed symmetrically and this is easy. Then you get an equation relating $R$, $G$ and $\omega$ which yields the relation that I give. $\endgroup$ – Steven Mathey Aug 18 '15 at 16:34

Physically, a negative total energy is a necessary and sufficient condition for avoiding all particles flying off to infinity separately.

However, it is always possible for some particles to be given enough energy to escape a system. In fact, this tends to happen in real life. Planets get ejected from their solar systems over millions of years, and stars "evaporate" from globular clusters. The escaping bodies get energy from the rest of the system. Note that two Newtonian point masses can have an arbitrarily negative total energy as their orbital separation shrinks, so one can in principle extract arbitrary energy from any pair of particles to give to other particles.

In my experience, though, this sort of behavior is almost certainly indicative of a bug. The three most common pitfalls are:

  1. Using the wrong integration algorithm. If you're using an unstable scheme like forward-Euler, it is expected that everything will go wrong. People who know better than to use forward-Euler will still often recommend 4th-order Runge-Kutta. This is better, but still a bad choice for long-term integrations, since it doesn't conserve total energy in the system. Leapfrog is a good choice to start with.
  2. Using too large a timestep. How exactly are you calculating the timestep? All too often people choose a number that looks small, without answering the necessary question of "small compared to what?"
  3. Not softening the potential. Remember from above how point particles can get arbitrarily close, seeing arbitrarily rapidly changing potentials and undergoing arbitrarily large velocities and accelerations? Arbitrarily large values are bad for simulations. At a minimum, timesteps should decrease as needed when particles get close to each other. To avoid timesteps going to zero, "softened" potentials are often employed. For example, change Newton's law of gravity to be $F = G m_1 m_2 / (r^2 + \epsilon^2)$ for some small length scale $\epsilon$ of your choosing.

For diagnosing the problem, try plotting the particles' positions simultaneously with the system's total energy. If you see the energy change every time two particles get close, you are suffering from hard potentials and large time steps. If you simply see a steady trend in the energy, you probably aren't using a symplectic integrator.

As for the initial conditions, again, everything should be fine with just a negative total energy. For small numbers of objects, you can try Steven Mathey's suggestion of particles equally spaced on a circle. Note this is physically unstable, so you shouldn't expect the formation to hold indefinitely.1

Once you have enough particles (at least 20, or better yet at least 100), you can leverage statistical averaging. Choose a velocity dispersion $\sigma$. Randomly place particles to sample the density distribution $\rho = \sigma^2 / 2\pi G r^2$. Give them velocities drawn from a Maxwell-Boltzmann distribution ($\sigma^2 = kT/m$), which will be isotropically distributed and independent of position. You have now constructed a singular isothermal sphere.

1In fact, no simulation should be expected to work arbitrarily far into the future, and you should always be aware of the limits of any simulations. Symplectic integrators preserve energy and momentum for long times, but they quickly scramble the internal state (you wouldn't necessarily use one to predict the location of Mars 10 years from now). Higher-order non-symplectic methods are often better at short-term accuracy (but unreliable for determining things like whether the solar system is stable on billion-year timescales).


Not the answer you're looking for? Browse other questions tagged or ask your own question.