.There have been many variations on the theme of solving Kepler's Equation for elliptic orbital motion. I thought I had found a new method, but someone else had already found it about 5 years ago. I am grateful to the editor of Celestial Mechanics and Dynamical Astronomy for sending me a photocopy of the original paper by F. Landis Markley where the same approximation was also employed. His version giving truly excellent initial estimates.
I think this new method deserves to be better known as it saves considerable time on iterative methods that start from poor initial choices. Kepler's equation looks innocent enough
E = M + e sin E
But it becomes troublesome when e ~ 1 and M ~ 0.
The new method is based on a Pade approximation to sin(x) and was published in a paper titled Kepler Equation Solver, by F. Landis Markley in Celestial Mechanics and Dynamical Astronomy vol 63, 101-111, 1995. The first part of the paper deals with getting a good initial estimate, and the second on refining it to full double precision.
The normal Taylor series for sin(x) is
sin(x) = x - x^3/3! + x^5/5! - x^7/7! + ...
sin(x) = x - x^3/6 + x^5/120
This is very accurate for small x, but rapidly deviates from the true sin(x) function as x increases.
The Pade approximation for sin(x) is
Pade{sin(x)} = x (60-7x^2)/(60+3x^2)
This exactly matches the coefficients of the truncated 5th order polynomial Taylor series, but has much better global convergence and mimics sin(x) rather well. This is best seen graphically.
.
Landis actually started with the polynomial representation:
sin(E) ~ E - E^3/6 + E^5/12a
Where a is a parameter which is about 10, but can be varied slightly to make sin(0) = sin(pi) = 0.
This yields a new Pade approximant for sin(E)
Pade{ sin(E) } = E (6a - (a-3)E^2)/(6a + 3E^2)
Substituting this into Kepler's Equation gives a cubic equation:
[ 3(1-e) + ae] E^3 - 3M E^2 + 6a(1-e)E - 6aM = 0
My own derivation was for a constant a=10 and gave a maximum relative error of 4% at M=180 degrees, and required an empirical correction to be applied for large values of M. But by adjusting the parameter a Landis was able to get the relative error everywhere less than 0.03%.
Some ingenious algebra in the original paper shows that by taking a(e,M) as a slowly varying function of the known values e, M an excellent level of agreement can be obtained between the true solution and the truncated cubic approximation over the full range of elliptic motion.
a( e, M ) = (3 pi^2 + 1.6 pi ( pi - |M| )/(1+e))/(pi^2-6)
The solution of a cubic equation for it's roots is well covered in standard texts, but there is a small trick needed for Keplers equation to avoid some small answers being the result of the difference of two similar numbers.
The final refinement step is possible because the answer is close enough to the true solution to solve by explicit differences. There is no need to iteratively solve the equation or work out sin(E) or cos(E) more than once!
Again given the function to solve for a root:
f(E) = E - e sin(E) - M = 0
One standard approach to solving it is via Newton Raphson iteration.
En= E - f(E)/f '(E)
Where f '(E) = 1 - e cos(E)
This works provided the initial guess is close to the required root. But it can go haywire if it is started from too far away from the root. A problem which tends to arise when e ~ 1 and M ~ 0. It is instructive to try and solve Kepler's equation this way for the case M = 0.1, e = 0.9999 with the initial guess E = M in a spreadsheet.
If you trust downloading a spreadsheet then here is an Excel version which allows you to try various values of e, M and see how simple iteration, Newton-Raphson and Halley's method converge on the answer. When e ~ 1 it is pretty obvious that E = M is a very poor starting approximation and will lead to instability. A starting guess that is computationally cheap and not too bad is :
E = M + ( sign(M) - M/pi) e/2
Halley's third order method is somewhat better behaved in these awkward cases and involves little extra effort for a considerable increase in robustness.
En= E - f(E)/ (f '(E) - f(E) f "(E)/(2 f (E)))
Where f "(E) = e sin(E)
This formulation converges slightly faster, and is much less inclined to bounce chaoticaly if started with an unsuitable initial guess. The detailed treatment in the original paper extends this difference approach to fifth order which generates full double machine precision results
