Indice | Previo: Problema2Cuerpos.OrbitaEspacio | Siguiente: Problema2Cuerpos.SolucionTiempo.AproximacionKepler

7.1. El problema de los dos cuerpos en el tiempo

Como vimos en la Sección Constantes de movimiento, ninguna de las constantes de movimiento identificadas en el problema relativo de los dos cuerpos, depende explícitamente del tiempo. Esto implica que con las cuadraturas disponibles lo único que podemos hacer es expresar la posición y velocidad en términos de un parámetro, la anomalía verdadera $f$ (ver Sección La ecuación de la trayectoria) sin tener todavía ninguna clave de cómo este parámetro depende del tiempo. En esta sección nos proponemos justamente eso, encontrar la manera de calcular la anomalía verdadera del vector relativo (o cualquiera de las otras anomalías definidas antes) para cualquier instante del tiempo. Llamamos a este el problema de Kepler y su solución no es otra cosa que el objetivo perseguido desde el principio de este libro.

Para ello debemos primero responder a la pregunta: ¿hay alguna propiedad geométrica en el movimiento del vector relativo sobre la cónica cuyo valor puede predecirse exactamente en el tiempo? Examinando los teoremas del movimiento orbital (ver Sección Teoremas del movimiento orbital) nos damos cuenta que la única cantidad con esta propiedad es el área $A$ barrida por el radio vector (area del sector de cónica) que, deacuerdo con el Teo. (ncuerpos_solucion), cambia uniformemente con el tiempo a razón de (Ec. velocidad_areal):

$$ \frac{\mathrm{d}A}{\mathrm{d}t}=\frac{h}{2} $$

Si integramos esta ecuación el área del sector de cónica barrida por el vector relativo entre el tiempo de paso por el periapsis $t_p$ y un tiempo arbitrario $t$ será:

\begin{equation} \label{eq:area_sector_tiempo} \Delta A = \frac{h}{2}(t-t_p) \end{equation}

Ahora bien, como vimos en la Sección Áreas de las cónicas el área del sector $\Delta A$ depende explícitamente de las anomalías (verdadera o excéntrica), de modo la Ec. (area_sector_tiempo) combinada con la expresión explícita de $\Delta A$ como función de las anomalías proveerá la respuestas final a la pregunta formulada aquí.

7.10.1. La ecuación de Halley

Aunque no es común presentarlo en este orden e históricamente el caso de la parábola fue el último en ser planteado como se muestra aquí, comenzaremos por resolver la pregunta de cómo calcular la anomalía verdadera del vector relativo en el problema de los dos cuerpos que se mueven sobre una parábola, $e=1$ (p.e. un cometa.)

En la Sección Área de las cónicas habíamos deducido geométricamente la relación entre el área de un sector de parábola y la anomalía verdadera correspondiente:

$$ \Delta A_\mathrm{parabola}=\frac{1}{4}p^2\left(\tan\frac{f}{2}+\frac{1}{3}\tan^3\frac{f}{2}\right) $$

Dado que por la Ec. ({eq:area_sector_tiempo) sabemos cuánto vale $\Delta A$ en cualquier tiempo, podemos escribir explícitamente:

$$ \frac{1}{2}h(t-t_p)=\frac{1}{4}p^2\left(\tan\frac{f}{2}+\frac{1}{3}\tan^3\frac{f}{2}\right) $$

Si usamos $h=\sqrt{\mu p}$ obtenemos después de algunas manipulaciones algebraicas elementales:

\begin{equation} \label{eq:ecuacion_barker} 6\sqrt{\frac{\mu}{p^3}}(t-t_p)=\tan^3\frac{f}{2}+3\tan\frac{f}{2} \end{equation}

Llamaremos a esta ecuación, la ecuación de Halley y su solución nos brinda finalmente la respuesta a la pregunta original.

Un poco de historia: ¿Halley o Barker?. Es común presentar la Ec. (ecuacion_halley) en la forma:

$$ 2\sqrt{\frac{\mu}{p^3}}(t-t_p)=\tan\frac{f}{2}+\frac{1}{3}\tan^3\frac{f}{2} $$

y referirse a ella como la Ecuación de Barker. Así lo hacen prácticamente la totalidad de los textos de mecánica celeste escritos en los últimos 250 años.

Una indagación histórica más juiciosa (Colwell, 1993) ha revelado, sin embargo, que el primero en deducir y utilizar esta ecuación fue Edmund Halley ("Edmond Hali".)

Halley (1656-1742, ver Figura (halley)) fue un virtuoso inglés contemporáneo de Newton, cuyos trabajos científicos abarcaron áreas tan diversas como la astronomía (en la que hizo sus principales contribuciones), la física, la meteorología, las matemáticas y la geofísica. Fue Halley en 1684 quién presentó a Newton el problema que conduciría a este a su teoría general de la gravitación: qué tipo de trayectoria seguirá un cuerpo que se mueva bajo la influencia de una fuerza que disminuya con el cuadrado de la distancia. Newton respondió a este problema con la monografía "De Motu corporum in gyrum", que no es otra cosa que la "semilla" de la Mecánica Celeste moderna. En 1687, Halley culminó el "trabajo editorial" de su vida al publicar la primera edición de los Principia de Newton.

En 1705, más de veinte años después de la publicación de los Principia, Halley publicá la monografía "Una sinopsis de la astronomía de los cometas" $^1$ en la que compila información observacional relacionada con los cometas observados en su tiempo y en siglos anteriores y de provee tablas de posiciones pasadas y futuras de esos mismos cometas. Es en este trabajo en el que Halley introduce por primera vez la Ec. (ecuacion_halley), la cual utiliza para el cálculo de sus tablas. Allí, además, Halley realizá el trascendental descubrimiento por el que sería recordado por siempre. Según sus análisis tres cometas que habían sido observado en 1531, 1607 y 1682 (esté último observado y registrado por él mismo) eran en realidad uno solo. Armado con la teoría de Newton Halley predice que el cometa se aproximaría nuevamente al Sol y a la Tierra en el año 1758. Lamentablemente no viviría lo suficiente para hacer realidad su predicción, que se cumplió sin falta. Hoy, el mencionado cometa lleva su nombre, 1P/Halley o "cometa Halley". Sus últimas apariciones (antes de la preparación de este libro) se produjeron en 1910 y 1986 y las próximas serán en 2061 y 2134 (ver Figura (Halley).)

¿Por qué entonces la mayoría de los autores llaman a la Ec. (ecuacion_halley), "ecuación de Barker"?

En 1757, apenas un año antes del regreso del cometa Halley, El meteorologo inglés Thomas Barker escribió una monografía titulada "Un relato sobre los descubrimientos relacionados con los cometas" en la que no solo realizaba predicciones precisas de la posición de los cometas sobre órbitas parabólicas sino también aplicaba los mismos métodos teóricos para estudiar el problema más práctico del movimiento de proyectiles. Barker uso en su trabajo la Ec. (ecuacion_halley) para sus predicciones y no sabemos si la dedujo por su cuenta o la tomó del trabajo de Halley (como seguramente debió ocurrir) sin hacer necesariamente referencia a él (una práctica común en la época, especialmente si se trataba de resultados que no erán muy difíciles de deducir.) En 1793, el astrónomo Henry Englefield publico su influyente libro "Sobre la determinación de las órbitas de los cometas" en el que uso, entre otros, los métodos de Barker. Para dar crédito a su autor, fue Englefield el que llamo por primera vez a la Ec. (ecuacion_halley) "ecuación de Barker" y de allí se nutrieron los autores sucesivos.

Halley.

Figura 7.87. Halley.

¿Tiene la ecuación de Halley una solución analítica? Si usamos la variable auxiliar $z\equiv\tan(f/2)$ la Ec. (ecuacion_barker) se puede escribir como:

$$ z^3+3z-2M_p=0 $$

donde,

\begin{equation} \label{eq:Mp} M_p\equiv3\sqrt{\frac{\mu}{p^3}}(t-t_p) \end{equation}

una cantidad que llamaremos la anomalía media parabólica.

Esta ecuación cúbica tiene solo una raiz real (Meire, 1985):

$$ z=\sqrt[3]{M_p+\sqrt{M_p^2+1}}+\sqrt[3]{M_p-\sqrt{M_p^2+1}} $$

Si definimos la variable auxiliar $y\equiv\sqrt[3]{M_p+\sqrt{M_p^2+1}}$, la solución a la ecuación de Halley se puede escribir finalmente como:

$$ \tan\frac{f}{2}=y-\frac{1}{y} $$

Una aproximación interesante y útil se obtiene cuando consideramos el caso de tiempos muy cercanos al paso por el periapsis $(t-t_p)\rightarrow 0$. En esta situación $M_p\ll 1$ y una aproximación a $y$ se puede obtener con el teorema del binomio:

\begin{eqnarray} \nonumber y & \approx & 1+\frac{1}{3}M_p\\ \frac{1}{y} & \approx & 1-\frac{1}{3}M_p \end{eqnarray}

Por otro lado como $f\ll 1$,

$$ \tan\frac{f}{2}\approx \frac{f}{2} $$

con lo que resulta:

\begin{equation} \label{eq:f_Mp} f\approx \frac{4}{3}M_p\;\mathrm{para}\;M_p\ll 1 \end{equation}

Como acostumbramos aquí, podemos "verificar" este resultado usando el siguiente algoritmo:

In [1]:
#Constantes del sistema
mu=1.0
h=3.0

#Tamaño de la parabola
p=h**2/mu

#Tiempo de paso por el periapsis
tp=0.0

#Tiempo en el que deseamos calcular f
t=1.0

#Anomalía media parabólica
from numpy import sqrt
Mp=3*sqrt(mu/p**3)*(t-tp)

#Variable auxiliar
y=(Mp+sqrt(1+Mp**2))**(1./3)

#Raíz de la ecuación de Halley
z=y-1/y

#Anomalía verdadera
from numpy import arctan
f=2*arctan(z)

#Aproximación de la anomalía media
faprox=(4./3)*Mp

#Polinomio cúbico en z
polinomio=z**3+3*z-2*Mp

Porpiedades del sistema: h=3, mu=1, tp=0
Solución al problema de Kepler en t = 1:
Anomalía parabólica: Mp = 6.3662 grados
Variable auxiliar: y = 1.0376528
Anomalía verdadera: z = tan(f/2) = 0.0739393
Anomalía verdadera: f = 8.4574333 grados
Aproximación: 4/3 Mp = 8.4882636 grados
Polinomio: z^3 + 3z - 2Mp = 7.494005416219807e-16

Nota: La anomalía media parabólica no es un ángulo. Aunque en el resultado del algoritmo anterior, hemos presentado el valor de la anomalía media parabólica $M_p$ en grados, no puede interpretarse esta cantidad como un verdadero ángulo.

7.10.2. La ecuación de Kepler

En la Sección Área de las cónicas habíamos deducido geométricamente el área de un sector de elipse como función de la anomalía excéntrica (Ec. area_sector_elipse). Usando ahora el teorema de áreas en la Ec. (area_sector_tiempo) podemos escribir la relación:

$$ \frac{h}{2}(t-t_p) = \frac{1}{2} ab (E - e \sin E) $$

Si tenemos en cuenta el teorema armónico (Teo. area_sector_tiempo), a saber que

$$a^3 n^2=\mu$$

y las relaciones que $p=h^2/\mu$ y $a=p/(1-e^2)$ obtenemos finalmente la ecuación:

\begin{equation} \label{eq:ecuacion_kepler} M = E - e \sin E \end{equation}

donde $M\equiv n(t-t_p)$ se conoce como la anomalía media elíptica o simplemente anomalía media.

A esta ecuación fundamental se la conoce universalmente como la Ecuación de Kepler.

El mismo procedimiento puede usarse para deducir una ecuación análoga en el caso de la hipérbola.

Usando el área del sector en la Ec. (area_sector_hiperbola) y aplicando el teorema de áreas obtenemos:

$$ \frac{h}{2}(t-t_p) = \frac{1}{2} |a|\beta (e \sinh F - F) $$

donde $\beta=|a|(e^2-1)$ y $F$ es la anomalía excéntrica sobre la hipérbola.

Una manipulación algebraica similar a la usada para deducir la Ecuación de Kepler en el caso elíptico conduce a:

\begin{equation} \label{eq:ecuacion_kepler_hiperbola} M_h = e \sin F - F \end{equation}

donde $M_h\equiv n_h(t-t_p)$ se conoce como la anomalía media hiperbólica y hemos introducido una nueva cantidad $n_h$, que guarda una relación con el valor absoluto del semieje mayor $|a|$, similar a la que guarda la velocidad angular promedio $n$ en el caso elíptico:

$$ |a|^3 n_h^2=\mu $$

Llamamos a la Ec. (ecuacion_kepler_hiperbola), la ecuación de Kepler hiperbólica.

7.10.3. La función generalizada de Kepler

Las ecuaciones de Kepler para elipses e hipérbolas son muy parecidas pero no pueden escribirse como una sola ecuación de términos de las mismas funciones explícitas. Sin embargo es posible usar una "parametrización" o notación que permite unificarlas, al menos para los propósitos de escribir por ejemplo algunos algoritmos de solución.

Para ello definamos las siguientes variables y funciones auxiliares:

$$ \sigma,\mathrm{c}(G),\mathrm{s}(G),\mathrm{t}(G)= \left\{ \begin{array}{ll} +1,\cos G,\sin G,\tan G & e<1\\ -1,\cosh G,\sinh G,\tanh G & e>1 \end{array} \right. $$

donde $G$ es la anomalía excéntrica generalizada ($E$ o $F$ dependiendo de la cónica).

Con esta parametrización podemos definir la función generalizada de Kepler:

\begin{equation} \label{eq:kepler_generalizada} k(G;M,e)=\sigma[G-e\;\mathrm{s}(G)]-M \end{equation}

donde $M=n(t-t_p)$, $n^2|a|^3=\mu$ y $a=p(1-e^2)$.

Las derivadas de esta función con respecto de $G$, que serán utilizadas más adelante, son por su parte:

\begin{eqnarray} \label{eq:kepler_generalizada_derivada1} k'(G;M,e)\equiv\frac{\mathrm{d}k}{\mathrm{d}G} & = & \sigma[1-e\;\mathrm{c}(G)]\\ \label{eq:kepler_generalizada_derivada2} k''(G;M,e)\equiv\frac{\mathrm{d}^2k}{\mathrm{d}G^2} & = & e\;\mathrm{s}(G)\\ \end{eqnarray}

La ecuación de Kepler generalizada será entonces:

\begin{equation} \label{eq:kepler_generalizada} k(G;M,e)=0 \end{equation}

y su solución nos da el valor de la anomalía excéntrica generalizada $G$, una vez provistos los valores de la anomalía media $M$ y la excentricidad $e$ de la cónica.

Con el valor de $G$, la anomalía verdadera $f$ se calcula finalmente como:

\begin{equation} \label{eq:f_G} \tan\frac{f}{2}=\sqrt{\frac{1+e}{\sigma(1-e)}}\;\mathrm{t}\left(\frac{G}{2}\right) \end{equation}

La siguiente rutina permite implementar la función generalizada de Kepler y calcular además de su valor, el de su primera y segunda derivada:

In [3]:
def funcion_kepler(G,M=0,e=0):
    #Parametro sigma
    sigma=+1 if e<1 else -1
    #Funciones cG, sG
    from numpy import cos,cosh,sin,sinh
    cG=cos(G) if e<1 else cosh(G)
    sG=sin(G) if e<1 else sinh(G)
    #Función de Kepler
    k=sigma*(G-e*sG)-M
    #Primera derivada
    kp=sigma*(1-e*cG)
    #Segunda derivada
    kpp=e*sG
    return k,kp,kpp

En la Figura (code:plot_funcion_kepler) se presentan curvas de la función generalizada de Kepler para distintos valores de la excéntricidad, tanto en el caso de elipse como en el de hipérbolas. En cada caso la intersección de la función generalizada de Kepler con el eje $G$ es igual al valor de la anomalía excéntrica correspondiente a la respectiva anomalía media y excentricidad.

Figura 7.88. Gráficos de la función generalizada de Kepler $k(G;M,e)$ para $M=\pi/2$ y para distintos valores de la excentricidad (elipses, líneas continuas e hipérbolas, lineas rayadas.) El intercepto de cada curva con el eje horizontal G (línea continua horizontal) provee el valor de la anomalía excéntrica $G$ correspondiente a los valores de $M$ y de $e$ respectivos.

En las libretas que encontrará con la versión electrónica del libro podrá ver una versión interactiva de esta gráfica.

Figura 7.89.

7.10.4. Interpretación geométrica de la anomalía media

La solución al problema relativo de los dos cuerpos en las secciones precedentes, nos ha revelado que además de las dos anomalías ya conocidas para indicar la posición relativa, la anomalía verdadera $f$ y la anomalía excéntrica $E$ (o $F$ en el caso de la hipérbola) existe una tercera anomalía, la anomalía media, $M$ (o $M_h$ en el caso de la hipérbola) que también es única para cada punto sobre la trayectoria.

La anomalía media es, además, la única cantidad que varia uniformemente con el tiempo y por la misma razón, puede ser predicha trivialmente conociendo la velocidad angular promedio $n$ del vector relativo en el caso de movimiento elíptico o el parámetro $n_h$ en el caso de la hipérbola.

¿Qué interpretación geométrica tiene la anomalía media $M$?

De la misma manera que la anomalía excéntrica $F$ en el caso de una hipérbola no tiene una interpretación geométrica elemental, como la que vimos para el caso de la anomalía excéntrica $E$ en una elipse$^2$, así mismo, la anomalía media hiperbólica $M_h$ carece también de dicha interpretación. Sin embargo, en el movimiento sobre una elipse, que es con mucho el más estudiado desde los tiempos de Kepler, se han encontrado multiples interpretaciones para $M$.

En la Figura (anomalia_media) adapatamos un diagrama del capítulo 60 del libro Astronomia Nova$^3$ de Johannes Kepler y en la que en 1609, el astrónomo Prusiano dedujo por primera vez la Ecuación que lleva su nombre.

Construcción geométrica original de Johannes Kepler para interpretar la anomalía media $M$ o su suplemento $M'

Figura 7.90. Construcción geométrica original de Johannes Kepler para interpretar la anomalía media $M$ o su suplemento $M'

Kepler interpreto geométricamente la Ec. (ecuacion_kepler) en términos de las áreas de sectores sobre la elipse. Para entender su argumento debemos primero definir los suplementos de la anomalía excéntrica $E'\equiv\pi-E$ y de la anomalía media $M'\equiv\pi-M$.

En términos de estos suplementos la ecuación de Kepler (Ec. ecuacion_kepler) se escribe:

$$ M'=E'+e\sin E' $$

A medida que el cuerpo se mueve sobre la elipse, el valor de la anomalía excéntrica $E$ y de su suplemento $E'$ cambian con una velocidad variable (recordemos que por el teorema de áreas el cuerpo se mueve más rápido cerca al periapsis O.)

Cuando el cuerpo se encuentra en el punto P, toda el área sombreada en la Figura (anomalia_media) (que es la que le falta barrer al radio dirigido del foco a la circunferencia circunscrita, antes de llegar al afelio), será igual al área del triángulo RCF más el área del sector de círculo RCA:

$$ A_\mathrm{RFA}=A_\mathrm{RCA}+A_\mathrm{RCF} $$

Pero $A_\mathrm{RCA}=a^2 E'/2$ (área de un sector circular) y $A_\mathrm{RCF}=a^2 e \sin E'/2$ (área de un triángulo), de modo que reemplazando el área total del sector sombreado en la Figura (anomalia_media)(anomalia_media) queda:

\begin{eqnarray} \nonumber A_\mathrm{RFA} & = & \frac{1}{2}a^2(E'+e\sin E')\\ \nonumber & = & \frac{1}{2}a^2 M' \end{eqnarray}

donde en el último paso hemos usado la ecuación de Kepler en términos de los suplementos.

Este último resultado llevo a Kepler a planetar el problema del cálculo de la posición planetaria sobre una elipse como aquel en el que el Astrónomo, que puede predecir con suma facilidad el valor de la anomalía media $M=n(t-t_p)$, debe encontrar el lugar geométrico del punto R sobre la circunferencia circunscrita (y que esta justo arriba del punto sobre la elipse buscado buscado), tal que el área RFA sea igual a $a^2M'/2$.

Dado que el punto F no está en el centro de la circunferencia circunscrita, $M$ no se puede interpretar como ningún ángulo en la construcción de la Figura (anomalia_media). Como consecuencia de este hecho, y como lo reconoció Kepler desde aquel entonces, el problema de encontrar $E'$ dado un valor de $M'$ es altamente no trivial. A la fecha no se conoce ninguna solución geométrica al problema de Kepler, es decir, no podemos encontrar la posición del punto R usando solamente una regla y un compás.

En términos algebraicos la conclusión de Kepler es equivalente a reconocer que en la ecuación:

$$ M=E-e\sin E $$

es imposible despejar algebraicamente $E$. A este tipo de ecuaciones se las conoce en matemáticas como ecuaciones trascendentales.

Esta conclusión es válida también en el caso de la ecuación de Kepler para el movimiento sobre una hipérbola:

$$ M_h=e\sinh F - F $$

Si queremos entonces resolver el problema de los dos cuerpos en el tiempo, al menos cuando $e\neq 1$ tenemos que encontrar una manera aproximada de resolver esta ecuación. Este será justamente el tema de la siguiente sección. No sobra, sin embargo, hacer aquí una última reflexión, con la que podemos cerrar el esfuerzo de este capítulo para encontrar una solución al problema de los dos cuerpos.

A pesar de que logramos encontrar un número suficiente de cuadraturas para el problema e incluso encontramos una cantidad geométrica que depende de la posición y que cambia de forma predecible con el tiempo (el área barrida por el radiovector), la solución completa del problema de los dos cuerpos dependerá en últimas de resolver numéricamente una ecuación trascendental. En conclusión, ni siquiera la versión más simple del problema general de los N cuerpos admite una solución algebraica cerrada. Sin embargo y como veremos en la Sección Solución en series de la ecuación de Kepler es posible encontrar una solución en términos de series de potencias uniformemente convergentes (funciones analíticas), satisfaciendo así las condiciones al problema de los N cuerpos formulado en el Capítulo El problema de los N cuerpos.

NOTAS AL PIE:

  1. Una copia digital del manuscrito esta disponible en línea en https://archive.org/details/synopsisofastron00hall/page/n4

  1. Aunque vale admitir que hay una interpretación geométrica curiosa para $F$ que puede encontrarse en Portilla, 2019

  1. Para una versión digital en latín, disponible en línea vea: https://www.e-rara.ch/zut/content/pageview/162861

Indice | Previo: Problema2Cuerpos.OrbitaEspacio | Siguiente: Problema2Cuerpos.SolucionTiempo.AproximacionKepler