<p>I can reproduce your JPL Z coordinate plot if I use Earth's heliocentric position relative to the <em>J2000</em> ecliptic.
Naturally the amplitude has a minimum around the year 2000:</p>

<p><img src="https://i.sstatic.net/JxxqZ.png" alt="EMB.z, J2000 ecliptic, vs. time"></p>

<p>Relative to the ecliptic <em>of date</em>, which accounts for precession, the heliocentric Z coordinate of the Earth-Moon barycenter is much smaller:</p>

<p><img src="https://i.sstatic.net/7EgUn.png" alt="EMB.z, ecliptic of date, vs. time"></p>

<p>The IAU 2006 precession model has two components, detailed in <a href="https://ui.adsabs.harvard.edu/abs/2003A%26A...412..567C/abstract" rel="nofollow noreferrer">Capitaine et al. 2003</a>:</p>

<ul>
<li>precession of the <em>ecliptic</em>, a 47"/century change in the plane of Earth's orbit</li>
<li>precession of the <em>equator</em>, a 1.4&deg;/century change in the plane and axis of Earth's rotation</li>
</ul>

<p>Fifty years of ecliptic precession amount to 23.5" or 0.000114 radian, and the first plot is consistent with that.
The classical precession of the equinoxes (the intersection of the two planes) is a combination of both components.</p>

<p><a href="https://ui.adsabs.harvard.edu/abs/1994AJ....108..711W/abstract" rel="nofollow noreferrer">Williams 1994</a> examines Earth's equatorial precession (dominated by the Sun and Moon) and finds that Venus's effect is larger than Jupiter's.
When looking at <a href="/a/34752">fluctuations in Earth's year length</a>, I noticed a cycle similar to Venus's synodic period.
I expect that adding Venus to your system would reproduce some ecliptic precession.</p>

CC BY-SA 4.0
2020-01-31T23:14:20.850Z
Thanks for the edit!
uhoh
CC BY-SA 4.0
2020-02-02T00:26:26.577Z
Thanks for the information! I had tried adding Venus before I rewrote my code to use AU instead of KM and it didn't make much difference then since it was quite inaccurate. Adding it now has made it much similar! I didn't realise the degree to which Venus has!
jameslfc19
CC BY-SA 4.0
2020-02-03T06:19:32.520Z
@jameslfc19 Calculate $m/r^2$ for the distance of closest approach and Venus turns out to be almost as "strong" as Jupiter. Then factoring in the higher frequency of its perturbation and the similarity with Earth's period.
uhoh
CC BY-SA 4.0
2020-01-31T08:44:14.200Z
2020-01-31T22:58:28.767Z
2020-01-31T22:58:28.767Z
Mike G
<p>I'm currently building an N-body simulator for my masters project and I'm using JPL's Horizons system to compare against.</p>

<p>I started with just a couple orbits and everything seemed to be working nicely (After fixing many mistakes implementing a Fourth Order Integrator) but after extending to 10+ orbits I notice that in the Z-axis the Earth's amplitude increases linearly, according to JPL, but my simulation has a constant amplitude. JPL's result is shown in red in the bottom left:</p>

<p><a href="https://i.sstatic.net/qcFwA.png" rel="noreferrer"><img src="https://i.sstatic.net/qcFwA.png" alt=""></a>  </p>

<p>What is causing this? Currently I only have the Sun, Earth, Moon and Jupiter in the system, but I didn't even realise the amplitude increased like this? What am I missing or got wrong? </p>

<p>(Also, the graph only shows 3 simulated orbits, and 80 JPL orbits)</p>

CC BY-SA 4.0
2020-01-29T20:09:20.013Z
What year is day 0?
Mike G
CC BY-SA 4.0
2020-01-29T21:57:18.940Z
@MikeG 2019-10-15. I think this is due to Orbital Inclination, I added Saturn to the simulation and it's closer now, but it's still slightly off using RK4 and a time step of 0.00001 days. Not sure if this is just the accuracy of the method or if there is something else. https://i.imgur.com/VaIjHUV.png
jameslfc19
CC BY-SA 4.0
2020-01-30T01:03:00.363Z
The Sun has a weak J2 but I don't know if that's it. See for example answers to [How to calculate the planets and moons beyond Newtons's gravitational force?](https://space.stackexchange.com/q/23408/12102)
uhoh
CC BY-SA 4.0
2020-02-01T02:37:14.293Z
BTW, although you can use RK4 in long timespan celestial mechanics sims, the errors can lead to non-conservation of energy (which causes orbits to spiral in or out unless you use a really tiny timestep). The cure for this is to use a [symplectic integrator](https://en.wikipedia.org/wiki/Symplectic_integrator). I quite like synchronized [Leapfrog integration](https://en.wikipedia.org/wiki/Leapfrog_integration); higher order integrators can be constructed using Yoshida coefficients. There's a nice discussion of this in http://www.artcompsci.org/kali/vol/two_body_problem_2/ch07.html#rdocsect46
PM 2Ring
CC BY-SA 4.0
2020-02-02T00:21:40.380Z
@PM2Ring Thanks! The Leapfrog Integrator works well!
jameslfc19
CC BY-SA 4.0
2020-02-02T04:49:55.023Z
Assuming you are using double precision rather than infinite precision arithmetic, your time step of 0.00001 days is far too small of a value for RK4, by a factor of about 5000 for Earth. Unfortunately, writing on ODE solvers make it appear that the only factors leading to error are technique order and step size. This is not the case. You might as well be using a second order technique with that tiny step size.
David Hammen
CC BY-SA 4.0
2020-02-02T04:53:48.600Z
@PM2Ring Unless you're using ridiculously small step sizes, leapfrog is rather awful compared to higher order techniques. Of course, with ridiculously small step sizes, every technique is rather awful. JPL, for example, does not use anything close to leapfrog to generate its ephemerides. It instead uses a non-symplectic Adams type integrator, and takes rather large steps.
David Hammen
CC BY-SA 4.0
2020-02-02T05:26:23.507Z
@David Ok, maybe plain Leapfrog isn't fantastic, but have you tried the 4th order or 8th order versions, using Yoshida coefficients? You might be surprised...
PM 2Ring
CC BY-SA 4.0
2020-02-02T05:49:25.037Z
@David I must confess that I'm somewhat surprised that JPL don't use a symplectic integrator. (True, symplectic integrators don't *exactly* conserve the energy, they conserve a closely related "shadow" Hamiltonian). OTOH, I have to concede that John Couch Adams has a good reputation with numerical integration, having predicted Neptune. ;)
PM 2Ring
CC BY-SA 4.0
2020-02-03T18:32:08.653Z
@DavidHammen Yeah turns out my RK4 implementation was considerably wrong (I was using the acceleration calculation from the beginning of the step for each intermediate calculation) I can now use time steps greater than 0.01 and achieve accurate results! Also I implemented the 4th order Yoshida Integrator and that produces almost exactly the same as RK4 from a small sample of runs. Next I need a way to benchmark each integrator at different time steps, any pointers?
jameslfc19
<orbit><solar-system><n-body-simulations>
Peculiar Amplitude of Earth's Orbit in Z Axis
jameslfc19
# Peculiar Amplitude of Earth's Orbit in Z Axis

(Asked by: username_0 on 2020-01-29)

I'm currently building an N-body simulator for my masters project and I'm using JPL's Horizons system to compare against.

I started with just a couple orbits and everything seemed to be working nicely (After fixing many mistakes implementing a Fourth Order Integrator) but after extending to 10+ orbits I notice that in the Z-axis the Earth's amplitude increases linearly, according to JPL, but my simulation has a constant amplitude. JPL's result is shown in red in the bottom left:

[![](https://i.sstatic.net/qcFwA.png)](https://i.sstatic.net/qcFwA.png)

What is causing this? Currently I only have the Sun, Earth, Moon and Jupiter in the system, but I didn't even realise the amplitude increased like this? What am I missing or got wrong?

(Also, the graph only shows 3 simulated orbits, and 80 JPL orbits)

**Tags:** orbit, solar-system, n-body-simulations

> **Comments:**
> **username_1**: What year is day 0?  
> **username_0**: @username_1 2019-10-15. I think this is due to Orbital Inclination, I added Saturn to the simulation and it's closer now, but it's still slightly off using RK4 and a time step of 0.00001 days. Not sure if this is just the accuracy of the method or if there is something else. https://i.imgur.com/VaIjHUV.png  
> **username_2**: The Sun has a weak J2 but I don't know if that's it. See for example answers to [How to calculate the planets and moons beyond Newtons's gravitational force?](https://space.stackexchange.com/q/23408/12102)  
> **username_3**: BTW, although you can use RK4 in long timespan celestial mechanics sims, the errors can lead to non-conservation of energy (which causes orbits to spiral in or out unless you use a really tiny timestep). The cure for this is to use a [symplectic integrator](https://en.wikipedia.org/wiki/Symplectic_integrator). I quite like synchronized [Leapfrog integration](https://en.wikipedia.org/wiki/Leapfrog_integration); higher order integrators can be constructed using Yoshida coefficients. There's a nice discussion of this in http://www.artcompsci.org/kali/vol/two_body_problem_2/ch07.html#rdocsect46  
> **username_0**: @username_3 Thanks! The Leapfrog Integrator works well!  
> **username_4**: Assuming you are using double precision rather than infinite precision arithmetic, your time step of 0.00001 days is far too small of a value for RK4, by a factor of about 5000 for Earth. Unfortunately, writing on ODE solvers make it appear that the only factors leading to error are technique order and step size. This is not the case. You might as well be using a second order technique with that tiny step size.  
> **username_4**: @username_3 Unless you're using ridiculously small step sizes, leapfrog is rather awful compared to higher order techniques. Of course, with ridiculously small step sizes, every technique is rather awful. JPL, for example, does not use anything close to leapfrog to generate its ephemerides. It instead uses a non-symplectic Adams type integrator, and takes rather large steps.  
> **username_3**: @username_4 Ok, maybe plain Leapfrog isn't fantastic, but have you tried the 4th order or 8th order versions, using Yoshida coefficients? You might be surprised...  
> **username_3**: @username_4 I must confess that I'm somewhat surprised that JPL don't use a symplectic integrator. (True, symplectic integrators don't *exactly* conserve the energy, they conserve a closely related "shadow" Hamiltonian). OTOH, I have to concede that John Couch Adams has a good reputation with numerical integration, having predicted Neptune. ;)  
> **username_0**: @username_4 Yeah turns out my RK4 implementation was considerably wrong (I was using the acceleration calculation from the beginning of the step for each intermediate calculation) I can now use time steps greater than 0.01 and achieve accurate results! Also I implemented the 4th order Yoshida Integrator and that produces almost exactly the same as RK4 from a small sample of runs. Next I need a way to benchmark each integrator at different time steps, any pointers?  

---

## 1 answers

### Answer by username_1 on 2020-01-31. Score: 4 | Accepted answer

I can reproduce your JPL Z coordinate plot if I use Earth's heliocentric position relative to the *J2000* ecliptic.
Naturally the amplitude has a minimum around the year 2000:

![EMB.z, J2000 ecliptic, vs. time](https://i.sstatic.net/JxxqZ.png)

Relative to the ecliptic *of date*, which accounts for precession, the heliocentric Z coordinate of the Earth-Moon barycenter is much smaller:

![EMB.z, ecliptic of date, vs. time](https://i.sstatic.net/7EgUn.png)

The IAU 2006 precession model has two components, detailed in [Capitaine et al. 2003](https://ui.adsabs.harvard.edu/abs/2003A%26A...412..567C/abstract):

* precession of the *ecliptic*, a 47"/century change in the plane of Earth's orbit
* precession of the *equator*, a 1.4°/century change in the plane and axis of Earth's rotation

Fifty years of ecliptic precession amount to 23.5" or 0.000114 radian, and the first plot is consistent with that.
The classical precession of the equinoxes (the intersection of the two planes) is a combination of both components.

[Williams 1994](https://ui.adsabs.harvard.edu/abs/1994AJ....108..711W/abstract) examines Earth's equatorial precession (dominated by the Sun and Moon) and finds that Venus's effect is larger than Jupiter's.
When looking at [fluctuations in Earth's year length](/a/34752), I noticed a cycle similar to Venus's synodic period.
I expect that adding Venus to your system would reproduce some ecliptic precession.

> **Comments:**
> **username_2**: Thanks for the edit!  
> **username_0**: Thanks for the information! I had tried adding Venus before I rewrote my code to use AU instead of KM and it didn't make much difference then since it was quite inaccurate. Adding it now has made it much similar! I didn't realise the degree to which Venus has!  
> **username_2**: @username_0 Calculate $m/r^2$ for the distance of closest approach and Venus turns out to be almost as "strong" as Jupiter. Then factoring in the higher frequency of its perturbation and the similarity with Earth's period.  

---

