<blockquote>
<p>I used this n-body algorithm.</p>
</blockquote>
<p>Do not use that algorithm, repeated many times over. As a starter, the algorithm uses <a href="https://en.wikipedia.org/wiki/Semi-implicit_Euler_method" rel="noreferrer">symplectic Euler</a>. Paraphrasing from the movie &quot;Jaws&quot;, &quot;You're going to need a better integrator.&quot; 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.</p>
<p>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: <span class="math-container">$(a+b)+c = a+(b+c)$</span>. 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.</p>
<p>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</p>
<ul>
<li>A few million steps per orbit for a good second order integrator such as leapfrog,</li>
<li>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).</li>
<li>A few hundred steps per orbit for a high order Adams family integrator, particularly one that distinguishes between position and velocity.</li>
<li>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.</li>
</ul>
<p>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.</p>
<p>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.</p>
<p>To make matters worse yet (again), the parameter the source calls &quot;mass&quot; is not mass. It is instead a gravitational parameter, expressed in goofy units. The units of a gravitational parameter are length cubed / time squared.</p>
<p>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 <span class="math-container">$\text{au}^3/\,\text{day}^2$</span>. 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 <code>header.441</code> from <a href="https://ssd.jpl.nasa.gov/ftp/eph/planets/ascii/de441/" rel="noreferrer">https://ssd.jpl.nasa.gov/ftp/eph/planets/ascii/de441/</a>. 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 <code>GMx</code>, where <code>x</code> 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:</p>
<ul>
<li>GMS: 0.295912208284119561e-03 (Sun)</li>
<li>GM1: 0.491250019488931818e-10 (Mercury)</li>
<li>GM2: 0.724345233264411869e-09 (Venus)</li>
<li>GMB: 0.899701139294734660e-09 (Earth system)</li>
<li>GM4: 0.954954882972581189e-10 (Mars system)</li>
<li>GM5: 0.282534582522579175e-06 (Jupiter system)</li>
<li>GM6: 0.845970599337629027e-07 (Saturn system)</li>
<li>GM7: 0.129202656496823994e-07 (Uranus system)</li>
<li>GM8: 0.152435734788519386e-07 (Neptune system)</li>
<li>GM9: 0.217509646489335811e-11 (Pluto system)</li>
</ul>
<p>All values are in <span class="math-container">$\text{au}^3/\,\text{day}^2$</span>, which is exactly what you want for an n-body solar system model with length in astronomical units and time in days.</p>
<p>Bottom line:</p>
<ul>
<li>Use a better integrator.</li>
<li>Use a smaller timestep.</li>
<li>Use published values (published on-line in a reliable source counts as &quot;published&quot;).</li>
</ul>
<blockquote>
<p>Is using this data in a 64-bit float n-body simulation a fools errand?</p>
</blockquote>
<p>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.</p>
<p>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.</p>

CC BY-SA 4.0
2022-09-07T17:48:10.970Z
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?
Prallax
CC BY-SA 4.0
2022-09-07T22:43:57.367Z
@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.
David Hammen
CC BY-SA 4.0
2022-09-07T22:57:06.213Z
@Prallax For an "art project", at least switch to [leapfrog](https://en.wikipedia.org/wiki/Leapfrog_integration), [velocity verlet](https://en.wikipedia.org/wiki/Verlet_integration#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.
David Hammen
CC BY-SA 4.0
2022-09-09T16:12:13.067Z
@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 :-)
Michael Cole
CC BY-SA 4.0
2022-09-07T14:59:01.690Z
2022-09-07T15:08:39.563Z
2022-09-07T15:08:39.563Z
David Hammen
<p>I wanted to make a n-body simulation for an art project.</p>
<p>I used this <a href="https://github.com/lemire/ComputerLanguageBenchmark/blob/master/nbody/nbody.typescript-3.typescript" rel="nofollow noreferrer">n-body algorithm</a>, which includes a stable system for the Sun + gas giants, but no inner planets.</p>
<p>I wanted to simulate inner planets as well, so took <a href="https://ssd.jpl.nasa.gov/horizons/app.html#/" rel="nofollow noreferrer">JPL data from Horizons in AU and days</a> and plugged it into the system.  But, it doesn't create stable orbits.</p>
<p>Is using this data in a 64-bit float n-body simulation a fools errand?</p>
<p>Is there a better way to get this data for all planets + Pluto + Chiron?  Thanks!</p>

CC BY-SA 4.0
2022-09-06T02:47:43.563Z
What does 'not create stable orbits' mean exactly? Over what time?
planetmaker
CC BY-SA 4.0
2022-09-06T04:30:27.100Z
"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: https://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.
uhoh
CC BY-SA 4.0
2022-09-06T04:34:14.697Z
Or you could modify `advance()` and implement a "split-step reverse Euler" method, again from a cookbook of your choice.
uhoh
CC BY-SA 4.0
2022-09-06T08:25:27.273Z
@uhoh instead of ranting how one could dare use this algorithm (quite reasonable when one is not familiar with different integrators) and expect good results, you might turn this slightly expanded into a good answer.
planetmaker
CC BY-SA 4.0
2022-09-06T09:27:15.987Z
@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`)
uhoh
CC BY-SA 4.0
2022-09-06T09:30:05.333Z
@planetmaker but if I start offering advice on the various flavors of Euler and RK then David Hammen will scold me for not getting it right. I'm secretly hoping he'll see this question and offer some more specific numerical advice or add an answer. *The question is only 7 hours old!* Give it some time.
uhoh
CC BY-SA 4.0
2022-09-06T09:30:59.693Z
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.
PM 2Ring
CC BY-SA 4.0
2022-09-06T10:11:12.923Z
Also, DAYS_PER_YEAR probably should use the [Gaussian year](https://en.wikipedia.org/wiki/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.
PM 2Ring
CC BY-SA 4.0
2022-09-07T15:24:46.147Z
@planetmaker yay! [happy ending!](https://astronomy.stackexchange.com/a/50395/7982)
uhoh
CC BY-SA 4.0
2022-09-09T16:13:16.107Z
Hey guys, thanks for the help.  I love how the astronomy community is small enough yall know each other.  Have an awesome day!
Michael Cole
<n-body-simulations><ephemeris>
How to get JPL Ephemeris data to work in n-body simulation or why is Mercury flying to Jupiter
Michael Cole
# How to get JPL Ephemeris data to work in n-body simulation or why is Mercury flying to Jupiter

(Asked by: username_0 on 2022-09-06)

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

I used this [n-body algorithm](https://github.com/lemire/ComputerLanguageBenchmark/blob/master/nbody/nbody.typescript-3.typescript), 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](https://ssd.jpl.nasa.gov/horizons/app.html#/) 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!

**Tags:** n-body-simulations, ephemeris
> **Comments:**
> **username_1**: What does 'not create stable orbits' mean exactly? Over what time?  
> **username_2**: "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: https://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.  
> **username_2**: Or you could modify `advance()` and implement a "split-step reverse Euler" method, again from a cookbook of your choice.  
> **username_1**: @uhoh instead of ranting how one could dare use this algorithm (quite reasonable when one is not familiar with different integrators) and expect good results, you might turn this slightly expanded into a good answer.  
> **username_2**: @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`)  
> **username_2**: @planetmaker but if I start offering advice on the various flavors of Euler and RK then David Hammen will scold me for not getting it right. I'm secretly hoping he'll see this question and offer some more specific numerical advice or add an answer. *The question is only 7 hours old!* Give it some time.  
> **username_3**: 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.  
> **username_3**: Also, DAYS_PER_YEAR probably should use the [Gaussian year](https://en.wikipedia.org/wiki/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.  
> **username_2**: @planetmaker yay! [happy ending!](https://astronomy.stackexchange.com/a/50395/7982)  
> **username_0**: Hey guys, thanks for the help.  I love how the astronomy community is small enough yall know each other.  Have an awesome day!  

---

## 1 answers

### Answer by username_4 on 2022-09-07. Score: 7

> I used this n-body algorithm.

Do not use that algorithm, repeated many times over. As a starter, the algorithm uses [symplectic Euler](https://en.wikipedia.org/wiki/Semi-implicit_Euler_method). 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.

> **Comments:**
> **username_5**: 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?  
> **username_4**: @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.  
> **username_4**: @Prallax For an "art project", at least switch to [leapfrog](https://en.wikipedia.org/wiki/Leapfrog_integration), [velocity verlet](https://en.wikipedia.org/wiki/Verlet_integration#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.  
> **username_0**: @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 :-)  

---

