I wanted to make a n-body simulation for an art project.

I used this n-body algorithm, which includes a stable system for the Sun + gas giants, but no inner planets.

I wanted to simulate inner planets as well, so took JPL data from Horizons in AU and days and plugged it into the system. But, it doesn't create stable orbits.

Is using this data in a 64-bit float n-body simulation a fools errand?

Is there a better way to get this data for all planets + Pluto + Chiron? Thanks!

  • 1
    $\begingroup$ "Is using this data in a 64-bit float n-body simulation a fools errand?" No, it's an excellent idea! But using a random bit of code from a benchmarking test which implements a really simple numerical integrator is! If you keep making your step sizes smaller and smaller, eventually you might get something reasonable, but just replace advance(dt: number) with something more reasonable. Here's a bit of Python: space.stackexchange.com/a/23409/12102 I used a canned integrator but you can look up "RK4(5)" with variable step size in any numerical cookbook for your favorite language. $\endgroup$
    – uhoh
    Sep 6, 2022 at 4:30
  • 1
    $\begingroup$ @planetmaker I disagree with your characterization of my helpful advice as a "rant", but here's my thinking. Based on the OP's profile in Stack Overflow and other software sites they're quite capable of grabbing a few lines of code from a numerical cookbook and implementing it as a function in the program they have already adopted, which they might find even fun. Ideally they'll get it working, then post their solutions as answer & pick up rep. In all my years I've only been able to learn and remember three languages; english, math, and python (and maybe a little FORTRAN, e.g. GOSUB 100) $\endgroup$
    – uhoh
    Sep 6, 2022 at 9:27
  • 2
    $\begingroup$ A better integrator will certainly help. But I suspect you might have an error with your velocity conversion. Make sure your inner planet velocities are in AU / year. $\endgroup$
    – PM 2Ring
    Sep 6, 2022 at 9:30
  • 1
    $\begingroup$ Also, DAYS_PER_YEAR probably should use the Gaussian year, but that's just a minor thing and not the cause of Mercury's escape. You definitely need a smaller time step: 0.01 years is ok for the outer planets, but not for Mercury. Using JPL data, the Gaussian year is ~365.2568983272363971 days. $\endgroup$
    – PM 2Ring
    Sep 6, 2022 at 10:11
  • 1
    $\begingroup$ Hey guys, thanks for the help. I love how the astronomy community is small enough yall know each other. Have an awesome day! $\endgroup$ Sep 9, 2022 at 16:13

1 Answer 1


I used this n-body algorithm.

Do not use that algorithm, repeated many times over. As a starter, the algorithm uses symplectic Euler. Paraphrasing from the movie "Jaws", "You're going to need a better integrator." Better integrators are more complex, but they allow you to take larger step sizes. 100 steps per orbit (e.g., the Earth in the test code) is far too coarse for integrators that are far better than symplectic Euler, and it is beyond too coarse for symplectic Euler. There are some integrators where 100 steps per orbit is okay. Symplectic Euler does not fall into that class.

To make matters worse, lines 51-57 are subject to a nasty problem with using typical computer representations of the real numbers. Addition is associative in the real numbers: $(a+b)+c = a+(b+c)$. That is not true with typical representations of the real numbers such as those in javascript, python, whatever. You should accumulate the accelerations (or acceleration * dt) before adding that delta v to the velocity.

To make matters even worse, the optimal time step for Mercury is very different from the optimal time step for Pluto. Ideally, you'll use different step sizes for each body. The optimal step size depends on the integrator, the floating point precision, and the magnitude of the perturbations. For a perfectly circular, unperturbed orbit I found that the optimal step size using 64 bit floating point is about

  • A few million steps per orbit for a good second order integrator such as leapfrog,
  • About five thousand steps per orbit for a good fourth order integrator such as RK4 (but note that RK4 does not conserve angular momentum or energy).
  • A few hundred steps per orbit for a high order Adams family integrator, particularly one that distinguishes between position and velocity.
  • A step per orbit for an extremely high order Runge-Kutta integrator. However, that very high order Runge-Kutta integrator will inevitably make several hundred derivative calls per step.

I couldn't find the optimal step size for symplectic Euler. It's in the hundreds of millions (or more) of steps per orbit. Don't use symplectic Euler.

To make matters worse yet, the initial state (position and velocity) in that source are goofy. Use reliable sources (Horizons is a reliable source), all taken at the same time. Calculate the position and velocity relative to the solar system barycenter, which can be done with Horizons. Note: You should have Horizons calculate the position and velocity of the Earth-Moon barycenter rather than the Earth.

To make matters worse yet (again), the parameter the source calls "mass" is not mass. It is instead a gravitational parameter, expressed in goofy units. The units of a gravitational parameter are length cubed / time squared.

You want to use astronomical units as the unit length of time and days (86400 seconds, no leap seconds please) as the unit of time, which means you want gravitational parameters in units of $\text{au}^3/\,\text{day}^2$. These are exactly the units that JPL uses internally to calculate their Development Ephemerides. Horizons currently uses DE441. You can get the underlying data from JPL by downloading header.441 from https://ssd.jpl.nasa.gov/ftp/eph/planets/ascii/de441/. You are looking for GROUP 1040 (a list of names) and GROUP 1041 (a list of values). You want the values for the parameters named GMx, where x is S for the Sun, 1 for Mercury, 2 for Venus, B for the Earth+Moon, 4 for Mars, etc. I've done this for you below:

  • GMS: 0.295912208284119561e-03 (Sun)
  • GM1: 0.491250019488931818e-10 (Mercury)
  • GM2: 0.724345233264411869e-09 (Venus)
  • GMB: 0.899701139294734660e-09 (Earth system)
  • GM4: 0.954954882972581189e-10 (Mars system)
  • GM5: 0.282534582522579175e-06 (Jupiter system)
  • GM6: 0.845970599337629027e-07 (Saturn system)
  • GM7: 0.129202656496823994e-07 (Uranus system)
  • GM8: 0.152435734788519386e-07 (Neptune system)
  • GM9: 0.217509646489335811e-11 (Pluto system)

All values are in $\text{au}^3/\,\text{day}^2$, which is exactly what you want for an n-body solar system model with length in astronomical units and time in days.

Bottom line:

  • Use a better integrator.
  • Use a smaller timestep.
  • Use published values (published on-line in a reliable source counts as "published").

Is using this data in a 64-bit float n-body simulation a fools errand?

No. A 64-bit floating point n-body simulation should be good enough for your purposes for a few hundred years, maybe even a few thousand years -- if you use a good integrator.

That said, JPL uses 128 bit floating point numbers with an adaptive step size, adaptive order Adams-style integrator to generate their Development Ephemerides. They also use a lot of corrections (general relatively, oblateness, ...) that are well beyond what you want to do.

  • 1
    $\begingroup$ Maybe an overkill for an art project, but otherwise, great answer! Just a question: how do you define the optimal time step? I mean, optimal is always relative to a goal. Which is the definition in this case? $\endgroup$
    – Prallax
    Sep 7, 2022 at 17:48
  • 1
    $\begingroup$ @Prallax I would define optimal time step as the time step that minimizes error. There are two primary sources of error: the integration technique itself and how the real numbers are represented. Errors from the integration technique dominate when the time step is overly large. Decreasing the step size in this regime reduces the error. Errors from the representation dominate when the time step is overly small. Increasing the step size in this regime reduces the error. This means there's a crossover point where the dominant error source switches. This crossover point is the optimal step size. $\endgroup$ Sep 7, 2022 at 22:43
  • 1
    $\begingroup$ @Prallax For an "art project", at least switch to leapfrog, velocity verlet, or something similar (these are easy to implement), and reduce the step size. 100 steps per orbit will result in garbage for all but the very best of integrators. For an "art project", a hundred thousand steps per Mercury's orbital period might suffice. Never use symplectic Euler. $\endgroup$ Sep 7, 2022 at 22:57
  • $\begingroup$ @DavidHammen Yes! Thanks for the help. Yeah, changing to Leapfrog or similar will work. I figured it was because the steps were too "big" for the inner planets, but I want the gas giants to show some hustle orbiting the sun :-) $\endgroup$ Sep 9, 2022 at 16:12

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .