Page 1 of 1

Newton's Method for solving the Equation of Kepler

Posted: 11.06.2010, 02:50
by bdm
I'll begin with a brief explanation:

The equation of Kepler is:
E = M + e sin E
It is used for computing the position of an object in its orbit.
This equation is transcendental (has no inverse). To solve for E, various methods exist. One good method is Newton's method:
E(i+1) = E(i) + (M + e sin E(i) - E(i)) / (1 + e cos E(i))

Many textbooks suggest using E[0] = M as a starting value, but then state that the equation has pathological behaviour with slow convergence if e is large and M is small. This includes Meeus' Astronomical Algorithms. Then the texts suggests various ways of preventing the pathological behaviour, such as limiting the value of Newton's Method, or passing the result through sin and then arcsin.

The pathological behaviour can be seen in this plot of M versus e, with number of iterations plotted:
Image

The top left corner magnified 4 times, showing the fractal nature of this function:
Image

However, all textbooks that state that the pathological behaviour must be worked around are incorrect.

It turns out that the fault lies with the choice of E(0) = M as a starting value. The pathological behaviour only occurs whenever E(i) is less than E, M < 1 (radian), and e > 0.95 (approximately).

To derive a better starting value, choose a value for E(0) that is greater than 1. When one chooses E(0) = pi/2, E(1) = M + e. We can then skip this step and just use E(0) = M + e.

It turns out that this small improvement eliminates the pathological behaviour entirely. With this improvement, Newton's Method is entirely reliable for all 0 <= M <= pi, and all 0 <= e < 1.

Here is a diagram showing this:
Image

The obvious diagonal line is also of interest. This corresponds to E = M + e, where the function converges very quickly.