How to get JPL Ephemeris data to work in n-body simulation or why is Mercury flying to Jupiter

Question

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!

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.