3

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?

math_lover
  • 5
  • 1
  • 6
  • 2
    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: https://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., https://scicomp.stackexchange.com/questions/206) – Kirill May 12 '17 at 21:12
  • 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. – math_lover May 12 '17 at 21:16
  • 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. – Kirill May 12 '17 at 21:22
  • 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? – math_lover May 12 '17 at 21:50
  • Ok. According to some papers, measuring the energy drift of the simulation is a good means of measuring the acccuracy. The problem is the papers which I have seen are not clear on how the total energy is calculated at each step. Do they mean, just compute the sum of the kinetic energies and the graviational potential energies of the system at each step (wouldn't this slow down the code alot?) – math_lover May 12 '17 at 22:06
  • 4
    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. – Richard Zhang May 13 '17 at 16:44
  • @RichardZhang : indeed. This is helpful, thank you. I will use this as a metric to test the accuracy of the various integrators, along with energy drift. – math_lover May 13 '17 at 23:09

0 Answers0