Iteration Methods for Solving Kepler's Equation

​
Set the eccentricity e in Kepler's equation M = E - e sin E, where M is the mean anomaly and E is the excentric anomaly.
eccentricity e
0.4
select the result
mode
timing
curves
show
curve
count
error
set the mean anomaly range
n
500
x
L
0.
d
1.05
shift
0.
method parameters
p
-9
i
max
500
steps
3
After Johannes Kepler had discovered the equal area law of planetary motion, he tried to apply this law to determine the position of a planet in an elliptical orbit as a function of time. By ingenious geometrical manipulation, he reduced this task to what has become known as Kepler's problem: for any angle
M
(the mean anomaly) and any number
e∈[0,1)
, determine an angle
E
(the eccentric anomaly) such that Kepler's equation
M=E-esinE
holds. Kepler was unhappy about not being able to find a closed-form solution and had to resort to numerical trial and error, with precomputed numerical tables providing a starting approximation.
This Demonstration shows a contemporary version of Kepler's approach, an iterative algorithm that can solve the problem to any desired accuracy.
The values of
M
form an equally spaced lattice on the
x
axis. In each of the iteration methods considered, a value of
E
is computed for each value of
M
. With "timing", the cumulative computation times for the respective iteration methods are represented as a bar chart. For each iteration method, selecting "curves" shows one of three quantities as a function of
M
: the value of
E
, the number of iterations performed, and the absolute error
|M-E+esinE|
.
Mouse over the control labels and graphics for full explanations.

Details

Snapshot 1: The pure Newton iteration (notice "steps" is 0) shows excessive residual error for
e
close to 1 (see peaks in the blue curve).
Snapshot 2: The Newton iteration preceded by three steps of averaged standard iteration yields tiny residual errors for all
e<1
. Standard iteration (red curve) fails to reach the accuracy goal near
M=π
.
Snapshot 3: This and the following snapshots deal with the larger eccentricity
e=0.9
. Here we show the computation time for 5000 evaluations of
E
for each of the three iteration methods. Due to the consistently low number of iterations, the Newton method (blue curve) is fastest.
Snapshot 4: All methods are seen to keep the residual error below the value
-9
10
(see
p=-9
) set by the termination condition. But the Newton iteration underestimates this by about five orders of magnitude.
Snapshot 5: The large number of required iterations of standard iteration (red curve) explains that this method, although involving the simplest formula, is slowest.
Snapshot 6: The computed values of
E
coincide for all iteration methods in graphical resolution.
Snapshot 7: This is no longer the case if we restrict the allowed iterations to a low value, such as 10 in the present case. Notice that by selecting
x
L
and
d
appropriately, we show the critical region around
M=π
where standard iteration (red) asks for many iterations. Therefore, the error caused by restricting the number of iterations becomes visible.
Consider three different but interrelated iteration schemes. To formulate them compactly, abbreviate Kepler's equation as
E=f(E)
, where
M
and
e
are implicit parameters. The most obvious method, which works very well for, say,
e<0.2
, is nesting function
f
:
(a)
E
i+1
=f(
E
i
)
.
The next one is (a) with a simple provision against oscillation by averaging over two adjacent simple iteration results. This method seems not to be in common use:
(b)
E
⋆
=f(
E
i
)
,
E
⋆⋆
=f(
E
⋆
)
,
E
i+1
=
1
2
(
E
⋆
+
E
⋆⋆
)
.
It always converges, even for
e=1
, but it needs many iterations if
E
is near to an integer multiple of
2π
. The most sophisticated method is defined as a few steps (e.g. three steps) of method (b) followed by the Newton iteration:
(c)
E
i+1
=
E
i
-
E
i
-f(
E
i
)
1-f'(
E
i
)
.
Running all methods with the same termination condition
|
E
i+1
-
E
i
|<ϵ
, the error
|M-E+esinE|
in method (c) is orders of magnitude smaller than in methods (a) and (b). The Demonstration shows that method (c) definitively outperforms (b). The role of (b) thus becomes a robust and simple initialization method for the Newton iteration.

Permanent Citation

Ulrich Mutze
​
​"Iteration Methods for Solving Kepler's Equation"​
​http://demonstrations.wolfram.com/IterationMethodsForSolvingKeplersEquation/​
​Wolfram Demonstrations Project​
​Published: June 23, 2015