En la página titulada El problema de Euler de los tres cuerpos, estudiamos una situación en la que dos de ellos están fijos en el espacio y una tercera partícula de pequeña masa se mueve en el espacio circundante. Se supone el movimiento de la partícula se confina a un plano que contiene a los cuerpos fijos.

En esta página, se plantea un problema algo más complicado. Supongamos que el sistema formado por la Tierra y la Luna (o por el Sol y Júpiter) se mueven en órbitas circulares alrededor de su centro de masas común.

Datos del sistema Tierra-Luna:

• Masa de la Tierra, M_T=5.98·10²⁴ kg
• Radio de la Tierra, R_T=6370 km=6.37·10⁶ m
• Masa de la Luna, M_L=7.34·10²² kg
• Radio de la Luna, R_L=1740 km=1.74·10⁶ m
• Distancia entre la Tierra y la Luna, d=384400 km=384.4·10⁶ m
• Velocidad angular de rotación del sistema Tierra-Luna, ω=2.76·10^-6 rad/s

Movimiento de una partícula bajo la influencia de la Tierra y la Luna

En la página titulada Puntos de Lagrange estudiamos el movimiento circular uniforme de cada uno de los cuerpos celestes alrededor del centro de masas.

Sea d=r_T+r_L la distancia entre los centros de la Tierra y la Luna. Situamos un Sistema de Referencia Inercial con origen en el centro de masas del sistema Tierra-Luna. En el instante t, el ángulo que forma la recta que une el centro de la Tierra con el centro de la Luna forma un ángulo ωt, con el eje X. Donde la velocidad angular ω se calcula a partir de la dinámica del movimiento circular tal como se describe en la página previa

${\displaystyle {\omega }^2=G\frac{M_L+M_T}{d^3}\omega =2.67·{10}^{-6}\text{rad/s}}$

La posición de la partícula en el instante t es (x, y). La posición de la Luna es (r_L·cos(ωt), r_L·sin(ωt)). La posición de la Tierra es (-r_T·cos(ωt), -r_T·sin(ωt))

${\displaystyle r_T=\frac{M_L}{M_L+M_T}d,r_L=\frac{M_T}{M_L+M_T}d}$

La distancia entre la Luna y la partícula y entre la Tierra y la partícula son, respectivamente

${\displaystyle \begin{array}{l}{d}_L=\sqrt{{(x-r_L\cos (\omega t))}^2+{(y-r_L\text{sin}(\omega t))}^2}\\ d_T=\sqrt{{(x+r_T\cos (\omega t))}^2+{(y+r_T\text{sin}(\omega t))}^2}\end{array}}$

Las componentes de la fuerza de atracción que ejerce la Luna y la Tierra sobre la partícula son

${\displaystyle \begin{array}{l}{F}_x=-G\frac{M_{T}m}{d_T^2}\cos {\theta }_T-G\frac{M_{L}m}{d_L^2}\cos {\theta }_L=-G\frac{M_{T}m}{d_T^2}\frac{x+r_T\cos (\omega t)}{d_T}-G\frac{M_{L}m}{d_L^2}\frac{x-r_L\cos (\omega t)}{d_T}\\ F_y=-G\frac{M_{T}m}{d_T^2}\text{sin}{\theta }_T-G\frac{M_{L}m}{d_L^2}\text{sin}{\theta }_L=-G\frac{M_{T}m}{d_T^2}\frac{y+r_T\text{sin}(\omega t)}{d_T}-G\frac{M_{L}m}{d_L^2}\frac{y-r_L\text{sin}(\omega t)}{d_T}\end{array}}$

Aplicamos la segunda ley de Newton md²x/dt²=F_x, y md²y/dt²=F_y

${\displaystyle \begin{array}{l}\frac{d^{2}x}{d{t}^2}=-G\frac{M_T}{d_T^3}\left(x+r_T\cos (\omega t)\right)-G\frac{M_L}{d_L^3}\left(x-r_L\cos (\omega t)\right)\\ \frac{d^{2}y}{d{t}^2}=-G\frac{M_T}{d_T^3}\left(y+r_T\text{sin}(\omega t)\right)-G\frac{M_L}{d_L^3}\left(y-r_L\text{sin}(\omega t)\right)\end{array}}$

Las fuerzas se pueden derivar, alternativamente, de la energía potencial

${\displaystyle \begin{array}{l}V(x,y)=-G\frac{M_{T}m}{d_T}-G\frac{M_{L}m}{d_L}\\ \begin{array}{l}{F}_x=-\frac{\partial V}{\partial x}=-G\frac{M_{T}m}{d_T^3}\left(x+r_T\cos (\omega t)\right)-G\frac{M_{L}m}{d_L^3}\left(x-r_L\cos (\omega t)\right)\hfill \\ F_y=-\frac{\partial V}{\partial y}=-G\frac{M_{T}m}{d_T^3}\left(y+r_T\sin (\omega t)\right)-G\frac{M_{L}m}{d_L^3}\left(y-r_L\sin (\omega t)\right)\hfill \end{array}\end{array}}$

Sistema de Referencia en rotación

Establecemos un nuevo Sistema de Referencia No Inercial con el mismo origen que se mueve con velocidad angular ω, con respecto del Sistema de Referencia Inercial. El eje X_R del nuevo Sistema de Referencia es la recta que pasa por el centro de la Tierra y el Centro de la Luna.

En la figura, se muestra la relación entre las posición de la partícula (x, y) en el Sistema de Referencia Inercial y la posición (x_R, y_R) en el Sistema de Referencia en rotación.

x=x_R·cos(ωt)-y_R·sin(ωt)
y=x_R·sin(ωt)+y_R·cos(ωt)

Las distancias de la partícula al centro de la Tierra y al centro de la Luna son, respectivamente

${\displaystyle \begin{array}{l}{d}_L=\sqrt{{(x_R-r_L)}^2+y_R^2}\\ d_T=\sqrt{{(x_R+r_T)}^2+y_R^2}\end{array}}$

Calculamos las derivadas segundas de x e y

${\displaystyle \begin{array}{l}\frac{dx}{dt}=\omega \left(-x_R\text{sin}(\omega t)-y_R\cos (\omega t)\right)+\frac{d{x}_R}{dt}\cos (\omega t)-\frac{d{y}_R}{dt}\text{sin}(\omega t)\\ \frac{dy}{dt}=\omega \left(x_R\cos (\omega t)-y_R\text{sin}(\omega t)\right)+\frac{d{x}_R}{dt}\text{sin}(\omega t)+\frac{d{y}_R}{dt}\cos (\omega t)\end{array}}$

${\displaystyle \begin{array}{l}\frac{d^{2}x}{d{t}^2}=\{\begin{array}{l}-{\omega }^2\left(x_R\cos (\omega t)-y_R\text{sin}(\omega t)\right)-2\omega \left(\frac{d{x}_R}{dt}\text{sin}(\omega t)+\frac{d{y}_R}{dt}\cos (\omega t)\right)\\ +\frac{d^2{x}_R}{d{t}^2}\cos (\omega t)-\frac{d^2{y}_R}{d{t}^2}\text{sin}(\omega t)\end{array}\\ \frac{d^{2}y}{d{t}^2}=\{\begin{array}{l}-{\omega }^2\left(x_R\text{sin}(\omega t)+y_R\cos (\omega t)\right)+2\omega \left(\frac{d{x}_R}{dt}\cos (\omega t)-\frac{d{y}_R}{dt}\text{sin}(\omega t)\right)\\ +\frac{d^2{x}_R}{d{t}^2}\text{sin}(\omega t)+\frac{d^2{y}_R}{d{t}^2}\cos (\omega t)\end{array}\end{array}}$

Introduciendo estas expresiones en las dos ecuaciones diferenciales de segundo orden, e igualando los coeficientes de cos(ωt) y de sin(ωt), obtenemos el siguiente sistema de dos educaciones diferenciales

${\displaystyle \begin{array}{l}\frac{d^2{x}_R}{d{t}^2}-2\omega \frac{d{y}_R}{dt}-{\omega }^2{x}_R=-G\frac{M_T}{d_T^3}\left(x_R+r_T\right)-G\frac{M_L}{d_L^3}\left(x_R-r_L\right)\\ \frac{d^2{y}_R}{d{t}^2}+2\omega \frac{d{x}_R}{dt}-{\omega }^2{y}_R=-G\frac{M_T}{d_T^3}{y}_R-G\frac{M_L}{d_L^3}{y}_R\end{array}}$

Los dos términos en el primer miembro representan las componentes de las fuerzas de Coriolis y centrífuga por unidad de masa.

Constante del movimiento

Empecemos por el caso más simple, el movimiento de una partícula de masa m bajo la fuerza de atracción de un cuerpo fijo en el espacio de masa M, las ecuaciones del movimiento eran

${\displaystyle \begin{array}{l}\frac{d^{2}x}{d{t}^2}=-\frac{GM}{r^3}x{r}^2=x^2+y^2\\ \frac{d^{2}y}{d{t}^2}=-\frac{GM}{r^3}y\end{array}}$

Multiplicamos la primera ecuación por dx/dt y la segunda ecuación por dy/dt y teniendo en cuenta que

${\displaystyle \begin{array}{l}{v}^2=v_x^2+v_y^2\\ \frac{d{v}^2}{dt}=2v_x\frac{d{v}_x}{dt}+2v_y\frac{d{v}_y}{dt}=2\frac{dx}{dt}\frac{d^{2}x}{d{t}^2}+2\frac{dy}{dt}\frac{d^{2}y}{d{t}^2}=2\dot{x}\frac{d\dot{x}}{dt}+2\dot{y}\frac{d\dot{y}}{dt}\\ \frac{d(x^2+y^2)}{dt}=2x\frac{dx}{dt}+2y\frac{dy}{dt}=2x\dot{x}+2y\dot{y}\end{array}}$

Donde el símbolo punto encima de la letra indica derivada respecto del tiempo. Llegamos a

${\displaystyle \frac{1}{2}\frac{d}{dt}({\dot{x}}^2+{\dot{y}}^2)=-\frac{GM}{r^3}\frac{1}{2}\frac{d}{dt}(x^2+y^2)}$

Integramos cada uno de los dos miembros

${\displaystyle \begin{array}{l}\frac{1}{2}{\displaystyle \int d{v}^2=}\frac{1}{2}{v}^2\\ -\frac{1}{2}GM{\displaystyle \int \frac{d(r^2)}{r^3}}=-GM{\displaystyle \int \frac{dr}{r^2}}=GM\frac{1}{r}\\ \frac{1}{2}{v}^2-\frac{GM}{r}=C\end{array}}$

Donde C es una constante de integración, que es la energía por unidad de masa que se mantiene constante en todos los puntos de la trayectoria

La constante J de Jacobi

Volviendo de nuevo, a las dos ecuaciones diferenciales que describen el movimiento de la partícula en el Sistema de Referencia en rotación. Multiplicamos la primera ecuación por dx_R/dt y la segunda ecuación por dy_R/dt y sumamos ambas ecuaciones, llegamos a

$$\begin{gathered} {\displaystyle \frac{1}{2}\frac{d}{dt}({\dot{x}}_R^2+{\dot{y}}_R^2)-\frac{1}{2}{\omega }^2\frac{d}{dt}(x_R^2+y_R^2)= \\ -\frac{1}{2}\frac{G{M}_T}{d_T^3}\frac{d}{dt}\left({(x_R+r_T)}^2+y_R^2\right)-\frac{1}{2}\frac{G{M}_L}{d_L^3}\frac{d}{dt}\left({(x_R-r_L)}^2+y_R^2\right)} \end{gathered}$$

Integrando de forma similar a caso de una sola fuerza de atracción

${\displaystyle \frac{1}{2}({\dot{x}}_R^2+{\dot{y}}_R^2)-\frac{1}{2}{\omega }^2(x_R^2+y_R^2)=\frac{G{M}_T}{d_T}+\frac{G{M}_L}{d_L}+J}$

J se denomina constante de Jacobi

${\displaystyle J=\frac{1}{2}\left({\left(\frac{d{x}_R}{dt}\right)}^2+{\left(\frac{d{y}_R}{dt}\right)}^2\right)-\frac{1}{2}{\omega }^2\left(x_R^2+y_R^2\right)-\frac{G{M}_T}{\sqrt{{\left(x_R+r_T\right)}^2+y_R^2}}-\frac{G{M}_L}{\sqrt{{\left(x_R-r_L\right)}^2+y_R^2}}}$

Ecuaciones del movimiento adimensionales

Establecemos un sistema de unidades tal que la distancia se mide en unidades de la distancia entre la Tierra y la Luna d=r_T+r_L y el tiempo se mide en unidades tal que un periodo P (tiempo que tarda en dar una vuelta completa la Tierra y la Luna alrededor del c.m.) es 2π

${\displaystyle \omega =\sqrt{G\frac{M_L+M_T}{d^3}},P=\frac{2\pi }{\omega }}$

Llamando

${\displaystyle X_R=\frac{x_R}{d},Y_R=\frac{y_R}{d},\tau =\omega t,\alpha =\frac{r_T}{d},1-\alpha =\frac{r_L}{d}}$

Las dos ecuaciones diferenciales se convierten en

${\displaystyle \begin{array}{l}{\omega }^{2}d\frac{d^2{X}_R}{d{\tau }^2}-2{\omega }^{2}d\frac{d{Y}_R}{d\tau }-{\omega }^{2}d·X_R=-G\frac{M_T}{d^3{\left({\left(X_R+\frac{r_T}{d}\right)}^2+Y_R^2\right)}^{3/2}}d\left(X_R+\frac{r_T}{d}\right)-G\frac{M_L}{d^3{\left({\left(X_R-\frac{r_L}{d}\right)}^2+Y_R^2\right)}^{3/2}}d\left(X_R-\frac{r_L}{d}\right)\hfill \\ {\omega }^{2}d\frac{d^2{Y}_R}{d{\tau }^2}+2{\omega }^{2}d\frac{d{X}_R}{dt}-{\omega }^{2}d·Y_R=-G\frac{M_T}{d^3{\left({\left(X_R+\frac{r_T}{d}\right)}^2+Y_R^2\right)}^{3/2}}d·Y_R-G\frac{M_L}{d^3{\left({\left(X_R-\frac{r_L}{d}\right)}^2+Y_R^2\right)}^{3/2}}d·Y_R\hfill \end{array}}$

Simplificando, teniendo en cuenta que α=M_L/(M_L+M_T) y 1-α=M_T/(M_L+M_T)

${\displaystyle \{\begin{array}{l}\frac{d^2{X}_R}{d{\tau }^2}-2\frac{d{Y}_R}{d\tau }=X_R-(1-\alpha )\frac{(X_R+\alpha )}{{\left({\left(X_R+\alpha \right)}^2+Y_R^2\right)}^{3/2}}-\alpha \frac{(X_R-1+\alpha )}{{\left({\left(X_R-1+\alpha \right)}^2+Y_R^2\right)}^{3/2}}\\ \frac{d^2{Y}_R}{d{\tau }^2}+2\frac{d{X}_R}{dt}=Y_R-(1-\alpha )\frac{Y_R}{{\left({\left(X_R+\alpha \right)}^2+Y_R^2\right)}^{3/2}}-\alpha \frac{Y_R}{{\left({\left(X_R-1+\alpha \right)}^2+Y_R^2\right)}^{3/2}}\end{array}}$

Resolveremos el sistema de dos ecuaciones diferenciales por el procedimiento numérico ode45 con las siguientes condiciones iniciales, en el instante τ=0, la partícula parte de la posición (X₀, Y₀) con velocidad (u₀, v₀)

La constante de Jacobi se expresa en el nuevo sistema de unidades de la siguiente forma

${\displaystyle \begin{array}{l}J=\frac{1}{2}{\omega }^2{d}^2\left({\left(\frac{d{X}_R}{d\tau }\right)}^2+{\left(\frac{d{Y}_R}{d\tau }\right)}^2\right)-\frac{1}{2}{\omega }^2{d}^2\left(X_R^2+Y_R^2\right)-\frac{G{M}_T}{d\sqrt{{\left(X_R+\alpha \right)}^2+Y_R^2}}-\frac{G{M}_L}{d\sqrt{{\left(X_R-1+\alpha \right)}^2+Y_R^2}}\\ J={\omega }^2{d}^2\left\{\frac{1}{2}\left({\left(\frac{d{X}_R}{d\tau }\right)}^2+{\left(\frac{d{Y}_R}{d\tau }\right)}^2\right)-\frac{1}{2}\left(x_R^2+y_R^2\right)-\frac{1-\alpha }{\sqrt{{\left(X_R+\alpha \right)}^2+Y_R^2}}-\frac{\alpha }{\sqrt{{\left(X_R-1+\alpha \right)}^2+Y_R^2}}\right\}\end{array}}$

Movimiento en las proximidades del punto L₄ del sistema Sol-Júpiter

Consideremos el movimiento de los asteroides troyanos en el sistema formado por el Sol y el planeta Júpiter. Los datos de las masas de los dos cuerpos celestes son:

• Masa de Júpiter, M_L=1.898·10²⁷ kg
• Masa del Sol, M_T=1.989·10³⁰ kg
• Distancia media entre el Sol y Júpiter, d=5.203 UA=7.7837·10¹¹ m

El parámetro α=M_L/(M_L+M_T)=9.5334·10^-4. El periodo P=2π/ω o tiempo que tarda en dar una vuelta completa el sistema formado por el Sol y Júpiter es de P=3.7443·10⁸ s=11.8731 años. Por ejemplo, un tiempo τ=118, equivale a 118/(2π)=18.78 vueltas (ya que en dar una vuelta completa, el sistema Sol-Júpiter emplea un tiempo de 2π en unidades adimensionales). El tiempo real es de 18.78·11.8731=223 años

Las coordenadas del punto L₄ de Lagrange son

${\displaystyle a=\frac{1}{2}-\alpha ,b=\frac{\sqrt{3}}{2}}$

El segundo miembro del sistema de dos ecuaciones diferenciales, se obtiene derivando el potencial V(X_R, Y_R) respecto de X_R y respecto de Y_R

${\displaystyle \begin{array}{l}\{\begin{array}{l}\frac{d^2{X}_R}{d{\tau }^2}-2\frac{d{Y}_R}{d\tau }=-\frac{\partial V}{\partial X_R}\\ \frac{d^2{Y}_R}{d{\tau }^2}+2\frac{d{X}_R}{dt}=-\frac{\partial V}{\partial Y_R}\end{array}\\ V\left(X_R,Y_R\right)=-\frac{1}{2}(X_R^2+Y_R^2)-\frac{(1-\alpha )}{\sqrt{{\left(X_R+\alpha \right)}^2+Y_R^2}}-\frac{\alpha }{\sqrt{{\left(X_R-1+\alpha \right)}^2+Y_R^2}}\end{array}}$

Tomamos los primeros términos del desarrollo en serie del gradiente de potencial en las proximidades del punto L₄. Téngase en cuenta que el gradiente de potencial V es nulo en la posición de equilibrio L₄

${\displaystyle \begin{array}{l}\{\begin{array}{l}\frac{\partial V}{\partial X_R}=\left(X_R-a\right){\frac{{\partial }^{2}V}{\partial X_R^2}|}_{\begin{array}{l}{X}_R=a\\ Y_R=b\end{array}}+\left(Y_R-b\right){\frac{{\partial }^{2}V}{\partial X_R\partial Y_R}|}_{\begin{array}{l}{X}_R=a\\ Y_R=b\end{array}}\\ \frac{\partial V}{\partial Y_R}=\left(X_R-a\right){\frac{{\partial }^{2}V}{\partial X_R\partial Y_R}|}_{\begin{array}{l}{X}_R=a\\ Y_R=b\end{array}}+\left(Y_R-b\right){\frac{{\partial }^{2}V}{\partial Y_R^2}|}_{\begin{array}{l}{X}_R=a\\ Y_R=b\end{array}}\end{array}\\ \{\begin{array}{l}\frac{\partial V}{\partial X_R}=c_1\left(X_R-a\right)+c_3\left(Y_R-b\right)\\ \frac{\partial V}{\partial X_R}=c_3\left(X_R-a\right)+c_2\left(Y_R-b\right)\end{array}\end{array}}$

Hemos calculado las derivadas segundas en la sección Estabilidad la página titulada Puntos de Lagrange. Utilizaremos Math Symbolic de MATLAB para calcular el valor de las derivadas segundas en la posición (a, b) de L₄ es decir, los coeficientes c₁, c₂ y c₃

En el código, el parámetro α se representa por la variable simbólica k

syms x y k; % k es alfa
V=-(1-k)/sqrt((x+k)^2+y^2)-k/sqrt((x-1+k)^2+y^2)-(x^2+y^2)/2;
Vxx=diff(V,x,2);
c1=subs(Vxx,{x,y},{1/2-k,sqrt(3)/2});
Vyy=diff(V,y,2);
c2=subs(Vyy,{x,y},{1/2-k,sqrt(3)/2});
Vxy=diff(diff(V,x),y);
c3=subs(Vxy,{x,y},{1/2-k,sqrt(3)/2});

>> c1,c2
c1 =-3/4
c2 =-9/4
>> simplify(c3)
ans =(3*3^(1/2)*(2*k - 1))/4

El sistema de ecuaciones diferenciales del movimiento se escribe de forma aproximada

${\displaystyle \{\begin{array}{l}\frac{d^2{X}_R}{d{\tau }^2}-2\frac{d{Y}_R}{d\tau }\approx \frac{3}{4}\left(X_R-a\right)+\frac{3\sqrt{3}}{4}(1-2\alpha )\left(Y_R-b\right)\\ \frac{d^2{Y}_R}{d{\tau }^2}+2\frac{d{X}_R}{dt}\approx \frac{3\sqrt{3}}{4}(1-2\alpha )\left(X_R-a\right)+\frac{9}{4}\left(Y_R-b\right)\end{array}}$

Para describir el movimiento aproximado en las cercanías del punto L₄(a,b), ensayamos una solución de la forma,

${\displaystyle \begin{array}{l}{X}_R=a+A\exp (\lambda \tau )\\ Y_R=b+B\exp (\lambda \tau )\end{array}}$

Obtenemos un sistema homogéneo de dos ecuaciones con dos incógnitas A y B

${\displaystyle \begin{array}{l}\{\begin{array}{l}{\lambda }^{2}A{e}^{\lambda \tau }-2\lambda B{e}^{\lambda \tau }=\frac{3}{4}A{e}^{\lambda \tau }+\frac{3\sqrt{3}}{4}(1-2\alpha )B{e}^{\lambda \tau }\\ {\lambda }^{2}B{e}^{\lambda \tau }+2\lambda A{e}^{\lambda \tau }=\frac{3\sqrt{3}}{4}(1-2\alpha )A{e}^{\lambda \tau }+\frac{9}{4}B{e}^{\lambda \tau }\end{array}\\ \{\begin{array}{l}\left({\lambda }^2-\frac{3}{4}\right)A-\left(2\lambda +\frac{3\sqrt{3}}{4}(1-2\alpha )\right)B=0\\ \left(2\lambda -\frac{3\sqrt{3}}{4}(1-2\alpha )\right)A+\left({\lambda }^2-\frac{9}{4}\right)B=0\end{array}\end{array}}$

El determinante de los coeficientes deberá ser nulo

${\displaystyle \begin{array}{l}\left({\lambda }^2-\frac{3}{4}\right)\left({\lambda }^2-\frac{9}{4}\right)+\left(2\lambda +\frac{3\sqrt{3}}{4}(1-2\alpha )\right)\left(2\lambda -\frac{3\sqrt{3}}{4}(1-2\alpha )\right)=0\\ {\lambda }^4+{\lambda }^2+\frac{27}{4}\alpha \left(1-\alpha \right)=0\\ {\lambda }^2=\frac{-1\pm \sqrt{1-27\alpha \left(1-\alpha \right)}}{2}\end{array}}$

Para que el movimiento sea oscilatorio λ² tiene que ser negativo, es decir, λ un número imaginario. Para que esto ocurra, el discrimiante Δ=1-27α(1-α)≥0

Representamos la función y=1-27x(1-x)

>> fplot(@(x) 1-27*x.*(1-x),[0,1])
>> grid on

Para α<α_c el discriminante Δ es positivo

${\displaystyle {\alpha }_c=\frac{1-\sqrt{1-\frac{4}{27}}}{2}=\text{0}\text{.0385}}$

Cuando se cumple esta condición, las dos raíces imaginarias de la ecuación de cuarto grado son

${\displaystyle i{\Omega }_1=\sqrt{\frac{-1-\sqrt{\Delta }}{2}},i{\Omega }_2=\sqrt{\frac{-1+\sqrt{\Delta }}{2}}}$

donde $i=\sqrt{-1}$

Como α es pequeño, hacemos la aproximación, $\sqrt{1-x}\approx 1-\frac{x}{2}$

>> syms x;
>> taylor(sqrt(1-x),x)
ans =- (7*x^5)/256 - (5*x^4)/128 - x^3/16 - x^2/8 - x/2 + 1

Las frecuencias angulares son

${\displaystyle {\Omega }_1\approx \sqrt{\frac{1+\left(1-\frac{27}{2}\alpha \right)}{2}}=\sqrt{1-\frac{27}{4}\alpha },{\Omega }_2\approx \sqrt{\frac{1-\left(1-\frac{27}{2}\alpha \right)}{2}}=\sqrt{\frac{27}{4}\alpha }}$

Resolvemos el sistema de dos ecuaciones diferenciales por el procedimiento numérico ode45 con las siguientes condiciones iniciales, en el instante τ=0, la partícula parte de la posición (X₀=a+0.001, Y₀=b+0.002) cercana a L₄ con velocidad (u₀=0, v₀=0), en reposo. Representamos la trayectoria durante un tiempo τ=80π, que equivale a 40 vueltas del sistema Sol-Júpiter

MT=1.989e30; %masa del Sol
ML=1.898e27; %masa de Júpiter
alfa=ML/(ML+MT);

hold on
x0=[1/2-alfa+0.001,0,sqrt(3)/2+0.002, 0]; %condiciones iniciales
J0=(x0(2)^2+x0(4)^2)/2-(x0(1)^2+x0(3)^2)/2-(1-alfa)/sqrt((x0(1)+alfa)^2+x0(3)^2)-
alfa/sqrt((x0(1)-1+alfa)^2+x0(3)^2);
f=@(t,x) [x(2); 2*x(4)+x(1)-(1-alfa)*(x(1)+alfa)/((x(1)+alfa)^2+x(3)^2)^(3/2)-
alfa*(x(1)-1+alfa)/((x(1)-1+alfa)^2+x(3)^2)^(3/2); x(4); -2*x(2)+x(3)-
(1-alfa)*x(3)/((x(1)+alfa)^2+x(3)^2)^(3/2)-
alfa*x(3)/((x(1)-1+alfa)^2+x(3)^2)^(3/2)];
[t,x]=ode45(f,[0,80*pi],x0);
plot(x(:,1),x(:,3))
J=(x(:,2).^2+x(:,4).^2)/2-(x(:,1).^2+x(:,3).^2)/2-(1-alfa).
/sqrt((x(:,1)+alfa).^2+x(:,3).^2)-alfa./sqrt((x(:,1)-1+alfa).^2+x(:,3).^2);
plot(1/2-alfa,sqrt(3)/2,'ro','markersize',3,'markeredgecolor','r',
'markerfacecolor','r') %L4
hold off
axis square
grid on
xlabel('x')
ylabel('y')
title('Asteriodes troyanos')

El punto de color rojo señala la posición de L₄

Comprobamos que la constante de jacobi, J en unidades adimensionales, se mantiene aproximadamente constante

>> J0
J0 =   -1.4995
>> J
J =
   -1.4995
   -1.4995
   ...
   -1.4995
   -1.4995

Ajuste a la suma de dos funciones armónicas

Utilizamos el script para representar la coordenda X_R en función del tiempo τ. Vamos a comprobar que se ajusta a la función

${\displaystyle X_R(\tau )\approx A_0+A_1\cos \left({\Omega }_1\tau +{\phi }_1\right)+A_2\cos \left({\Omega }_2\tau +{\phi }_2\right)}$

Los datos (τ, X_R) proporcionados por el procemiento numérico ode45 se ajustan a dicha función utilizando fminsearch, o alternativamente, nlinfit (procedimiento de mínimos cuadrados)

MT=1.989e30; %masa del Sol
ML=1.898e27; %masa de Júpiter
alfa=ML/(ML+MT);

hold on
x0=[1/2-alfa+0.001,0,sqrt(3)/2+0.002, 0];
f=@(t,x) [x(2); 2*x(4)+x(1)-(1-alfa)*(x(1)+alfa)/((x(1)+alfa)^2+x(3)^2)^(3/2)-
alfa*(x(1)-1+alfa)/((x(1)-1+alfa)^2+x(3)^2)^(3/2); x(4); -2*x(2)+x(3)-
(1-alfa)*x(3)/((x(1)+alfa)^2+x(3)^2)^(3/2)-alfa*x(3)
/((x(1)-1+alfa)^2+x(3)^2)^(3/2)];
[t,x]=ode45(f,[0,80*pi],x0);
plot(t,x(:,1))

%ajuste
w1=sqrt(27*alfa/4);
w2=1-27*alfa/8;
f=@(a,t) a(1)+a(2)*cos(w1*t+a(3))+a(4)*cos(w2*t+a(4));      
error=@(a) sum((x(:,1)-f(a,t)).^2);
a0=[0.5, -0.1, 1.6, 0.01, 1.8];  %valor inicial
af=fminsearch(error,a0);
g=@(t) f(af,t);
fplot(g,[0,80*pi])
legend('e. d.','ajuste','location','southeast')

hold off
grid on
xlabel('x')
ylabel('y')
title('Asteriodes troyanos')

En color azul, la trayectoria, solución del sistema de dos ecuaciones diferenciales. En rosa, el ajuste a la función suma de dos armónicos de frecuencias angulares Ω₁ y Ω₂ y un término constante

>> af
af =    0.4790   -0.1474    1.7213   -0.0034    2.9514

Los valores de los parámetrros de ajuste son: A₀=0.4790, A₁=-0.1474, φ₁=1.7213, A₂=-0.0034 y φ₂= 2.9514

Trayectoria en el Sistema de Referencia en Rotación y en el Sistema de Referencia Inercial

Resolvemos el sistema de dos ecuaciones diferenciales por el procedimiento numérico ode45 con las siguientes condiciones iniciales, en el instante τ=0, la partícula parte de la posición (X₀=a+0.024, Y₀=b+0.055) cercana a L₄ con velocidad (u₀=0.0778, v₀=-0.0429). Representamos la trayectoria durante un tiempo τ=80π, que equivale a 40 vueltas del sistema Sol-Júpiter

MT=1.989e30; %masa del Sol
ML=1.898e27; %masa de Júpiter
alfa=ML/(ML+MT);

hold on
x0=[1/2-alfa+0.024,0.0778,sqrt(3)/2+0.055, -0.0429];
J0=(x0(2)^2+x0(4)^2)/2-(x0(1)^2+x0(3)^2)/2-(1-alfa)/sqrt((x0(1)+alfa)^2+x0(3)^2)-
alfa/sqrt((x0(1)-1+alfa)^2+x0(3)^2);
f=@(t,x) [x(2); 2*x(4)+x(1)-(1-alfa)*(x(1)+alfa)/((x(1)+alfa)^2+x(3)^2)^(3/2)-
alfa*(x(1)-1+alfa)/((x(1)-1+alfa)^2+x(3)^2)^(3/2); x(4); -2*x(2)+x(3)-
(1-alfa)*x(3)/((x(1)+alfa)^2+x(3)^2)^(3/2)-alfa*x(3)/((x(1)-1+
alfa)^2+x(3)^2)^(3/2)];
[t,x]=ode45(f,[0,80*pi],x0);
plot(x(:,1),x(:,3))
J=(x(:,2).^2+x(:,4).^2)/2-(x(:,1).^2+x(:,3).^2)/2-(1-alfa).
/sqrt((x(:,1)+alfa).^2+x(:,3).^2)-alfa./sqrt((x(:,1)-1+alfa).^2+x(:,3).^2);
plot(1/2-alfa,sqrt(3)/2,'ro','markersize',3,'markeredgecolor','r',
'markerfacecolor','r') %L4
plot(1/2-alfa,-sqrt(3)/2,'ro','markersize',3,'markeredgecolor','r',
'markerfacecolor','r') %L5

hold off
axis square
grid on
xlabel('x')
ylabel('y')
title('Asteriodes troyanos')

Los puntos de color rojo señalan las posiciones de L₄ (arriba) y L₅ (abajo)

Comprobamos que la constante de Jacobi J en unidades adimensionales, se mantiene aproximadamente constante

J0 =
   -1.5007
>> 
J =
   -1.5007
   -1.5007
   ....
   -1.5007
   -1.5007

Representamos esta trayectoria en el Sistema de Referencia Inercial. Como hemos demostrado al principio de esta página, la relación entre la posición de la partícula (X,Y) en el sistema inercial y la posición de dicha partícula (X_R,Y_R) en el sistema en rotación, en unidades adimensionales, es

${\displaystyle \{\begin{array}{l}X=X_R\cos \tau -Y_R\sin \tau \\ Y=X_R\sin \tau +Y_R\cos \tau \end{array}}$

MT=1.989e30; %masa del Sol
ML=1.898e27; %masa de Júpiter
alfa=ML/(ML+MT);

x0=[1/2-alfa+0.024,0.0778,sqrt(3)/2+0.055, -0.0429];
f=@(t,x) [x(2); 2*x(4)+x(1)-(1-alfa)*(x(1)+alfa)/((x(1)+alfa)^2+x(3)^2)^(3/2)
-alfa*(x(1)-1+alfa)/((x(1)-1+alfa)^2+x(3)^2)^(3/2); x(4); -2*x(2)+x(3)-
(1-alfa)*x(3)/((x(1)+alfa)^2+x(3)^2)^(3/2)-alfa*x(3)/((x(1)-1+
alfa)^2+x(3)^2)^(3/2)];
[t,x]=ode45(f,[0,80*pi],x0);
X=x(:,1).*cos(t)-x(:,3).*sin(t);
Y=x(:,1).*sin(t)+x(:,3).*cos(t);
plot(X,Y)
axis square
grid on
xlabel('x')
ylabel('y')
title('Asteriodes troyanos')

La trayectoria que sigue el asteroide es circular, pasando del exterior al interior en ciertos intervalos de tiempo. El lector pueda apreciar el comportamiento del asteroide en la siguiente animación

Actividades

Se elige

• Representar la trayectoria del asteroide en el Sistema de Referencia en Rotación, activando el botón de radio S. R. Rotación
• Representar la trayectoria del asteroide en el Sistema de Referencia Inercial, activando el botón de radio S. R. Inercial

En la parte inferior, se proporciona el tanto por ciento de error, de la constante J del movimiento. Cuando es mayor que la unidad, el programa se detiene.

Cuando se elige S.R. Inercial, observamos que hacia τ=10, el asteroide se mueve desde la circunferencia exterior hacia la interior, en la que permenece hasta aproximadamente hasta τ=102, que vuelve a la circunferencia exterior

Viaje de la Tierra a la Luna

Supongamos que la Tierra y la Luna están fijos en el espacio. Vamos a calcular la velocidad mínima con la que se debe disparar una bala, desde la superficie de la Tierra, a lo largo en la dirección que une los centros de los dos cuerpos celestes, para que llegue a la Luna.

Entre la Tierra y la Luna hay un punto de equilibrio situado a una distancia x del centro de la Tierra. Un objeto de masa m situado en dicho punto experimenta dos fuerzas de atracción iguales y de sentido contrario. Como la masa de la Tierra es mayor que la masa de la Luna, dicho punto x_e estará cerca de la Luna.

${\displaystyle G\frac{M_{T}m}{x^2}=G\frac{M_{L}m}{{(d-x)}^2}{x}_e=\frac{d}{1+\sqrt{M_L/M_T}}}$

La posición de equilibrio en la línea que une el centro de la Tierra con el centro de la Luna,  se encuentra a x_e=345.7 ·10⁶ m del centro de la Tierra.

G=6.67e-11; %constante G
MT=5.98e24; %Masa de la Tierra
ML=7.34e22; %masa de la Luna
RT=6.37e6; %radio de la Tierra
RL=7.34e6; %radio de la Luna
d=384.4e6; %distancia entre sus centros
x=(1:0.5:59)*RT;
F=-G*MT./x.^2+G*ML./(d-x).^2;
plot(x/RT,F/9.8)
grid on
ylim([-0.1,0.05])
xlabel('r/R_T')
ylabel('F/(mg)')
title('Fuerza que ejerce la Tierra y la Luna')

La energía potencial de la partícula bajo la acción de la fuerza de atracción de la Tierra y de la Luna es

${\displaystyle E_p(x)=-G\frac{M_{T}m}{x}-G\frac{M_{L}m}{d-x}}$

Vamos a calcular la velocidad v con la que tenemos que disparar la bala desde la superficie de la Tierra para que llegue al punto de equilibrio x_e con una velocidad nula, tal como se aprecia en la figura

La energía de la bala en la superficie de la Tierra x=R_T es

${\displaystyle E=\frac{1}{2}m{v}_T^2-G\frac{M_{T}m}{R_T}-G\frac{M_{L}m}{d-R_T}}$

La energía de la bala en la posición de equilibrio, es

 ${\displaystyle E=\frac{1}{2}m{v}_e^2-G\frac{M_{T}m}{x_e}-G\frac{M_{L}m}{d-x_e}}$

Poniendo v_e=0, y aplicando el principio de conservación de la energía, obtenemos la velocidad mínima de disparo v_T=11076.8 m/s

Comparamos esta velocidad, con la velocidad de escape de la superficie de la Tierra únicamente

${\displaystyle 0=\frac{1}{2}m{v}^2-G\frac{M_{T}m}{R_T}}$

v=11190.7 m/s

La Tierra y la Luna en movimiento de rotación

En el sistema Tierra-Luna, el eje X pasa por el centro de los dos cuerpos, y el origen se sitúa en el centro de masa del sistema, el centro de la Tierra a la izquierda a una distancia r_T=αd, y el centro de la Luna a una distancia r_L=(1-α)d

Una nave espacial describe una órbita circular alrededor de la Tierra a una altura h de su superficie a una distancia r de su centro. Aplicando la dinámica del movimiento circular uniforme obtenemos la velocidad v del satélite.

${\displaystyle m\frac{v^2}{r}=\frac{GMm}{r^2}v=\sqrt{\frac{GM}{r}}}$

Por ejemplo, un satélite artificial que circunda la Tierra a una altura de h=1000 km o bien r=7.37·10⁶ m lleva una velocidad de v=7356.6 m/s y tarda 1.75 horas en dar una vuelta completa. Un satélite geostacionario que tarda un día en dar una vuelta completa, se encuentra a una altura h=42250.5-3670=38580 km de altura y su velocidad es v=3072.5 m/s

Cuando la nave espacial se encuentra en la posición tal que hace un ángulo θ, con la dirección que une el centro de la Tierra y el centro de la Luna se encienden los cohetes por un breve intervalo de tiempo, dándole una velocidad adicional Δv en dirección tangente a la trayectoria

En el instante t=0, la nave espacial parte de la posición

x₀=-r_T+r·cosθ
y₀= r·sinθ

Con velocidad inicial

v_0x=-(v+Δv)·sinθ
v_0y=(v+Δv)·cosθ

Donde Δv es el incremento de velocidad proporcionado por los cohetes de la nave de forma casi instantánea.

La posición inicial de la nave se convierte en unidades adimensionales dividiéndola entre d (distancia Tierra-Luna). La velocidad inicial de la nave se convierte en unidades adimensionales dividiéndola entre ωd

Representamos la trayectoria de la nave espacial. La posición del centro de la Tierra en color azul y de la Luna en color rojo

MT=5.98e24; %masa de la Tierra
ML=7.34e22; %masa de la Luna
alfa=ML/(ML+MT);
d=384.4e6; %distancia Tierra-Luna
RT=6.37e6; %radio de la Tierra
w=sqrt(6.67e-11*(ML+MT)/d^3); %velocidad angular de rotación

ang=250*pi/180; %angulo de disparo
h=25480*1000; %altura en metros
inc=1190; %incremento de velocidad
x0=zeros(1,4); %condiciones iniciales
x0(1)=(-alfa*d+(h+RT)*cos(ang))/d; 
x0(3)=(h+RT)*sin(ang)/d;
v=sqrt(6.67e-11*MT/(RT+h)); %velocidad de la nave
x0(2)=-(v+inc)*sin(ang)/(w*d);
x0(4)=(v+inc)*cos(ang)/(w*d);
tspan=w*10*24*3600;
f=@(t,x) [x(2); 2*x(4)+x(1)-(1-alfa)*(x(1)+alfa)/((x(1)+alfa)^2+x(3)^2)^(3/2)-
alfa*(x(1)-1+alfa)/((x(1)-1+alfa)^2+x(3)^2)^(3/2); x(4); -2*x(2)+x(3)-
(1-alfa)*x(3)/((x(1)+alfa)^2+x(3)^2)^(3/2)-alfa*x(3)/((x(1)-1+alfa)^2+
x(3)^2)^(3/2)];
J0=(x0(2)^2+x0(4)^2)/2-(x0(1)^2+x0(3)^2)/2-(1-alfa)/sqrt((x0(1)+alfa)^2+x0(3)^2)
-alfa/sqrt((x0(1)-1+alfa)^2+x0(3)^2);
[~,x]=ode45(f,[0,tspan],x0);

hold on
plot(x(:,1),x(:,3))
J=(x(:,2).^2+x(:,4).^2)/2-(x(:,1).^2+x(:,3).^2)/2-(1-alfa)./sqrt((x(:,1)+alfa).^2+
x(:,3).^2)-alfa./sqrt((x(:,1)-1+alfa).^2+x(:,3).^2);
plot(-alfa,0,'bo','markersize',4,'markeredgecolor','b','markerfacecolor','b')
plot(1-alfa,0,'ro','markersize',4,'markeredgecolor','r','markerfacecolor','r')
hold off
grid on
xlabel('x')
ylabel('y')
title('Viaje de la Tierra a la Luna')

Comprobamos que el valor J se mantiene aproximadamente constante

>> J0
J0 =
   -1.2936
>> J
J =
   -1.2936
   -1.2936
   ...
   -1.2930
   -1.2930

Actividades

Se introduce

• La altura h en km de la órbita circular de la nave espacial sobre la superficie de la Tierra, en el control titulado Altura nave.

Se pulsa el botón titulado Nuevo y a continuación, ►

Se observa el movimiento del sistema Tierra-Luna respecto del Sistema de Referencia Inercial con origen en el centro de masa de ambos cuerpos.

Se observa el movimiento de la nave espacial alrededor de la Tierra y a la vez, la Tierra girando alrededor del c.m. del sistema.

• La posición angular θ, en grados, de la nave espacial respecto de la recta que une el centro de la Tierra y de la Luna, en el control titulado Angulo.

• La velocidad adicional Δv, en m/s que en dirección tangente a la trayectoria circular proporcionan los cohetes, en el control titulado Incremento velocidad.

Se pulsa el botón titulado Lanzar

En la parte superior derecha, se proporcionan los datos del tiempo en días desde el momento del lanzamiento, la posición x_R e y_R de la nave espacial en fracciones del radio de la Tierra, respecto del Sistema de Referencia en rotación

En la parte inferior, se proporciona el tanto por ciento de error, de la constante J del movimiento. Cuando es mayor que la unidad, el programa se detiene.

Nota: Se ha de advertir al lector, que como el paso de integración de las ecuaciones diferenciales del movimiento es variable, la velocidad del punto que representa la partícula, no se corresponde con la velocidad real de la partícula.

Referencias

O. L. de Lange, J. Pierrus. Solved Problems in Classical Mechanics. Analytical and numerical solutions with comments. Oxford University Press (2010). Question 11.9, pp. 348-352

Harmon N. J. Leidel C., Lindner J. F. Optimal exit: Solar escape as a restricted three-body problem. Am. J. Phys. 71(9) September 2003, pp. 871-877