2.1. El problema de los tres cuerpos. El movimiento de una nave espacial.

Modelo:

Nuestro sistema se compone de tres cuerpos: Tierra, Luna y nave espacial que consideramos que no tienen estructura interna y que están caracterizados de forma única por sus masas: MT, ML y m respectivamente. La única interacción presente es la gravitatoria. Vamos a suponer por simplicidad que todos cuerpos se mueven en el mismo plano y además la Luna gira a velocidad angular constante w alrededor de la Tierra en un orbita circular de radio fijo dTL. La nave parte desde la superficie terrestre con velocidad v0 y ángulo q0 con respecto a unos ejes arbitrarios centrados en la Tierra que podemos suponer como fija. Finalmente la nave puede gastar una cantidad de energía E en cambiar su velocidad mediante la acción de su motor. Supondremos que el motor de la nave funciona como un pulso instantáneo.

Problema:

Teoría:

Tomamos como origen del sistema de referencia a la Tierra. La Luna, puesto que se mueve a velocidad angular w constante alrededor de la Tierra, tendrá como coordenadas en el tiempo t: xL(t) = dTLcos(wt) e yL(t) = dTLsin(wt). Nos falta encontrar las ecuaciones del movimiento de la nave espacial.

Sean xo(t) e yo(t) las coordenadas de la nave en el instante t. Puesto que la interacción tiene simetría polar es conveniente trabajar con coordenadas polares: xo(t) = r(t)cos(f(t)), yo(t) = r(t)sin(f(t)), donde r(t) = [xo(t)2+yo(t)2]1/2 no es más que la distancia de la nave a la Tierra. En estas coordenadas, la distancia de la nave a la Luna es rL(t) = [r(t)2+dTL2-2r(t)cos(f-wt)]1/2. La energía cinética de la nave es:

T = 1
2
m é
ê
ë
æ
ç
è
d xo
dt
ö
÷
ø
2

 
+ ( d yo
dt
)2 ù
ú
û
= 1
2
m[(t r)2+r2(t f)2 ]
(1)
y la energía potencial:
V = -G mMT
r
-G m ML
rL
(2)

De esta forma el Lagrangiano no es más que L = T-V. Para construir el Hamiltoniano hemos de obtener los momentos conjugados:

pr = L
r
= mt r       pf = L
(t f)
= mr2tf
(3)

Asi, el Hamiltoniano queda:

H = pr t r +pft f-L = pr2
2m
+ pf2
2mr2
-G mMT
r
-G m ML
rL
(4)

Finalmente, las ecuaciones del movimiento de la nave son:

t r = H
pr
= pr
m
(5)
t f = H
pf
= pf
mr2
(6)
t pr = - H
r
= pf2
mr3
-G mMT
r2
-G mML
rL3
[r-dTLcos(f-wt) ]
(7)
t pf = - H
f
= -G mML
rL3
rdTLsin(f-wt)
(8)

Método numérico:

Hemos de resolver numéricamente el conjunto de cuatro ecuaciones diferenciales no lineales. Para ello utilizamos el algoritmo de Runge-Kutta de cuarto orden que en resumen funciona como sigue: Sea la ecuación diferencial a resolver:

t y(t) = f(y(t),t)
(9)
donde y(t) puede ser un vector n-dimensional. Sea y(t0) = y0 dado. Entonces y(t0+h) se puede aproximar por
y(t0+h) = y(t0)+ 1
6
[ k(1)+2k(2)+2k(3)+k(4) ]+O(h5)
(10)
donde
k(1) = hf(y0,t0)
(11)
k(2) = hf(y0+ k(1)
2
,t0+ h
2
)
(12)
k(3) = hf(y0+ k(2)
2
,t0+ h
2
)
(13)
k(4) = hf(y0+k(3),t0+h)
(14)

En general, el algoritmo de iteración para un sistema de ecuaciones diferenciales:

t yn(t) = fn(y1,y2,....,yn;t)    n = 1,¼, N
(15)
es

El error cometido con este algoritmo es de orden h5, asi pues, cuanto más pequeño se tome h más preciso será el cálculo aunque más tardaremos en alcanzar tiempos largos. Antes de elegir un valor de h comentemos las escalas de tiempo típicas en nuestro problema. Una nave espacial en una órbita terrestre tiene un periodo de unos 90 minutos, asi para describir con precisión una órbita terrestre tendríamos que elegir h < 1min. Por otra parte, si la nave se encuentra a mitad de camino entre la Tierra y la Luna, los cambios de velocidad aparecen en el rango de las horas. De esta forma vemos que si queremos describir conjuntamente las órbitas y los desplazamientos entre la Luna y la Tierra de la nave debieramos utilizar la h más pequeña que en nuestro caso son h @ 1min. Esta elección supone, en una misión típica de una semana, el realizar unas 104 iteraciones del anterior algoritmo lo que hace que el cálculo sea lento. Una forma de evitar esto es el intentar adaptar el valor de h a cada momento de la evolución. Un algoritmo que realiza esta función es el siguiente:

Comentemos el anterior algoritmo. La idea principal es estudiar la variación relativa del error que obtenemos si relizamos una iteración de tamaño h o si realizamos dos iteraciones de tamaño h/2. Sea y(t+h1;h2) el resultado iterar varias veces con h = h2 hasta evolucionar el sistema un tiempo t+h1. En particular, sabemos que el error del algoritmo de Runge-Kutta es de orden eh = C h5, cuando realizamos una iteración de tamaño h. Si realizamos dos iteraciones de tamaño h/2, el error acumulado será del orden eh/2 = C (h/2)5 2 = eh/16. Esto es, eh-eh/2 = eh(1-1/16), asi eh @ 16 |eh-eh/2|/15. Puesto que y(t+h;h) = yexacta(t+h)+eh y y(t+h;h/2) = yexacta(t+h)+eh/2, tenemos que eh @ 16|y(t0+h;h/2)-y(t0+h;h) |/15. Conociendo eh y el emax que toleramos, podemos obtener el hmax correspondiente a nuestro error tolerado. Asi sabiendo que eh/emax = h5/hmax5 obtenemos hmax = h(emax/e)1/5. Si hmax < h/2 quiere decir que hemos evaluado las trayectorias utilizando h/2 que dan errores más grandes que nuestra tolerancia y por lo tanto hemos de recalcular todo con una h menor, en particular h = hmax. Si la h inicial es menor que la hmax quiere decir que estamos realizando el cálculo con excesiva precisión (respecto a la tolerancia que nos hemos impuesto) y por lo tanto podemos incrementar h, en particular h = 2h.

Un último comentario se refiere a la influencia de los valores de las variables en el cálculo. Como bien sabemos, el ordenador introduce una fuente de error intrinseca debido al redondeo que ejecuta en cada operación aritmética. Estos errores se magnifican cuando trabajamos con variables que tienen valores muy dispares entre si. Para minimizar este efecto debemos de garantizar que nuestras variables r, f, pr y pf son del mismo orden de magnitud. Una forma de garantizar esto es reescalar las variables, lo cual se puede hacer de varias formas. En particular podemos usar: [r\tilde] = r/dTL, f, [p\tilde]r = pr/mdTL y [p\tilde]f = pf/mdTL2. En este caso las ecuaciones del movimiento quedan reducidas a

t ~
r
 
= ~
p
 

r 
(16)
tf =
~
p

f

~
r
2
(17)
t ~
p
 

r 
=
~
p
2
f

~
r
3
-D é
ê
ê
ê
ë
1
~
r
2
+ m
~
r
 
¢3
æ
è
~
r
 
- cos(f-wt) ö
ø
ù
ú
ú
ú
û
(18)
t ~
p
 

f 
= -
mD ~
r
 

~
r
 
¢3
sin(f-wt)
(19)

donde D = GMT/dTL3, m = ML/MT y [r\tilde]¢ = [1+[r\tilde]2-2[r\tilde]cos(f-wt)]1/2. Los valores numéricos reales que se deben de utilizar son: G = 6.67×10-11, MT = 5.9736 ×1024, ML = 0.07349 ×1024, dTL = 3.844 ×108 y w = 2.3606208 ×10-6. Además, como la nave debe de despegar de la superficie terrestre y no tocar la superficie de la Luna, es necesario utilizar los radios terrestre: 6.378160×106 y lunar: 1.737400×106 (todos los datos se dan en unidades del sistema internacional).

Ayudas: