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 \le 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). 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?

  • 2
    $\begingroup$ First, your question isn't really specific to n-body problems, it's more like a general question of how to do verification and validation. E.g., there's this question already: scicomp.stackexchange.com/questions/3451. Second, why would it matter that the model doesn't neatly fit on one display? You should automate your tests and compare the results to a reference solution, not look at them to see if they look roughly right (e.g., scicomp.stackexchange.com/questions/206) $\endgroup$
    – Kirill
    May 12, 2017 at 21:12
  • $\begingroup$ Well my simulation uses a visual display so I can see real time what is happening. My end goal is to search for periodic solutions to the three body problem but for now I'm just trying to obtain an accurate and fast simulation. It matters that everything fits on the screen because then I can see real time what is happening in the simulation, because the time of the simulation is displayed on the screen. Also there are some specificities here to the n-body problem, namely the choice of units. $\endgroup$
    – math_lover
    May 12, 2017 at 21:16
  • $\begingroup$ What I'm saying is that inspecting the output by hand is an unreasonably hard way to ensure correctness and accuracy (in general, not just n-body). You could just partially rewrite it to automate it and so make the task easier for yourself. $\endgroup$
    – Kirill
    May 12, 2017 at 21:22
  • $\begingroup$ Well I'm not planning on doing it by hand. I want to place some measurement on the display screen which updates real time. For example I took the simple case of a planet orbiting in a circle around a much more massive sun which essentially remains fixed in space and updated the error between the distance between the two bodies on the screen and the theoretically constant radius. I could then visualize the error vs time and if I wanted to I could plot it on a graph. My question is what general metrics should I use to evaluate accuracy? $\endgroup$
    – math_lover
    May 12, 2017 at 21:50
  • 4
    $\begingroup$ I believe the n-body problem has analytical solutions for the 2-body case, and even for certain special instances of the 3-body cases. It may be helpful for you to simulate these numerically, and to compare your results against the analytical solution. You should observe, for example, that the numerical converges towards the analytic as the timesteps are tightened. $\endgroup$ May 13, 2017 at 16:44


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.