• Teoría: Mecánica clásica, ecuaciones de Hamilton, campos centrales.
• Ténica numérica: Resolución de ecuaciones diferenciales por el método de Runge-Kutta con paso de integración variable.
• Programación: Fortran, C, introducción a X-windows, gnuplot.

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: M_T, M_L 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 d_TL. La nave parte desde la superficie terrestre con velocidad v₀ y ángulo q₀ 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:

• Obligatorio:  Diseñar y escribir el programa que resuelve las ecuaciones del movimiento de una nave moviéndose bajo la acción del campo gravitatorio terrestre y lunar.

• Voluntario: El problema que se nos plantea es el de diseñar la misión del Apolo XI. Para ello suponemos que la nave tiene unos motores muy potentes que le permiten instantáneamente cambiar el módulo de la velocidad (no asi su dirección). Se pregunta:

  • (1) Dar un plan de vuelo preciso para que una nave despegue de la superficie terrestre, orbite la luna durante dos dias y luego regrese a la Tierra. La velocidad máxima de reentrada (cuando la nave contacta con la Tierra) debe de estar en la ventana de 500-1000 km/h. Esto es, se piden, los valores q₀, v₀ y los instantes en los que se deben de encender los motores. Notar que los motores siempre dan el mismo pulso instantáneo, de igual magnitud, la cual también ha de ser dada como parte del disenño de la misión. También se pide dar la duración de la misión y el gasto energético total en la misma.

  • (2) Comparar la anterior misión con lo que ocurriría si en lugar de una misión directa Tierra-Luna-Tierra hiciésemos: (i) orbitar la Tierra dos veces, (ii) Ir a la Luna, (iii) orbitar la Luna dos dias, (iv) volver a la Tierra (v) orbitar la Tierra dos veces y (vi) caer sobre la Tierra. Diseñar esta misión un comparar las tiempos de duración de la misma y su gasto energético total con los de la anterior misión directa.

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: x_L(t) = d_TLcos(wt) e y_L(t) = d_TLsin(wt). Nos falta encontrar las ecuaciones del movimiento de la nave espacial.

Sean x_o(t) e y_o(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: x_o(t) = r(t)cos(f(t)), y_o(t) = r(t)sin(f(t)), donde r(t) = [x_o(t)²+y_o(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 r_L(t) = [r(t)²+d_TL²-2r(t)cos(f-wt)]^1/2. La energía cinética de la nave es:

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

Asi, el Hamiltoniano queda:

Finalmente, las ecuaciones del movimiento de la nave son:

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:

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

• (0) Dar y_n = y_n⁽⁰⁾ para t = 0
• (1) Evaluar k_n⁽¹⁾ = hf_n(y₁,y₂,...,y_n;t) para n = 1,...,N.
• (2) Evaluar k_n⁽²⁾ = hf_n(y₁+k₁⁽¹⁾/2, y₂+k₂⁽¹⁾/2,...,y_n+k_n⁽¹⁾/2;t+h/2) para n = 1,...,N.
• (3) Evaluar k_n⁽³⁾ = hf_n(y₁+k₁⁽²⁾/2, y₂+k₂⁽²⁾/2,...,y_n+k_n⁽²⁾/2;t+h/2) para n = 1,...,N.
• (4) Evaluar k_n⁽⁴⁾ = hf_n(y₁+k₁⁽³⁾, y₂+k₂⁽³⁾,...,y_n+k_n⁽³⁾;t+h) para n = 1,...,N.
• (5) y_n(t+h) = y_n(t)+[1/6][k_n⁽¹⁾+2k_n⁽²⁾+2k_n⁽³⁾+k_n⁽⁴⁾ ]
• (6) t = t+h. Ir a (1)

El error cometido con este algoritmo es de orden h⁵, 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 10⁴ 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:

• (0) Dar una h inicial, t = 0, el error máximo que toleramos e_max @ h⁵, t₀ y y(t₀).
• (1) Evaluar e = 16|y(t₀+h;h/2)-y(t₀+h;h)|/15 (si y tiene varias componentes, calculamos la r de cada una de ellas y tomamos la de máximo valor).
• (2) Evaluar h_max = h/s, donde s = max((e/e_max)^0.2,10^-8).
• (3) Si s > 2 entonces h = h_max Ir a (1), Si s < 2 entonces t = t+h, y = y(t₀+h;h/2).
• (4) Si h < h_max entonces h = 2h. Ir a (1).

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+h₁;h₂) el resultado iterar varias veces con h = h₂ hasta evolucionar el sistema un tiempo t+h₁. En particular, sabemos que el error del algoritmo de Runge-Kutta es de orden e_h = C h⁵, cuando realizamos una iteración de tamaño h. Si realizamos dos iteraciones de tamaño h/2, el error acumulado será del orden e_h/2 = C (h/2)⁵ 2 = e_h/16. Esto es, e_h-e_h/2 = e_h(1-1/16), asi e_h @ 16 |e_h-e_h/2|/15. Puesto que y(t+h;h) = y_exacta(t+h)+e_h y y(t+h;h/2) = y_exacta(t+h)+e_h/2, tenemos que e_h @ 16|y(t₀+h;h/2)-y(t₀+h;h) |/15. Conociendo e_h y el e_max que toleramos, podemos obtener el h_max correspondiente a nuestro error tolerado. Asi sabiendo que e_h/e_max = h⁵/h_max⁵ obtenemos h_max = h(e_max/e)^1/5. Si h_max < 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 = h_max. Si la h inicial es menor que la h_max 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, p_r y p_f 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/d_TL, f, [p\tilde]_r = p_r/md_TL y [p\tilde]_f = p_f/md_TL². En este caso las ecuaciones del movimiento quedan reducidas a

donde D = GM_T/d_TL³, m = M_L/M_T y [r\tilde]¢ = [1+[r\tilde]²-2[r\tilde]cos(f-wt)]^1/2. Los valores numéricos reales que se deben de utilizar son: G = 6.67×10^-11, M_T = 5.9736 ×10²⁴, M_L = 0.07349 ×10²⁴, d_TL = 3.844 ×10⁸ 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×10⁶ y lunar: 1.737400×10⁶ (todos los datos se dan en unidades del sistema internacional).

Ayudas:

• Comparar trayectorias y tiempos de cálculo obtenidos cuando se utiliza h fija o h-adaptada.
• Demostrar que H¢ = H+wp_f es una constante del movimiento, i.e. dH¢/dt = 0. Utilizar este hecho para conocer el error acumulado por la utilización del algoritmo de Runge-Kutta adaptativo.
• Información y curiosidades sobre las misiones Apolo puede encontrarse en el web: http://www.nasa.gov/kids.html.