Kepler’s Equation of Elliptical Motion Solver

Kepler’s Equation of Elliptical Motion is a famous example of an equation for which an exact, analytical, solution is not possible. In its most succinct form, it is stated as follows:

E - e sin(E) = M

where e is the eccentricity of the ellipse described by the motion
(0 < e < 1);
E is the eccentric anomaly; and
M is the mean anomaly, 2 π /T (t is the time since periapse, and T is the period of the motion).
Note that there is no way to express E explicitly (i.e. – E = a simple expression). Solutions are found either (i) by making approximations, or (ii) computing numerical solutions to a desired degree of accuracy. This post considers some aspects of the second option: computing numerical solutions.

The task of computing numerical solutions to Kepler’s Equation of Elliptical Motion (KE) falls into the general category of root-finding. The equation is arranged as follows

E - e sin(E) – M = 0

which is of the general form f(x) = 0.
It is now in the form applicable to the general task of root-finding (i.e. – finding values of x that make the function f(x) equal to 0, called the roots).

The task of root-finding is a well-developed field of mathematics and computer science. As far as I know, all root-finding algorithms take an iterative approach to computing a solution to a desired degree of accuracy: first, an initial guess for the root is made and checked, then a closer solution is approximated and checked, and this process is repeated until the desired accuracy is obtained.

A number of numerical techniques exist for computing solutions to KE. Some are rather simple approaches, such as the bisection method; others are more sophisticated. In fact, the computing power of today’s computers provides the capability for even less-than-optimum techniques to compute solutions in the blink of an eye. However, motivated by a desire to be more pedantic, a few aspects of motion defined by KE can be highlighted, and their influence on the approaches taken to KE’s solution can be observed. First, consider the two figures below, which plot the motion of two variables defined by KE over one period, T. Figure 1 is a graph of E and M versus t; E is the solid red line and M is the dashed blue line. Figure 2 is a graph of sin(E) and sin(M) versus t; sin(E) is the solid red line and sin(M) is the dashed blue line.

 

 

Three factors come into play for optimizing an approach for solving KE:
1) the size of the interval;
2) the algorithm;
3) the initial guess for a root.

The Size of the Interval

The interval in which a root exists is usually input to a root-finding algorithm. If the interval over which a solution will be sought is small, the computer program has less work to do. The smaller the interval is to begin with, the faster it will be narrowed to an acceptable result.

Rearranging KE, we get

E = M + e sin(E)

Because sin(E) alternates between 1 and -1, the magnitude of the second term on the right side of the equation is restricted to e. In other words, E will never be larger than

E = M + e

and never be smaller than

E = M - e

In fact, in most real-world situations, e is quite small; the eccentricities of the orbits of most of the planets around the Sun are much smaller than 1. For example, the eccentricity of the orbit of Mars, which has a comparatively large eccentricity, is about 0.093. In other words, to compute a point in the orbit of Mars using Kepler’s Equation, the first guess would already be accurate to one decimal point because the interval itself is already less than 0.1. Not much more can be done to reduce the size of the interval for the purpose of hastening a solution of KE.

The Algorithm

Picking an appropriate algorithm is another factor that determines how fast a solution to KE is computed. As can be seen from the graphs above, the motion is smooth and there are no irregularities; first--and higher-order--derivatives exist and are continuous. This characteristic makes algorithms that use these derivatives applicable to the task of solving
KE--qualifying it for a wide range of fast-converging root-finding algorithms. The Newton-Raphson method, which uses the function and its first derivative, is often used for this purpose. However, care must be taken to ensure the algorithm does not diverge away from the solution; otherwise, it converges to a solution at a quadratic rate. Furthermore, other algorithms offer quick convergence, for example, Halley's Method converges at a cubic rate and yet other algorithms converge at a quartic rate. These faster-converging algorithms offer these capabilities by including higher-order terms of the Taylor series when forming the approximation for successive iterations.

(Note that the utility to which the title of this post links uses a combination of the bisection and secant method instead of one of these more sophisticated methods. When the utility was first posted, I planned to write other utilities and wanted an algorithm general enough to use for other root-finding problems; not one that was particularly well-suited to a specific problem.)

The Initial Guess for a Solution

The initial estimate for a root could be anywhere within the interval; however, the closer the first guess is to the actual root, the fewer iterations the algorithm has to make to find the actual root (i.e. - the quicker the program will execute).
Consider the first half of the period:     0 < t < T/2.
For this portion of the motion, sin(E) > 0 and M <= E <= M + e. So the initial estimate for E, E_o, should be in this interval.

One often-recommended initial value for E_o is the following:

Another approach for picking E_o would be to compute some points of E and save them in a table. An interpolating function could then be created and used to provide E_o values for future computations of E. This approach would be especially applicable if E values for a particular ellipse, or over a certain range only, were to be computed.

Other schemes for picking for E_o exist, but many of them get more complicated than I would like to examine in this post.

To summarize, consideration of the characteristics of motion defined by KE yields approaches for its solution that are much more efficient than a simple approach that, on the contrary, takes no account of KE’s characteristics. By (i) minimizing the initial size of the interval over which a solution is to be sought, (ii) applying an appropriate algorithm, and (iii) making an excellent initial estimation for a root, solutions can be computed extremely fast. A variety of methods--the simple Bisection Method, the Newton-Raphson Method, Halley’s Method, and even faster converging iterative approaches--are available to suit a user’s needs, or simply to satisfy a desire to be more pedantic. In conclusion, although an exact solution of Kepler’s Equation does not yet exist, many powerful tools and technique are available that allow numerical solutions to be computed quickly and accurately.

This post presented a brief overview of existing information on the topic of solving KE. For the sake of intellectual curiosity, I do enjoy exploring KE; any work that could be considered research-ish (i.e. – any work done toward the goal of discovering new information) will be posted on my other blog: RocketScientists.ca.

Edit: 20 October 2012
My WordPress blog on rocketscientists.ca was hacked so I deleted the entire site, had my host reset the account, and started over. The main site is back up, but I do not know when--or if--the blog side of the site will ever be re-posted, so I will copy some of those posts over to this blog. And the hotlink given in the paragraph above has been edited to point to the home page.