Home » Space Flight/Orbital Mechanics » Kepler’s equation, Iterative and Non-Iterative Solver Comparison

Kepler’s equation, Iterative and Non-Iterative Solver Comparison

Kepler’s equation M = E-e*sin(E) M – Mean anomaly [0..pi] rad e – Eccentricity [0..1] E – Eccentric anomaly [rad]

```clear all;
clc;```

Non-Iterative Solver

```tic;
j = 1;
for M = 0:0.01:pi
i = 1;
for e = 0:0.01:1
a   =  (1 - e)*3/(4*e+0.5);
b   =  -M/(4*e+0.5);
y   =  (b^2/4 + a^3/27)^0.5;
x   =  (-b/2 + y)^(1/3) - (b/2 + y)^(1/3);
w   =  x - 0.078*x^5/(1 + e);
E   =  M + e*(3*w - 4*w^3);
% Two times applying 5th-order Newton correction
for k = 1:2
f   =  (E - e*sin(E)- M);
fd  =  1 - e*cos(E);
f2d =  e*sin(E);
f3d = -e*cos(E);
f4d =  e*sin(E);
E = E -f/fd*(1 + f*f2d/(2*fd^2) + f^2*(3*f2d^2 - fd*f3d)/(6*fd^4)+...
(10*fd*f2d*f3d - 15*f2d^3-fd^2*f4d)*f^3/(24*fd^6));
end
E_nit(i,j) = E;
i = i + 1;
end
j = j + 1;
end
t_nit = toc;
fprintf('Calculation time = %4.6f [s] \n',t_nit);
Calculation time = 0.149860 [s]```

Newton Iterative Solver

```tic;
j = 1;
eps = 1E-15;        % Tolerance
for M = 0:0.01:pi
i = 1;
for e = 0:0.01:1
En  = M;
Ens = En - (En-e*sin(En)- M)/(1 - e*cos(En));
while ( abs(Ens-En) > eps )
En = Ens;
Ens = En - (En - e*sin(En) - M)/(1 - e*cos(En));
end;
E_it(i,j) = Ens;
i = i + 1;
end
j = j + 1;
end
t_it = toc;
fprintf('Calculation time = %4.6f [s] \n',t_it);```
`Calculation time = 0.893673 [s]`
```Error = log10(abs(E_it - E_nit));
e = 0:0.01:1;
M = (0:0.01:pi);
contourf (M,e,Error);