Puntos de Lagrange

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

La posición del centro de masas del sistema Tierra–Luna se encuentra entre el centro de la Tierra y el centro de la Luna a una distancia r_T de la Tierra y r_L de la Luna, tal como se muestra en la figura. Como la masa de la Tierra es mayor que la masa de la Luna M_T>M_L luego, r_T<r_L. Situando el origen en el centro de masas.

$0 = \frac{- M_{T} r_{T} + M_{L} r_{L}}{M_{T} + M_{L}} d = r_{T} + r_{L}$

Despejamos r_T y r_L

$r_{T} = \frac{d \cdot M_{L}}{M_{L} + M_{T}} = 4 {.656·10}^{6} m r_{L} = \frac{d \cdot M_{T}}{M_{L} + M_{T}} = 379.344 {·10}^{6} m$

La posición del centro de masas del sistema Tierra-Luna está en el interior de la Tierra, más cerca de la superficie que del centro.

Supondremos que el centro de la Tierra y el centro de la Luna se mueven en órbitas circulares de radios r_T y r_L, respectivamente, alrededor de su centro de masas común.

El centro de la Luna describe una trayectoria circular de radio r_L bajo la acción de la fuerza de atracción de la Tierra, que dista d de su centro. Si ω es la velocidad angular constante. Aplicando la dinámica del movimiento circular uniforme a la Luna

$M_{L} ω^{2} r_{L} = G \frac{M_{T} M_{L}}{d^{2}}$

El centro de la Tierra describe una trayectoria circular de radio r_T bajo la acción de la fuerza de atracción de la Luna, que dista d de su centro. Si ω es la velocidad angular constante. Aplicando la dinámica del movimiento circular uniforme a la Tierra

$M_{T} ω^{2} r_{T} = G \frac{M_{T} M_{L}}{d^{2}}$

Despejamos la velocidad angular ω de una u otra ecuación

$ω^{2} = G \frac{M_{L} + M_{T}}{d^{3}} ω = 2.67 \cdot 10^{- 6} rad/s$

El periodo P=2π/ω=27.2 días

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

$\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}$

Energía potencial

En la página titulada Trabajo, energía cinética y energía potencial, vimos que la fuerza que ejerce un muelle elástico sobre una partícula, F=-kx es conservativa y la energía potencial es kx²/2, tomando como nivel cero de energía potencial cuando el muelle está sin deformar, x=0

Cuando una partícula de masa m está situada a una distancia r del eje de rotación de una plataforma que gira con velocidad angular constante ω, experimenta una fuerza centrífuga mω²r. Esta fuerza es conservativa y su energía potencial es -mω²r²/2, tomando como nivel cero de energía potencial cuando la partícula está en el eje r=0.

La energía potencial de la partícula de masa m es

$V (x_{R}, y_{R}) = - G \frac{M_{T} m}{\sqrt{{(x_{R} + r_{T})}^{2} + y_{R}^{2}}} - G \frac{M_{L} m}{\sqrt{{(x_{R} - r_{L})}^{2} + y_{R}^{2}}} - \frac{1}{2} m ω^{2} (x_{R}^{2} + y_{R}^{2})$

Para representar esta función, establecemos un sistema de unidades en el que la distancia se mide en unidades de la separación entre la Tierra y la luna d=r_T+r_L. Llamaremos

$\begin{array}{l} X_{R} = \frac{x_{R}}{r_{T} + r_{L}}, Y_{R} = \frac{y_{R}}{r_{T} + r_{L}} \\ α = \frac{r_{T}}{r_{T} + r_{L}} = \frac{M_{L}}{M_{L} + M_{T}}, 1 - α = \frac{r_{L}}{r_{T} + r_{L}} = \frac{M_{T}}{M_{L} + M_{T}} \end{array}$

La energía potencial es

$V (X_{R}, Y_{R}) = - m G {\frac{1}{r_{T} + r_{L}} \frac{M_{T}}{\sqrt{{(X_{R} + α)}^{2} + Y_{R}^{2}}} + \frac{1}{r_{T} + r_{L}} \frac{M_{L}}{\sqrt{{(X_{R} - 1 + α)}^{2} + Y_{R}^{2}}} + \frac{1}{2} \frac{M_{L}}{r_{T}} (X_{R}^{2} + Y_{R}^{2})}$

Teniendo en cuenta que r_TM_T=r_LM_L

$\begin{array}{l} V (X_{R}, Y_{R}) = - m G {\frac{1}{r_{T} + r_{L}} \frac{r_{L}}{r_{T}} \frac{M_{L}}{\sqrt{{(X_{R} + α)}^{2} + Y_{R}^{2}}} + \frac{1}{r_{T} + r_{L}} \frac{M_{L}}{\sqrt{{(X_{R} - 1 + α)}^{2} + Y_{R}^{2}}} + \frac{1}{2} \frac{M_{L}}{r_{T}} (X_{R}^{2} + Y_{R}^{2})} = \\ - m G \frac{M_{L}}{r_{T}} {\frac{(1 - α)}{\sqrt{{(X_{R} + α)}^{2} + Y_{R}^{2}}} + \frac{α}{\sqrt{{(X_{R} - 1 + α)}^{2} + Y_{R}^{2}}} + \frac{1}{2} (X_{R}^{2} + Y_{R}^{2})} \\ - m G \frac{M_{L} + M_{T}}{d} {\frac{(1 - α)}{\sqrt{{(X_{R} + α)}^{2} + Y_{R}^{2}}} + \frac{α}{\sqrt{{(X_{R} - 1 + α)}^{2} + Y_{R}^{2}}} + \frac{1}{2} (X_{R}^{2} + Y_{R}^{2})} \end{array}$

V(X_R, Y_R)→-∞ en las posiciones de los centros de la Tierra (-α, 0) y de la Luna (1-α, 0)
V(X_R, Y_R)→-∞ cuando X_R e Y_R son grandes como consecuencia del último término, debido a la rotación.

Para el sistema Tierra-Luna el parámetro α=M_L/(M_T+M_L)=0.0121. Representamos la función V(X_R, Y_R) que es proporcional a

$- \frac{(1 - α)}{\sqrt{{(X_{R} + α)}^{2} + Y_{R}^{2}}} - \frac{α}{\sqrt{{(X_{R} - 1 + α)}^{2} + Y_{R}^{2}}} - \frac{1}{2} (X_{R}^{2} + Y_{R}^{2})$

MT=5.98e24; %masa de la Tierra
ML=7.349e22; %masa de la Luna
alfa=ML/(ML+MT);

[x,y] = meshgrid(-1.5:0.05:1.5);
V=-(1-alfa)./sqrt((x+alfa).^2+y.^2)-alfa./sqrt((x-1+alfa).^2+y.^2)-(x.^2+y.^2)/2;
surf(x,y,V)
zlim([-2.5,-1])
xlabel('x'); 
ylabel('y'); 
zlabel('V(x,y)')
title('Energía potencial')
view(30,45)

Posiciones de equilibrio

Hay cinco posiciones de equilibrio o la fuerza es nula, son los denominados puntos de Lagrange

$\vec{F} = - \frac{\partial V}{\partial X_{R}} \hat{i} - \frac{\partial V}{\partial Y_{R}} \hat{j}$

Las dos derivadas parciales deben de ser nulas simultáneamente

$\begin{array}{l} \frac{\partial V}{\partial X_{R}} = (1 - α) \frac{(X_{R} + α)}{{({(X_{R} + α)}^{2} + Y_{R}^{2})}^{3 / 2}} + α \frac{(X_{R} - 1 + α)}{{({(X_{R} - 1 + α)}^{2} + Y_{R}^{2})}^{3 / 2}} - X_{R} = 0 \\ \frac{\partial V}{\partial Y_{R}} = (1 - α) \frac{Y_{R}}{{({(X_{R} + α)}^{2} + Y_{R}^{2})}^{3 / 2}} + α \frac{Y_{R}}{{({(X_{R} - 1 + α)}^{2} + Y_{R}^{2})}^{3 / 2}} - Y_{R} = 0 \end{array}$

Puntos L₄ y L₅

Las dos posiciones más sencillas de encontrar son las denominadas L₄ y L₅, son simétricas respecto del eje X (dirección que une la Tierra y la Luna). Cuando Y_R≠0, dividimos la segunda ecuación por Y_R y tenemos el sistema de ecuaciones

$\begin{array}{l} (1 - α) \frac{(X_{R} + α)}{{({(X_{R} + α)}^{2} + Y_{R}^{2})}^{3 / 2}} + α \frac{(X_{R} - 1 + α)}{{({(X_{R} - 1 + α)}^{2} + Y_{R}^{2})}^{3 / 2}} - X_{R} = 0 \\ (1 - α) \frac{1}{{({(X_{R} + α)}^{2} + Y_{R}^{2})}^{3 / 2}} + α \frac{1}{{({(X_{R} - 1 + α)}^{2} + Y_{R}^{2})}^{3 / 2}} - 1 = 0 \end{array}$

Multiplicamos la segunda por -(X_R+α) y la sumamos a la primera

$\begin{array}{l} - α \frac{1}{{({(X_{R} - 1 + α)}^{2} + Y_{R}^{2})}^{3 / 2}} + α = 0 \\ {(X_{R} - 1 + α)}^{2} + Y_{R}^{2} = 1 \end{array}$

Multiplicamos la segunda por -(X_R-1+α) y la sumamos a la primera

$\begin{array}{l} \frac{1 - α}{{({(X_{R} + α)}^{2} + Y_{R}^{2})}^{3 / 2}} - (1 - α) = 0 \\ {(X_{R} + α)}^{2} + Y_{R}^{2} = 1 \end{array}$

Los dos puntos de Lagrange L₄ y L₅ son las intersecciones de las circunferencias de centros: la Tierra (0,-α) y la Luna (0, 1-α) y radio unidad. Resolviendo el sistema de dos ecuaciones con dos incógnitas obtenemos,

$X_{R} = \frac{1}{2} - α, Y_{R} = \pm \frac{\sqrt{3}}{2}$

Representamos la energía potencial en las proximidades de unos de los dos puntos

MT=5.98e24; %masa de la Tierra
ML=7.349e22; %masa de la Luna
alfa=ML/(ML+MT);

%punto de Lagrange (0.8517,0)
X=0.3:0.01:0.7;
Y=0.7:0.01:1.1;
[x,y] = meshgrid(X,Y);
V=@(x,y) -(1-alfa)./sqrt((x+alfa).^2+y.^2)-alfa./sqrt((x-1+alfa).^2+y.^2)
-(x.^2+y.^2)/2;
hold on
mesh(x,y,V(x,y))
x1=1/2-alfa; y1=sqrt(3)/2; %punto de Lagrange
yy=ones(1,length(X))*y1;
plot3(X, yy,V(X,y1), 'color','r')
xx=ones(1,length(Y))*x1;
plot3(xx, Y,V(x1,Y), 'color','r')
plot3(x1,y1,V(x1,y1),'ko','markersize',4,'markeredgecolor','k',
'markerfacecolor','k')
hold off
xlabel('x'); 
ylabel('y'); 
zlabel('V(x,y)')
title('Potencial V(x,y)')
view(25,15)

Puntos L₁, L₂ y L₃

Los otros tres puntos de Lagrange L₁, L₂ y L₃, se encuentran en el eje X, (línea que une el centro de la Tierra y de la Luna)

Representamos la energía potencial V(X_R, 0)

MT=5.98e24; %masa de la Tierra
ML=7.349e22; %masa de la Luna
alfa=ML/(ML+MT);

V=@(x,y) -(1-alfa)./sqrt((x+alfa).^2+y.^2)-alfa./
sqrt((x-1+alfa).^2+y.^2)-(x.^2+y.^2)/2;
f=@(x) V(x,0);
fplot(f,[-1.5,1.5])
grid on
ylim([-2.5,-1])
xlabel('x'); 
ylabel('V(x,0)')
title('Energía potencial, Y=0')

Observamos que hay tres posiciones X_R en la que V(X_R, 0) presenta un máximo

Cuando Y_R=0, se cumple que $\frac{\partial V}{\partial Y_{R}} = 0$ y buscamos las raíces con la primera ecuación

$(1 - α) \frac{(X_{R} + α)}{{| X_{R} + α |}^{3}} + α \frac{(X_{R} - 1 + α)}{{| X_{R} - 1 + α |}^{3}} - X_{R} = 0$

Para X_R<-α, los dos numeradores son negativos, simplificamos

$- \frac{(1 - α)}{{(X_{R} + α)}^{2}} - \frac{α}{{(X_{R} - 1 + α)}^{2}} - X_{R} = 0$

>> syms x k;
>> z=(1-k)*(x-1-k)^2+k*(x+k)^2+x*(x-1-k)^2*(x+k)^2;
>> zz=collect(z)
zz =x^5 - 2*x^4 + ((k + 1)^2 + k^2 - 2*k*(2*k + 2))*x^3 + (2*k*(k + 1)^2 - 
k^2*(2*k + 2) + 1)*x^2 + ((2*k + 2)*(k - 1) + k^2*(k + 1)^2 + 2*k^2)*x +
 k^3 - (k - 1)*(k + 1)^2
>> z=subs(zz,k,alfa);
>> p=sym2poly(z);
>> roots(p)
ans =
  -1.0050 + 0.0000i
   0.4917 + 0.8589i
   0.4917 - 0.8589i
   1.0108 + 0.0791i
   1.0108 - 0.0791i

Vemos que se obtiene un polinomio zz de quinto grado, cuyos coefcientes dependen del valor de k=α. Se obtienen cinco raíces mediante roots una de las cuales es real, la primera

Para -α<X_R<1-α, (ente los centros de la Tierra y la Luna), el numerador de la izquierda es positivo y el de la derecha es negativo, simplificamos

$\frac{1 - α}{{(X_{R} + α)}^{2}} - \frac{α}{{(X_{R} - 1 + α)}^{2}} - X_{R} = 0$

>>  z=-(1-k)*(x-1-k)^2+k*(x+k)^2+x*(x-1-k)^2*(x+k)^2;
>> zz=subs(z,k,alfa);
>> p=sym2poly(zz)
>> roots(p)
ans =
  -0.5068 + 0.8614i
  -0.5068 - 0.8614i
   1.0809 + 0.1460i
   1.0809 - 0.1460i
   0.8517 + 0.0000i

Se obtienen cinco raíces, una es real, la última

Para X_R>1-α, el numerador de la izquierda y de la derecha son positivos, simplificamos

$\frac{1 - α}{{(X_{R} + α)}^{2}} + \frac{α}{{(X_{R} - 1 + α)}^{2}} - X_{R} = 0$

>> z=-(1-k)*(x-1-k)^2-k*(x+k)^2+x*(x-1-k)^2*(x+k)^2;
>> zz=subs(z,k,alfa);
>> p=sym2poly(zz);
>> roots(p)
ans =
  -0.5054 + 0.8636i
  -0.5054 - 0.8636i
   1.1729 + 0.0000i
   0.9189 + 0.1315i
   0.9189 - 0.1315i

Se obtienen cinco raíces, una es real, la tercera

Juntamos todas las porciones de código en un único script

MT=5.98e24; %masa de la Tierra
ML=7.349e22; %masa de la Luna
alfa=ML/(ML+MT);

syms k x;
z=(1-k)*(x-1-k)^2+k*(x+k)^2+x*(x-1-k)^2*(x+k)^2;
zz=subs(z,k,alfa);
p=sym2poly(zz);
raices=roots(p);
for i=1:length(raices)
    if isreal(raices(i))
        disp(raices(i))
    end
end

z=-(1-k)*(x-1-k)^2+k*(x+k)^2+x*(x-1-k)^2*(x+k)^2;
zz=subs(z,k,alfa);
p=sym2poly(zz);
raices=roots(p);
for i=1:length(raices)
    if isreal(raices(i))
        disp(raices(i))
    end
end

z=-(1-k)*(x-1-k)^2-k*(x+k)^2+x*(x-1-k)^2*(x+k)^2;
zz=subs(z,k,alfa);
p=sym2poly(zz);
raices=roots(p);
for i=1:length(raices)
    if isreal(raices(i))
        disp(raices(i))
    end
end

Las tres posiciones (X_R, 0) de los puntos de Lagrange L₃, L₁ y L₂ son las siguientes

   -1.0050
    0.8517
    1.1729

Observamos más de cerca la energía potencial V(X_R, Y_R) en las proximidades de L₁

MT=5.98e24; %masa de la Tierra
ML=7.349e22; %masa de la Luna
alfa=1/(1+MT/ML);

%punto de Lagrange (0.8517,0)
X=0.8:0.005:0.9;
Y=-0.1:0.005:0.1;
[x,y] = meshgrid(X,Y);
V=@(x,y) -(1-alfa)./sqrt((x+alfa).^2+y.^2)-alfa./sqrt((x-1+alfa).^2+y.^2)
-(x.^2+y.^2)/2;
hold on
mesh(x,y,V(x,y))
x1=0.8517;
plot3(x1,0,V(x1,0),'ko','markersize',4,'markeredgecolor','k',
'markerfacecolor','k')
yy=zeros(1,length(X));
plot3(X, yy,V(X,0), 'color','r')
xx=ones(1,length(Y))*x1;
plot3(xx, Y,V(x1,Y), 'color','r')
hold off
%zlim([-2.5,-1])
xlabel('x'); 
ylabel('y'); 
zlabel('V(x,y)')
view(30,45)

Se trata de un punto de ensilladura, es máximo a lo largo del eje X (línea que une el centro de la Tierra y de la Luna) y mínimo a lo largo de la dirección perpendicular que pasa por el punto

Estabilidad

Calculamos las derivadas segundas de V(X_R, Y_R)

$\begin{array}{l} \frac{\partial^{2} V}{\partial X_{R}^{2}} = - 3 \frac{(1 - α) {(X_{R} + α)}^{2}}{{({(X_{R} + α)}^{2} + Y_{R}^{2})}^{5 / 2}} + \frac{(1 - α)}{{({(X_{R} + α)}^{2} + Y_{R}^{2})}^{3 / 2}} - 3 \frac{α {(X_{R} - 1 + α)}^{2}}{{({(X_{R} - 1 + α)}^{2} + Y_{R}^{2})}^{5 / 2}} + \frac{α}{{({(X_{R} - 1 + α)}^{2} + Y_{R}^{2})}^{3 / 2}} - 1 \\ \frac{\partial^{2} V}{\partial Y_{R}^{2}} = - 3 \frac{(1 - α) Y_{R}^{2}}{{({(X_{R} + α)}^{2} + Y_{R}^{2})}^{5 / 2}} + \frac{(1 - α)}{{({(X_{R} + α)}^{2} + Y_{R}^{2})}^{3 / 2}} - 3 \frac{α Y_{R}^{2}}{{({(X_{R} - 1 + α)}^{2} + Y_{R}^{2})}^{5 / 2}} + \frac{α}{{({(X_{R} - 1 + α)}^{2} + Y_{R}^{2})}^{3 / 2}} - 1 \\ \frac{\partial}{\partial Y_{R}} (\frac{\partial V}{\partial X_{R}}) = \frac{\partial}{\partial X_{R}} (\frac{\partial V}{\partial Y_{R}}) = - 3 \frac{(1 - α) (X_{R} + α) Y_{R}}{{({(X_{R} + α)}^{2} + Y_{R}^{2})}^{5 / 2}} - 3 \frac{α (X_{R} - 1 + α) Y_{R}}{{({(X_{R} - 1 + α)}^{2} + Y_{R}^{2})}^{5 / 2}} \end{array}$

En este script, Math Symbolic de MATLAB realiza el cálculo de las derivadas segundas

MT=5.98e24; %masa de la Tierra
ML=7.349e22; %masa de la Luna
alfa=ML/(ML+MT);

%puntos de Lagrange
xL=[-1.0050, 0.8517, 1.1729, 1/2-alfa, 1/2-alfa];
yL=[0,0,0,sqrt(3)/2,-sqrt(3)/2];

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);
Vyy=diff(V,y,2);
Vxy=diff(diff(V,x),y);
for i=1:length(xL)
    disp([double(subs(Vxx,{k,x,y},{alfa, xL(i),yL(i)})), 
double(subs(Vyy,{k,x,y},{alfa, xL(i),yL(i)})),
double(subs(Vxy,{k,x,y},{alfa, xL(i),yL(i)}))])
end

   -3.0217    0.0109         0
  -13.6834    5.3417         0
   -6.0195    1.5097         0
   -0.7500   -2.2500   -1.2675
   -0.7500   -2.2500    1.2675

Los tres primeros puntos de Lagrange L₁, L₂ y L₃ cumplen

$\frac{\partial^{2} V}{\partial X_{R}^{2}} < 0, \frac{\partial^{2} V}{\partial Y_{R}^{2}} > 0, \frac{\partial^{2} V}{\partial X_{R} \cdot \partial Y_{R}} = 0$

Los dos restantes puntos de Lagrange L₄ y L₅ cumplen la condicción de máximo y por tanto, son inestables

$\frac{\partial^{2} V}{\partial X_{R}^{2}} < 0, \frac{\partial^{2} V}{\partial Y_{R}^{2}} < 0, \frac{\partial^{2} V}{\partial X_{R} \cdot \partial Y_{R}} \neq 0$

Si estuvieramos en un sistema de referencia fijo, concluiríamos que los cinco puntos de equilibrio son inestables, la fuerza $\vec{F} = - \frac{\partial V}{\partial X_{R}} \hat{i} - \frac{\partial V}{\partial Y_{R}} \hat{j}$ , sobre una partícula de masa m estaría dirigida hacia afuera desde cada uno de dichos puntos cuando la partícula está cerca de los mismos. Sin embargo, el sistema de referencia está en rotación y en particular, para los puntos L₄ y L₅ hay que tener en cuenta la fuerza de Coriolis que como es perpendicular al vector velocidad no realiza trabajo alguno. Es posible, por tanto, que partículas describan trayectorias alrededor de estos dos puntos, siguiendo curvas de nivel (potencial constante). La energía cinética y potencial de dichas partículas se mantendría aproximadamente constantes