La ecuación de Kepler y su resolución numérica

1. Motivación

La lección anterior demostró que la órbita es una cónica, pero no dijo dónde está el cuerpo en cada instante. La ley de las áreas, $r^2\dot\nu = h$, es una ecuación diferencial en $\nu(t)$ sin solución elemental: Kepler la resolvió en 1609 introduciendo un ángulo auxiliar, y el resultado es una ecuación trascendente que llevó tres siglos de métodos numéricos (Colwell 1993). Toda efeméride, desde las tablas de Kepler hasta las de JPL, pasa por resolverla.

2. Planteamiento

Para el problema de dos cuerpos reducido con $\mu = G(m_1+m_2)$, la conservación del momento angular $h = r^2\dot\nu$ y la ecuación de la órbita conducen a la anomalía excéntrica $E$, definida por

Geométricamente, $E$ es el ángulo, visto desde el centro de la elipse, del punto de la circunferencia auxiliar de radio $a$ que se proyecta perpendicularmente al eje mayor sobre la posición del cuerpo. En coordenadas perifocales (origen en el foco, $x$ hacia el periapsis) eso da

y de ahí $r = \sqrt{x^2+y^2} = a(1-e\cos E)$. Calculemos $h = x\dot y - y\dot x$ con $\dot x = -a\sin E\,\dot E$, $\dot y = a\sqrt{1-e^2}\cos E\,\dot E$:

Como $h = \sqrt{\mu a(1-e^2)} = na^2\sqrt{1-e^2}$ con $n = \sqrt{\mu/a^3}$, queda $(1-e\cos E)\,dE = n\,dt$ e, integrando desde el paso por el periapsis $\tau$:

(equivalentemente: sustituir en $h\,dt = r^2d\nu$ e integrar). $M$ es la anomalía media: el ángulo que habría recorrido un cuerpo ficticio con movimiento uniforme $n$. La ecuación de Kepler es trascendente en $E$: dado el tiempo, hay que resolverla numéricamente.

3. Newton–Raphson

Sea $f(E) = E - e\sin E - M$. Entonces $f'(E) = 1 - e\cos E \ge 1-e > 0$, así que $f$ es estrictamente creciente y la raíz es única. La iteración

converge cuadráticamente cerca de la raíz. Para $e \to 1$ y $M \to 0$ la semilla importa: usamos la de Danby, $E_0 = M + 0.85\,e\,\operatorname{sgn}(\sin M)$. Implementación en src/lib/math/kepler.ts.

3.1 Convergencia cuadrática con cota explícita

Proposición. Si $E^*$ es la raíz, $|E_{k+1} - E^*| \le C\,|E_k - E^*|^2$ con $C = \dfrac{e}{2(1-e)}$.

Demostración. Por Taylor con resto de Lagrange alrededor de $E_k$: $0 = f(E^*) = f(E_k) + f'(E_k)(E^*-E_k) + \tfrac12 f''(\xi)(E^*-E_k)^2$ para algún $\xi$ entre $E_k$ y $E^*$. Dividiendo por $f'(E_k)$ y usando la definición de $E_{k+1}$: $E^* - E_{k+1} = -\dfrac{f''(\xi)}{2f'(E_k)}(E^*-E_k)^2$. Con $f''(E) = e\sin E$, $|f''| \le e$ y $f' \ge 1-e$ se obtiene la cota. $\blacksquare$

La convergencia es global además de cuadrática: como $f$ es creciente y convexa en $(0,\pi)$ (allí $f'' = e\sin E > 0$), Newton desde cualquier $E_0 \ge E^*$ produce una sucesión monótona decreciente hacia la raíz. La cota $C$ explica el problema con $e\to1$: para $e = 0.99$, $C = 49.5$, y la convergencia cuadrática solo empieza cuando $|E_k - E^*| \lesssim 1/C = 0.02$ rad.

3.2 El papel de la semilla

Conteo de iteraciones hasta $|\Delta E| < 10^{-14}$, barriendo $M\in(0,\pi)$ (calculado con el algoritmo de kepler.ts):

El caso patológico es $e = 0.99$, $M = 0.01$ con $E_0 = M$: $f'(E_0) = 1 - 0.99\cos0.01 \approx 0.010$ y el primer paso salta a $E_1 = 0.995$, lejísimos de la raíz $E^* = 0.342$; después Newton baja monótonamente ($0.660$, $0.462$, $0.368$, $0.344$, $0.3423$) y necesita 9 iteraciones. La semilla de Danby, $E_0 = M + 0.85e = 0.852$, está ya del lado convexo y converge en 8; para $M = 0.1$ la diferencia es 21 frente a 5. La regla general: arrancar por encima de la raíz garantiza monotonía.

4. Las tres anomalías

De $x = a(\cos E - e) = r\cos\nu$ e $y = a\sqrt{1-e^2}\sin E = r\sin\nu$, con $r = a(1-e\cos E)$:

y la fórmula del ángulo mitad del planteamiento; en la práctica se usa $\nu = 2\operatorname{atan2}\big(\sqrt{1+e}\sin\tfrac E2,\ \sqrt{1-e}\cos\tfrac E2\big)$, que da el cuadrante correcto sin ramas (trueAnomaly en kepler.ts). Para $e$ pequeña conviene la ecuación del centro, el desarrollo de $\nu - M$ en serie de Fourier de $M$ (Problema 5):

Verificación con Marte (J2000). Horizons da $M = 19.35648275°$ y $e = 0.0933151015$. Newton–Raphson converge en 4 iteraciones a $E = 21.2984971°$, y $\nu = 23.3331966°$, exactamente la anomalía verdadera TA $= 23.33319664°$ publicada por Horizons. La ecuación del centro truncada a orden $e^2$ da $\nu - M = 3.934°$ frente al valor exacto $3.977°$: error de $0.043°$, que es justo el término siguiente, $\tfrac{13}{12}e^3\sin3M = 0.043°$.

5. Órbitas hiperbólicas y parabólicas

Para $e > 1$ se define la anomalía hiperbólica $F$ por $r = a(e\cosh F - 1)$ con $a > 0$ el semieje de la hipérbola, y el mismo cálculo (Problema 3) conduce a

Aquí $f(F) = e\sinh F - F - M_h$ tiene $f' = e\cosh F - 1 \ge e - 1 > 0$ y Newton converge sin dificultad; una buena semilla es $F_0 = \ln\big(2M_h/e + 1.8\big)$ para $M_h$ grande. Para $e = 1$ (cometas casi parabólicos, el caso mal condicionado de la lección 0) no hay anomalía excéntrica: con $D = \tan(\nu/2)$ y la distancia al periapsis $q$, la ecuación de Barker es cúbica y se resuelve en forma cerrada:

Para evitar la discontinuidad entre los tres regímenes, los propagadores modernos usan variables universales (Curtis §3.7; Vallado §2.2.5), una sola ecuación $\chi$ válida para toda $e$.

6. Análisis numérico y de condicionamiento

Criterio de parada. Parar cuando $|\Delta E_k| < \epsilon$ es adecuado porque, con convergencia cuadrática, el error tras el último paso es $\sim C|\Delta E_k|^2$, mucho menor que $\Delta E_k$. Parar por $|f(E_k)| < \epsilon$ es engañoso para $e\to1$: $|f| = |f'|\,|E - E^*| \approx (1-e)|E-E^*|$, así que un residuo de $10^{-14}$ puede esconder un error de $10^{-12}$ en $E$.

Condicionamiento. La solución depende de los datos con $\dfrac{\partial E}{\partial M} = \dfrac1{1-e\cos E}$: cerca del periapsis de una órbita con $e = 0.99$ un error $\delta M$ se amplifica por $100$. No es un defecto del método sino del problema: el cuerpo se mueve tan deprisa en el periapsis que un pequeño error de tiempo equivale a un gran error angular. Para $\nu$ la amplificación se compone con $d\nu/dE = \sqrt{1-e^2}/(1-e\cos E)$, de nuevo $\sim 1/(1-e)$ en el periapsis.

Coste. Con la semilla de Danby, 4–5 iteraciones bastan para cualquier $e < 0.99$ a precisión de máquina; cada una cuesta un seno y un coseno. Markley (1995) y otros dan solvers sin iteración (aproximación inicial cúbica más una corrección de Halley de orden alto) con error $< 10^{-15}$ en una sola evaluación, útiles cuando hay que resolver millones de veces (vectorización, simulaciones de poblaciones). El simulador usa Newton con Danby, que para una sola órbita es más que suficiente.

7. Simulador

Observa, con $e = 0.9$, cómo los puntos equiespaciados en $M$ (tiempo) se apelotonan en el apoapsis y se separan en el periapsis: es la figura de $dE/dM = 1/(1-e\cos E)$. Comprueba que el tiempo entre periapsis y $\nu = 90°$ es una fracción del periodo mucho menor que $1/4$ (Problema 6).

Problemas

1. (Demostración) Demuestra que $f''(E) = e\sin E$ y completa la demostración de la cota $|E_{k+1}-E^*| \le C|E_k - E^*|^2$ con $C = e/(2(1-e))$. Calcula $C$ para $e = 0.5$, $0.9$, $0.99$ y el radio de convergencia cuadrática $\sim 1/C$.

2. (Computacional) Compara el número de iteraciones con semilla $E_0 = M$ y con la de Danby para $e\in\{0.5, 0.9, 0.99\}$ y $M$ en una malla de $(0,\pi)$. Reproduce la tabla de §3.2 y localiza el $M$ que más iteraciones exige en cada caso. Representa, para $e = 0.99$ y $M = 0.01$, las sucesivas $E_k$ de ambas semillas.

import numpy as np
def newton(M, e, E0, tol=1e-14, kmax=100):
    E = E0
    for k in range(1, kmax+1):
        dE = (E - e*np.sin(E) - M)/(1 - e*np.cos(E)); E -= dE
        if abs(dE) < tol: return k
    return kmax
for e in [0.5, 0.9, 0.99]:
    Ms = np.arange(0.01, np.pi, 0.01)
    a = [newton(M, e, M) for M in Ms]; b = [newton(M, e, M + 0.85*e) for M in Ms]
    print(e, max(a), np.mean(a), max(b), np.mean(b))

3. (Demostración) Deriva la ecuación de Kepler hiperbólica $M_h = e\sinh F - F$ a partir de $x = a(e - \cosh F)$, $y = a\sqrt{e^2-1}\sinh F$ (hipérbola con foco en el origen y $a > 0$).

4. (Demostración) Demuestra que la iteración de Newton es monótona decreciente si $E_0 \in [E^*, \pi]$ (para $M\in(0,\pi)$), y que la semilla de Danby cumple $E_0 \ge E^*$.

5. (Ecuación del centro) Deduce los dos primeros términos $\nu - M = 2e\sin M + \tfrac54e^2\sin2M + O(e^3)$. Sugerencia: invierte la ecuación de Kepler por iteración ($E = M + e\sin M + \tfrac12e^2\sin2M + \dots$) y desarrolla $\tan(\nu/2)$.

6. (Tiempos) Calcula la fracción del periodo que transcurre entre el periapsis y $\nu = 90°$ para $e = 0.0933$ (Marte), $0.5$ y $0.9$.

Autoevaluación

Lecturas

• Murray & Dermott (1999), §2.4–2.5: deducción de la ecuación de Kepler y las series en $e$ (ecuación del centro).
• Danby (1988), §6.6–6.7: la semilla de Danby, métodos de orden superior (Halley, Danby–Burkardt) y el caso hiperbólico.
• Curtis (2020), cap. 3: algoritmos 3.1–3.3 y variables universales, con ejemplos numéricos resueltos.
• Vallado (2013), §2.2: algoritmos 2–4 y discusión de la convergencia y de las semillas.
• Colwell (1993), Solving Kepler’s Equation over Three Centuries: historia completa de los métodos, de Kepler a la era de los computadores.
• Markley (1995), DOI 10.1007/BF00691917: solver sin iteración con error de $10^{-15}$.