Astronomy:Kepler's 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 derived by Johannes Kepler in 1609 in Chapter 60 of his Astronomia nova, and in book V of his Epitome of Copernican Astronomy (1621) Kepler proposed an iterative solution to the equation. This equation and its solution, however, first appeared in 9th century work of Habash al-Hasib al-Marwazi related to problems of parallax. The equation has played an important role in the history of both physics and mathematics, particularly classical celestial mechanics.

Equation

Kepler's equation is

[math]\displaystyle{ M = E - e \sin E }[/math]

where [math]\displaystyle{ M }[/math] is the mean anomaly, [math]\displaystyle{ E }[/math] is the eccentric anomaly, and [math]\displaystyle{ e }[/math] is the eccentricity.

The 'eccentric anomaly' [math]\displaystyle{ E }[/math] 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 [math]\displaystyle{ x = a(1 - e) }[/math], [math]\displaystyle{ y = 0 }[/math], at time [math]\displaystyle{ t = t_0 }[/math], then to find out the position of the body at any time, you first calculate the mean anomaly [math]\displaystyle{ M }[/math] from the time and the mean motion [math]\displaystyle{ n }[/math] by the formula [math]\displaystyle{ M = n(t - t_0) }[/math], then solve the Kepler equation above to get [math]\displaystyle{ E }[/math], then get the coordinates from:

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

where [math]\displaystyle{ a }[/math] is the semi-major axis, [math]\displaystyle{ b }[/math] the semi-minor axis.

Kepler's equation is a transcendental equation because sine is a transcendental function, meaning it cannot be solved for [math]\displaystyle{ E }[/math] algebraically. Numerical analysis and series expansions are generally required to evaluate [math]\displaystyle{ E }[/math].

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 ([math]\displaystyle{ 0 \le e \lt 1 }[/math]). The hyperbolic Kepler equation is used for hyperbolic trajectories ([math]\displaystyle{ e \gt 1 }[/math]). The radial Kepler equation is used for linear (radial) trajectories ([math]\displaystyle{ e = 1 }[/math]). Barker's equation is used for parabolic trajectories ([math]\displaystyle{ e = 1 }[/math]).

When [math]\displaystyle{ e = 0 }[/math], the orbit is circular. Increasing [math]\displaystyle{ e }[/math] causes the circle to become elliptical. When [math]\displaystyle{ e = 1 }[/math], there are three possibilities:

• a parabolic trajectory,
• a trajectory going in or out along an infinite ray emanating from the centre of attraction,
• or a trajectory that goes back and forth along a line segment from the centre of attraction to a point at some distance away.

A slight increase in [math]\displaystyle{ e }[/math] above 1 results in a hyperbolic orbit with a turning angle of just under 180 degrees. Further increases reduce the turning angle, and as [math]\displaystyle{ e }[/math] goes to infinity, the orbit becomes a straight line of infinite length.

Hyperbolic Kepler equation

The Hyperbolic Kepler equation is:

[math]\displaystyle{ M = e \sinh H - H }[/math]

where [math]\displaystyle{ H }[/math] is the hyperbolic eccentric anomaly. This equation is derived by redefining M to be the square root of −1 times the right-hand side of the elliptical equation:

[math]\displaystyle{ M = i \left( E - e \sin E \right) }[/math]

(in which [math]\displaystyle{ E }[/math] is now imaginary) and then replacing [math]\displaystyle{ E }[/math] by [math]\displaystyle{ iH }[/math].

Radial Kepler equation

The Radial Kepler equation is:

[math]\displaystyle{ t(x) = \sin^{-1}( \sqrt{ x } ) - \sqrt{ x ( 1 - x ) } }[/math]

where [math]\displaystyle{ t }[/math] is proportional to time and [math]\displaystyle{ x }[/math] is proportional to the distance from the centre of attraction along the ray. This equation is derived by multiplying Kepler's equation by 1/2 and setting [math]\displaystyle{ e }[/math] to 1:

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

and then making the substitution

[math]\displaystyle{ E = 2 \sin^{-1}(\sqrt{ x }). }[/math]

Inverse problem

Calculating [math]\displaystyle{ M }[/math] for a given value of [math]\displaystyle{ E }[/math] is straightforward. However, solving for [math]\displaystyle{ E }[/math] when [math]\displaystyle{ M }[/math] is given can be considerably more challenging. There is no closed-form solution.

One can write an infinite series expression for the solution to Kepler's equation using Lagrange inversion, but the series does not converge for all combinations of [math]\displaystyle{ e }[/math] and [math]\displaystyle{ M }[/math] (see below).

Confusion over the solvability of Kepler's equation has persisted in the literature for four centuries. Kepler himself expressed doubt at the possibility of finding a general solution:

I am sufficiently satisfied that it [Kepler's equation] cannot be solved a priori, on account of the different nature of the arc and the sine. But if I am mistaken, and any one shall point out the way to me, he will be in my eyes the great Apollonius.

—Johannes Kepler

Fourier series expansion (with respect to [math]\displaystyle{ M }[/math]) using Bessel functions is

[math]\displaystyle{ E = M + \sum_{m=1}^\infty \frac{2}{m} J_m(me) \sin(mM), \quad e\le 1, \quad M \in [-\pi,\pi]. }[/math]

With respect to [math]\displaystyle{ e }[/math], it is a Kapteyn series.

Inverse Kepler equation

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

[math]\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 \ne 1 \end{cases} }[/math]

Evaluating this yields:

[math]\displaystyle{ E = \begin{cases} \displaystyle s + \frac{1}{60} s^3 + \frac{1}{1400}s^5 + \frac{1}{25200}s^7 + \frac{43}{17248000}s^9 + \frac{ 1213}{7207200000 }s^{11} + \frac{151439}{12713500800000 }s^{13}+ \cdots \text{ with }s = ( 6 M )^{1/3} , & e = 1\\ \\ \displaystyle \frac{1}{1-e} M - \frac{e}{(1-e)^4 } \frac{M^3}{3!} + \frac{(9 e^2 + e)}{(1-e)^7 } \frac{M^5}{5!} - \frac{(225 e^3 + 54 e^2 + e) }{(1-e)^{10} } \frac{M^7}{7!} + \frac{ (11025 e^4 + 4131 e^3 + 243 e^2 + e) }{(1-e)^{13} } \frac{M^9}{9!}+ \cdots , & e \ne 1 \end{cases} }[/math]

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 Maclaurin series. Such 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. However, [math]\displaystyle{ E }[/math] is not an entire function of [math]\displaystyle{ M }[/math] at a given non-zero [math]\displaystyle{ e }[/math]. Indeed, the derivative

[math]\displaystyle{ \mathrm{dM}/\mathrm{d}E=1-e\cos E }[/math]

goes to zero at an infinite set of complex numbers when [math]\displaystyle{ e \lt 1, }[/math] the nearest to zero being at [math]\displaystyle{ E=\pm i\cosh^{-1}(1/e), }[/math] and at these two points

[math]\displaystyle{ M=E-e\sin E=\pm i\left(\cosh^{-1}(1/e)-\sqrt{1-e^2}\right) }[/math]

(where inverse cosh is taken to be positive), and [math]\displaystyle{ \mathrm{d}E/\mathrm{d}M }[/math] goes to infinity at these values of [math]\displaystyle{ M }[/math]. This means that the radius of convergence of the Maclaurin series is [math]\displaystyle{ \cosh^{-1}(1/e)-\sqrt{1-e^2} }[/math] and the series will not converge for values of [math]\displaystyle{ M }[/math] larger than this. The series can also be used for the hyperbolic case, in which case the radius of convergence is [math]\displaystyle{ \cos^{-1}(1/e)-\sqrt{e^2-1}. }[/math] The series for when [math]\displaystyle{ e = 1 }[/math] converges when [math]\displaystyle{ M \lt 2\pi }[/math].

While this solution is the simplest in a certain mathematical sense,, other solutions are preferable for most applications. Alternatively, Kepler's equation can be solved numerically.

The solution for [math]\displaystyle{ e \ne 1 }[/math] was found by Karl Stumpff in 1968, but its significance wasn't recognized.

One can also write a Maclaurin series in [math]\displaystyle{ e }[/math]. This series does not converge when [math]\displaystyle{ e }[/math] is larger than the Laplace limit (about 0.66), regardless of the value of [math]\displaystyle{ M }[/math] (unless [math]\displaystyle{ M }[/math] is a multiple of 2π), but it converges for all [math]\displaystyle{ M }[/math] if [math]\displaystyle{ e }[/math] is less than the Laplace limit. The coefficients in the series, other than the first (which is simply [math]\displaystyle{ M }[/math]), depend on [math]\displaystyle{ M }[/math] in a periodic way with period 2π.

Inverse radial Kepler equation

The inverse radial Kepler equation ([math]\displaystyle{ e = 1 }[/math]) can also be written as:

[math]\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] }[/math]

Evaluating this yields:

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

To obtain this result using Mathematica:

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

Numerical approximation of inverse problem

Newton's method

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

[math]\displaystyle{ f(E) = E - e \sin(E) - M(t) }[/math]

This can be done iteratively via Newton's method:

[math]\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})} }[/math]

Note that [math]\displaystyle{ E }[/math] and [math]\displaystyle{ M }[/math] are in units of radians in this computation. This iteration is repeated until desired accuracy is obtained (e.g. when [math]\displaystyle{ f(E) }[/math] < desired accuracy). For most elliptical orbits an initial value of [math]\displaystyle{ E_0 = M(t) }[/math] is sufficient. For orbits with [math]\displaystyle{ e \gt 0.8 }[/math], a initial value of [math]\displaystyle{ E_0 = \pi }[/math] can be used. Numerous works developed accurate (but also more complex) start guesses. If [math]\displaystyle{ e }[/math] is identically 1, then the derivative of [math]\displaystyle{ f }[/math], which is in the denominator of Newton's method, can get close to zero, making derivative-based methods such as Newton-Raphson, secant, or regula falsi numerically unstable. In that case, the bisection method will provide guaranteed convergence, particularly since the solution can be bounded in a small initial interval. On modern computers, it is possible to achieve 4 or 5 digits of accuracy in 17 to 18 iterations. 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.

Fixed-point iteration

A related method starts by noting that [math]\displaystyle{ E = M + e \sin{E} }[/math]. Repeatedly substituting the expression on the right for the [math]\displaystyle{ E }[/math] on the right yields a simple fixed-point iteration algorithm for evaluating [math]\displaystyle{ E(e,M) }[/math]. This method is identical to Kepler's 1621 solution.

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, [math]\displaystyle{ n }[/math], depends on the value of [math]\displaystyle{ e }[/math]. The hyperbolic form similarly has [math]\displaystyle{ H = e \sinh H - M }[/math].

This method is related to the Newton's method solution above in that

[math]\displaystyle{ E_{n+1} = E_{n} - \frac{E_{n} - e \sin(E_{n}) - M(t)}{ 1 - e \cos(E_{n})} = E_{n} + \frac {(M + e \sin{E_{n}} - E_{n})(1 + e \cos{E_{n}})}{1 - e^2 (\cos{E_{n}})^2} }[/math]

To first order in the small quantities [math]\displaystyle{ M-E_{n} }[/math] and [math]\displaystyle{ e }[/math],

[math]\displaystyle{ E_{n+1} \approx M + e \sin{E_{n}} }[/math].

See also

• Kepler's laws of planetary motion
• Kepler problem
• Kepler problem in general relativity
• Radial trajectory

External links

• Danby, John M.; Burkardt, Thomas M. (1983). "The solution of Kepler's equation. I". Celestial Mechanics 31 (2): 95–107. doi:10.1007/BF01686811. Bibcode: 1983CeMec..31...95D. 
• Conway, Bruce A. (1986). "24th Aerospace Sciences Meeting". doi:10.2514/6.1986-84. 
• Mikkola, Seppo (1987). "A cubic approximation for Kepler's equation". Celestial Mechanics 40 (3): 329–334. doi:10.1007/BF01235850. Bibcode: 1987CeMec..40..329M. https://articles.adsabs.harvard.edu/pdf/1987CeMec..40..329M.pdf. 
• Nijenhuis, Albert (1991). "Solving Kepler's equation with high efficiency and accuracy". Celestial Mechanics and Dynamical Astronomy 51 (4): 319–330. doi:10.1007/BF00052925. Bibcode: 1991CeMDA..51..319N. 
• Markley, F. Landis (1995). "Kepler equation solver". Celestial Mechanics and Dynamical Astronomy 63 (1): 101–111. doi:10.1007/BF00691917. Bibcode: 1995CeMDA..63..101M. 
• Fukushima, Toshio (1996). "A method solving kepler's equation without transcendental function evaluations". Celestial Mechanics and Dynamical Astronomy 66 (3): 309–319. doi:10.1007/BF00049384. Bibcode: 1996CeMDA..66..309F. 
• Charles, Edgar D.; Tatum, Jeremy B. (1997). "The convergence of Newton-Raphson iteration with Kepler's equation". Celestial Mechanics and Dynamical Astronomy 69 (4): 357–372. doi:10.1023/A:1008200607490. Bibcode: 1997CeMDA..69..357C. 
• Stumpf, Laura (1999). "Chaotic behaviour in the Newton iterative function associated with Kepler's equation". Celestial Mechanics and Dynamical Astronomy 74 (2): 95–109. doi:10.1023/A:1008339416143. Bibcode: 1999CeMDA..74...95S. 
• Palacios, Manuel (2002). "Kepler equation and accelerated Newton method". Journal of Computational and Applied Mathematics 138 (2): 335–346. doi:10.1016/S0377-0427(01)00369-7. Bibcode: 2002JCoAM.138..335P. 
• Boyd, John P. (2007). "Rootfinding for a transcendental equation without a first guess: Polynomialization of Kepler's equation through Chebyshev polynomial equation of the sine". Applied Numerical Mathematics 57 (1): 12–18. doi:10.1016/j.apnum.2005.11.010. 
• Pál, András (2009). "An analytical solution for Kepler's problem". Monthly Notices of the Royal Astronomical Society 396 (3): 1737–1742. doi:10.1111/j.1365-2966.2009.14853.x. Bibcode: 2009MNRAS.396.1737P. 
• Esmaelzadeh, Reza; Ghadiri, Hossein (2014). "Appropriate starter for solving the Kepler's equation". International Journal of Computer Applications 89 (7): 31–38. doi:10.5120/15517-4394. Bibcode: 2014IJCA...89g..31E. 
• Zechmeister, Mathias (2018). "CORDIC-like method for solving Kepler's equation". Astronomy and Astrophysics 619: A128. doi:10.1051/0004-6361/201833162. Bibcode: 2018A&A...619A.128Z. 

• Kepler's Equation at Wolfram Mathworld