Indice | Previo: Problema2Cuerpos.SolucionTiempo.Sintesis | Siguiente: Problema2Cuerpos.AproximacionJerarquico
Como hemos visto, el cálculo de la anomalía verdadera del vector relativo como función del tiempo en el problema de los dos cuerpos, se consigue después de resolver la ecuación de Halley o la ecuación de Kepler. En términos algorítmicos, como hemos visto en esta sección, el procedimiento puede hacerse bastante engorroso en tanto se hace necesario evaluar en distintos pasos sobre que cónica exactamente ocurre el movimiento.
En 1912, Karl Sundman descubrió un cambio de variable muy útil en el problema de los dos cuerpos que permite obtener lo que podríamos llamar una ecuación de kepler unificada o universal que es válida sin importar la cónica sobre la que se mueva el vector relativo. La teoría de Sundman, fue desarrollada con muchas y diversas variantes durante el resto de los 1900, hasta convertirse en la que se conoce hoy en día como la formulación universal del problema de Kepler.
La formulación presentada aquí es original pero no necesariamente es la más general o rigurosa de las que se han concebido. Para algunas versiones alternativas se invita al lector a explorar los ahora textos clásicos de Bate y otros (Bate et al., 1971) y Danby (Danby, 1992) y las referencias incluídas en ellos.
Considere por ejemplo la ecuación de Halley escrita en la forma:
$$ \sqrt{\frac{\mu}{p^3}}(t-t_p)=\frac{1}{2}\tan\frac{f}{2}+\frac{1}{6}\tan^3\frac{f}{2} $$Es fácil ver mostrar que si introducimos la variable auxiliar $x=\sqrt{p}\tan(f/2)$, la ecuación adopta la forma:
\begin{equation} \label{eq:parabola_universal} \sqrt{\mu}(t-t_p)=qx+\frac{x^3}{3!} \end{equation}donde $q=p/2$.
¿Qué puede tener de especial esta forma particular de la ecuación de Halley?
Escribamos ahora la ecuación de Kepler para órbitas elípticas:
$$ \sqrt{\frac{\mu}{a^3}}(t-t_p)=E-e\sin E $$Si usamos la expansión en series de potencias de la función seno:
$$ \sin E=E-\frac{E^3}{3!}+\frac{E^5}{5!}+\ldots $$la ecuación de Kepler queda:
$$ \sqrt{\frac{\mu}{a^3}}(t-t_p)=(1-e)E+\frac{E^3}{3!}-\frac{E^5}{5!}+\ldots $$Ahora bien, reconociendo que $q=a(1-e)$ y haciendo $x\equiv\sqrt{a}E$ y $\alpha\equiv1/a$, la ecuación de Kepler para órbitas elípticas se puede escribir de la forma:
\begin{equation} \label{eq:kepler_universal} \sqrt{\mu}(t-t_p)=qx+ex^3\left(\frac{1}{3!}-\alpha\frac{x^2}{5!}+\alpha^2\frac{x^4}{7!}\ldots\right) \end{equation}El parecido entre esta versión de la ecuación de Kepler y la ecuación de Halley escrita en la forma de la Ec. (parabola_universal) es simplemente asombroso. Basta notar que cuando se hace $e=1$ y $a\rightarrow\infty$ en la ecuación de Kepler escrita en la forma de la Ec. (kepler_universal) se obtiene justamente la ecuación de Halley. Dos ecuaciones que fueron obtenidas con procedimientos geométricos y analíticos completamente diferentes, al expresarse en la forma de series infinitas revelan su parentesco de forma increible.
Pero las sorpresas continúan cuando consideramos la ecuación de Kepler en el caso hiperbólico:
$$ \sqrt{\frac{\mu}{a^3}}(t-t_p)=e\sinh F-F $$Usando, de forma análoga a como lo hicimos en el caso de la elipse, la expansión en series de Taylor de la función seno hiperbólico:
$$ \sinh F=F+\frac{F^3}{3!}+\frac{F^5}{5!}+\ldots $$la ecuación de Kepler en este caso se escribe como:
$$ \sqrt{\frac{\mu}{|a|^3}}(t-t_p)=(e-1)F+\frac{F^3}{3!}+\frac{F^5}{5!}+\ldots $$Si reconocemos que $a<0$, $q=|a|(1-e)$ y hacemos $x=\sqrt{|a|}F$, obtenemos nuevamente la Ec. (kepler_universal):
$$ \sqrt{\mu}(t-t_p)=qx+ex^3\left(\frac{1}{3!}-\alpha\frac{x^2}{5!}+\alpha^2\frac{x^4}{7!}\ldots\right) $$En el paréntesis del lado derecho de la última ecuación se reconocen claramente las series de Stumpff que introdujimos en la Sección Series infinitas. En términos de estas series la ecuación universal de Kepler queda:
\begin{equation} \label{eq:kepler_universal_stumpff} \sqrt{\mu}(t-t_p)=qx+ex^3c_3(\alpha x^2) \end{equation}donde la variable universal $x$ se puede expresar en términos de las anomalías como:
$$ x= \left\{ \begin{array}{cc} \sqrt{p}\tan(f/2) & \mathrm{Si}\; e=1\\ \sqrt{|a|}G & \mathrm{Si}\; e=1\\ \end{array} \right. $$y $G$ es la anomalía excéntrica generalizada ($G=E$ en el caso de la elipse y $G=F$ en el caso de la hipérbola.)
La solución a la ecuación universal de Kepler se puede obtener usando algunos de los métodos que vimos en la Sección Solución numérica a la ecuación de Kepler.
Para ello, sin embargo es necesario primero escribir una rutina que permita calcular la función universal de Kepler:
$$ k_x(x;M,e,q)=qx+ex^3c_3(\alpha x^2)-M_u $$donde $M_u=\sqrt{\mu}(t-t_p)$ es una anomalía media generalizada (valida para cualquier cónica) y sus derivadas, para lo cual podemos usar las relaciones de recurrencia de las derivadas de las series de Stumpff vistas en la Sección Series infinitas (Ecs. stumpff_derivadas_recurrencia):
\begin{eqnarray} \nonumber \frac{\mathrm{d}k_u}{\mathrm{d}x} & = & q+e x^2c_2(\alpha x^2)\\ \frac{\mathrm{d}^2k_u}{\mathrm{d}x^2} & = & e x c_1(\alpha x^2)\\ \end{eqnarray}La rutina que implementa la ecuación universal de Kepler será:
def funcion_universal_kepler(x,M,e,q):
#Parametro alga
alfa=(1-e)/q
#Funcion universal de Kepler
from pymcel.export import serie_stumpff
k=q*x+e*x**3*serie_stumpff(alfa*x**2,3)-M
kp=q+e*x**2*serie_stumpff(alfa*x**2,2)
kpp=q+e*x*serie_stumpff(alfa*x**2,1)
return k,kp,kpp
Usando esta rutina podemos ahora aplicar el método de Newton o el más general método de Laguerre para resolver la ecuación:
#Propiedades del sistema
mu=1.0
p=2.0
tp=0.0
t=10.0
e=1.5
a=p/(1-e**2)
q=p/(1+e)
#Anomalía media generalizada
from numpy import sqrt
M=sqrt(mu)*(t-tp)
#Solución usando variables universales
from pymcel.export import metodo_laguerre
x,errorx,ni=metodo_laguerre(funcion_universal_kepler,
x0=M,args=(M,e,q),delta=1e-14)
Euni=x/sqrt(abs(a))
errorEuni=errorx/sqrt(abs(a))
#Solución usando la ecuación de Kepler
from pymcel.export import kepler_newton
n=sqrt(mu/abs(a)**3)
M=n*(t-tp)
E,errorE,niE=kepler_newton(M,e,delta=1e-14)
Si la aproximación de variables universales que vimos en la sección anterior tuviera el único propósito de unificar la ecuación de Kepler para distintas cónicas, simplificar su solución o escribir algoritmos más efectivos para encontrarla, este formalismo sería poco más que una curiosidad matemática o numérica sin mayor trascendencia. Sin embargo, hay otros aspecto de esta manera de parametrizar el problema que hace de las variables universales un tema de primera importancia en la solución al problema de los dos cuerpos en mecánica celeste.
Considere por ejemplo la siguiente cuestión. En el conjunto pasos descritos en la Sección Síntesis del problema de los dos cuerpos, vimos que para predecir el estado en cualquier instante del tiempo del vector relativo, es necesario primero convertir el vector de estado inicial ${\vec x}(t_0):[{\vec r}(t_0)\;\dot{\vec r}(t_0)]$ en elementos orbitales, de allí, encontrar, usando la ecuación de Kepler o de Halley, la anomalía verdadera en el instante deseado y finalmente convertir los elementos orbitales, incluyendo la nueva anomalía, en el vector de estado buscado ${\vec x}(t):[{\vec r}(t)\;\dot{\vec r}(t)]$. ¿Habrá una manera de encontrar el vector de estado en el tiempo $t$ usando solamente el vector de estado en el instante inicial, sin pasar por la conversión hacia y desde los elementos orbitales?
Si nos restringimos al plano natural de la cónica, sabemos que los vectores de ${\vec r}(t)$, $\dot{\vec r}(t)$ son coplanares con los vectores ${\vec r}(t_0)$, $\dot{\vec r}(t_0)$. Esto implica que podemos expresar, por ejemplo, el vector posición en el tiempo $t$ como una combinación lineal de la posición y la velocidad inicial:
\begin{equation} \label{eq:r_fg} {\vec r}(t) = f {\vec r}(t_0)+g \dot{\vec r}(t_0) \end{equation}donde $f$ y $g$ son cantidades desconocidas que dependeran naturalmente del tiempo.
Dado que ${\vec r}(t_0)$, $\dot{\vec r}(t_0)$ son constantes, la velocidad en el tiempo $t$ será por su parte:
\begin{equation} \label{eq:v_fg} {\vec v}(t) = \dot{f} {\vec r}(t_0)+\dot{g} \dot{\vec r}(t_0) \end{equation}¿Cómo son las funciones $f$ y $g$? ¿se pueden escribir de modo que no dependan de los elementos orbitales y solo lo hagan del tiempo o de otras cantidades estrictamente relacionadas con la posición y velocidad inicial? Aunque no nos detendremos aquí a ahondar en esta importante cuestión, nos bastará, para los propósitos de este libro, con decir que la respuesta a estas preguntas es afirmativa y que las expresiones más generales de las funciones $f$ y $g$ se obtienen precisamente en el formalismo de las variables universales. El lector curioso puede conocer la deducción de las fórmulas provistas a continuación en cualquier texto avanzado de mecánica celeste (Danby, 1992). Algunas deducciones relacionadas se incluyen en los problemas al final de este capítulo.
Las funciones $f$ y $g$ como función del tiempo y de su estado inicial se pueden calcular usando las ecuaciones:
\begin{eqnarray} \label{eq:f_s} f(s) & = & 1-\left(\frac{\mu}{r_0}\right)s^2 c_2(\beta s^2)\\ \label{eq:g_s} g(s) & = & t-t_0-\mu s^3 c_3(\beta s^2)\\ \label{eq:dotf_s} \dot{f}(s) & = & -\left(\frac{\mu}{rr_0}\right)s c_1(\beta s^2)\\ \label{eq:dotg_s} \dot{g}(s) & = & 1-\left(\frac{\mu}{r}\right)s^2 c_2(\beta s^2)\\ \end{eqnarray}donde $\beta=\mu/a$ y la variable universal $s$ (una variable similar a la $x$ usada en la sección anterior) satisface la ecuación:
\begin{equation} \label{eq:kepler_universal_s} k_s(s;r,r_0,\beta,\mu,t-t_0)\equiv r_0 s c_1(\beta s^2)+r_0 \dot{r}_0 s^2 c_2(\beta s^2)+\mu s^3 c_3(\beta s^2)-(t-t_0)=0 \end{equation}que llamaremos la ecuación de Kepler en la variable universal s.
La rutina a continuación calcula la función de Kepler en la variable universal s, $k_s$, y sus derivadas $k'_s$ y $k''_s$, usando para esto último las relaciones de recurrencia de las derivadas de las series de Stumpff vistas en la Sección Series infinitas (Ecs. stumpff_derivadas_recurrencia):
def funcion_universal_kepler_s(s,r0,rdot0,beta,mu,M):
#Variable auxiliar
u=beta*s**2
#Series de Stumpff requeridas
from pymcel.export import serie_stumpff
c0=serie_stumpff(u,0)
s1c1=s*serie_stumpff(u,1)
s2c2=s**2*serie_stumpff(u,2)
s3c3=s**3*serie_stumpff(u,3)
#Ecuación universal de Kepler en s y sus derivadas
k=r0*s1c1+r0*rdot0*s2c2+mu*s3c3-M
kp=r0*c0+r0*rdot0*s1c1+mu*s2c2
kpp=(mu-r0*beta)*s1c1+r0*rdot0*c0
return k,kp,kpp
El poder de las Ecs. (r_fg) y (v_fg) además de la forma explícita de las funciones $f$ y $g$ y sus derivadas, provistas en las Ecs. (f_s)-(dotg_s), no puede menospreciars. Si bien en este caso la variable universal $s$ carece de la interpretación geométrica que encontramos en la sección anterior para la variable $x$, es claro que esta cantidad nos permite obtener, independientemente del tipo de cónica, la posición relativa de los cuerpos en cualquier instante del tiempo, conociendo unicamente el estado relativo inicial y sin requerir el cálculo intermedio de elementos orbitales o anomalías.
La siguiente rutina implemente las Ecs. (f_s)-(kepler_universal_s) y resuelve el problema de Kepler encontrando la posición y velocidad relativa de un sistema de dos cuerpos en un tiempo $t$ cuando se conoce su posición y velocidad relativa en un tiempo de referencia $t_0$.
def propaga_f_g(mu,rvec0,vvec0,t0,t,delta=1e-14,verbose=False):
from numpy.linalg import norm
from numpy import dot,cross
#Calcular r0, rdot0
r0=norm(rvec0)
rdot0=dot(rvec0,vvec0)/r0
#Calcula el valor del parámetro beta
hvec=cross(rvec0,vvec0)
h=norm(hvec)
e=norm(cross(vvec0,hvec)/mu-rvec0/norm(rvec0))
p=h**2/mu
q=p/(1+e)
beta=mu*(1-e)/q
#Equivalente a la anomalía media
M=t-t0
#Resuelve la ecuación universal de Kepler en s
sn=M/r0
from pymcel.export import metodo_laguerre
s,error,ni=metodo_laguerre(funcion_universal_kepler_s,
x0=sn,args=(r0,rdot0,beta,mu,M),delta=1e-15)
#Variable auxiliar
u=beta*s**2
#Series de Stumpff requeridas
from pymcel.export import serie_stumpff
s1c1=s*serie_stumpff(u,1)
s2c2=s**2*serie_stumpff(u,2)
s3c3=s**3*serie_stumpff(u,3)
#Calcula las funciones f,g
f=1-(mu/r0)*s2c2
g=M-mu*s3c3
#Calcula r
rvec=rvec0*f+vvec0*g
r=norm(rvec)
#Calcula las funciones f',g'
dotf=-(mu/(r*r0))*s1c1
dotg=1-(mu/r)*s2c2
#Calcula v
vvec=rvec0*dotf+vvec0*dotg
return s,f,g,dotf,dotg,rvec,vvec
Que podemos probar usando:
#Sistema
from numpy import array
t0=0
sistema=[
dict(m=1.0,
r=array([0.0,0.0,+0.3]),
v=array([+1.0,0.0,0.5])),
dict(m=0.5,
r=array([+1.0,0.0,0.0]),
v=array([0.0,+1.0,0.0])),
]
#Condiciones iniciales
m1=sistema[0]["m"]
r1_0=sistema[0]["r"]
v1_0=sistema[0]["v"]
m2=sistema[1]["m"]
r2_0=sistema[1]["r"]
v2_0=sistema[1]["v"]
#Posición y velocidad relativa inicial
rvec0=r1_0-r2_0
vvec0=v1_0-v2_0
s,f,g,dotf,dotg,rvec,vvec=propaga_f_g(mu,rvec0,vvec0,t0,t,verbose=True)
Que coincide con los resultados que habíamos obtenido con el Alg. (ejemplo_propaga_estado).
Indice | Previo: Problema2Cuerpos.SolucionTiempo.Sintesis | Siguiente: Problema2Cuerpos.AproximacionJerarquico