How can I evaluate the accuracy of my n-body simulation?

I am making an n-body simulation in python. There are many different methods to numerically solve the system of differential equations governing the gravitational interactions between the n particles; the challenge of n-body simulations is to find an integration scheme which is both accurate and fast.

For the moment, I am working with n small (say n ≤ 5). I began with the most naive method of Euler integration, which is first order. Even with two bodies, this method can produce pretty bad errors if the bodies come very close (since then the force between them becomes very large, so their accelerations become very large). Of course, one can always reduce the time-step to get more accuracy, but at the expense of computation time.

Next I would like to implement a higher order method which has better stability properties etc like leapfrog integration or Runge-Kutta.

I expect these methods to give much better results, but the problem is I am unsure how to measure quantitatively the accuracy of my simulation. What metric should I use?

One of the issues is that I can't tell if the simulation is giving realistic results, because I haven't been using the real value of the gravitational constant G and none of the distances/times/masses in my simulation have units.

The reason why I haven't been paying attention to units is because, say if I wanted to do a scale sun-earth-moon simulation; since the sun-earth distance is about 370 times the earth-moon distance, it's not possible for everything to fit inside the display window and also be able to distinguish the earth and the moon.

The other problem is using the real value of G and real masses etc leads to computations involving big numbers. But, keeping everything else fixed, choosing lower values of G reduces the errors (because the forces are then smaller). But when I reduce G, I don't know if I'm cheating by making it unrealistically small and getting good results, or if in fact my program is very accurate and it is the higher values of $G$ which are unrealistic. So I have no real way of telling if the simulation is accurate.

So how can I rigorously test the accuracy of the simulation, and what can I do concerning units?

• Just checking energy conservation over time will show you how badly Euler Integration is broken for physics simulation. The check with Computational Science for symplectic integrators (in fact that site is likely a better source of help than Physics). Once you've got that up and running the Virial theorem gives you another useful benchmark. – dmckee --- ex-moderator kitten May 12 '17 at 22:11
• @dmckee: Yes I read about energy drift being a useful metric for accuracy. But what I don't understand is isn't checking energy conservation going to slow down the simulation? If at each step I have to calculate the kinetic energies and all the potential energies, thats a lot of added computations... – Joshua Benabou May 12 '17 at 23:14
• You don't need to check conservation on every step. Even every hundred steps is probably excessive.. – dmckee --- ex-moderator kitten May 12 '17 at 23:49
• Against @dmckee's advice, computational science is not necessarily a good place to go for advice on this kind of problem. If you are planning to send a vehicle to Mars or Pluto, the symplectic integrators they will inevitably recommend are something you do not want to use. If you want to analyze the stability of the solar system over millions of years -- that also is not what you want to use. You want to use a technique that is highly accurate over the short term for the first problem, and for the latter, you'll want to use a geometric integrator (a generalization of symplectic integrators). – David Hammen May 12 '17 at 23:52
The first thing you should want to do in an orbital simulation is to toss $G$ and mass. Instead, use the product $GM$, or $\mu$ for short. This is observable; the gravitational parameter for our Sun is $1.32712440018\times10^{20}\pm9\times10^9\text{m}^3/s^2\$. Note the incredible precision compared to that of the gravitational constant. If you ever want to have a chance of non-embarrassing results on comparing your performance against the solar system. It's also easier. The acceleration of body #2 toward body #1 is $\frac{\mu_1}{||\vec r_1-\vec r_2||^3}(\vec r_1-\vec r_2)$, where $\mu_1$ is body $1's gravitational parameter. There are a number of things you can do to assess the quality of your simulation. One is to look at its performance in the two body problem. The two body problem is analytically solvable. You need to replicate this knowable behavior to some reasonable degree of accuracy over some reasonable span of time before you go onward. You will not be able to do this with Euler. Compare the worst-case performance of your simulation in terms of$\frac {||\vec r_{true} - \vec r_{sim}||}{r_{true}}$over the course of a smallish number of orbits. (I use 1, 3, 10, and 30 orbits for this test.) With regard to Euler, you'll need to reduce the step size so that you take millions of steps per orbit to have anything close to reasonable accuracy, and even then, your accuracy is shot after ten or so orbits. When playing with different techniques and different step sizes, you should eventually start seeing a general trend. Accuracy is poor across the board for extremely small step sizes (tens of million steps per orbit). The best technique won't fare much better than will Euler at these extremely small step sizes. This is the regime of machine errors. Which technique you use doesn't matter much in this regime of small step sizes. The errors that result from from using finite precision arithmetic dominate over the errors due to the integration techniques themselves in this small step size regime. Regardless of technique, accuracy in this regime should increase as you increase the step size. Use a log-log graph. The error should initially decrease linearly with increasing step size. With ever larger steps, you'll start seeing techniques break away from this decreasing error trend. The error instead increases with increasing step size increases, once again linearly on a log-log graph. Both the location of where the error starts increasing and the slope of that upward trend depend very much on the technique. With advanced techniques such as Gauss-Jackson or a geometrically correct Adams-Bashforth-Moulton integrator, you can make your step size a small fraction of an orbital period and see phenomenal accuracy. Some techniques claim success with step sizes larger than an orbit. Another test is to compare how well you do in various three body configurations. Create a two body system with one of the objects 1000 times more massive than the other, and then add a third body of negligible mass orbiting the larger at various locations. If you make the small object orbit elliptically very close to the large object, you should see it's periapsis precess. This precession is something that can be modeled via various approximation techniques. Does your simulation replicate those? (This, by the way, is a key win for general relativity. Mercury's Newtonian precession differs a bit from its observed precession.) Now make that small object orbit far from the two larger bodies. The key feature you should see here is that the orbital rate of that remote, small body is essentially that of a body orbiting a single body with the combined mass of the two inner bodies. These are short term tests. What if you want to see how your integrator fares over thousands of orbits, millions of orbits, or even more? Regardless of which technique you use, your accuracy is shot. However, your system should obey the conservation laws. Over the course of a few hundred orbits, how well does your n-body simulation conserve linear momentum, angular momentum, and energy? Some techniques that are highly accurate in the short term are rather lousy in the long term, and vice versa. • Ok this is a really great answer and very helpful. I will have some follow up questions tomorrow. But as for assessing the accuracy by looking at the special case$n=2\$, I already considered comparing the errors for integration schemes after the same number of orbits. Unfortunately I guess I have to count the number of orbits manually, because I can't think of any easy way to do this in the code. – Joshua Benabou May 12 '17 at 23:17