Orbit lab

Units and the model

The lab integrates Newton's equations for N point masses,

\ddot{\mathbf r}_i = G \sum_{j \ne i} m_j \frac{\mathbf r_j - \mathbf r_i}{\lvert \mathbf r_j - \mathbf r_i \rvert^3}, \qquad i = 1, \dots, N,

in astronomical units (AU), years and solar masses. In these units the gravitational constant is

G = 4\pi^2 \ \mathrm{AU^3\, M_\odot^{-1}\, yr^{-2}},

which just says that a massless body on a circular orbit of radius 1 AU around one solar mass takes one year to go round. That's the Earth, more or less. Kepler's third law then reads T^2 = a^3/M with T in years, a in AU and M in solar masses [1, Ch. 2]. The year defined this way differs from the calendar year by parts per million, which doesn't matter here.

The exact dynamics conserve two quantities, and the lab watches both. One is the total energy

E = \sum_i \tfrac12 m_i \lvert \dot{\mathbf r}_i \rvert^2 \; - \; \sum_{i<j} \frac{G m_i m_j}{\lvert \mathbf r_i - \mathbf r_j \rvert},

and the other is the total angular momentum \mathbf L = \sum_i m_i\, \mathbf r_i \times \dot{\mathbf r}_i. Every scenario starts in the centre of mass frame with zero total momentum, so the centre of mass (the small cross on the canvas) stays at the origin.

Kepler's laws and Newton's derivation

Kepler pulled three rules out of Tycho Brahe's observations of Mars. Planets move on ellipses with the Sun at one focus. The line from the Sun to a planet sweeps out equal areas in equal times. And the square of the period is proportional to the cube of the semimajor axis [1, Ch. 1]. The first two appeared in Astronomia Nova (1609), the third in Harmonices Mundi (1619).

Newton got all three from a single force law. Book 1, Proposition 1 of the Principia shows that any force directed toward a fixed centre makes the radius sweep equal areas in equal times; Proposition 11 shows that motion on an ellipse with the force directed to a focus requires a force inversely proportional to the square of the distance [2]. In modern language, the equal areas law is conservation of angular momentum: the areal velocity is

\frac{dA}{dt} = \frac{\lvert \mathbf r \times \dot{\mathbf r} \rvert}{2} = \frac{h}{2},

which is constant for any central force. The third law follows from integrating that rate over one orbit: the area of the ellipse is \pi a b, so T = 2\pi a b/h, and with h^2 = G M a (1-e^2) one gets

T = 2\pi \sqrt{\frac{a^3}{G(m_1 + m_2)}}.

Notice the total mass. When two bodies have comparable mass, both orbit their common centre of mass, and the period depends on m_1 + m_2 [1, Ch. 2].

The two-body problem and orbital elements

Two bodies reduce exactly to one. The relative vector \mathbf r = \mathbf r_2 - \mathbf r_1 obeys \ddot{\mathbf r} = -\mu \mathbf r / r^3 with \mu = G(m_1+m_2). Its solution is a conic. A bound orbit is described by six numbers, the orbital elements: semimajor axis a, eccentricity e, inclination I, longitude of the ascending node \Omega, argument of perihelion \omega and an angle that fixes where the body is along the orbit, usually the mean anomaly M [1, Ch. 2].

Position as a function of time comes from Kepler's equation,

M = n(t - \tau) = E - e \sin E, \qquad n = \sqrt{\mu / a^3},

where E is the eccentric anomaly. Given E, the position in the orbital plane is x' = a(\cos E - e), y' = a\sqrt{1-e^2}\,\sin E, and three rotations through \omega, I and \Omega carry it into the ecliptic frame. The lab solves Kepler's equation by Newton iteration and differentiates the result for the velocity.

That's how the lab places the inner planets. JPL publishes Keplerian elements fitted to its ephemerides, by E. M. Standish and J. G. Williams, together with the recipe above [3]. The lab takes Table 1 (valid 1800 to 2050) at J2000 and converts the published mean longitude L and longitude of perihelion \varpi to M = L - \varpi and \omega = \varpi - \Omega. The planets' masses come from the Sun to planet mass ratios of the IAU Working Group on Numerical Standards [4]. JPL's own accuracy table gives errors of order 15 to 40 arcseconds for these planets over the fitted interval. After that, the integrated orbits drift away from the real ones, because the lab leaves out the outer planets that tug on them.

The lab also runs this backwards. From the relative position and velocity of two bodies it computes the orbital energy \varepsilon = v^2/2 - \mu/r, hence the osculating semimajor axis a = -\mu/(2\varepsilon), and prints the Kepler period of that ellipse next to the period it measures.

The three-body problem, Poincaré and chaos

Add a third body and the closed form is gone. In his 1890 memoir in Acta Mathematica, Poincaré developed integral invariants and the theory of periodic solutions for the problem, and in the second part studied asymptotic surfaces, the sets of orbits that approach a periodic orbit in the far past or far future [5]. Where such surfaces cross, orbits wander in tangled ways that hinge on the tiniest change in the starting point. These crossings, homoclinic tangles in modern language, are one of the basic ways chaos arises in Hamiltonian systems.

The cleanest way to see it is Burrau's problem of 1913 [6]: masses 3, 4 and 5 at the vertices of a 3-4-5 right triangle, each opposite the side of matching length, released from rest. Burrau's integration ended at t \approx 3.35 (in units with G = 1). Szebehely and Peters completed it in 1967 with regularised variables and an adaptive fifth order Runge-Kutta scheme [7]. They found close approach after close approach, the closest about 4 \times 10^{-4} apart, and an ending in which masses 4 and 5 pair off for good while mass 3 is thrown out. Their comment on a near return of the initial configuration at t \approx 31.66: "The sensitivity of this dynamical system to initial conditions is well demonstrated in Fig. 5." [7] In the lab's units time is divided by 2\pi, so the ejection happens after about 11 years.

Lagrange points and Trojans

In 1772 Lagrange found that three bodies can turn together, rigidly, at the corners of an equilateral triangle. In the restricted problem (a massless particle moving in the field of two masses on circular orbits), this gives two equilibria in the rotating frame, L_4 and L_5, 60 degrees ahead of and behind the smaller mass. They are linearly stable only if the mass ratio \mu = m_2/(m_1+m_2) satisfies 27\mu(1-\mu) < 1, that is \mu < 0.0385 [1, Ch. 3]. Jupiter's \mu \approx 9.5 \times 10^{-4} is far below that, and sure enough its Trojan asteroids cluster around L_4 and L_5 [1, Ch. 3]. A body displaced from L_4 librates around it on a tadpole-shaped path in the rotating frame; larger displacements give horseshoe orbits that enclose L_3, L_4 and L_5 [1, Ch. 3].

Choreographies

Picture three bodies chasing each other around one closed curve, forever. That's a choreography: a periodic solution in which all the bodies share a single path. Moore found the figure-eight for three equal masses numerically in 1993, using a gradient flow on the action [8; 9]. Chenciner and Montgomery proved it exists in 2000 with a variational argument [9]. In their abstract the three bodies "chase each other around a fixed eight-shaped curve", with zero angular momentum and a high degree of symmetry. The paper prints the initial conditions computed by Carles Simó, with G = m_i = 1: \mathbf x_1 = -\mathbf x_2 = (0.97000436, -0.24308753), \mathbf x_3 = 0, \dot{\mathbf x}_3 = -2\dot{\mathbf x}_1 = -2\dot{\mathbf x}_2 = (-0.93240737, -0.86473146), and period 6.32591398 [9]. The lab uses these with velocities scaled by 2\pi for G = 4\pi^2, so one period lasts 1.0068 years. Simó's computations also indicate that the orbit is stable [9]. Systematic numerical searches have since found many more families of planar three-body periodic orbits [10].

Integrating Newton's equations

You can pick from four integrators, all advancing positions \mathbf q and velocities \mathbf v by a step h.

Explicit Euler sets \mathbf q_{n+1} = \mathbf q_n + h \mathbf v_n and \mathbf v_{n+1} = \mathbf v_n + h\, \mathbf a(\mathbf q_n). It's first order: the global error scales as h.

Classical Runge-Kutta (RK4) evaluates the force four times per step and combines the slopes; its global error scales as h^4 [14, Ch. II].

Leapfrog, in its velocity Verlet or Störmer-Verlet form, kicks the velocity by half a step, drifts the position a full step and kicks again:

\mathbf v_{n+1/2} = \mathbf v_n + \tfrac h2 \mathbf a(\mathbf q_n), \quad \mathbf q_{n+1} = \mathbf q_n + h\, \mathbf v_{n+1/2}, \quad \mathbf v_{n+1} = \mathbf v_{n+1/2} + \tfrac h2 \mathbf a(\mathbf q_{n+1}).

Verlet used it for molecular dynamics in 1967 [11]. It's second order and needs one new force evaluation per step. It runs the same backwards as forwards. And it's symplectic, meaning it preserves the phase space area form exactly, as the true flow does [14, Ch. VI; 15].

Yoshida's fourth order method composes three leapfrog steps of lengths w_1 h, w_0 h, w_1 h with

w_1 = \frac{1}{2 - 2^{1/3}}, \qquad w_0 = -\frac{2^{1/3}}{2 - 2^{1/3}}.

The weights sum to one and 2 w_1^3 + w_0^3 = 0, which cancels the leading error term of the symmetric leapfrog. Yes, the middle step runs backward in time. Yoshida showed how such compositions build symmetric symplectic integrators of any even order [13]; Ruth had earlier constructed explicit canonical integrators for Hamiltonian systems [12]. A composition of symplectic maps is symplectic, so Yoshida's method inherits leapfrog's structure.

Why symplectic integrators keep energy

A non-symplectic method makes a small energy error at each step, and nothing stops those errors from piling up. For RK4 on a Kepler orbit the energy error grows roughly linearly in time. For Euler it grows so fast you can watch the orbit spiral outward.

A symplectic method doesn't, and backward error analysis explains why. The numerical map of a symplectic integrator of order p, applied with step h to a Hamiltonian H, is very nearly the exact time-h flow of a modified Hamiltonian

\tilde H = H + h^p H_{p+1} + h^{p+1} H_{p+2} + \cdots

The series does not converge in general, but truncated optimally it describes the numerical solution up to an error exponentially small in 1/h over exponentially long times [16; 14, Ch. IX]. The numerical solution therefore conserves \tilde H almost exactly, and since \tilde H differs from H by O(h^p), the energy error stays bounded at O(h^p) with no drift, for times far longer than any simulation runs [14, Ch. IX and X]. The error wobbles, peaking at close approaches where the force changes fastest, but it doesn't grow.

Angular momentum is even better off. It is a quadratic invariant of the form \mathbf q^{\mathsf T} C \mathbf p, and symplectic partitioned Runge-Kutta methods such as leapfrog conserve such invariants exactly, up to round-off [14, Ch. IV]. Euler and RK4 do not.

There's a catch. The modified Hamiltonian exists only if the step is fixed. Changing h from step to step according to the state, as the lab's adaptive mode does, breaks the argument, and a naive variable step symplectic method loses its long-time energy behaviour [14, Ch. VIII]. The lab uses adaptive steps only to survive close encounters, where a fixed step would be hopeless. Production codes for planetary systems use symplectic maps built around the Kepler problem itself [17].

Circumbinary planets

Kepler-16b orbits both stars of an eclipsing binary. Doyle and coworkers measured stellar masses of 0.6897 and 0.20255 solar masses, a binary semimajor axis of 0.22431 AU with eccentricity 0.15944, and a planet of 0.333 Jupiter masses on a nearly circular orbit at 0.7048 AU [19]. Holman and Wiegert fitted the inner edge of the stable region around a binary from long numerical integrations,

\frac{a_c}{a_b} = 1.60 + 5.10 e - 2.22 e^2 + 4.12 \mu - 4.27 e\mu - 5.09 \mu^2 + 4.61 e^2 \mu^2,

good to a few percent [18]. For Kepler-16 this gives a_c \approx 0.65 AU. The planet sits only about 9 percent outside the limit, close to the edge of where a planet can survive at all.

Things to try in the lab

Euler versus leapfrog. The Earth on a circular orbit with a step of 0.005 yr, run twice. The red Euler orbit spirals outward and its energy error rises steadily on the log plot; the leapfrog orbit stays on the circle with an error below 10^{-6}. Try other steps: Euler's error per orbit scales as h, leapfrog's bounded error as h^2.

RK4 versus Yoshida, 300 years. Both are fourth order, on a comet with e = 0.8. RK4's error is smaller at first, but it climbs with every orbit; Yoshida's spikes at each perihelion and comes back down. Within a few decades the curves cross.

Kepler's second law. Wedges are drawn every twelfth of a period. Near perihelion they are short and wide, near aphelion long and thin. The wedge spread readout stays near 10^{-5}. Turn the step up to 0.02 yr and watch it grow.

Inner planets from JPL elements. Mercury, Venus, the Earth-Moon barycentre and Mars, placed from Standish's J2000 elements. Mercury's 7 degree inclination is in the model but invisible in this top view.

A Trojan near L4. In the corotating frame the Sun and Jupiter stand still and a massless Trojan, started 66 degrees ahead of Jupiter, librates around L_4; the lab measures a libration period of roughly 150 years. Throw it a little faster to make a horseshoe, or much faster to lose it.

Figure-eight choreography. Run it with Yoshida, then with Euler at the same step and watch the eight unravel. Nudge one body by dragging it and the choreography breaks.

Pythagorean three-body problem. With adaptive steps the lab reproduces Szebehely and Peters' outcome: 4 and 5 pair off and 3 leaves toward the upper right. Switch to fixed steps or another integrator and a different body may escape, or none at all. A chaotic outcome is only as good as the integration that produced it.

Circumbinary planet. Watch the binary wobble the planet's orbit. Then press on the planet and throw it inward, inside 0.65 AU, and see how fast the binary flings it out or swallows it.

Reading the lab's instruments

Time is in years of the G = 4\pi^2 system. Energy is the total E above, in \mathrm{M_\odot\,AU^2\,yr^{-2}}. \lvert\Delta E / E_0\rvert compares it with its value at the last restart or throw; the comparison row shows the same for the second integrator, which runs on exactly the same sequence of steps. Values below 10^{-16} are shown as 10^{-16}, the round-off floor.

Angular momentum L_z is the component perpendicular to the screen. \lvert\Delta \mathbf L\rvert / L scale is the change in the full vector, divided by \lvert \mathbf L_0 \rvert or, when that vanishes as in the figure-eight and the Pythagorean problem, by M\sqrt{G M R} with M the total mass and R the scenario's view radius.

Period, measured appears for two-body scenarios. The lab follows the angle of the second body around the first, unwrapped, and records the moments it completes each full turn, interpolating between steps; the readout is the time between the last two. Period, Kepler III is 2\pi\sqrt{a^3/G(m_1+m_2)} with a from the current energy. Their difference is the integrator's phase error.

Wedge area spread is (A_{\max} - A_{\min})/\bar A over the last twelve wedges, each closed exactly at its boundary time by interpolation. Kepler's second law says it is zero.

The plots show both relative errors against time on logarithmic axes, where linear growth appears as a curve that keeps climbing and bounded oscillation as a band. The canvas draws trails in the chosen frame, bodies with radii growing with the logarithm of their mass, the comparison run in red and the centre of mass as a cross. Press on empty sky to add a body of one Earth mass, or on a body to grab it; drag out its velocity, read the speed against the local circular speed and let go. The throw is measured relative to the rotating frame when that frame is on.

References

[1] C. D. Murray and S. F. Dermott, Solar System Dynamics (Cambridge University Press, 1999). Ch. 1 "Structure of the Solar System"; Ch. 2 "The Two-Body Problem"; Ch. 3 "The Restricted Three-Body Problem". DOI: 10.1017/CBO9781139174817

[2] I. Newton, The Principia: Mathematical Principles of Natural Philosophy, translated by I. B. Cohen and A. Whitman (University of California Press, 1999). Book 1, Propositions 1 and 11.

[3] E. M. Standish and J. G. Williams, "Approximate Positions of the Planets," JPL Solar System Dynamics, ssd.jpl.nasa.gov/planets/approx_pos.html, Table 1. The page reproduces their 1992 article with permission.

[4] IAU Working Group on Numerical Standards for Fundamental Astronomy, "Current Best Estimates," iau-a3.gitlab.io/NSFA/NSFA_cbe.html; B. Luzum et al., "The IAU 2009 system of astronomical constants," Celestial Mechanics and Dynamical Astronomy 110, 293-304 (2011). DOI: 10.1007/s10569-011-9352-4

[5] H. Poincaré, "Sur le problème des trois corps et les équations de la dynamique," Acta Mathematica 13, 1-270 (1890). DOI of the introduction: 10.1007/BF02392506; second part, ch. II "Etude des surfaces asymptotiques": 10.1007/BF02392511

[6] C. Burrau, "Numerische Berechnung eines Spezialfalles des Dreikörperproblems," Astronomische Nachrichten 195, 113-118 (1913). DOI: 10.1002/asna.19131950602

[7] V. Szebehely and C. F. Peters, "Complete solution of a general problem of three bodies," The Astronomical Journal 72, 876-883 (1967). DOI: 10.1086/110355

[8] C. Moore, "Braids in classical dynamics," Physical Review Letters 70, 3675-3679 (1993). DOI: 10.1103/PhysRevLett.70.3675

[9] A. Chenciner and R. Montgomery, "A Remarkable Periodic Solution of the Three-Body Problem in the Case of Equal Masses," Annals of Mathematics 152, 881-901 (2000). DOI: 10.2307/2661357. arXiv: math/0011268

[10] M. Šuvakov and V. Dmitrašinović, "Three Classes of Newtonian Three-Body Planar Periodic Orbits," Physical Review Letters 110, 114301 (2013). DOI: 10.1103/PhysRevLett.110.114301. arXiv: 1303.0181

[11] L. Verlet, "Computer 'Experiments' on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules," Physical Review 159, 98-103 (1967). DOI: 10.1103/PhysRev.159.98

[12] R. D. Ruth, "A Canonical Integration Technique," IEEE Transactions on Nuclear Science 30, 2669-2671 (1983). DOI: 10.1109/TNS.1983.4332919

[13] H. Yoshida, "Construction of higher order symplectic integrators," Physics Letters A 150, 262-268 (1990). DOI: 10.1016/0375-9601(90)90092-3

[14] E. Hairer, C. Lubich and G. Wanner, Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations, 2nd ed. (Springer, 2006). Ch. I "Examples and Numerical Experiments"; Ch. II "Numerical Integrators"; Ch. IV "Conservation of First Integrals and Methods on Manifolds"; Ch. VI "Symplectic Integration of Hamiltonian Systems"; Ch. VIII "Structure-Preserving Implementation"; Ch. IX "Backward Error Analysis and Structure Preservation"; Ch. X "Hamiltonian Perturbation Theory and Symplectic Integrators". DOI: 10.1007/3-540-30666-8

[15] E. Hairer, C. Lubich and G. Wanner, "Geometric numerical integration illustrated by the Störmer-Verlet method," Acta Numerica 12, 399-450 (2003). DOI: 10.1017/S0962492902000144

[16] G. Benettin and A. Giorgilli, "On the Hamiltonian interpolation of near-to-the identity symplectic mappings with application to symplectic integration algorithms," Journal of Statistical Physics 74, 1117-1143 (1994). DOI: 10.1007/BF02188219

[17] J. Wisdom and M. Holman, "Symplectic maps for the n-body problem," The Astronomical Journal 102, 1528 (1991). DOI: 10.1086/115978

[18] M. J. Holman and P. A. Wiegert, "Long-Term Stability of Planets in Binary Systems," The Astronomical Journal 117, 621-628 (1999). DOI: 10.1086/300695. arXiv: astro-ph/9809315

[19] L. R. Doyle et al., "Kepler-16: A Transiting Circumbinary Planet," Science 333, 1602-1606 (2011). DOI: 10.1126/science.1210923. arXiv: 1109.3432