7.11. Una síntesis del problema de los dos cuerpos¶
En las secciones anteriores hemos demostrado que es posible resolver analíticamente el problema de los dos cuerpos, que podemos formular en términos generales como:
Dadas las posiciones y velocidades de dos partículas en un tiempo de referencia dado $t_0$, ${\vec x}_1(t_0):[{\vec r}_1(t_0)\;\dot{\vec r}_1(t_0)]$, ${\vec x}_2(t_0):[{\vec r}_2(t_0)\;\dot{\vec r}_2(t_0)]$, calcular las posiciones y velocidades de las partículas en cualquier instante de tiempo pasado o futuro $t$.
Podríamos resumir la solución encontrada en este capítulo, con el siguiente conjunto de operaciones o pasos, que partiendo de los vectores de estado inicial ${\vec x}_1(t_0)$, ${\vec x}_2(t_0)$, nos permiten obtener los vectores de estado ${\vec x}_1(t)$, ${\vec x}_2(t)$ en cualquier tiempo $t$:
Encontrar la posición y velocidad inicial del centro centro de masa del sistema y expresar el vector de estado en ese sistema de referencia:
$$ {\vec x}_\mathrm{CM}(t_0)=\frac{m_1{\vec x}_1(t_0)+m_2{\vec x}_2(t_0)}{M} $$ donde $M=m_1+m_2$ y
\begin{eqnarray} \nonumber {\vec x}_1'(t_0) & = & {\vec x}_1(t_0)-{\vec x}_\mathrm{CM}(t_0)\\ \nonumber {\vec x}_2'(t_0) & = & {\vec x}_2(t_0)-{\vec x}_\mathrm{CM}(t_0)\\ \end{eqnarray}
- Calcular el vector de estado relativo ${\vec x}:({\vec r}\;\dot{\vec r})$ para el estado inicial $${\vec x}(t_0)={\vec x}_1'(t_0)-{\vec x}_2'(t_0)$$ y el parámetro gravitacional del sistema $\mu=GM$.
Calcular las constantes de movimiento:
\begin{eqnarray} \nonumber \vec h & = & \vec r\times\dot{\vec r} \\ \nonumber \vec{e} & = & \frac{\dot{\vec{r}} \times \vec{h}}{\mu} - \frac{\vec{r}}{r} \end{eqnarray}
Usando las constantes de movimiento, calcular los parámetros de tamaño $p$ (semilatus rectum) y forma $e$ (excentricidad) de la cónica descrita por el vector relativo:
\begin{eqnarray} \nonumber p & = & \frac{h^2}{\mu}\\ \nonumber e & = & |\vec e| \end{eqnarray}
- Usando las constantes de movimiento calcular los elementos orbitales $(p,e,i,\Omega,\omega,f_0)$. Para ello se utilizan las Ecs. (det_p)-(f_cuadrante).
Calcula la anomalía media $M_0$ correspondiente a $f_0$. Si la órbita es una parábola $e=1$, usamos para ello simplemente la ecuación de Halley:
$$ M_0=\frac{1}{2}\left(\tan^3\frac{f_0}{2}+3\tan\frac{f_0}{2}\right) $$ en caso de tratarse de una órbita elíptica o hiperbólica debemos primer obtener la anomalía excéntrica general $G_0$ usando la inversa de la Ec. (\ref{eq}) :
$$ \mathrm{t}\left(\frac{G_0}{2}\right)=\sqrt{\frac{\sigma(1-e)}{1+e}}\tan\frac{f_0}{2} $$
donde:
$$ \sigma,G,\mathrm{s}(G),\mathrm{t}(G)= \left\{ \begin{array}{ll} +1,E,\sin G,\tan G & e<1\\ -1,F,\sinh G,\tanh G & e>1 \end{array} \right. $$ para luego obtener $M_0$ de la ecuación de Kepler:
$$ M_0=\sigma[G_0-e\;\mathrm{s}(G_0)] $$
Calculamos la anomalía en el tiempo $t$:
$$M(t)=M_0+n(t-t_0)$$ donde,
$$ n= \left\{ \begin{array}{cc} 3\sqrt{\mu/p^3} & \mathrm{Si}\;e=1\\ \sqrt{\mu/|a|^3} & \mathrm{Si}\;e\neq 1\\ \end{array} \right. $$
En el caso de elipse o parábola, con la anomalía media calculamos la anomalía excéntrica $G$ para el tiempo en cuestión, resolviendo la ecuación de Kepler:
$$ M(t)=\sigma[G-e\mathrm{s}(G)] $$ usando los métodos vistos en este capítulo. Una vez obtenida $G$ despejamos la anomalía verdadera $f$ usando:
$$ \tan\frac{f}{2}=\sqrt{\frac{\sigma(1-e)}{1+e}}\;\mathrm{t}\left(\frac{G}{2}\right) $$
En el caso de la parábola, la anomalía verdadera se obtiene directamente despejando de la ecuación de Halley:
$$ \tan\frac{f}{2}=y-\frac{1}{y} $$ con $y\equiv\sqrt[3]{M+\sqrt{M^2+1}}$.
- Una vez calculada $f$ el vector de estado relativo se puede obtener aplicando las ecuaciones (elementos_estado_f) y (elementos_dotx)-(elementos_dotz).
Para calcular la posición y velocidad de las partículas en el sistema de referencia inercial original, calculamos primero la posición del centro de masa en $t$:
$$ {\vec r}_\mathrm{CM}(t)={\vec r}_\mathrm{CM}(t_0)+(t-t_0){\vec v}_\mathrm{CM} $$ y obtenemos el vector de estado de las partículas con:
\begin{eqnarray} {\vec x}_1(t) & = & {\vec x}_\mathrm{CM}(t)+\frac{m_2}{M}{\vec x}(t)\\ {\vec x}_2(t) & = & {\vec x}_\mathrm{CM}(t)-\frac{m_1}{M}{\vec x}(t)\\ \end{eqnarray}
Y el problema queda finalmente resuelto.
7.11.1. Un ejemplo numérico¶
Como siempre ha sido la filosofía de este libro, la mejor manera de verificar si lo visto en este capítulo y que hemos resumido en la "receta" general que presentamos en párrafos precedentes, ha sido realmente asimilado, es ponerlo en práctica en una situación concreta.
Para ello supongamos que queremos predecir de forma exacta el movimiento del siguiente sistema de dos cuerpos:
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])),
]
Para describir las propiedades de las partículas y las condiciones iniciales, hemos usado la notación que introdujimos en el Capítulo El problema de n cuerpos. La razón para hacerlo se verá más adelante.
Con esta definición podemos calcular las condiciones iniciales del sistema de dos cuerpos:
#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
Los vectores rvec0 y vvec0 serán usados frecuentemente en esta y las próximas secciones, cuando resolvamos el problema con distintas aproximaciones.
Supongamos que queremos calcular el estado de las partículas en el tiempo $t=10.0$.
En el algoritmo provisto a continuación se detallan los cálculos correspondientes a cada uno de los pasos enumerados anteriormente para realizar esta predicción con la teoría vista en este capítulo. En el algoritmo, hemos especificado usando comentarios y de la manera más clara posible, cada una de las etapas del cálculo. Así mismo y para ahorrar espacio, usamos algunas de las rutinas que escribimos en secciones precedentes para realizar tareas específicas (p.e. resolver la ecuación de Kepler, convertir del vector de estado a los elementos orbitales, etc.)
def propaga_estado(sistema,t0,t,verbose=0):
########################################################
# Preparación del cálculo
########################################################
#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"]
if verbose:
print(f"r1_0 = {r1_0}, v1_0 = {v1_0}")
print(f"r2_0 = {r2_0}, v2_0 = {v2_0}")
Mtot=m1+m2
#En unidades canónicas G=1
mu=Mtot
#Paso 1: estado del centro de masa
r_CM_0=(m1*r1_0+m2*r2_0)/Mtot
v_CM_0=(m1*v1_0+m2*v2_0)/Mtot
if verbose:print(f"r_CM_0 = {r_CM_0}, v_CM_0 = {v_CM_0}")
#Paso 2: Condiciones iniciales relativas
r_0=r1_0-r2_0
v_0=v1_0-v2_0
if verbose:print(f"r_0 = {r_0}, v_0 = {v_0}")
#Paso 3: Constantes de movimiento
from numpy import cross
from numpy.linalg import norm
hvec=cross(r_0,v_0)
evec=cross(v_0,hvec)/mu-r_0/norm(r_0)
if verbose:print(f"hvec = {hvec}, evec = {evec}")
#Paso 4 y 5: Elementos orbitales
from pymcel.export import estado_a_elementos
from numpy import hstack
p,e,i,W,w,f0=estado_a_elementos(mu,hstack((r_0,v_0)))
from numpy import pi
if verbose:
print(f"Elementos: {p}, {e}, {i*180/pi}, {W*180/pi}, {w*180/pi}, {f0*180/pi}")
#Paso 6: Anomalía media inicial
if e==1:
from numpy import tan
tanf02=tan(f0/2)
#Ecuación de Halley
M0=0.5*(tanf02**3+3*tanf02)
else:
from numpy import sin,cos,sinh,cosh,tan,tanh
from numpy import sqrt,arctan,arctanh
sigma=+1 if e<1 else -1
s=sin if e<1 else sinh
c=cos if e<1 else cosh
ta=tan if e<1 else tanh
at=arctan if e<1 else arctanh
#Anomalía excéntrica
G0=2*at(sqrt(sigma*(1-e)/(1+e))*tan(f0/2))
#Ecuación de Kepler
M0=sigma*(G0-e*s(G0))
if verbose:print(f"M0 = {M0*180/pi}")
########################################################
# Aquí viene la predicción
########################################################
#Paso 7: Anomalía media en t
if e==1:
n=3*sqrt(mu/p**3)
else:
a=p/(1-e**2)
n=sqrt(mu/abs(a)**3)
M=M0+n*(t-t0)
if verbose:print(f"n = {n}, M = {M*180/pi}")
#Paso 8: Anomalía verdadera en t:
from numpy import arctan
if e==1:
y=(M+sqrt(M**2+1))**(1./3)
f=2*arctan(y-1/y)
else:
from pymcel.export import kepler_newton
G,error,ni=kepler_newton(M,e,M,1e-14)
f=2*arctan(sqrt((1+e)/(sigma*(1-e)))*ta(G/2))
if verbose:print(f"f = {f*180/pi}")
#Paso 9: de elementos a estado
from pymcel.export import elementos_a_estado
from numpy import array
x=elementos_a_estado(mu,array([p,e,i,W,w,f]))
r=x[:3]
v=x[3:]
if verbose:
print(f"r = {r}, v = {v}")
print(f"h = {cross(r,v)}")
#Paso 10: estado en el sistema de referencia original
v_CM=v_CM_0
r_CM=r_CM_0+v_CM_0*(t-t0)
if verbose:print(f"r_CM = {r_CM}, v_CM = {v_CM}")
r1=r_CM+(m2/Mtot)*r
v1=v_CM+(m2/Mtot)*v
r2=r_CM-(m1/Mtot)*r
v2=v_CM-(m1/Mtot)*v
#Variables requeridas para comparaciones
if verbose:
from numpy import dot
print(f"f0={f0};f={f};r={norm(r)};r0={norm(r_0)};rdot0={dot(r_0,v_0)/norm(r_0)}")
return r1,v1,r2,v2,r,v
Usando la rutina podemos calcular entonces la posición y velocidad de las partículas en el tiempo deseado $t=10$:
t=10.0
r1,v1,r2,v2,rvec,vvec=propaga_estado(sistema,t0,t)
Para comprobar la validez de los resultados teóricos que implementamos en esta rutina, y como hicimos en la Sección La órbita en el espacio, podemos comparar nuestra predicción analítica con aquella obtenida numéricamente:
from numpy import linspace
ts=linspace(0,30.0,20)
from pymcel.export import ncuerpos_solucion
rs,vs,rps,vps,constantes=ncuerpos_solucion(sistema,ts)
Calculemos ahora el vector de estado de cada partícula, para cada uno de los tiempos en los que calculamos con la rutina numérica la evolución del sistema:
from numpy import zeros_like
rs_teo=zeros_like(rs)
vs_teo=zeros_like(vs)
for i,t in enumerate(ts):
r1,v1,r2,v2,rvec,vvec=propaga_estado(sistema,t0,t)
rs_teo[0,i]=r1
rs_teo[1,i]=r2
vs_teo[0,i]=v1
vs_teo[1,i]=v2
Podemos hacer un gráfico de la trayectoria de las partículas en el espacio:
from pymcel.export import plot_ncuerpos_3d
fig=plot_ncuerpos_3d(rs,vs,lw=1,marker='o',mfc='None')
ax=fig.gca()
ax.plot(rs_teo[0,:,0],rs_teo[0,:,1],rs_teo[0,:,2],'k+');
ax.plot(rs_teo[1,:,0],rs_teo[1,:,1],rs_teo[1,:,2],'r+');
Con este gráfico tenemos una primera comprobación de la coincidencia casi perfecta entre la solución del problema de los dos cuerpos desarrollada en este capítulo y lo que esperaríamos en un sistema real.
Una comprobación más rigurosa se realiza comparando una a una todas las componentes de la posición y la velocidad calculadas con nuestra teoría analítica, con aquellas obtenidas por la integración numérica de las ecuaciones de movimiento:
error_rs=abs(rs-rs_teo)
error_vs=abs(vs-vs_teo)
import matplotlib.pyplot as plt
fig=plt.figure()
ax=fig.gca()
for i in range(2):
for j in range(3):
ax.plot(ts[1:],error_rs[i,1:,j],label="")
ax.plot(ts[1:],error_vs[i,1:,j],label="")
#Decoración
ax.set_xlabel("$t$ (u.c.)")
ax.set_ylabel("|Teoría - Método Numérico|")
ax.set_yscale("log")
Como puede apreciarse en la Figura (code:error_teoria), la comparación entre las componentes del vector de estado calculado con la teoría y el obtenido con la solución numérica de la e.d.m., tiene dos características bastante notables:
La diferencia entre ambos es muy pequeño. Este resultado permite concluir que la teoría es capaz de predecir con gran precisión la posición y velocidad de las partículas, que es justamente lo que buscabamos.
La diferencia no es nula y al parecer crece con el tiempo. Este importante hecho, en lugar de revelar una limitación intrínseca de nuestros resultados analíticos, en realidad es prueba de la precisión limitada de los métodos numéricos utilizados para resolver la e.d.m. Adicionalmente, el hecho de que el error de estos métodos tenga una tendencia creciente con el tiempo, es una limitación bien conocida de los métodos numéricos aplicados en la solución de la e.d.m. de n cuerpos en mecánica celeste, un problema sobre el que volveremos en el [Capítulo El formalismo hamiltoniano].