Solving Kepler's equation — Newton-Raphson iteration

WHY do we even need this?

Given a time $t$, we want to know where a satellite is on its ellipse. The clean chain is:

$\text{time } t \;\longrightarrow\; M \;\longrightarrow\; E \;\longrightarrow\; \nu \;\longrightarrow\; (r,\theta)$

• $M$ = mean anomaly — a fictitious angle that grows uniformly with time.
• $E$ = eccentric anomaly — the geometric angle on the auxiliary circle.
• $\nu$ = true anomaly — the actual angle of the body from the focus.

The trouble sits between $M$ and $E$.

HOW is Kepler's equation derived (from scratch)?

We need to earn $M = E - e\sin E$, not just quote it. It comes from Kepler's 2nd law (equal areas in equal times).

Setup. Draw the ellipse (semi-axes $a,b$) and its circumscribing auxiliary circle of radius $a$. The eccentric anomaly $E$ is the angle at the centre to the point on the circle above the body.

Step 1 — Area swept scales with time. By the 2nd law, area swept from perihelion is proportional to elapsed time: $\frac{A_{\text{ellipse sector}}}{A_{\text{ellipse total}}} = \frac{t-t_p}{T} = \frac{M}{2\pi}.$ Why? Equal areas / equal times means fraction of area = fraction of period, and $M/2\pi$ is exactly the fraction of the period elapsed (since $M=n(t-t_p)$ and $nT=2\pi$).

Step 2 — Squash the circle into the ellipse. Every $y$-coordinate on the circle is scaled by $b/a$ to land on the ellipse. So every area on the circle is scaled by $b/a$: $A_{\text{ellipse sector}} = \frac{b}{a}\,A_{\text{circle sector}}.$

Step 3 — Area of the circle sector. The circle "sector" swept from perihelion, measured from the focus, splits into a pie-slice about the centre minus a triangle: $A_{\text{circle sector}} = \underbrace{\tfrac12 a^2 E}_{\text{pie slice, angle }E} - \underbrace{\tfrac12 a^2 e\sin E}_{\text{triangle to focus}}.$ Why the triangle? The focus is offset from the centre by $c = ae$. The triangle formed by centre, focus and the circle point has base $ae$ and height $a\sin E$, area $\tfrac12(ae)(a\sin E)$.

Step 4 — Assemble. Total ellipse area $=\pi a b$. Plug in: $\frac{\tfrac{b}{a}\!\left(\tfrac12 a^2 E - \tfrac12 a^2 e\sin E\right)}{\pi a b} = \frac{M}{2\pi}.$ The $a$'s and $b$'s cancel beautifully: $\frac{\tfrac12 a b (E - e\sin E)}{\pi a b} = \frac{E-e\sin E}{2\pi} = \frac{M}{2\pi} \;\Rightarrow\; \boxed{M = E - e\sin E}.$

WHY Newton-Raphson? (the 80/20 method)

We must solve $M = E - e\sin E$ for $E$ given $M,e$. No closed form exists. Define the root function:

$f(E) = E - e\sin E - M.$

We want $f(E)=0$. Newton-Raphson is the workhorse: fast (quadratic convergence), needs only $f$ and $f'$, and here $f'$ is trivial.

Why $f'=1-e\cos E$? Differentiate $f(E)=E - e\sin E - M$ term by term: $\frac{d}{dE}E=1$, $\frac{d}{dE}(-e\sin E)=-e\cos E$, $M$ is constant. Since $e<1$, $f'>0$ always — the function is strictly increasing, so the root is unique and Newton is well-behaved.

Worked Example 1 — moderate eccentricity

Given $M = 0.8$ rad, $e = 0.3$. Find $E$.

Seed: $E_0 = M + e\sin M = 0.8 + 0.3(0.7174) = 1.0152$. Why this seed? It pre-corrects for the $-e\sin E$ term, landing us close to the root so fewer iterations are needed.

Iteration 1: $f = 1.0152 - 0.3\sin(1.0152) - 0.8 = 1.0152 - 0.2555 - 0.8 = -0.0403$ $f' = 1 - 0.3\cos(1.0152) = 1 - 0.1585 = 0.8415$ $E_1 = 1.0152 - \frac{-0.0403}{0.8415} = 1.0631$ Why this step? We move $E$ up because $f<0$ meant our guess was too small.

Iteration 2: $f = 1.0631 - 0.3\sin(1.0631) - 0.8 = 1.0631 - 0.2626 - 0.8 = 0.0005$ $E_2 = 1.0631 - \frac{0.0005}{0.8558} \approx 1.06266$

Already $|f|\sim 10^{-4}$; one more step gives machine precision. $E \approx 1.0627$ rad. Notice $E>M$ — the body has sped past the uniform angle, as expected.

Worked Example 2 — high eccentricity (where seeds matter)

Given $M = 0.1$ rad, $e = 0.9$.

Seed choice — why care? For high $e$, $E_0 = M$ can converge slowly or oscillate. Use $E_0 = \pi$ only if $M$ near $\pi$; here $M$ is small so seed $E_0 = M + e\sin M = 0.1 + 0.9(0.0998) = 0.1898$.

Iter 1: $f = 0.1898 - 0.9\sin(0.1898) - 0.1 = 0.1898 - 0.1700 - 0.1 = -0.0802$ $f' = 1 - 0.9\cos(0.1898) = 1 - 0.8838 = 0.1162$ $E_1 = 0.1898 - (-0.0802/0.1162) = 0.8801$

Why the big jump? Near high $e$, $f'$ is small (denominator $\approx 0.12$), so Newton takes a large step. This is exactly why high-$e$ orbits need care — the correction is violent.

Iter 2: $f = 0.8801 - 0.9\sin(0.8801) - 0.1 = 0.8801 - 0.6909 - 0.1 = 0.0892$ $f' = 1 - 0.9\cos(0.8801) = 1 - 0.5731 = 0.4269$ $E_2 = 0.8801 - 0.0892/0.4269 = 0.6712$

Iter 3: $f = 0.6712 - 0.9\sin(0.6712) - 0.1 = 0.6712 - 0.5595 - 0.1 = 0.0117$; $f'=0.2914$; $E_3=0.6310$ Continuing converges to $E \approx 0.6249$ rad. Slower — high $e$ is genuinely harder.

From $E$ to true anomaly $\nu$

Once $E$ is found: $\tan\frac{\nu}{2} = \sqrt{\frac{1+e}{1-e}}\,\tan\frac{E}{2}, \qquad r = a(1-e\cos E).$ Why $r=a(1-e\cos E)$? Project the geometry: the body's distance from the focus, expressed through the eccentric anomaly, gives this compact form (max at aphelion $E=\pi$: $r=a(1+e)$; min at perihelion $E=0$: $r=a(1-e)$ — check it).

Flashcards

Connections

• Kepler's Laws of Planetary Motion — the 2nd law births this equation.
• Eccentric Anomaly and the Auxiliary Circle — geometric meaning of $E$.
• True Anomaly and the Orbit Equation — next step after solving for $E$.
• Newton-Raphson Method — general root-finding, convergence conditions.
• Mean Motion and Orbital Period — where $n$ and $M$ come from.
• Fixed-Point Iteration for Kepler — the simpler (slower) alternative $E_{k+1}=M+e\sin E_k$.

Concept Map

Hinglish (regional understanding)