# Which timestep should I use for a $N$-body simulation of the Solar system?

I am trying to implement a $N$-body simulation of the Solar system and I am stuck on the issue of the simulation timestep size. At first I would like to use a simple Euler stepping scheme. Knowing that the error is dependent on the size of the timestep, I am not sure what to use for the model.

Would using a timestep size of one day be enough to get valid results? And also I would need to input the timestep in the model using seconds (60*60*24), correct?

Finally, where can I find values for the initial conditions that I should input in the model (initial positions & velocities vectors).

Sorry if these questions are kinda basic but my physics is about 15 years rusty... :)

• The time step is not the only source of error. The integration scheme has inherent problems. The simplest Euler method does not conserve energy, and will present errors no matter how small the time step. (This depends on you tolerance for errors, of course). There are better integrators. You might search for "symplectic integrator" to get started. – garyp May 28 '16 at 16:19
• Would Computational Science be a better home for this question? – Qmechanic May 28 '16 at 20:55
• @Qmechanic I think my question is about computational physics. Computational Science is too broad in that sense. – BigONotation May 29 '16 at 13:55
• @garyp Yes I am aware Euler has issues :) But I just want to get a "quick and dirty" model working. Then I will switch to something more sophisticated such as Runge-Kutta 4. Will also definitely look into the "symplectic integrator", thanks for the pointer! I am using C++ as an implementation language, so the idea is to make the different integration algorithms "switchable". – BigONotation May 29 '16 at 21:17

The standard way to choose a time-step is to run a test simulation with multiple bodies and plot the total energy of the system versus time. The total energy should remain (roughly) constant. If your step size is too large then you will get energy drift. So simply find the largest time-step that does not produce energy drift. In the case of modelling the solar system, my hunch is that a time-step of 1 day would be too large. I'd suggest trialling something closer to 1 hour.

By the way, more advanced codes tend to use variable time-steps, e.g..

And also I would need to input the timestep in the model using seconds (60*60*24), correct?

That entirely depends on the units that the code uses. It's generally regarded as bad practice to use SI units though for modelling astronomical systems due to the huge numbers involved.

You may find the Wiki article on Numerical model of the Solar System to be informative.

• thanks a lot for all of the clarifications and pointers. This is very helpful! There is still one point I am not sure about: how should I express the problem using the astronomical system of units? From my understanding, this would be about scaling all quantities. For example the Sun would have a mass of 1, Jupiter 1/1048, and Earth 1/332950. Same thing for distances, which would be expressed in astronomical units (so distance Sun-Earth = 1 AU). For time, I would use the Day. So the 1 hour time-step would be 1/24 Day. Otherwise I would apply Newton's gravitational law as is. Correct? – BigONotation May 28 '16 at 17:39
• For the initial conditions, here is the answer given in a linked question astronomy.stackexchange.com/questions/2416/… – BigONotation May 28 '16 at 17:49
• @BigLudinski Choose any units you like, and then work out what the gravitational constant $G$ is in those particular units. As long as you're consistent then that's all there is to it. – lemon May 28 '16 at 19:23
• thanks a lot for all of your feedback! Also found this great answer regarding scaling the SI units: astronomy.stackexchange.com/questions/11879/… . So in essence if I use the astronomical system of units (time = Day, mass = Solar mass, length = AU), the value of G = 2.95912208286e-4 – BigONotation May 29 '16 at 13:35

Direct integration schemes give bad results. You can do much better by using the exact solution in the absence of gravitational interactions between the planets. You can then set up a variation of constants approach where you take the integration constants (which are the orbital parameters) as dynamical variables and write the differential equations (where you now include the mutual gravitational attraction between the planets) in terms of these variables.

Because you now only have slowly varying functions to deal with, you can take much larger time steps. The downside is that the math is much more complicated involving functions that need to be evaluated using numerical integrations.