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?