Kepler's equation

(Redirected from Kepler equation)
Kepler's equation solutions for five different eccentricities between 0 and 1

In orbital mechanics, Kepler's equation relates various geometric properties of the orbit of a body subject to a central force.

It was first derived by Johannes Kepler in 1609 in Chapter 60 of his Astronomia nova,[1][2] and in book V of his Epitome of Copernican Astronomy (1621) Kepler proposed an iterative solution to the equation.[3][4] The equation has played an important role in the history of both physics and mathematics, particularly classical celestial mechanics.

Equation

Kepler's equation is

 ${\displaystyle M=E-e\sin E}$

where M is the mean anomaly, E is the eccentric anomaly, and e is the eccentricity.

The 'eccentric anomaly' E is useful to compute the position of a point moving in a Keplerian orbit. As for instance, if the body passes the periastron at coordinates x = a(1 − e), y = 0, at time t = t0, then to find out the position of the body at any time, you first calculate the mean anomaly M from the time and the mean motion n by the formula M = n(tt0), then solve the Kepler equation above to get E, then get the coordinates from:

 ${\displaystyle {\begin{array}{lcl}x&=&a(\cos E-e)\\y&=&b\sin E\end{array}}}$

Kepler's equation is a transcendental equation because sine is a transcendental function, meaning it cannot be solved for E algebraically. Numerical analysis and series expansions are generally required to evaluate E.

Alternate forms

There are several forms of Kepler's equation. Each form is associated with a specific type of orbit. The standard Kepler equation is used for elliptic orbits (0 ≤ e < 1). The hyperbolic Kepler equation is used for hyperbolic orbits (e ≫ 1). The radial Kepler equation is used for linear (radial) orbits (e = 1). Barker's equation is used for parabolic orbits (e = 1). When e = 1, Kepler's equation is not associated with an orbit.

When e = 0, the orbit is circular. Increasing e causes the circle to flatten into an ellipse. When e = 1, the orbit is completely flat, and it appears to be either a segment if the orbit is closed, or a ray if the orbit is open. An infinitesimal increase to e results in a hyperbolic orbit with a turning angle of 180 degrees, and the orbit appears to be a ray. Further increases reduce the turning angle, and as e goes to infinity, the orbit becomes a straight line of infinite length.

Hyperbolic Kepler equation

The Hyperbolic Kepler equation is:

 ${\displaystyle M=e\sinh H-H}$

where H is the hyperbolic eccentric anomaly. This equation is derived by multiplying Kepler's equation by the square root of −1; i = √−1 for imaginary unit, and replacing

${\displaystyle E=iH}$

to obtain

${\displaystyle M=i\left(E-e\sin E\right)}$

 ${\displaystyle t(x)=\sin ^{-1}({\sqrt {x}})-{\sqrt {x(1-x)}}}$

where t is time, and x is the distance along an x-axis. This equation is derived by multiplying Kepler's equation by 1/2 making the replacement

${\displaystyle E=2\sin ^{-1}({\sqrt {x}})}$

and setting e = 1 gives

${\displaystyle t(x)={\frac {1}{2}}\left[E-\sin(E)\right].}$

Inverse problem

Calculating M for a given value of E is straightforward. However, solving for E when M is given can be considerably more challenging.

Kepler's equation can be solved for E analytically by Lagrange inversion. The solution of Kepler's equation given by two Taylor series below.

Confusion over the solvability of Kepler's equation has persisted in the literature for four centuries.[5] Kepler himself expressed doubt at the possibility of ﬁnding a general solution.

Inverse Kepler equation

The inverse Kepler equation is the solution of Kepler's equation for all real values of ${\displaystyle e}$:

${\displaystyle E={\begin{cases}\displaystyle \sum _{n=1}^{\infty }{\frac {M^{\frac {n}{3}}}{n!}}\lim _{\theta \to 0^{+}}\!{\Bigg (}{\frac {\mathrm {d} ^{\,n-1}}{\mathrm {d} \theta ^{\,n-1}}}{\bigg (}{\bigg (}{\frac {\theta }{\sqrt[{3}]{\theta -\sin(\theta )}}}{\bigg )}^{\!\!\!n}{\bigg )}{\Bigg )},&e=1\\\displaystyle \sum _{n=1}^{\infty }{\frac {M^{n}}{n!}}\lim _{\theta \to 0^{+}}\!{\Bigg (}{\frac {\mathrm {d} ^{\,n-1}}{\mathrm {d} \theta ^{\,n-1}}}{\bigg (}{\Big (}{\frac {\theta }{\theta -e\sin(\theta )}}{\Big )}^{\!n}{\bigg )}{\Bigg )},&e\neq 1\end{cases}}}$

Evaluating this yields:

${\displaystyle E={\begin{cases}\displaystyle x+{\frac {1}{60}}x^{3}+{\frac {1}{1400}}x^{5}+{\frac {1}{25200}}x^{7}+{\frac {43}{17248000}}x^{9}+{\frac {1213}{7207200000}}x^{11}+{\frac {151439}{12713500800000}}x^{13}+\cdots \ |\ x=(6M)^{\frac {1}{3}},&e=1\\\\\displaystyle {\frac {1}{1-e}}M-{\frac {e}{(1-e)^{4}}}{\frac {M^{3}}{3!}}+{\frac {(9e^{2}+e)}{(1-e)^{7}}}{\frac {M^{5}}{5!}}-{\frac {(225e^{3}+54e^{2}+e)}{(1-e)^{10}}}{\frac {M^{7}}{7!}}+{\frac {(11025e^{4}+4131e^{3}+243e^{2}+e)}{(1-e)^{13}}}{\frac {M^{9}}{9!}}+\cdots ,&e\neq 1\end{cases}}}$

These series can be reproduced in Mathematica with the InverseSeries operation.

InverseSeries[Series[M - Sin[M], {M, 0, 10}]]
InverseSeries[Series[M - e Sin[M], {M, 0, 10}]]

These functions are simple Taylor series. Taylor series representations of transcendental functions are considered to be definitions of those functions. Therefore, this solution is a formal definition of the inverse Kepler equation. While this solution is the simplest in a certain mathematical sense, for values of ${\displaystyle e}$ near 1 the convergence is very poor, other solutions are preferable for most applications. Alternatively, Kepler's equation can be solved numerically.

The solution for e ≠ 1 was discovered by Karl Stumpff in 1968,[7] but its significance wasn't recognized.[8]

Algorithm for the inverse Kepler equation

Note that ${\displaystyle E=M+e\sin {E}}$. Repeatedly substituting the expression on the right for the ${\displaystyle E}$ on the right yields a simple fixed point iteration algorithm for evaluating ${\displaystyle E(e,M)}$.

function E(e,M,n)
E = M
for k = 1 to n
E = M + e*sin E
next k
return E


The number of iterations, ${\displaystyle n}$, depends on the value of ${\displaystyle e}$.

Using Newton's method to find the zero of ${\displaystyle f=E-e\sin {E}-M}$ one gets,

${\displaystyle \Delta E={\frac {M+e\cdot \sin {E}-E}{1-e\cdot \cos {E}}}={\frac {(M+e\cdot \sin {E}-E)(1+e\cdot \cos {E})}{1-e^{2}\cdot (\cos {E})^{2}}}\approx M-E+e\cdot \sin {E}}$

to first order in the small quantities ${\displaystyle M-E}$ and ${\displaystyle e}$, or,

${\displaystyle E'=E+\Delta E=M+e\cdot \sin {E}}$.

The algorithm actually uses a simplified version of Newton's method to improve on a value for E. Although Newton's method usually converges faster to the correct value for E there are situations where it requires division by 0 and the method above is then preferred.

The inverse radial Kepler equation is:

${\displaystyle x(t)=\sum _{n=1}^{\infty }\left[\lim _{r\to 0^{+}}\left({\frac {t^{{\frac {2}{3}}n}}{n!}}{\frac {\mathrm {d} ^{\,n-1}}{\mathrm {d} r^{\,n-1}}}\!\left(r^{n}\left({\frac {3}{2}}{\Big (}\sin ^{-1}({\sqrt {r}})-{\sqrt {r-r^{2}}}{\Big )}\right)^{\!-{\frac {2}{3}}n}\right)\right)\right]}$

Evaluating this yields:

${\displaystyle x(t)=p-{\frac {1}{5}}p^{2}-{\frac {3}{175}}p^{3}-{\frac {23}{7875}}p^{4}-{\frac {1894}{3931875}}p^{5}-{\frac {3293}{21896875}}p^{6}-{\frac {2418092}{62077640625}}p^{7}-\ \cdots \ {\bigg |}{p=\left({\tfrac {3}{2}}t\right)^{2/3}}}$

To obtain this result using Mathematica:

InverseSeries[Series[ArcSin[Sqrt[t]] - Sqrt[(1 - t) t], {t, 0, 15}]]

Numerical approximation of inverse problem

For most applications, the inverse problem can be computed numerically by finding the root of the function:

${\displaystyle f(E)=E-e\sin(E)-M(t)}$

This can be done iteratively via Newton's method:

${\displaystyle E_{n+1}=E_{n}-{\frac {f(E_{n})}{f'(E_{n})}}=E_{n}-{\frac {E_{n}-e\sin(E_{n})-M(t)}{1-e\cos(E_{n})}}}$

Note that E and M are in units of radians in this computation. This iteration is repeated until desired accuracy is obtained (e.g. when f(E) < desired accuracy). For most elliptical orbits an initial value of E0 = M(t) is sufficient. For orbits with e > 0.8, an initial value of E0 = π should be used.[9] A similar approach can be used for the hyperbolic form of Kepler's equation. In the case of a parabolic trajectory, Barker's equation is used.