4
$\begingroup$

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)

Improve this question
$\endgroup$
10
  • 1
    $\begingroup$ What year is day 0? $\endgroup$ – Mike G Jan 29 '20 at 20:09
  • 1
    $\begingroup$ @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. i.imgur.com/VaIjHUV.png $\endgroup$ – jameslfc19 Jan 29 '20 at 21:57
  • $\begingroup$ 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? $\endgroup$ – uhoh Jan 30 '20 at 1:03
  • 1
    $\begingroup$ 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 artcompsci.org/kali/vol/two_body_problem_2/ch07.html#rdocsect46 $\endgroup$ – PM 2Ring Feb 1 '20 at 2:37
  • 1
    $\begingroup$ @PM2Ring Thanks! The Leapfrog Integrator works well! $\endgroup$ – jameslfc19 Feb 2 '20 at 0:21
4
$\begingroup$

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.

Improve this answer
$\endgroup$
3
  • $\begingroup$ Thanks for the edit! $\endgroup$ – uhoh Jan 31 '20 at 23:14
  • 2
    $\begingroup$ 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! $\endgroup$ – jameslfc19 Feb 2 '20 at 0:26
  • $\begingroup$ @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. $\endgroup$ – uhoh Feb 3 '20 at 6:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.