How do you calculate the effects of precession on elliptical orbits?

Question

Kepler's first law states that planets (and all celestial bodies orbiting another body) travel in elliptical orbits, which have well-known formulas that make it relatively easy to calculate the orbital elements and associated behavior. However, ongoing precession means that the orbit is constantly changing - and so the planet isn't actually traveling in the ellipse it originally set out on! You can calculate precession and its related effects (this question and answer are helpful), but is there any way to calculate how the elliptical orbit will be "deformed" by precession?

Answer 1

A good starting point would be <insert name of some scientist from long ago> planetary equations of motion. For example, there are Lagrange's planetary equations (sometimes called the Lagrange-Laplace planetary equations), Gauss' planetary equations, Delaunay's planetary equations, Hill's planetary equations, and several more. The common theme amongst these various planetary equations is that they yield the time derivatives of various orbital elements as a function of the partial derivatives of the perturbing force / perturbing potential with respect to some generalized position.

In general, the only words that can describe the result of this process at first is "hot mess." A hot mess did not deter those brilliant minds of old. Via various simplifying assumptions and long term time averaging, they came up with fairly simple descriptions of, for example, $\left \langle \frac{d\omega}{dt} \right\rangle$ (apsidal precession) and $\left \langle \frac{d\Omega}{dt} \right\rangle$ (planar precession). You can see some of this in the cited 1900 work by Hill below.

While these techniques are old, these planetary equations are still used today. That sometimes you do get a "hot mess" is okay now that we have computers. People are using planetary equations coupled with geometric integration techniques to yield integrators that are fast, accurate, stable, and conserve angular momentum and energy over long spans of time. (Normally, you can't have all of these. You're lucky if you get just two or three.) Another nice feature of these planetary equations is that they let you see features such as resonances that are otherwise obscured by the truly "hot mess" of the cartesian equations of motion.

Selected reference material, sorted by date:

Hill (1900), "On the Extension of Delaunay's Method in the Lunar Theory to the General Problem of Planetary Motion," Transactions of the American Mathematical Society, 1.2:205-242.

Vallado (1997 and later), "Fundamentals of Astrodynamics and Applications", various publishers. Other than the hole it punches through your wallet, you can't go wrong with this book.

Efroimsky (2002), "Equations for the keplerian elements: hidden symmetry," Institute for Mathematics and its Applications

Efroimsky and Goldreich (2003), "Gauge symmetry of the N-body problem in the Hamilton–Jacobi approach." Journal of Mathematical Physics, 44.12:5958-5977.

Wyatt (2006-2009), Graduate lecture course on planetary systems, Institute of Astronomy, Cambridge.
The results of the Lagrange planetary equations are presented on slide 6.

Ketchum et al. (2013), "Mean Motion Resonances in Exoplanet Systems: An Investigation into Nodding Behavior." The Astrophysical Journal 762.2.

Answer 2

The only truly confocal elliptical orbit is that of a bound test particle in the central potential $-k/r$ or, equivalently, that of two point-like (with spherically symmetric internal mass distributions) masses attracting each other with Newtonian gravity (and having negative total energy, i.e. being bound to each other).

Everything else is non-elliptic (unbound orbits are parabolic or hyperbolic), but most deviations are small. Small deviations can arise from a number of sources, including quadrupole terms in the mass distribution of the bodies (in particualr the Sun), non-gravitational forces (radiation pressure and gas drag on dust grains), non-Newtonian (GR) effects, perturbations from other objects (all the other planets). Newton himself was well aware of this last effect.

If the deviations are small, then the traditional way to estimate them is perturbation theory, where one integrates the perturbing force along the unperturbed (elliptical) orbit. For example, in order to obtain the precession of the periapse, one could integrate the changes to the eccentricity vector. A rotation of that vector corresponds to periapse precession. See my answer to this question, for an example of exactly that.

Answer 3

David Hammen wrote

People are using planetary equations coupled with geometric integration techniques...

You could also try (what I call) a simple finite-step simulation using Newton's laws to operate on object masses, positions, velocities and accelerations. I'm not sure if this falls within what David calls "geometric integration techniques". My point is that you can do it without incorporating the planetary equations. Disadvantage = the simulator "cuts corners" by using approximations and this leads to behaviours in the model which are artefacts. These disadvantages can be overcome by using other techniques. Advantage = it makes the code design easier, it avoids the suspicion that the planetary equations (and their assumptions) are driving the show.

You do not need to be an expert in numerical methods to use the simple Leapfrog Integration technique (described in detail in Feynman Lectures vol I) to model Newtonian Precession in Solar system orbits over periods of up to a few centuries. By running simulations at various time steps (e.g. $dt=1200s, 600s, 300s, 100s$) plotting the results in Excel, fitting a curve and extrapolating to $dt=0$ you can obtain results for long-term average Newtonian precession which are within 1% of the accepted figures. Another advantage over analytical methods which produce long-term average results is that you can examine behaviours at shorter time-scales. For example if you graph perihelion direction vs. time for a certain planet (e.g. Mercury) you can see the $\approx 11.9 $ year periodic fluctuations in precession rate resulting from the movement of Jupiter around the Sun. It is also a lot of fun (and very easy once you have written the basic code) to play "what if?" simulations by varying the number and properties of the bodies in the system and even by adding additional Non-Newtonian forces.

To quote Feymnan:-

It may be that in one cycle of calculation, depending on the problem, we may have 30 multiplications, or something like that, so one cycle will take 300 microseconds. That means that we can do 3000 cycles of computation per second. In order to get an accuracy, of, say, one part in a billion, we would need 4×10^5 cycles to correspond to one revolution of a planet around the sun. That corresponds to a computation time of 130 seconds or about two minutes. Thus it takes only two minutes to follow Jupiter around the sun, with all the perturbations of all the planets correct to one part in a billion, by this method!

But you do have to think carefully about what you can reliably infer from the simulations - for example if your time-step is longer than a few hundred seconds the simulation will indicate precession in the opposite direction to that which really occurs (i.e. retrograde when it should be prograde).