4

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:

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)

astrosnapper
  • 8,357
  • 1
  • 21
  • 40
jameslfc19
  • 143
  • 3
  • 1
    What year is day 0? – Mike G Jan 29 '20 at 20:09
  • 1
    @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 Jan 29 '20 at 21:57
  • 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? – uhoh Jan 30 '20 at 01:03
  • 1
    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. I quite like synchronized 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 Feb 01 '20 at 02:37
  • 1
    @PM2Ring Thanks! The Leapfrog Integrator works well! – jameslfc19 Feb 02 '20 at 00:21
  • 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 Feb 02 '20 at 04:49
  • @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 Feb 02 '20 at 04:53
  • @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 Feb 02 '20 at 05:26
  • @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 Feb 02 '20 at 05:49
  • @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 Feb 03 '20 at 18:32

1 Answers1

4

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

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

The IAU 2006 precession model has two components, detailed in Capitaine et al. 2003:

  • 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 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, I noticed a cycle similar to Venus's synodic period. I expect that adding Venus to your system would reproduce some ecliptic precession.

Mike G
  • 18,640
  • 1
  • 26
  • 65
  • Thanks for the edit! – uhoh Jan 31 '20 at 23:14
  • 2
    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 Feb 02 '20 at 00:26
  • @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 Feb 03 '20 at 06:19