Indice | Previo: Problema2Cuerpos.SolucionTiempo.VariablesUniversales | Siguiente: Problema2Cuerpos.Perturbaciones
7.13. Aproximación de dos cuerpos a sistemas jerárquicos¶
En la Sección Motivación habíamos mencionado que muchos sistemas de N cuerpos en el Universo son en realidad sistema jerarquicos de N cuerpos (Ver Definición \ref{box:def:jerarquicos), es decir la dinámica de estos sistemas puede modelarse como la de $N-1$ sistemas de dos cuerpos interactuantes. En las secciones anteriores desarrollamos en detalle la solución analítica para el problema de los dos cuerpos y es tiempo de que pongamos a prueba esa afirmación original.
Para ello consideremos nuevamente el sistema que habíamos estudiado numéricamente al principio del capítulo:
sistema=[
dict(
m=10.0,
r=[1,0,0],
v=[0,1,0.5]),
dict(
m=1.0,
r=[1.5,0,0],
v=[0,-3,1],
),
dict(
m=0.01,
r=[-1,0,0],
v=[0,3,1],
)
]
Como podemos comprobar examinando las trayectorias de las partículas representadas en la Figura (code:ncuerpos_jerarquico1) y en la Figura (code:ncuerpos_jerarquico1_CM), este sistema puede describirse como un sistema anidado de tres cuerpos:
- Las partículas 0 y 1 forman un sistema de dos cuerpos central (lo llamaremos el sistema A).
- Las partícula 2 y el sistema A forma un segundo sistema de dos cuerpos (lo llamaremos el sistema B).
Las propiedades iniciales del sistema A pueden describirse simplemente usando:
from numpy import array
sistemaA=[
dict(
m=10.0,
r=array([1,0,0]),
v=array([0,1,0.5])
),
dict(
m=1.0,
r=array([1.5,0,0]),
v=array([0,-3,1]),
)
]
Para definir las propiedades iniciales del sistema B necesitamos conocer primero la masa total, posición y velocidad del centro de masa del sistema A:
masaA=sistemaA[0]["m"]+sistemaA[1]["m"]
r_CM_A=(sistemaA[0]["m"]*sistemaA[0]["r"]+\
sistemaA[1]["m"]*sistemaA[1]["r"])/masaA
v_CM_A=(sistemaA[0]["m"]*sistemaA[0]["v"]+\
sistemaA[1]["m"]*sistemaA[1]["v"])/masaA
Con estas propiedades calculadas podemos ahora definir el sistema B:
from numpy import array
sistemaB=[
dict(
m=masaA,
r=r_CM_A,
v=v_CM_A,
),
dict(
m=0.01,
r=array([-1,0,0]),
v=array([0,3,1]),
)
]
masaB=sistemaB[0]["m"]+sistemaB[1]["m"]
Podemos ahora utilizar las rutinas desarrolladas en secciones anteriores para predecir la posición de las componentes de cada sistema:
from pymcel.export import propaga_estado
t0=0
t=1.0
#Sistema A
r1,v1,r2,v2,rvecA,vvecA=propaga_estado(sistemaA,t0,t)
#Sistema B
rA,vA,r3,v3,rvecB,vvecB=propaga_estado(sistemaB,t0,t)
Para compararla con las posiciones "reales" resolvamos numéricamente el problema de los N cuerpos:
from pymcel.export import ncuerpos_solucion
ts=[t0,t]
rs,vs,rps,vps,constantes=ncuerpos_solucion(sistema,ts)
Ahora podemos comparar el valor de r1 y v1 con el de rs[0,1,:] y vs[0,1,:]:
La coincidencia no es perfecta pero no hay que minimizar este "logro" increíble. Recordemos que estamos comparando aquí una solución analítica (la de dos problemas de dos cuerpos) con una solución completamente numérica (la del sistema completo.) Ambas soluciones son completamente independientes, y aún así, al menos para este sistema conseguimos calcular posiciones y velocidades con una precisión relativa del orden del 1%.
Una comparación gráfica entre las dos soluciones hace más notable lo conseguido:
from pymcel.export import propaga_estado
from numpy import linspace
#Ventana de integración
t0=0.0
T=5.0
#Solución al problema de los N cuerpos
from pymcel.export import ncuerpos_solucion
ts=linspace(t0,T,200)
rs_num,vs_num,rps_num,vps_num,constantes=ncuerpos_solucion(sistema,ts)
from numpy import zeros_like
rs_aprox=zeros_like(rs_num)
Nt=20
ts=linspace(t0,T,Nt)
for i,t in enumerate(ts):
#Sistema A
ra1,va1,ra2,va2,ravecA,vavecA=propaga_estado(sistemaA,t0,t)
#Sistema B
raA,vaA,ra3,va3,ravecB,vavecB=propaga_estado(sistemaB,t0,t)
rs_aprox[0,i]=ra1
rs_aprox[1,i]=ra2
rs_aprox[2,i]=ra3
#Grafíco
from pymcel.export import plot_ncuerpos_3d
fig=plot_ncuerpos_3d(rs_num,vs_num)
ax=fig.gca()
ax.plot(rs_aprox[0,:Nt,0],rs_aprox[0,:Nt,1],rs_aprox[0,:Nt,2],'b+')
ax.plot(rs_aprox[1,:Nt,0],rs_aprox[1,:Nt,1],rs_aprox[1,:Nt,2],'r+')
ax.plot(rs_aprox[2,:Nt,0],rs_aprox[2,:Nt,1],rs_aprox[2,:Nt,2],'g+')
¿Podemos mejorar de alguna manera la aproximación obtenida aquí? Hay una manera de hacerlo para este sistema, al menos en el caso de las predicciones de los vectores de estado de los cuerpos 1 y 2. Para ello debemos primero reconocer que en el sistema B la posición del centro de masa del sistema A describe una trayectoria elíptica. Sin embargo, en el sistema A (en el que se ha asumido que la tercera partícula no existe), ese mismo centro de masa se mueve siguiendo una trayectoria rectilinea. Es obvio que en la realidad la trayectoria del que llamamos aquí el "centro de masa" de las partículas 1 y 2, no sigue una trayectporia tan simple, debido naturalmente a la presencia perturbadora de la partícula 3.
Una mejor aproximación a la posición de las partículas 1 y 2 se puede obtener entonces si usamos los vectores relativos rvecA y vvecA que obtuvimos con la solución al problema de los dos cuerpos en el sistema A, pero en lugar de asumir que están referidos a un centro de masa que se mueve inercialmente, lo hacemos al centro de masa cuyo movimiento es más complejo y que resulta de la interacción con la partícula 3 (vectores rA y vA):
r1=rA+sistemaA[1]["m"]/masaA*rvecA
r2=rA-sistemaA[0]["m"]/masaA*rvecA
v1=vA+sistemaA[1]["m"]/masaA*vvecA
v2=vA-sistemaA[0]["m"]/masaA*vvecA
Obtenemos inmediatamente un incremento sustancial en la precisión de nuestras predicciones para este sitema: pasamos de predicciones con precisiones cercanas al 1% a unas nuevas que coinciden con la solución numérica en un nivel igual o inferior al 0.01%.
Como vemos este tipo de correcciones de "segundo orden" pueden mejorar sustancialmente las predicciones de los vectores de estado, al menos de una parte de las partículas del sistema (en nuestro ejemplo, la precisión del vector de estado de la partícula 3 no fue mejorada en lo absoluto.) Sin embargo, no existe una regla general que indique qué correcciones exactamente debemos aplicar y cada sistema debe analizarse separadamente.
Un poco de historia: Kepler y la mecánica celeste de precisión. La mecánica celeste de precisión comenzó en el año 1600 cuando Johannes Kepler se unió al equipo de Tycho Brahe con la tarea expresa de usar los datos acumulados por este último (o al menos una parte de ellos) para poner a prueba modelos matemáticos alternativos al modelo Ptolemaico.
A lo largo de cerca de 20 años, Tycho había acumulado observaciones de los planetas, las estrellas e incluso de algunos cometas, usando instrumentos enormes de alta precisión (ver Figura (tycho_instrumentos).) Mientras que en los siglos que precedieron a Tycho, la precisión de las observaciones astronómicas (realizadas normalmente con pequeños cuadrantes, sextantes, astrolabios y otros instrumentos) era en el mejor de los casos de varios minutos de arco, la mayoría de las observaciones de Tycho alcanzaba una precisión cercana al minuto de arco y en algunos casos llegaba a ser de unas decenas de segundos de arco (al menos para las posiciones de las estrellas.)
Ni el modelo Ptolemaico del universo, que había sido usado por más de 1500 años para las predicciones astronómicas, ni el más reciente modelo Copernicano podían predecir con la precisión conseguida por los instrumentos de Tycho, las posiciones observadas de los planetas. Fue precisamente buscando mejorar la precisión de estas predicciones, como Kepler descubrió finalmente sus leyes del movimiento planetario y sembro con ellas la semilla de la mecánica celeste en el fertil suelo intelectual de los 1600.
No debemos sin embargo olvidar que la descripción de las órbitas planetarias como lo hizo Kepler, asumiendo que son elipses con foco en el Sol (o en la Tierra para el caso de la Luna), no es otra cosa que una aproximación de dos cuerpos a un sistema jerarquico, como lo hemos descrito en esta sección (en particular el sistema solar puede catalogarse como un sistema central en las categorías introducidas en la Sección Motivación.)
En esta aproximación Kepler asumió, sin saberlo, que el sistema solar, que en aquel tiempo estaba formado por solo 8 cuerpos, el Sol, Mercurio, Venus, la Tierra, la Luna, Marte, Júpiter y Saturno, era en realidad un sistema jerarquico de 7 pares de cuerpos: Sol-Mercurio, Sol-Venus, Sol-Tierra, Tierra-Luna, Sol-Marte, Sol-Júpiter y Sol-Saturno.
El éxito de la teoría de Kepler al predecir las posiciones planetarias (incluyendo la de la Luna) con una precisión inferior al minuto de arco, fue la primera demostración en la historia del poder que tiene el problema de los dos cuerpos en la descripción de sistemas de N cuerpos.
7.13.1. Predicciones en el Sistema Solar¶
Podemos usar todo lo visto en este capítulo, para poner a prueba la afirmación de que es posible predecir con una precisión aceptable para los estándares astronómicos del pasado (cercana o inferior a 1 minuto de arco), la posición de los planetas en el Sistema Solar. Para ello podemos asumir que nuestro sistema planetario, que en realidad esta formado por decenas de cuerpos de masas no despreciables, puede describirse como un sistema jerárquico de pares de dos cuerpos: Sol-Mercurio, Sol-Venus, Sol-Tierra, Sol-Marte, Sol-Asteroide, etc.
Podemos asumir (como lo hicieron implícitamente Kepler y sus contemporáneos) que la masa del Sol en estos pares de partículas es infinitamente mayor que la de los cuerpos más pequeño, $m_1=m_\odot\approx M$. Con esta aproximación, el Sol se puede asumir fijo en el centro de masa del sistema (ver Ecs. r1_r2_r y v1_v2_v):
$$ \begin{array}{rcl} \vec r_\odot & \approx & 0\\ \vec v_\odot & \approx & 0\\ \end{array} $$mientras que el estado de los otros cuerpos es simplemente igual al vector de estado relativo en el sistema:
$$ \begin{array}{rcl} \vec r_\mathrm{planeta} & \approx & \vec r\\ \vec v_\mathrm{planeta} & \approx & \vec v \end{array} $$Para poner a prueba la aproximación calcularemos a continuación los vectores de estado de la Tierra y Marte y deduciremos a partir de ellos la posición en el cielo del último planeta que después compararemos con su posición real obtenida con predicciones modernas de alta precisión.
Comencemos entonce por determinar el vector de estado de la Tierra y Marte respecto al Sol en una fecha específica. Para ello podemos valernos de la rutina spkezr que introdujimos en el Capítulo El problema de los N cuerpos.
from spiceypy import furnsh
furnsh("pymcel/data/de430.tpc")
furnsh("pymcel/data/de430.bsp")
#Fecha de referencia J2000.0
t0=0
#Parametro gravitacional del sistema
from spiceypy import bodvrd
mu=bodvrd("SUN","GM",1)[1][0]
#Vectores de estado de la Tierra y Marte
from spiceypy import spkezr
tierra,tluz=spkezr("EARTH_BARYCENTER",t0,"ECLIPJ2000","None","SUN")
rtierra0=tierra[:3]
vtierra0=tierra[3:]
marte,tluz=spkezr("MARS_BARYCENTER",t0,"ECLIPJ2000","None","SUN")
rmarte0=marte[:3]
vmarte0=marte[3:]
Ahora usaremos la rutina propaga_f_g que introdujimos en la Sección Las funciones $f$ y $g$ para propagar las posiciones de la Tierra y Marte, 180 días (aproximadamente 6 meses) en el futuro:
#30 días más tarde
t=180*86400
#Propaga la posición de la Tierra
from pymcel.export import propaga_f_g
prediccion_tierra=propaga_f_g(mu,rtierra0,vtierra0,t0,t,verbose=True)
rtierra=prediccion_tierra[5]
vtierra=prediccion_tierra[6]
#Propaga la posición de Marte
prediccion_marte=propaga_f_g(mu,rmarte0,vmarte0,t0,t,verbose=True)
rmarte=prediccion_marte[5]
vmarte=prediccion_marte[6]
Podemos comparar las predicciones realizadas con la teoría de dos cuerpos con aquellas muy precisas calculadas con las integraciones numéricas implícitas en los datos de SPICE:
#Vectores de estado de la Tierra y Marte
tierra,tluz=spkezr("EARTH_BARYCENTER",t,"ECLIPJ2000","None","SUN")
marte,tluz=spkezr("MARS_BARYCENTER",t,"ECLIPJ2000","None","SUN")
Las coincidencias en este caso con los resultados obtenidos con la teoría de los dos cuerpos (Alg. prediccion_ss) son sencillamente notables. Pero vayamos más lejos y hagamos la comparación que nos propusimos hacer desde el principio, es decir, calculemos la diferencia entre la posición de Marte en el cielo predicha con la teoría de dos cuerpos y la real (que asumiremos es la que provee SPICE.)
Para ello nos valdremos de una nueva rutina reclat que permite convertir un vector posición en sus correspondientes coordenadas de latitud y longitud referidas al sistema de referencia usado (en este caso longitud y latitud eclíptica):
from spiceypy import reclat
#Distancia a Marte, longitud y latitud eclíptica del planeta
dist,long,lat=reclat(marte[:3]-tierra[:3])
#Valores usando la teoría aproximada
dist_aprox,long_aprox,lat_aprox=reclat(rmarte-rtierra)
Como lo observo Kepler, el error en la posición en el cielo se mantiene, por lo menos durante 6 meses por debajo de 1 minuto de arco. En la Figura (code:error_posicion_marte)) se muestra el error calculado en la posición en el cielo, vista desde la Tierra, predicha para el planeta rojo con la teoría de dos cuerpos.
7.13.2. Evolución de los elementos orbitales osculatrices¶
Hay un "fenómeno" muy interesante que permite entender en parte porque las predicciones de la teoría de los dos cuerpos usadas aquí para predecir el estado de las partículas en sistemas jerarquicos de N cuerpos, no son más precisas. De acuerdo con lo que vimos en la Sección La órbita osculatriz a cada vector de estado (sea este el vector de estado de una partícula o el del sistema relativo) podemos asociar un conjunto de elementos orbitales osculatrices. ¿Cómo son los elementos osculatrices de las partículas de los sistemas jerarquicos estudiados en esta sección?
Consideremos por ejemplo el sistema ficticio del principio. Obtengamos primero la solución numérica a las ecuaciones de movimiento de lo describen:
from numpy import linspace
ts=linspace(0.0,10.0,200)
from pymcel.export import ncuerpos_solucion
rs,vs,rps,vps,constantes=ncuerpos_solucion(sistema,ts)
Para calcular los elementos osculatrices, obtengamos primero la posición del centro de masa de los sistemas A y B en los que descompusimos el sistema completo:
#Sistema A
r_CM_A=(sistema[0]["m"]*rs[0]+sistema[1]["m"]*rs[1])/masaA
v_CM_A=(sistema[0]["m"]*vs[0]+sistema[1]["m"]*vs[1])/masaA
#Sistema B
r_CM_B=(sistemaB[0]["m"]*r_CM_A+sistemaB[1]["m"]*rs[2])/masaB
v_CM_B=(sistemaB[0]["m"]*v_CM_A+sistemaB[1]["m"]*vs[2])/masaB
Para calcular los elementos osculatrices de la partícula 3, por ejemplo, debemos primero expresar su vector de estado respecto al centro de masa de su respectivo subsistema (sistema B) y de allí usar los procedimientos matemáticos que vimos en la Sección Determinación de la órbita y que se implementan (como vimos en la Sección El problema de los dos cuerpos en SPICE) con la rutina oscelt:
r3_CM=rs[2]-r_CM_B
v3_CM=vs[2]-v_CM_B
Es=[]
for i,t in enumerate(ts):
from spiceypy import oscelt
Es+=[oscelt(list(r3_CM[i])+list(v3_CM[i]),t,masaB)]
from numpy import array
Es=array(Es)
Podemos ahora graficar los elementos orbitales $q$, $e$, $i$, $\Omega$, $\omega$ como función del tiempo:
import matplotlib.pyplot as plt
fig,axs=plt.subplots(5,1,figsize=(5,8),sharex=True)
#Periapsis
axs[0].plot(ts,Es[:,0])
axs[0].set_ylabel("$q$ (u.c.)")
#Excentricidad
axs[1].plot(ts,Es[:,1])
axs[1].set_ylabel("$e$ (u.c.)")
#Inclinación
axs[2].plot(ts,Es[:,2]*180/pi)
axs[2].set_ylabel("$i$ (grados)")
#Longitud del nodo ascendente
axs[3].plot(ts,Es[:,3]*180/pi)
axs[3].set_ylabel("$\Omega$ (grados)")
#Argumento del periapsis
cond=Es[:,4]>pi
Es[cond,4]=Es[cond,4]-2*pi
axs[4].plot(ts,Es[:,4]*180/pi)
axs[4].set_ylabel("$\omega$ (grados)")
#Decoración
axs[-1].set_xlabel("$t$ (u.c.)")
fig.tight_layout()
Como vemos en la Figura (code:variacion_elementos) el valor de los elementos orbitales osculatrices de la partícula 3 no es constante por lo que no es de extrañar que las predicciones que hicimos de la posición de esa partícula con la teoría de los dos cuerpos, y que esencialmente asumen que la órbita tiene elementos constantes e iguales a aquellos en $t=t_0$, no sean muy precisos.
Pero existe un elemento aún más interesante: las variaciones de los elementos no parecen ser "arbitrarias" o "aleatorias. En el caso de $q$, $e$, $i$ y $\omega$ los elementos "oscilan" de formas aparenetemente predecibles alrededor de un valor promedio, con amplitudes muy pequeñas y frecuencias similares. Por su parte, el valor de $\Omega$ parece aumentar monótonamente, con algunos cambios en la tasa de variación que ocurren también a intervalos regulares de tiempo.
7.13.3. Un ejemplo real: los elementos osculatrices de la Luna¶
El curioso comportamiento observado en el sistema ficticio estudiado en esta sección, es también observado en sistemas astronómicos reales. El caso más notable en el Sistema Solar es el de la Luna. En el algoritmo a continuación calculamos los elementos osculatrices de nuestro satélites, relativos al centro de masa del sistema Tierra-Luna, durante un año completo.
#Masa del sistema
from spiceypy import bodvrd
mu=bodvrd("EARTH","GM",1)[1][0]+bodvrd("MOON","GM",1)[1][0]
from numpy import linspace
ts=linspace(0,180*86400,100)
Es=[]
for t in ts:
#Vector de estado de la Luna
from spiceypy import spkezr
luna,tluz=spkezr("MOON",t,"J2000","None","EARTH_BARYCENTER")
#Elementos osculatrices
from spiceypy import oscelt
Es+=[oscelt(luna,t,mu)]
from numpy import array
Es=array(Es)
Una gráfica de los elementos osculatirces calculados aquí como función del tiempo (medido en periodos sinódicos, es decir el tiempo de Luna llena a Luna llena) se muestra en la Figura (code:variacion_elementos). Como vemos allí y como observamos también en el sistema ficticio, existen variaciones periódicas en los elementos orbitales de la Luna. Adicionalmente, el período de variación de la excentricidad y la distancia al periapsis de la Luna es exactamente igual al período sinódico: la excentricidad, por ejemplo siempre es mínima cuando la Luna esta cerca a la luna nueva.
¿De qué dependen estas variaciones en los elementos osculatrices? ¿cómo puede predecirse y de qué manera su predicción puede ayudarnos a describir sistemas jerárquicos de N cuerpos? La solución a este problema obsesionó a Newton hasta el final de su vida y fue necesario una nueva generación de matemáticos en los 1700, entre ellos Joseph Louis Lagrange, para resolver satisfactoriamente el problema. En la siguiente sección abordaremos algunas nociones básicas del problema. Un desarrollo a fondo del mismo es parte de textos más avanzados de Mecánica Celeste.
Indice | Previo: Problema2Cuerpos.SolucionTiempo.VariablesUniversales | Siguiente: Problema2Cuerpos.Perturbaciones