Restringido De Los Tres Cuerpos De Problema

Question

El movimiento de una nave espacial entre la Tierra y la Luna es un ejemplo de la infame Tres Cuerpo Problema. Se dice que en general una solución analítica para el TBP no es conocido debido a la complejidad de resolver el efecto de tres cuerpos que todas tirar el uno al otro mientras se mueve, un total de seis interacciones. El matemático Richard Arenstorf mientras que en la NASA resuelto un caso especial de este problema, mediante la simplificación de las interacciones de los cuatro, ya que, el efecto de la nave de la gravedad en el movimiento de la mucho más masivo de la Tierra y la Luna es prácticamente inexistente. Arenstorf encontrado una órbita estable para una nave espacial en órbita entre la Tierra y la Luna, en forma de '8'

http://en.wikipedia.org/wiki/Richard_Arenstorf

Fue Arenstorf la solución analítica, o hizo uso numérica mecanismos? Es el '8' de la forma de un camino óptimo, es decir, la ruta en la que la nave iba a ampliar la menor cantidad de energía? Si sí, ¿cómo fue este requisito incluido en la derivación en forma matemática?

Si alguien tiene una clara derivación para este problema, que sería genial, o cualquiera de los enlaces a los libros, otros papeles, etc.

Nota: al Parecer hubo una anterior relacionados con la mathoverflow pregunta sobre esta así:

http://mathoverflow.net/questions/52489/on-the-non-rigorous-calculations-of-the-trajectories-in-the-moon-landings

Arenstorf del informe técnico aquí

http://hdl.handle.net/2060/19630005545

La más limpia, la escueta descripción de las ecuaciones que pude encontrar fueron aquí

http://farside.ph.utexas.edu/teaching/celestial/Celestial.pdf

http://www.uibk.ac.at/mathematik/personal/csomos/arenstorf_en.pdf

Restringido 3 Cuerpo de los problemas mencionados en este libro (capítulo 8)

http://ads.harvard.edu/books/1989fcm..book/

Algunos vpython (visual python) código puede simular Apolo 13 viaje,

http://maths.anu.edu.au/comptlsci/Tutorial-Gravity/tutorial_apollo_link2.html

el resultado es el 8 de órbita como se ve aquí

Alguna explicación sobre las matemáticas, y la superposición de gravedades también se menciona en este enlace.

http://maths.anu.edu.au/comptlsci/Tutorial-Gravity/tutorial_apollo.html

No tengo idea de cómo esto se relaciona Arenstorf del método, creo que en este código la Tierra y la Luna no se mueve, pero pensé que sería bueno compartir.

Saludos,

Answer 1

En los comentarios, que ya estaban informados de que los datos numéricos se tenía que utilizar, así que voy a explicar a partir de ese punto.

La figura 8 diseño no es porque es una ruta óptima. Esto ocurre debido a la gravedad de la luna y la Tierra. Cuando la nave entra dentro de la esfera de influencia (SOI) de la luna, la nave espacial se tira hacia él. Si la nave espacial se mueve a la velocidad de escape, la luna va a perturbar el vuelo, pero la nave no una mosca. Con la velocidad actual, la gravedad de la luna es suficiente para causar que un orbital volar. Al salir de la luna, del SOL, la nave es de ser tragados por la Tierra. Dado que la trayectoria de la mosca fue lanzar la nave espacial lejos de la luna, la luna se cruza en su camino original, pero esta es de corta duración ya que la Tierra, a continuación, tira de nuevo. Si la nave hubiera recogido bastante la velocidad de la maniobra orbital a estar en una parabólica o hiperbólica trayectoria, podría haber escapado de la atracción de la Tierra y ha sido enviado al espacio.

Una forma de determinar las velocidades de los test en el diseño de un vuelo para encontrar la constante de Jacobi, $C$. Para un determinado $C$, el cero curvas de velocidad se determinan. Puesto que quería llegar a la luna, $C\geq -1.6649$ que corresponde a una velocidad inicial de, al menos,$10.85762$, pero con una velocidad de $11.01$ es la velocidad de escape de la Tierra de modo que la velocidad inicial tiene que ser menos de $v_{esc}$.

---

Podemos derivar las ecuaciones de movimiento sólo hasta un punto. Estoy saltando las básicas de derivación de los dos cuerpo a cuerpo problema por lo que nos podemos mover en el quid de la cuestión. $$\begin{alignat*}{4} r_{12} & = \sqrt{(x_1 - x_2)^2} & \qquad & x_1 &&{}= \text{Is the } x \text{ location of } m_1\\ x_2 & = x_1 + r_{12} & & &&{}\phantom{=} \text{relative to the center of gravity.}\\ x_1 & = \frac{-m_2}{m_1 + m_2}r_{12} & & \pi_1 &&{}= \frac{m_1}{m_1 + m_2}\\ x_2 & = \frac{m_1}{m_1 + m_2}r_{12} & & \pi_2 &&{}= \frac{m_2}{m_1 + m_2}\\ 0 & = m_1x_1 + m_2x_2 \end{alignat*}$$ Podemos describir la posición de $m$ $\mathbf{r} = x\hat{\mathbf{i}} + y\hat{\mathbf{j}} + z\hat{\mathbf{k}}$ en relación con el centro de gravedad, es decir, el origen. $$\begin{align} \mathbf{r}_1 & = (x - x_1)\hat{\mathbf{i}} + y\hat{\mathbf{j}} + z\hat{\mathbf{k}}\\ & = (x + \pi_2r_{12})\hat{\mathbf{i}} + y\hat{\mathbf{j}} + z\hat{\mathbf{k}}\tag{1}\\ \mathbf{r}_2 & = (x - x_2)\hat{\mathbf{i}} + y\hat{\mathbf{j}} + z\hat{\mathbf{k}}\\ & = (x - \pi_1r_{12})\hat{\mathbf{i}} + y\hat{\mathbf{j}} + z\hat{\mathbf{k}}\tag{2} \end{align}$$ Vamos a definir la aceleración absoluta donde $\omega$ es la inicial de angular velocidad es constante. A continuación,$\omega = \frac{2\pi}{T}$. $$\ddot{\mathbf{r}}_{\text{abs}} = \mathbf{a}_{\text{rel}} + \mathbf{a}_{\text{CG}} + \mathbf{\varOmega}\times(\mathbf{\varOmega}\times\mathbf{r}) + \dot{\mathbf{\varOmega}}\times\mathbf{r} + 2\mathbf{\varOmega}\times\mathbf{v}_{\text{rel}} \etiqueta{3}$$ donde $$\begin{alignat*}{4} \mathbf{a}_{\text{rel}} & = \text{Rectilinear acceleration relative to the frame} & \quad & \mathbf{\varOmega}\times(\mathbf{\varOmega}\times\mathbf{r}) &&{}= \text{Centripetal acceleration}\\ 2\mathbf{\varOmega}\times\mathbf{v}_{\text{rel}} & = \text{Coriolis acceleration} \end{alignat*}$$ Puesto que la velocidad del centro de gravedad es constante, $\mathbf{a}_{\text{CG}} = 0$, e $\dot{\mathbf{\varOmega}} = 0$ desde el la velocidad angular de una órbita circular es constante. Por lo tanto, $(3)$ se convierte en: $$\ddot{\mathbf{r}} = \mathbf{a}_{\text{rel}} + \mathbf{\varOmega}\times(\mathbf{\varOmega}\times\mathbf{r}) + 2\mathbf{\varOmega}\times\mathbf{v}_{\text{rel}} \etiqueta{4}$$ donde $$\begin{align} \mathbf{\varOmega} & = \varOmega\hat{\mathbf{k}}\tag{5}\\ \mathbf{r} & = x\hat{\mathbf{i}} + y\hat{\mathbf{j}} + z\hat{\mathbf{k}}\tag{6}\\ \dot{\mathbf{r}} & = \mathbf{v}_{\text{CG}} + \mathbf{\varOmega}\times\mathbf{r} + \mathbf{v}_{\text{rel}}\tag{7}\\ \mathbf{v}_{\text{rel}} & = \dot{x}\hat{\mathbf{i}} + \dot{y}\hat{\mathbf{j}} + \dot{z}\hat{\mathbf{k}}\tag{8}\\ \mathbf{a}_{\text{rel}} & = \ddot{x}\hat{\mathbf{i}} + \ddot{y}\hat{\mathbf{j}} + \ddot{z}\hat{\mathbf{k}}\tag{9} \end{align}$$ Después de la sustitución de $(5)$, $(6)$, $(8)$, y $(9)$ a $(3)$, obtenemos $$\ddot{\mathbf{r}} = \left(\ddot{x} - 2\varOmega\dot{y} - \varOmega^2x\right)\hat{\mathbf{i}} + \left(\ddot{y} + 2\varOmega\dot{x} - \varOmega^2y\right)\hat{\mathbf{j}} + \ddot{z}\hat{\mathbf{k}}. \etiqueta{10}$$ Newton $2^{\text{nd}}$ Ley del Movimiento es $m\mathbf{a} = \mathbf{F}_1 + \mathbf{F}_2$ donde $\mathbf{F}_1 = -\frac{Gm_1m}{r_1^3}\mathbf{r}_1$ y $\mathbf{F}_2 = -\frac{Gm_2m}{r_2^3}\mathbf{r}_2$. Deje $\mu_1 = Gm_1$$\mu_2 = Gm_2$. $$\begin{align} m\mathbf{a} & = \mathbf{F}_1 + \mathbf{F}_2\\ m\mathbf{a} & = -\frac{m\mu_1}{r_1^3}\mathbf{r}_1 - \frac{m\mu_2}{r_2^3}\mathbf{r}_2\\ \mathbf{a} & = -\frac{\mu_1}{r_1^3}\mathbf{r}_1 - \frac{\mu_2}{r_2^3}\mathbf{r}_2\\ \bigl(\ddot{x} - 2\varOmega\dot{y} - \varOmega^2x\bigr)\hat{\mathbf{i}} + \bigl(\ddot{y} + 2\varOmega\dot{x} - \varOmega^2y\bigr)\hat{\mathbf{j}} + \ddot{z}\hat{\mathbf{k}} & = -\frac{\mu_1}{r_1^3}\mathbf{r}_1 - \frac{\mu_2}{r_2^3}\mathbf{r}_2\\ \bigl(\ddot{x} - 2\varOmega\dot{y} - \varOmega^2x\bigr)\hat{\mathbf{i}} + \bigl(\ddot{y} + 2\varOmega\dot{x} - \varOmega^2y\bigr)\hat{\mathbf{j}} + \ddot{z}\hat{\mathbf{k}} & = \begin{aligned} - & \frac{\mu_1}{r_1^3}\Bigl[(x + \pi_2r_{12})\hat{\mathbf{i}} + \hat{\mathbf{j}} + \hat{\mathbf{k}}\Bigr]\\ - & \frac{\mu_2}{r_2^3}\Bigl[(x - \pi_1r_{12})\hat{\mathbf{i}} + \hat{\mathbf{j}} + \hat{\mathbf{k}}\Bigr] \end{aligned}\etiqueta{12} \end{align}$$ Ahora todo lo que tienes que hacer es igualar los coeficientes. $$\begin{align} \ddot{x} - 2\varOmega\dot{y} - \varOmega^2x & = -\frac{\mu_1}{r_1^3}(x + \pi_2r_{12}) - \frac{\mu_2}{r_2^3}(x - \pi_1r_{12})\\ \ddot{y} + 2\varOmega\dot{x} - \varOmega^2y & = -\frac{\mu_1}{r_1^3}y - \frac{\mu_2}{r_2^3}y\\ \ddot{z} & = -\frac{\mu_1}{r_1^3}z - \frac{\mu_2}{r_2^3}z \end{align}$$ Ahora tenemos el sistema no lineal de ecuaciones diferenciales ordinarias. Podemos suponer que la trayectoria es en avión y lo hacemos dejando $z = 0$, por lo que sólo tenemos dos ecuaciones restantes: $$\begin{align} \ddot{x} - 2\varOmega\dot{y} - \varOmega^2x & = -\frac{\mu_1}{r_1^3}(x + \pi_2r_{12}) - \frac{\mu_2}{r_2^3}(x - \pi_1r_{12})\\ \ddot{y} + 2\varOmega\dot{x} - \varOmega^2y & = -\frac{\mu_1}{r_1^3}y - \frac{\mu_2}{r_2^3}y \end{align}$$ Ahora para lograr la requieren trayectoria podemos utilizar la verdadera anomalía de $-90^{\circ}$, un ángulo de trayectoria de vuelo de $20^{\circ}$, y una velocidad inicial de $v_0 = 10.9148$ km/s. En este punto, no tenemos otra opción sino para integrar numéricamente. Aquí es un código que va a alcanzar el nivel deseado de la trama. Sin embargo, en la realidad apolo vuelo, que tenía a mediados de correcciones de rumbo y llegaron exactamente en el otro lado de la Tierra.

Python

#!/usr/bin/env ipython

import numpy as np
import scipy
from scipy.integrate import odeint
from mpl_toolkits.mplot3d import Axes3D
from scipy.optimize import brentq
from scipy.optimize import fsolve
import pylab

me = 5.974 * 10 ** 24  #  mass of the earth
mm = 7.348 * 10 ** 22  #  mass of the moon
G = 6.67259 * 10 ** -20  #  gravitational parameter
re = 6378.0  #  radius of the earth in km
rm = 1737.0  #  radius of the moon in km
r12 = 384400.0  #  distance between the CoM of the earth and moon
M = me + mm
d = 200  #  distance the spacecraft is above the Earth
pi2 = mm / M
mue = 398600.0  #  gravitational parameter of earth km^3/sec^2
mum = G * mm  #  grav param of the moon
omega = np.sqrt(mu / r12 ** 3)
vbo = 10.9148
nu = -np.pi*0.5
gamma = 20*np.pi/180.0  #  angle in radians of the flight path

vx = vbo * (np.sin(gamma) * np.cos(nu) - np.cos(gamma) * np.sin(nu))
#  velocity of the bo in the x direction
vy = vbo * (np.sin(gamma) * np.sin(nu) + np.cos(gamma) * np.cos(nu))
#  velocity of the bo in the y direction

xrel = (re + d)*np.cos(nu)-pi2*r12
#  spacecraft x location relative to the earth
yrel = (re + 200.0) * np.sin(nu)

u0 = [xrel, yrel, 0, vx, vy, 0]

def deriv(u, dt):
    return [u[3],  #  dotu[0] = u[3]
            u[4],  #  dotu[1] = u[4]
            u[5],  #  dotu[2] = u[5]
            (2 * omega * u[4] + omega ** 2 * u[0] - mue * (u[0] + pi2 * r12) /
             np.sqrt(((u[0] + pi2 * r12) ** 2 + u[1] ** 2) ** 3) - mum *
             (u[0] - pi1 * r12) /
             np.sqrt(((u[0] - pi1 * r12) ** 2 + u[1] ** 2) ** 3)),
            #  dotu[3] = that
            (-2 * omega * u[3] + omega ** 2 * u[1] - mue * u[1] /
             np.sqrt(((u[0] + pi2 * r12) ** 2 + u[1] ** 2) ** 3) - mum * u[1] /
             np.sqrt(((u[0] - pi1 * r12) ** 2 + u[1] ** 2) ** 3)),
            #  dotu[4] = that
            0]  #  dotu[5] = 0


dt = np.linspace(0.0, 535000.0, 535000.0)  #  secs to run the simulation
u = odeint(deriv, u0, dt)

x, y, z, x2, y2, z2 = u.T

# 3d plot
fig = pylab.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(x, y, z, color = 'r')

# 2d plot if you want
#fig2 = pylab.figure(2)
#ax2 = fig2.add_subplot(111)
#ax2.plot(x, y, color = 'r')

#  adding the moon
phi = np.linspace(0, 2 * np.pi, 100)
theta = np.linspace(0, np.pi, 100)
xm = rm * np.outer(np.cos(phi), np.sin(theta)) + r12 - pi2 * r12
ym = rm * np.outer(np.sin(phi), np.sin(theta))
zm = rm * np.outer(np.ones(np.size(phi)), np.cos(theta))
ax.plot_surface(xm, ym, zm, color = '#696969', linewidth = 0)
ax.auto_scale_xyz([-8000, 385000], [-8000, 385000], [-8000, 385000])
#  adding the earth
xe = re * np.outer(np.cos(phi), np.sin(theta)) - pi2 * r12
ye = re * np.outer(np.sin(phi), np.sin(theta))
ze = re * np.outer(np.ones(np.size(phi)), np.cos(theta))
ax.plot_surface(xe, ye, ze, color = '#4169E1', linewidth = 0)
ax.auto_scale_xyz([-8000, 385000], [-8000, 385000], [-8000, 385000])

ax.set_xlim3d(-20000, 385000)
ax.set_ylim3d(-20000, 80000)
ax.set_zlim3d(-50000, 50000)

pylab.show()

Matlab

secuencia de comandos:

days = 24*3600;
G = 6.6742e-20;
rmoon = 1737;
rearth = 6378;
r12 = 384400;
m1 = 5974e21;
m2 = 7348e19;
M = m1 + m2;
pi_1 = m1/M;
pi_2 = m2/M;
mu1 = 398600;
mu2 = 4903.02;
mu = mu1 + mu2;
W = sqrt(mu/r12^3);
x1 = -pi_2*r12;
x2 = pi_1*r12;

d0 = 200;
phi = -90;
v0 = 10.9148;
gamma = 20;
t0 = 0;
tf = 6.2*days;
r0 = rearth + d0;
x = r0*cosd(phi) + x1;
y = r0*sind(phi);

vx = v0*(sind(gamma)*cosd(phi) - cosd(gamma)*sind(phi));
vy = v0*(sind(gamma)*sind(phi) + cosd(gamma)*cosd(phi));
f0 = [x; y; vx; vy];

[t,f] = rkf45(@rates, [t0 tf], f0);
x = f(:,1);
y = f(:,2);
vx = f(:,3);
vy = f(:,4);
xf = x(end);
yf = y(end);
vxf = vx(end);
vyf = vy(end);

df = norm([xf - x2, yf - 0]) - rmoon;
vf = norm([vxf, vyf]);

plot(x,y)

funciones: 2 por separado

function [tout, yout] = rkf45(ode_function, tspan, y0, tolerance)
a = [0 1/4 3/8 12/13 1 1/2];
b = [0 0 0 0 0
     1/4 0 0 0 0
     3/32 9/32 0 0 0
     1932/2197 -7200/2197 7296/2197 0 0
     439/216 -8 3680/513 -845/4104 0
     -8/27 2 -3544/2565 1859/4104 -11/40];
c4 = [25/216 0 1408/2565 2197/4104 -1/5 0];
c5 = [16/135 0 6656/12825 28561/56430 -9/50 2/55];
if nargin < 4
    tol = 1.e-8;
else
    tol = tolerance;
end
t0 = tspan(1);
tf = tspan(2);
t = t0;
y = y0;
tout = t;
yout = y';
h = (tf - t0)/100; % Assumed initial time step.

while t < tf
    hmin = 16*eps(t);
    ti = t;
    yi = y;
    for i = 1:6
        t_inner = ti + a(i)*h;
        y_inner = yi;
        for j = 1:i-1
            y_inner = y_inner + h*b(i,j)*f(:,j);
        end
        f(:,i) = feval(ode_function, t_inner, y_inner);
    end
    te = h*f*(c4' - c5');
    te_max = max(abs(te));
    ymax = max(abs(y));
    te_allowed = tol*max(ymax,1.0);
    delta = (te_allowed/(te_max + eps))^(1/5);
    if te_max <= te_allowed
        h = min(h, tf-t);
        t = t + h;
        y = yi + h*f*c5';
        tout = [tout;t];
        yout = [yout;y'];
    end
    h = min(delta*h, 4*h);
    if h < hmin
        fprintf(['\n\n Warning: Step size fell below its minimum\n'...
            ' allowable value (%g) at time %g.\n\n'], hmin, t)
        return
    end
end

segunda función:

function dfdt = rates(t,f)
r12 = 384400;
m1 = 5974e21;
m2 = 7348e19;
M = m1 + m2;
pi_1 = m1/M;
pi_2 = m2/M;
mu1 = 398600;
mu2 = 4903.02;
mu = mu1 + mu2;
W = sqrt(mu/r12^3);
x1 = -pi_2*r12;
x2 = pi_1*r12;
x = f(1);
y = f(2);
vx = f(3);
vy = f(4);
r1 = norm([x + pi_2*r12, y]);
r2 = norm([x - pi_1*r12, y]);
ax = 2*W*vy + W^2*x - mu1*(x - x1)/r1^3 - mu2*(x - x2)/r2^3;
ay = -2*W*vx + W^2*y - (mu1/r1^3 + mu2/r2^3)*y;
dfdt = [vx; vy; ax; ay];
end