Viaje a la Luna: el problema de los tres cuerpos

1 Introducción

Enviar una nave a la Luna, el cuerpo celeste más próximo a la Tierra, a una distancia media aproximada de 384.400 km, implica considerar el problema de tres cuerpos que se mueven bajo la acción de fuerzas gravitatorias mutuas. El problema de dos cuerpos es resoluble en general -se puede reducir al de un solo cuerpo- pero la adición de un tercero lo hace irresoluble. Así, pues, el problema del movimiento de tres cuerpos que interactúan entre sí a través de fuerzas gravitatorias sigue sin resolverse por métodos analíticos al cabo de más de 200 años de estudio.

Como la masa de la nave es en cualquier caso muy pequeña en comparación con las de la Tierra y la Luna, cabe considerar el caso simplificado, conocido como el problema restringido de los tres cuerpos, en el que la masa de uno de ellos es despreciable frente cualquiera de las otras dos. Bajo estas condiciones la masa pequeña no perturba el movimiento de las mayores, moviéndose éstas sobre órbitas elípticas en torno a su centro de masa común.

No obstante estas simplificaciones, el problema tampoco es resoluble en general y las trayectorias pueden ser "caóticas". En la figura de la izquierda se muestra el posible movimiento de una masa M₃ en el campo campo gravitatorio que crean dos masas grandes e iguales M₁ y M₂. Para simplificar, se ha supuesto que M₁ y M₂ se mueven en órbita circular alrededor de su centro de masas, y se ha escogido un sistema de referencia centrado en dicho círculo y que rota de forma que M₁ y M₂ parecen estar en reposo. En el diagrama a continuación se representan las líneas de fuerza y superficies equipotenciales del campo gravitatorio creado por el sistema Tierra-Luna, suponiendo ambas fijas y despreciando el efecto del Sol. Obsérvese que las líneas de fuerza no son radiales. En el punto A el campo resultante es nulo (Am J. Phys. 33, 712 (1965)).

Aquí vamos a suponer la Tierra inmóvil y que la Luna sigue una órbita circular a su alrededor, como se representa en la figura a la derecha.

La posición de la Luna en el momento del lanzamiento es C₁ y la de llegada C₂. La trayectoria de la nave es por tanto una línea curva cuyo cálculo es, como se ha dicho, muy complicado. Puede ocurrir que la nave, en vez de caer en la Luna, pase cerca de su superficie. En este caso, dependiendo de su velocidad y de la distancia de aproximación máxima, la nave seguirá diferentes trayectorias. Puede quedar retenida en una órbita estable en torno la Luna, pasar cerca de ella y regresar a la Tierra o bien desviarse y perderse en el Sistema Solar. Este fue el caso de la sonda Luna 1, lanzada por la antigua Unión Soviética en enero de 1959. Después de acercarse a 8.000 km, Luna 1 entró en órbita elíptica alrededor del Sol con un periodo de 443 días. La primera sonda en llegar a la Luna fue la Luna 2, lanzada en septiembre de 1959 y cuyo viaje duró 35 horas. La siguiente figura muestra las fases de la misión Apolo que llevó a Armstrong y Aldrin a la superficie de la Luna en julio de 1969. En ella se observa que después del lanzamiento (1) la nave se coloca en una órbita terrestre (2) de la que se separa en el punto (3), para seguir hacia la Luna a lo largo de la trayectoria (4) hasta entrar entrar en órbita lunar (5). Posteriormente, el módulo lunar descendió hasa la superficie de la Luna (6) y, después de completar la misión, regresó a los módulos de comando y de servicio, que habían permanecieron en órbita (7), para inciar el viaje de regreso (8).

2 Modelo

Nuestro sistema se compone de tres cuerpos -Tierra, Luna y nave espacial- que consideramos sin estructura interna y caracterizados de forma única por sus masas, M_T, M_L y m, respectivamente. La única interacción presente es la gravitatoria. Vamos a suponer, por simplicidad, que todos los cuerpos se mueven en el mismo plano y, además, que la Luna gira con velocidad angular constante alrededor de la Tierra en una órbita de radio fijo d_TL. La nave parte de la superficie terrestre con volocidad v₀ y ángulo θ₀ con respecto a unos ejes arbitrarios centrados en la Tierra que suponemos fija. Finalmente, la nave puede gastar una cierta cantidad de energía E en cambiar de velocidad mediante la acción de un motor. Supondremos que este motor funciona como un pulso instantáneo.

Problemas:

Obligatorio: diseñar y escribir el programa que resuelve las ecuaciones del movimiento de una nave moviéndose bajo la acción de los campos gravitatorios terrestre y lunar.
Voluntario: Diseñar la misión del Apolo XI. Para ello se supone que la nave dispone de motores que le permiten cambiar instantáneamente el módulo de la velocidad (no así la dirección). Se pide:
Dar un plan de vuelo preciso para que una nave despegue de la superficie terrestre, orbite la Luna durante dos días y luego regrese a la Tierra. La velocidad máxima de reentrada debe estar en la ventana 500-1000 km/h. Esto es, se piden los valores θ₀ y v₀ y los instantes en los que se deben encender los motores. Notar que éstos siempre dan el mismo pulso instantáneo, cuya magnitud también se ha de proporcionar como parte del diseño de la misión. También se pide dar la duaración total de la misión y el gasto energético de la misma.
Comparar la misión anterior con lo que ocurriría si en lugar de una misión directa Tierra-Luna-Tierra se hiciera: (i) orbitar la Tierra dos veces, (ii) ir a la Luna, (iii) orbitar la Luna dos días, (iv) volver a la Tierra, (v) orbitar la Tierra dos veces y (vi) caer sobre la Tierra.

2.1 Teoría

Tomamos como origen del sistema de referencia a la Tierra. La Luna, puesto que se mueve con velocidad angular ω constante alrededor de la Tierra, tendrá las coordenadas temporales x_L(t) = d_TLcos(ωt) e y_L(t)=d_TLsin(ωt). Falta encontrar las ecuaciones del movimiento de la nave espacial. Sean x(t) e y(t) las coordenadas de la nave en el instante t. Puesto que la interacción tiene simetría polar, conveniene trabajar con coordenadas polares: x(t)=r(t) cos(ϕt), y(t)=r(t) sin(ϕt), donde r(t)=[x(t)²+y(t)²]^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, por aplicación del teorema del coseno, r_L=[r(t)²+d²_TL−2r(t)d_TLcos(ϕ−ωt)]^1/2. La energía cinética de la nave es

⎛
⎝

⋅

⎞
⎠

⎛
⎝

⋅

+ r²

⋅

⎞
⎠

(1)

y la energía potencial

V=−G

mM_T

−G

mM_L

r_L

(2)

A partir de ellas se construye el lagrangiano, L=T−V. Para obtener el hamiltoniano se han de calcular los momentos conjugados:

p_r=

∂L

∂

⋅

, p_ϕ=

∂L

∂

⋅

=mr²

⋅

(3)

Así, el hamiltoniano queda

H=p_r

⋅

+p_ϕ

⋅

−L=

p_r²

p²_ϕ

2mr²

−G

mM_T

−G

mM_L

r_L

(4)

Finalmente, las ecuaciones del movimiento de la nave son

⋅

∂H

∂p_r

p_r

(5)

⋅

∂H

∂p_ϕ

p_ϕ

mr²

(6)

⋅

−

∂H

∂r

p²_ϕ

mr³

−G

mM_T

r²

−G

mM_L

r_L³

[ r−d_TL cos(ϕ−ωt)],

(7)

⋅

−

∂H

∂ϕ

=−G

mM_L

r_L ³

r d_TLsin(ϕ−ωt).

(8)

2.2 Método numérico

Hemos de resolver numéricamente el anterior 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)=f(y(t),t),

(9)

donde y(t) es un vector de N dimensiones. Sea y(t₀)=y₀ dado. Entonces, y(t₀+h) se puede aproximar por

y(t₀+h)=y(t₀)+

[k⁽¹⁾ +2k⁽²⁾ +2k⁽³⁾ +k⁽⁴⁾ ] + O(h⁵),

(10)

donde

k⁽¹⁾

h f(y₀,t₀),

(11)

k⁽²⁾

h f

⎛
⎝

y₀+

k⁽¹⁾

,t₀+

⎞
⎠

(12)

k⁽³⁾

h f

⎛
⎝

y₀+

k⁽²⁾

,t₀+

⎞
⎠

(13)

k⁽⁴⁾

h f(y₀+ k⁽³⁾,t₀+h).

(14)

En general, el algoritmo de iteración para una sistema arbitrario de ecuaciones diferenciales · y_n(t)=f_n(y₁, y₂, …, y_n;t) (n=1,2, …,N) es:

Dar y_n=y_n(0) para t=0 y n=1,2, …,N.
Evaluar k⁽¹⁾_n=hf_n(y₁, y₂, …, y_n;t) para n=1,2, …,N
A Evaluar k⁽²⁾_n = h f_n(y₁+ [(k⁽¹⁾₁)/2],y₂+ [(k⁽¹⁾₂)/2], …, y_n+ [(k⁽¹⁾_n)/2]; t+[h/2]) para n=1,2, …,N.
Evaluar k⁽³⁾_n = h f_n(y₁+ [(k⁽²⁾₁)/2],y₂+ [(k⁽²⁾₂)/2], …, y_n+ [(k⁽²⁾_n)/2]; t+[h/2]) para n=1,2, …,N.
Evaluar k⁽⁴⁾_n = h f_n(y₁+ k⁽³⁾₁,y₂+ k⁽³⁾₂, …, y_n+ k⁽³⁾_n; t+h) para n=1,2, …,N.
y_n(t+h)=y_n(t)+ [1/6] [k⁽¹⁾ +2k⁽²⁾ +2k⁽³⁾ +k⁽⁴⁾ ].
t=t+h. Ir a (2).

El porqué de esta estructura recursiva se explica en el apéndice A, donde hace la derivación detallada del algoritmo de orden 2, se esboza la derivación del de orden 3 y se justifica el uso generalizado del de orden 4.

El error cometido con este algoritmo es de orden h⁵, luego cuanto más pequeño sea h tanto más preciso será el cálculo, aunque más se tardará en alcanzar tiempo largos. Antes de elegir un valor de h comentamos las escalas típicas de tiempo del problema. Una nave espacial en órbita terrestre tiene un periodo de unos 90 minutos. Por tanto, para describir con precisión una de tales órbitas se debería escoger h < 1 minuto. Por otra parte, si la nave se encuentra a mitad de camino entre la Tierra y la Luna, los tiempos característicos de los cambios de velocidad son del orden de las horas. De esta forma, se ve que si se quiere describir conjuntamente las órbitas y los desplazamientos entre la Tierra y la Luna se debe usar el paso más pequeño, esto es, h ≅ 1 minuto. Esta elección supone realizar unas 10⁴ del anterior algoritmo en una misión típica de una semana, lo que hace que el cálculo se lento. Se puede evitar esto adaptando el valor de h a cada momento de la evolución. El siguiente algoritmo realiza esta función:

Dar una h inicial, t=0, y el error máximo que se tolera ϵ ≅ h⁵, t₀ e y₀(t₀).
Evaluar ϵ = 16|y(t₀+h;h/2)−y(t₀+h;h)|/15. En caso de que y tenga varias componentes se ha de calcular el ϵ de cada una ellas y tomar el máximo.
Evaluar h_max=h/s, donde s=max{(ϵ/ϵ_max)^0.2,10⁻⁸}.
Si s > 2, entonces h=h_max. Ir a (2). Si s < 2, entoces t=t+h e y=y(t₀+h;h/2)
Si h < h_max, entonces h=2h. Ir a (2).

La idea principal del algoritmo anterior es estudiar la variación relativa del error que se obtiene cuando se realiza un iteración de paso h en comparación con dos iteraciones de paso h/2. Sea y(t+h₁;h₂) el resultado de iterar varias veces con h=h₂ hasta hacer evolucionar el sistema un tiempo t+h₁. En particular, se sabe que el error del algoritmo de Runge-Kutta es del orden ϵ_h=C h⁵ cuando la iteración tiene un paso h. El error acumulado tras dos iteraciones de paso h/2 será del orden de ϵ_h/2=2C (h/2)⁵=ϵ/16. Esto es, ϵ_h−ϵ_h/2=15 ϵ/16, luego ϵ_h = 16(ϵ_h−ϵ_h/2)/15. Puesto que y(t+h;h)=y_exacta(t+h)+ϵ_h e y(t+h;h/2)=y_exacta(t+h)+ϵ_h/2, se tiene que ϵ_h ≅ 16 [y(t₀+h;h/2)−y(t₀+h;h)]/15. Conociendo ϵ_h y la tolerancia ϵ_max se obtiene el paso h_max correspondiente al máximo error tolerado. Así, sabiendo que ϵ_h/ϵ_max=h⁵/h⁵_max se obtiene h_max = h (ϵ_max/ϵ_h)^1/5. Si h_max < h/2 quiere decir que se ha evaluado las trayectorias utilizando h/2, que da errores mayores que la tolerancia y, por tanto, hemos de recalcular todo con un h menor, en particular h=h_max. Si la h inicial es menor que h_max, quiere decir que estamos realizando el cálculo con excesiva precisión y, por tanto, se puede incrementar h, en particular h=2h.

comete ϵ_h−ϵ_h/2=y(t;h)−y(t;h/2). correspondiente

Por último, comentamos la influencia de los valores de las variables en el cálculo. Como es bien sabido, el ordenador introduce una fuente de error intrínseco dedido al redondeo que efectúa en cada operación aritmética, error que se amplifica cuando se trabaja con variables con valores muy dispares entre sí. Para minimizar este efecto se debe intentar que las variables r, ϕ, p_r y p_ϕ sean del mismo orden de magnitud. Una forma de conseguirlo es reescalarlas, lo cual se puede hacer de varias formas. En particular, se puede usar: ~r=r/d_TL,ϕ, ~p_r=p_r/md_TL y ~p_ϕ=p_ϕ/md²_TL. En este caso, las ecuaciones de movimiento son

⋅

(15)

⋅

(16)

⋅

2
ϕ

−∆

⎧
⎨
⎩

′³

⎡
⎣

−cos(ϕ−ωt)

⎤
⎦

⎫
⎬
⎭

(17)

⋅

−

∆μ

′³

sin(ϕ−ωt),

(18)

donde ∆ ≡ GM_T / d_TL³, μ ≡ M_L/M_T y ~r′ ≡ [1+~r²−2~r cos(ϕ−ωt) ]^1/2. Los valores numéricos que deben usarse son: G=6.67 ×10⁻¹¹, M_T=6.9736×10²⁴, M_L=0.07349 ×10²⁴, d_TL=3.844×10⁸ y ω = 2.6617×10⁻⁶. Además, como la nave despega desde la superficie terrestre y llegar a la lunar, son necesarios los radios de la Tierra, R_T=6.378160 ×10⁶, y de la Luna, R_L=1.7374 ×10⁶ (todos los datos se dan en unidades del sistema internacional).

2.3 Ayudas

Comparar las trayectorias y tiempos de cálculo obtenidos cuando se utiliza h fija o la h adaptada.
Demostrar que H′=H+ωp_ϕ es una constante del movimiento, esto es, · H′=0. Utilizar este hecho para conocer el error acumulado por la utilización del algoritmo de Runge-Kutta adaptativo.

A Sobre las condicones iniciales

Resolver las 4 ecuaciones de movimiento requiere proporcionar los valores iniciales de r, ϕ, p_r y p_ϕ. Los dos primeros determinan el punto de lanzamiento del cohete. En cuanto a los valores iniciales de los momentos p_r y p_ϕ, magnitudes de las que por lo general no se tiene intuición, pueden obtenerse a partir de la velocidad de lanzamiento v=(· x, · y)=(v cos(θ),v sin(θ)) sin más que recurrir a la definiciones. Así, se tienen las relaciones

p_r

m d_TL

d_TL

√

x²+y²

⋅

r d_TL

xvcos(θ)+yvsin(θ)

rd_TL

rvcos(θ)cos(ϕ)+rvsinθsin(ϕ)

rd_TL

cos(θ−ϕ),

p_ϕ

m d_TL²

dϕ

arctan

⎛
⎝

⎞
⎠

1+y²/x²

⎛
⎝

⎞
⎠

r²

(

⋅

x −y

⋅

) =

sin(θ− ϕ).

(19)

Salvo que el punto de lanzamiento está sobre la superficie de la Tierra, r=R_T, no existe un método general para asignar valores iniciales al resto de variables, ϕ, θ y v. Es aconsejable, no obstante, que v sea próxima a la velocidad de escape, puesto que éste es el mínimo valor que garantiza que la nave tiene la energía suficiente para escapar del campo gravitatorio de la Tierra. Finalmente, hay mucha libertad para escoger ϕ y θ, que dependen de dónde esté situada la Luna en el momento del lanzamiento.

B El método de Runge-Kutta

Sea la ecuación diferencial de primer orden

=f(x,y),

(20)

con la condición inicial y(x₀)=y₀, que pretendemos resolver en el intervalo [a,b]. Para ello, consideramos un conjunto de puntos de [a,b], equiespaciados, entre los cuales está el punto x₀¹. Sea y(x) una solución exacta de la ecución e y_j ≈ y(x_j) los valores aproximados resultado del cálculo numérico. En los métodos monopaso, el valor aproximado de y(x_j+1) viene dado mediante

y_j+1 = y_j+hg(x_j,y_j,f,h).

(21)

Por ejemplo, desarrollando en serie de Taylor,

y_j+1=y_j+h y′_j+

h²

y"_j+

h³

y"′_j + O(h⁴).

(22)

Usando la aproximación más simple se obtiene el método de Euler, y_j+1=y_j+h f(x_j,y_j), en el que se observa que g es la propia función f de la ecuación diferencial, es decir, la pendiente en (x_j,y_j). En este caso el error es del orden de h². Veamos cómo, con muy poco esfuerzo, se pueden obtener resultados mejores.

En el método de Runge-Kutta (RK), g es una media ponderada de valores de f(x,y) en el intervalo [x_j,x_j+1], y se dice que es de orden m si alcanza una aproximación equiparable a la del desarrollo de Taylor de ese orden. Una característica del método RK es que hace uso de la función f pero no de sus derivadas. Exponemos a continuación el método de segundo orden:

y_j+1=y_j+h[ak₁+bk₂].

(23)

El primer término de la media ponderada es siempre la pendiente al principio del intervalo, k₁=f(x_j,y_j)=y′(x_j). En cuanto al segundo,

k₂=f(x_j+λh,y_j+μh k₁),

(24)

donde 0 < λ ≤ 1. Los pasos a,b y los números λ,μ se fijan imponiendo que el algoritmo sea compatible con un desarrollo de Taylor de orden 2, y se ha impuesto la forma μh k₁ del segundo incremento para facilitar dicha comparación. Desarrollando k₂,

f(x_j +λh,y_j+μh f_j) = k₁+∂_x f(x_j,y_j)λh+∂_y f(x_j,y_j) μh k₁+ O(h²),

(25)

resulta

y_j+1=y_j +h (a+b)k₁ +h² b [ ∂_x f(x_j,y_j) λ+∂_y f (x_j,y_j) μk₁ ] + O(h³).

(26)

Comparamos ahora esta última expresión con el desarrollo de Taylor de y(x_j+1)=y(x_j+h),

y(x_j+1)=y(x_j)+hy′(x_j) +

h²

y"(x_j) + ….

(27)

Sustituyendo y(x_j) por su valor aproximado y_j y notando que

y"(x_j)=∂_x f(x_j,y_j)+ ∂_y f (x_j,y_j) f(x_j,y_j)

(28)

resulta

y_j+1 ≈ y_j+h f(x_j,y_j)+

h²

[ ∂_x f(x_j,y_j)λ+∂_y f(x_j,y_j) μf(x_j,y_j) ].

(29)

Comparando, se llega al sistema a+b=1, bλ = 1/2, bμ = 1/2, que es indeterminado. Dejando libre b, resulta a=1−b, λ = μ = 1/2b. Si b=1, lo que proporciona un método particular RK, conduce a a=0 y λ = μ = 1/2. En definitiva,

y_j+1=y_j +

⎛
⎝

x_j+

,y_j+

k₁

⎞
⎠

(30)

Un RK de tercer orden viene dado por

y_j+1=y_j+h[ak₁+bk₂+ck₃],

(31)

con

k₁

f(x_j,y_j),

k₂

f(x_j+λh,y_j+μh k1),

k₃

f(x_j+λ₂ h, y_j+μ₂ k₂+(λ₂−μ₂)hk₁ ).

(32)

Los tres pasos, a,b,c y los cuatro coeficientes incrementales λ,μ,λ₂ y μ₂ se calculan desarrollando k₂ y k₃ en serie de Taylor de dos variables hasta orden h², e identificando los factores que multiplican a h y h² con los correspondientes del desarrollo de Taylor de una variable de y(x_j+h). El sistema de ecuaciones así obtenido es indeterminado, y sus distintas soluciones corresponden a diferentes esquemas RK.

El RK de cuarto orden conjuga bien la precisión con el esfuerzo de computación. Uno de ellos es

y_j+1=y_j+h[

k₁+

k₂+

k₃+

k₄],

k₁=f(x_j,y_j), k₂ = f(x_j+

,y_j +

k₁),

k₂=f(x_j+

,y_j+

k₂), k₄=f(x_j+h,y_j + h k₃).

(33)

Footnotes:

¹En general, dado un paso h arbitrario, [x₀−mh,x₀+mh] no coincide con [a,b], es decir, los extremos a y b no coincidirán con puntos de división. Esto carece de importancia y escogeremos el intervalo de máxima longitud [x₀−mh,x₀+mh] contenido en [a,b].

File translated from T_EX by T_TH, version 4.03.
On 3 Feb 2014, 12:58.