El problema de dos cuerpos

1. Motivación

Dos masas puntuales que se atraen con la ley de Newton forman el único problema gravitatorio que se resuelve exactamente. Toda la mecánica celeste posterior (perturbaciones, tres cuerpos, maniobras espaciales) se construye sobre él: cada planeta, cometa o sonda se describe, instante a instante, como si siguiera la órbita de dos cuerpos que mejor ajusta su estado. En esta lección demostramos que esa órbita es una cónica, identificamos las cantidades que se conservan y comprobamos todo con el vector de estado real de Marte que publica el JPL.

2. Reducción a un cuerpo

Sean dos masas $m_1$, $m_2$ con posiciones inerciales $\mathbf R_1$, $\mathbf R_2$ y $G$ la constante de gravitación. Las ecuaciones de Newton son

Sumándolas, el centro de masas $\mathbf R_{cm} = (m_1\mathbf R_1 + m_2\mathbf R_2)/(m_1+m_2)$ cumple $\ddot{\mathbf R}_{cm} = 0$: se mueve en línea recta con velocidad constante (seis constantes de integración). Restando las ecuaciones divididas por su masa, el vector relativo $\mathbf r = \mathbf R_2 - \mathbf R_1$ obedece

que es el movimiento de una partícula en un campo central fijo de parámetro gravitatorio $\mu$ (m³/s²). La energía y el momento angular del sistema en el centro de masas se escriben con la masa reducida $m^* = m_1 m_2/(m_1+m_2)$: $E = \tfrac12 m^* v^2 - G m_1 m_2/r$ y $\mathbf L = m^*\,\mathbf r\times\mathbf v$, con $\mathbf v = \dot{\mathbf r}$. Dividiendo por $m^*$ obtenemos las cantidades específicas (por unidad de masa reducida) que usaremos en adelante:

Para un planeta alrededor del Sol, $\mu = GM_\odot(1 + m_p/M_\odot)$. El JPL lo refleja en su «Keplerian GM»: para Marte, $\mu = 1.32712482870\times10^{11}$ km³/s², que es $GM_\odot = 1.32712440041\times10^{11}$ más $GM_{\text{Marte}} = 4.2828\times10^{4}$ km³/s² (DE441).

3. Integrales primeras

3.1 Momento angular: el movimiento es plano

Como $\mathbf r\perp\mathbf h$ en todo instante, la órbita está contenida en el plano perpendicular a $\mathbf h$. En coordenadas polares $(r,\nu)$ de ese plano, $\mathbf v = \dot r\,\hat{\mathbf r} + r\dot\nu\,\hat{\boldsymbol\nu}$ y

que es la segunda ley de Kepler (ley de las áreas).

3.2 Energía

usando $\mathbf r\cdot\mathbf v = r\dot r$. Luego $\varepsilon = v^2/2 - \mu/r$ es constante.

3.3 El vector de Laplace–Runge–Lenz

Proposición. El vector $\mathbf A = \mathbf v\times\mathbf h - \mu\,\dfrac{\mathbf r}{r}$ se conserva.

Demostración. Como $\mathbf h$ es constante, $\dfrac{d}{dt}(\mathbf v\times\mathbf h) = \ddot{\mathbf r}\times\mathbf h$. Con la identidad $\mathbf a\times(\mathbf b\times\mathbf c) = \mathbf b(\mathbf a\cdot\mathbf c) - \mathbf c(\mathbf a\cdot\mathbf b)$:

Por otro lado, $\dfrac{d}{dt}\dfrac{\mathbf r}{r} = \dfrac{\mathbf v}{r} - \dfrac{\dot r}{r^2}\mathbf r = \dfrac{\mathbf v}{r} - \dfrac{(\mathbf r\cdot\mathbf v)}{r^3}\mathbf r$. Ambas expresiones coinciden, así que $\dot{\mathbf A} = \ddot{\mathbf r}\times\mathbf h - \mu\,\frac{d}{dt}(\mathbf r/r) = \mathbf 0$. $\blacksquare$

Dos propiedades inmediatas:

• $\mathbf A\cdot\mathbf h = 0$ (ambos sumandos son perpendiculares a $\mathbf h$): $\mathbf A$ está en el plano orbital.
• $A^2 = \mu^2 + 2\varepsilon h^2$. En efecto, $\mathbf v\perp\mathbf h$ da $\lVert\mathbf v\times\mathbf h\rVert^2 = v^2h^2$, y $(\mathbf v\times\mathbf h)\cdot\mathbf r = \mathbf h\cdot(\mathbf r\times\mathbf v) = h^2$, luego $A^2 = v^2h^2 - 2\mu h^2/r + \mu^2 = \mu^2 + 2h^2\left(\tfrac{v^2}{2}-\tfrac{\mu}{r}\right)$.

Contemos constantes: $\mathbf h$ (3), $\varepsilon$ (1) y $\mathbf A$ (3) dan 7, pero las dos relaciones anteriores reducen el número de integrales independientes a 5. La sexta constante que fija la solución es el instante de paso por el periapsis $\tau$, que aparecerá en la ecuación de Kepler (lección 1). El vector $\mathbf A$ fue hallado por Laplace en 1799 y redescubierto varias veces (Goldstein 1975); en mecánica cuántica es el responsable de la degeneración «accidental» del átomo de hidrógeno.

4. La ecuación de la órbita

Multiplicamos escalarmente $\mathbf A$ por $\mathbf r$:

Si $\nu$ es el ángulo entre $\mathbf A$ y $\mathbf r$ (la anomalía verdadera), $\mathbf A\cdot\mathbf r = A r\cos\nu$ y, despejando,

Es la ecuación polar de una cónica con un foco en el origen, de semilado recto $p$ y excentricidad $e$; el vector $\mathbf A$ apunta al periapsis ($\nu = 0$, $r$ mínimo). Esta es la primera ley de Kepler, generalizada a todas las cónicas. Dividiendo $A^2 = \mu^2 + 2\varepsilon h^2$ por $\mu^2$:

Para la elipse, $p = a(1-e^2)$ con $a$ el semieje mayor, y de la relación anterior $h^2 = \mu a(1-e^2)$ se obtiene

Los extremos de $r$ son $r_p = a(1-e)$ (periapsis) y $r_a = a(1+e)$ (apoapsis). Finalmente, el área de la elipse, $\pi ab$ con $b = a\sqrt{1-e^2}$, se barre a ritmo $h/2$, de modo que el periodo es $T = 2\pi ab/h$ y, sustituyendo $h$:

5. La ecuación de Binet

Hay una ruta alternativa, sin $\mathbf A$, que además sirve para cualquier fuerza central. Con $u = 1/r$ y $h = r^2\dot\nu$:

La componente radial de $\ddot{\mathbf r} = -\mu\mathbf r/r^3$ en polares es $\ddot r - r\dot\nu^2 = -\mu/r^2$, es decir, $-h^2u^2u'' - h^2u^3 = -\mu u^2$. Dividiendo por $-h^2u^2$:

Es la ecuación de un oscilador armónico forzado por una constante; su solución general es $u = \mu/h^2 + C\cos(\nu-\nu_0)$, que con $\nu_0 = 0$ y $C = \mu e/h^2$ reproduce exactamente $r = p/(1+e\cos\nu)$. Para una fuerza central arbitraria por unidad de masa $f(r)$ la ecuación es $u'' + u = -f(1/u)/(h^2u^2)$; el hecho de que solo las fuerzas $\propto r$ y $\propto r^{-2}$ den órbitas cerradas para todas las condiciones iniciales es el teorema de Bertrand (1873), y una perturbación $\propto r^{-3}$ produce una precesión del periapsis (Problema 4), que es el mecanismo detrás del avance relativista del perihelio de Mercurio.

6. Verificación con un cuerpo real: Marte en J2000

El JPL Horizons (efemérides DE441, Park et al. 2021) da para Marte, respecto al centro del Sol, en el plano de la eclíptica J2000 y en la época JD 2451545.0 TDB (2000-ene-01 12:00):

con $\mu = 1.3271248287\times10^{11}$ km³/s². Calculamos (en km, s; los valores están redondeados a 6 cifras):

Horizons publica, para la misma época, $a = 2.279391329739\times10^{8}$ km y $e = 0.09331510146$: nuestro cálculo coincide en todas las cifras mostradas (el test elementos.test.ts lo comprueba a $10^{-9}$ relativo). Dos comprobaciones más: la inclinación del plano orbital respecto a la eclíptica es $\arccos(h_z/h) = 1.84988°$ (Horizons: $1.849876°$), y $p/r = 1.08568 = 1 + e\cos\nu$ con $\nu = 23.333°$, que es exactamente la anomalía verdadera que da Horizons. La ecuación de la órbita funciona con diez cifras sobre un planeta real: el residuo lo explican las perturbaciones de los otros planetas, que cambian lentamente estos elementos osculadores (lección 4).

7. Análisis numérico

Cancelación en $e$. $\mathbf A/\mu$ es la diferencia de dos vectores de módulo próximo a 1 ($\lVert\mathbf v\times\mathbf h\rVert/\mu \approx 1 + e\cos\nu$ y $\lVert\mathbf r/r\rVert = 1$). Si la órbita es casi circular, la resta pierde cifras: el error relativo en $e$ es del orden de $\epsilon_{\text{máq}}/e$, con $\epsilon_{\text{máq}} \approx 2.2\times10^{-16}$ en doble precisión. Para Marte ($e \approx 0.09$) se pierde una cifra; para una órbita con $e \sim 10^{-6}$ se pierden seis. La alternativa $e^2 = 1 + 2\varepsilon h^2/\mu^2$ no ayuda (también es una resta $1 - (1 - e^2)$). Lo robusto es calcular directamente las componentes pequeñas $e\cos\nu = p/r - 1$ y $e\sin\nu = (\mathbf r\cdot\mathbf v)\,h/(\mu r)$, que se obtienen sin restar cantidades grandes casi iguales (Problema 6).

Sensibilidad de $a$ a la velocidad. De $a = \mu r/(2\mu - rv^2)$,

El factor $v^2/|\varepsilon| = 2(2a/r - 1)$ crece sin límite al acercarse a la parábola ($\varepsilon\to0$): la determinación orbital de un cometa casi parabólico a partir de pocas observaciones es intrínsecamente mal condicionada en $a$ y en $T$ ($\delta T/T = \tfrac32\,\delta a/a$), aunque $p$ y $e$ se obtengan bien.

Diagnóstico de integradores. Como $\mathbf h$, $\varepsilon$ y $\mathbf A$ son exactamente constantes, su deriva en una integración numérica mide el error del método. En RK4, la dirección de $\mathbf A$ gira lentamente (precesión numérica del periapsis) y $\varepsilon$ deriva de forma secular $\propto t\,h^4$; un integrador simpléctico (leapfrog) mantiene $\varepsilon$ acotada pero también hace precesar $\mathbf A$ (véase metodos-numericos/01). El simulador muestra la energía en cada instante.

8. Simulador

Preguntas de análisis:

1. Con $e = 0.9$, compara la rapidez en el periapsis y el apoapsis con la ecuación vis-viva: debe cumplirse $v_p/v_a = (1+e)/(1-e) = 19$.
2. Mantén $a$ fijo y aumenta $e$: el periodo no cambia (tercera ley) aunque la trayectoria sea muy distinta.
3. Observa que el vector velocidad nunca apunta hacia el foco salvo en los ápsides, donde es perpendicular a $\mathbf r$.

Problemas

1. (Demostración) Prueba que $\mathbf v\times\mathbf h = \mu\,\big[(1+e\cos\nu)\,\hat{\mathbf r} - e\sin\nu\,\hat{\boldsymbol\nu}\big]$ en la base polar $(\hat{\mathbf r},\hat{\boldsymbol\nu})$ del plano orbital, y deduce de nuevo que $\mathbf A = \mu e\,\hat{\mathbf p}$, con $\hat{\mathbf p}$ la dirección del periapsis.

2. (Demostración) Deduce la ecuación vis-viva $v^2 = \mu(2/r - 1/a)$ para la elipse y escribe su versión hiperbólica. ¿Qué vale $v$ cuando $r\to\infty$ en la hipérbola?

3. (Hodógrafa) Demuestra que el extremo del vector velocidad describe una circunferencia de radio $\mu/h$ centrada en $(0,\ \mu e/h)$ en el sistema perifocal (eje $x$ hacia el periapsis).

4. (Binet y precesión) Para la fuerza por unidad de masa $f(r) = -\mu/r^2 - \delta/r^3$ con $\delta \ll h^2$, resuelve la ecuación de Binet y demuestra que el periapsis avanza $\Delta\varpi \approx \pi\delta/h^2$ por revolución.

5. (Hipérbola) Un cuerpo llega del infinito con exceso $v_\infty$ y pasa a distancia mínima $r_p$ del centro. Demuestra que $e = 1 + r_p v_\infty^2/\mu$ y que las asíntotas forman con el eje un ángulo tal que el ángulo de deflexión total es $\delta = 2\arcsin(1/e)$.

6. (Computacional) Horizons da para la Tierra (centro Sol, eclíptica J2000, JD 2451545.0 TDB): $\mathbf r = (-2.649903368\times10^{7},\ 1.446972968\times10^{8},\ -6.111494260\times10^{2})$ km, $\mathbf v = (-29.79426007,\ -5.469294940,\ 1.817836785\times10^{-4})$ km/s, $\mu = 1.3271283864\times10^{11}$ km³/s². Escribe un programa que calcule $h$, $\varepsilon$, $a$, $e$, la inclinación $i = \arccos(h_z/h)$ y $\nu$ (usa las fórmulas robustas de la sección 7) y compáralo con Horizons: $a = 1.496650146\times10^{8}$ km, $e = 0.01711862907$, $i = 4.181\times10^{-4}$ grados, $\nu = 358.5689°$. ¿Por qué $a$ no es exactamente 1 UA ni $i$ exactamente 0?

import numpy as np
mu = 1.3271283864171489e11
r = np.array([-2.649903367743050e7, 1.446972967925493e8, -6.111494259536266e2])
v = np.array([-29.79426007043741, -5.469294939770602, 1.817836785027449e-4])
rn, vn = np.linalg.norm(r), np.linalg.norm(v)
h = np.cross(r, v); hn = np.linalg.norm(h)
eps = vn**2/2 - mu/rn;  a = -mu/(2*eps);  p = hn**2/mu
ecos = p/rn - 1;  esin = np.dot(r, v)*hn/(mu*rn)       # sin cancelación
e = np.hypot(ecos, esin);  nu = np.degrees(np.arctan2(esin, ecos)) % 360
i = np.degrees(np.arccos(h[2]/hn))
print(a, e, i, nu)   # 1.4966501455e8  0.017118629  4.181e-4  358.5689

7. (Virial) Demuestra que en una órbita elíptica los promedios temporales cumplen $\langle 1/r\rangle_t = 1/a$ y $\langle v^2\rangle_t = \mu/a$, es decir, $2\langle T\rangle = -\langle V\rangle$.

Autoevaluación

Lecturas

• Murray & Dermott (1999), Solar System Dynamics, cap. 2: el problema de dos cuerpos con la notación que usamos aquí (§2.3–2.5).
• Danby (1988), Fundamentals of Celestial Mechanics, cap. 6: tratamiento clásico con las tres anomalías y la hodógrafa.
• Curtis (2020), Orbital Mechanics for Engineering Students, cap. 2: enfoque ingenieril, con los ejemplos numéricos 2.x en km y km/s.
• Vallado (2013), Fundamentals of Astrodynamics and Applications, cap. 1 y §2.2: constantes, integrales y casos degenerados.
• Goldstein, Poole & Safko (2002), Classical Mechanics, §3.7–3.9: teorema de Bertrand y vector de Laplace–Runge–Lenz; Goldstein (1975), DOI 10.1119/1.9745, sobre la historia del vector.
• Park et al. (2021), DE440/DE441, DOI 10.3847/1538-3881/abd414: las efemérides de las que proceden los vectores de Horizons.