Solving Kepler equation for true or eccentric anomaly

Question

Is there any reason to always solve the Kepler equation for the eccentric anomaly, $E$, instead of the more meaningful (at least to me) true anomaly, $\theta$?

Solving for the eccentric anomaly usually would mean calculating the true anomaly from its solution and maybe even the initial guess of the eccentric anomaly from an initial guess of the true anomaly. So this would require more calculations, however the Kepler equation expressed in the eccentric anomaly is simpler and thus less expensive to calculate. When I search for solving Kepler's equation I only have some approaches which use the eccentric anomaly, however I can imagine that when the error tolerance of the solution is not very tight or the initial guess is good it might be cheaper to solve for the true anomaly. Or are there other reasons for solving for the eccentric anomaly?

The relevant equations, using Newton's method, are: $$ E = 2 \tan^{-1}\left(\sqrt{\frac{1-e}{1+e}}\tan\frac{\theta}{2}\right), $$ $$ M = E - e \sin{E}, $$ where $M$ is the mean anomaly and $e$ is the eccentricity. Such that using Newton's method results in $$ E_{n+1} = E_n - \frac{E_n - e \sin{E_n} - M}{1 - e \cos{E_n}}. $$ And for the true anomaly, $$ M = 2\tan^{-1}\left(\sqrt{\frac{1-e}{1+e}}\tan{\frac{\theta}{2}}\right)-\frac{e\sqrt{1-e^2}\sin{\theta}}{1+e\cos{\theta}}, $$ $$ \theta_{n+1} = \theta_n - \left(1 + e \cos{\theta_n}\right)^2 \frac{2\tan^{-1}\left(\sqrt{\frac{1-e}{1+e}}\tan{\frac{\theta_n}{2}}\right)-\frac{e\sqrt{1-e^2}\sin{\theta_n}}{1+e\cos{\theta_n}} - M}{\sqrt{(1 - e^2)^3}}. $$

Answer 1

The most robust way of answering this is to benchmark it. Failing that, there are several things to note (roughly in order of importance).

First, the most cheap floating-point operations on a modern CPU are addition and multiplication (both are equally fast; same as fused multiply-adds when available). Division is much slower (by a factor of ~20), trigonometric functions are also slower (~200), square roots and logs similar. In fact, sometimes (depends on architecture) trig functions, roots, logs are implemented in a library. The best reference for this is usually the optimization guide published by the whichever company made your CPU. So as a rule of thumb you want to minimize special functions first, then divisions, then multiplications and additions.

So by this measure the second formula is much worse: the there is one inverse trig function, one square root (it would be wasteful to implement your formula as written, with three square roots), and one sincos (since they are expensive it is better to evaluate sincos once and then write things like $\tan\frac\theta2$ in terms of those). Your first formula involves the other special functions, apart from sincos, only once at the end.

Second, you should not stop at Newton's method. For example, if you compare Newton's method with Halley's method, $$ \frac{e \sin (x)-e x \cos (x)+M}{1-e \cos (x)}, \qquad x-\frac{(e \cos (x)-1) (e \sin (x)+M-x)}{e^2+e (M-x) \sin (x)-2 e \cos (x)+1}, $$ both evaluate the expensive trigonometric functions at the same arguments, but with a little more algebra Halley's method can converge faster; the extra algebra might not outweigh the savings of evaluating trigonometric fewer times. So Halley's method (and other iterative methods) also need to be checked.

Third, you can precompute some things. For example, if you start by reducing the argument to the range $0<E<2\pi$, you can experimentally, in advance, find the maximum number of iterations taken by the method. Since the number of iterations is fairly small (depending on precision) and mispredicted branches can be expensive, it may also make sense to unroll the iteration loop by hand to a fixed sufficient number of iterations. Unrolled code can also more easily benefit from vectorization.

Fourth, it is difficult to intuitively predict which optimizations will do the best. And since this is just one fairly simple equation, it is probably best to benchmark many different approaches and find the best. When doing this, also consider what optimization flags your compiler supports for floating-point arithmetic (some of them are interesting and important, many people know about -ffast-math, but it actually decomposes into different helpful or harmful optimization flags; here is gcc's list for example). Another thing to do is to look at the assembly output of your compiler to see which CPU instructions it actually ends up using.

Fifth, if you need to solve this equation many times for the same eccentricity $e$, it is possible to rewrite the problem. If you consider the function $E = E(M)$, on the range $[0,2\pi]$, given by the solutions of the equation for fixed known $e$, you can approximate the function $E(M)$ using, for example, Chebyshev series, which takes only a small number of evaluations of $E(M)$ (which can be done with any root-finding method). Once you have a sufficiently close approximation, which might be a 20-term Chebyshev series or something like that, you can evaluate that later without needing to solve the equation again.

Answer 2

This question now has quite some age, but its subject crops up so repeatedly, that perhaps the following answer and references may still be helpful and not out of place.

(a) 'Is there a reason to ... solve for ... the eccentric anomaly' (when the true anomaly is more meaningful)?'

Practically all of the numerical difficulty is in 'Kepler's equation', i.e. the part of the calculation for an undisturbed elliptical orbit that starts with the mean anomaly and computes the eccentric anomaly. The onward calculation of true anomaly from eccentric, by contrast, can be simple enough and numerically well-behaved when properly arranged. (A better-behaved formula than the traditional 'tangent half-angle' formula quoted in the question is given by R Broucke & P Cefola (1973), [http://adsabs.harvard.edu/cgi-bin/nph-data_query?bibcode=1973CeMec...7..388B].) Direct series evaluations of true anomaly do exist, of course, but the coefficients are usually laborious to compute, except where that can be done once for all, e.g. when only a few predetermined eccentricity values will be used. There are also a few iterative procedures that seem to go direct from mean to true anomaly, but they don't entirely avoid the eccentric anomaly, and it's arguable they only mix and mask difficulties without gaining simplicity or efficiency.) So it's usual to tackle separately the Kepler-equation part of the problem.

(b) There are many published good procedures for solving Kepler's equation, thus arguably little point now in devising one afresh. The subject is entirely 400 years old, but it received renewed interest in the 1960s on, when astrodynamics developed increased demand for procedures for orbit-tracking and control. These demand (as intermediates) unattended automatic computation of many thousands of instances of Kepler's equation. What emerged then were needs not only for computational efficiency, but also for the avoidance of added kinds of numerical misbehavior, such as misconvergence or failure to converge from a minority of pathological input values. Here is a review and a brief selection of solutions from that recent era.

** There's a good explanatory survey of methods up to the mid-80s, and of the problems to be met, (plus some coded Fortran functions), in: A W Odell, R H Gooding, 'Procedures for Solving Kepler's Equation' [http://adsabs.harvard.edu/cgi-bin/nph-data_query?bibcode=1986CeMec..38..307O].

** A very useful paper is that of A Nijenhuis (1991), 'Solving Kepler's equation with high efficiency and accuracy' [substitute the following bibcode in the URL given above: bibcode=1991CeMDA..51..319N], it's accompanied by a complete and workable implementation (in Pascal, from 1991, but easily portable e.g. to C). Methods of generating a starting- approximation are selected according to the input values of M (radians) and e (#), and a recently-devised cubic-approximation (Mikkola, 1987) is used for starting values in the specially-problematic area near to the singularity at (M,e) = (0, 1).

** Another modern method, perhaps even potentially better than the foregoing, is described (without code example) by F L Markley (1995), 'Kepler Equation Solver' [substitute the following bibcode in the URL given above: bibcode=1995CeMDA..63..101M].

** There are various methods that claim to avoid the calculation of trig functions. One class avoids them at run-time, but the cost is a precomputed lookup table (e.g. S A Feinstein, C A McLaughlin, 'Dynamic discretization method for solving Kepler's equation', [substitute the following bibcode in the URL given above: bibcode=2006CeMDA..96...49F]. Another method ingeniously does avoid computing trig functions; it achieves that by evaluating sin E (via sin E/3) rather than E itself:
S Mikkola, 'A cubic approximation for Kepler's equation', [substitute the following bibcode in the URL given above: bibcode=1987CeMec..40..329M].

Hopefully there's something here to respond to most variants of this oft-recurring question.

Answer 3

This may be of interest to you: a paper recently (December 2014) published presents a method for computing solutions to Kepler's Equation analytically:

http://www.scirp.org/journal/PaperInformation.aspx?PaperID=52772

I haven't worked through the paper myself (it involves a two-dimensional Laplace Transform), but it may be relevant to your work.

Answer 4

You need to solve numerically the eccentric anomaly using several iterations, while the true anomaly is calculated straight from the eccentric anomaly, without iterations using : theta = atan(sqrt((1+e)/(1-e)*tan(E/2)^2))*2

If you try to solve the true anomaly directly you will need the derivative of a much more complex function to apply the newton-raphson method.

Maybe I didn't understand the question properly.