# Does the Newton-Rhapson solution to Kepler's equation diverge?

I have the following simple solution to the Kepler's equation which uses newton-rhapson iterative technique. I have read various papers about the initial guess and divergence near an eccentricty of 1, but I run the following test code and never get any divergent solutions. Am I missing something or does this solution always converge?

The following is a matlab script with a helper function:

n = 1e6;
FM = rand(n,1)*2*pi;
ECC = rand(n,1)*.9999999;

ecan = zeros(n,1);
E = zeros(n,1);
tic
for i = 1:n
[E(i)] = kepler(FM(i),ECC(i));
end
toc
close all
plot(E - FM - ECC.*sin(E));

%kepler solves kepler equation for conversion from mean to eccentric
%anomaly using newton-rhapson itteration

function [E] = kepler(M,ecc)

%Tolerance for function
eps = 1e-8;

%Set initial guess
En  = M;

Ens =    En - (En - ecc*sin(En)- M)/(1 - ecc*cos(En));

% begin itteration (newton-rhapson)
while (abs(Ens-En) > eps)

En = Ens;

Ens = En - (En - ecc*sin(En) - M)/(1 - ecc*cos(En));

end

E = Ens;

end

• Try an eccentricity of 0.999 and a mean anomaly of 0.15 or $2\pi$-0.15, and look at the values of Ens as the algorithm bounces around to a solution. You might have to adjust the 0.15 a bit to see truly bad behavior; I used python rather than Matlab. The intermediate steps bounce all over the place (up to $10^{18}$ !!) before the algorithm gets lucky and hits a window of convergence. Starting with an initial guess of $\pi$ avoids the bouncing around. – David Hammen Aug 19 '20 at 11:30
• You are thinking of diverging in terms of running off to infinity. Another way to look at divergence is that the algorithm takes a step far, far away from the solution. Newton-Raphson is notorious for doing so. That the algorithm oftentimes eventually does converge to the solution when starting with an initial guess of E=M for solving Kepler's equation is more a case of luck and the vagaries of representing the reals with IEEE floating point. – David Hammen Aug 20 '20 at 12:01

This question is really more of a general numerical analysis question, rather than a space exploration question (as there's no reason Kepler's equation is special in this regard), but I'll provide this guidance anyway:

The divergence you are expecting comes from the fact that the $$1-e \cos En$$ term gets close to zero as the eccentricity approaches 1, meaning the corrective term in Newton's method is large relative to the initial value.

Given that you have not found a divergence case yet, I suspect it's because you are using random values for your eccentricity, and you are just not getting close enough to trigger divergence -- not also that saying Newton's method diverges does not necessarily mean it diverges for every choice of start, it just means that it will do so for some starting values.

In addition to mathematical divergence, there's also the aspect of ill-conditionedness, where machine precision errors compound upon themselves to overwhelm the actual true values.

Try, instead of choosing random values of eccentricity, to force it to get progressively closer to 1 and see what happens.

There are some combinations of a large eccentricity and a bad guess for the initial value of the eccentric anomaly where a Newton-Raphson iterator of Kepler's equation appears to diverge. An initial guess of eccentric anomaly being equal to the mean anomaly seems completely reasonable and works quite nicely when eccentricity is not large. But when eccentricity is close to one, that is turns out to not be the case. There are places where this choice as an initial guess results in divergence.

There is an easy solution: Don't make a bad guess for the initial value of the eccentric anomaly. Quoting the abstract of The Convergence of Newton–Raphson Iteration with Kepler's Equation,

Conway (Celest. Mech. 39, 199–211, 1986) drew attention to the circumstance that when the Newton–Raphson algorithm is applied to Kepler's equation for very high eccentricities there are certain apparently capricious and random values of the eccentricity and mean anomaly for which convergence seems not to be easily reached when the starting guess for the eccentric anomaly is taken to be equal to the mean anomaly. We examine this chaotic behavior and show that rapid convergence is always reached if the first guess for the eccentric anomaly is π. We present graphs and an empirical formula for obtaining an even better first guess. We also examine an unstable situation where iterations oscillate between two incorrect results until the instability results in sudden convergence to the unique correct solution.