Simulating an orbit - numerically solving $M(E) = E + \sin(E)$

Question

Well for a given kepler orbit (which is a ellipse) $0 \leq e < 1$. There are several functions to describe the motion of an object.

$$r(\nu) = \frac{a (1 - e^2)}{1 + e \cos(\nu)}$$ Where $a$ is the semi major axis, $e$ the eccentricity, $\nu$ the true anomaly. And $r$ the distance from the focus point.

One of the basic functions is to use the mean anomaly $M$. The mean anomaly is an object is basically defined as an angle that grows linearly with time ($n$, dependent on the orbital period) $$M(t) = M_0 + nt$$

The mean anomaly can be expressed as a function of the eccentric anomaly ($E$) as: $$M = E - e \sin(E)$$ The problem is quite obvious now: the mean anomaly doesn't represent any geometric value, and to convert it to one, a function has to be solved numerically. (Simple geometry allows conversion of $E$ to and from $\nu$).

Now using a method such as Newton's this can be solved - and luckily the mean anomaly doesn't depend on the eccentric anomaly so the truncation error isn't that important.

However I wonder if there is a more clever way, considering the derivatives are very well defined:

$$M'(E) = 1-e\cos(E)$$ $$M''(E) = e\sin(E)$$ $$M'''(E) = e\cos(E)$$ $$M^4(E) = -e\sin(E)$$

And this repeats onwards then. So can a method with a better convergence rate than newton's method be described for this problem? And further more ($E_0$ means start of the root finding algorithm) - is initializing as follow good (will it guaranteed converge): $$E_0 = M$$ -Starting the root finding at by stating the mean anomaly is equal to the eccentric anomaly, which is the case for the apo/periapsis.

Answer 1

It seems to me that your post addresses two points

• convergence of Newton method
• methods converging faster than Newton method

Concerning the first point, assuming the starting point to be "reasonable", Newton method would converge without any overshoot of the solution if $$M(E_0) M''(E_0) >0$$ (if I properly remember, this was established by Darboux).

Concerning methods converging faster, the simplest would be Halley (cubic convergence) or Householder (quartic convergence).

Higher order schemes have been proposed (have a look at http://numbers.computation.free.fr/Constants/Algorithms/newton.html)

For illustration purposes, let us use $e=5$ and $E_0=3$.

Newton iterates will be $$E_1=2.61438412994111$$ $$E_2=2.59582212992534$$ $$E_3=2.59573908134712$$ $$E_4=2.59573907964980$$ which is the solution for fifteen significant figures.

Halley iterates will be $$E_1=2.60556706231111$$ $$E_2=2.59573931852909$$ $$E_3=2.59573907964980$$

Householder iterates will be $$E_1=2.59642488378724$$ $$E_2=2.59573907964981$$

Edit

By the way, for a good estimate of the solution, you could use the splendid 1400 year old approximation $$\sin(x) \simeq \frac{16 (\pi -x) x}{5 \pi ^2-4 (\pi -x) x}$$ proposed by Mahabhaskariya of Bhaskara I, a seventh-century Indian mathematician. Then, the estimate is just obtained at the price of a quadratic equation. Applied to $e=5$, this would give $E_0=2.5956$ from which convergence will be almost immediate.

Answer 2

You can express $E$ as a series in powers of $e$ using the Lagrange reversion formula:

$$ E = M + \sum_{k=1}^\infty \dfrac{e^k}{k!} \left(\dfrac{\partial}{\partial M}\right)^{k-1} \sin^k(M) $$