El problema de Euler de los tres cuerpos

Los dos cuerpos fijos en el espacio tiene masas M1 y M2 y están separados una distancia d. Estableceremos un Sistema de Referencia, con el origen en el primer cuerpo y el eje X es la recta que une sus centros, tal como se muestra en la figura.

Cuando la partícula se encuentra en la posición (x, y), las fuerzas de atracción que ejerce cada uno de los cuerpos sobre la pequeña partícula se muestran en la figura, sus módulos valen

F 1 =G M 1 m r 1 2 F 2 =G M 2 m r 2 2 r 1 = x 2 + y 2 r 2 = (xd) 2 + y 2   

Descomponemos las fuerzas. Las ecuaciones del movimiento de la partícula de masa m son:

m d 2 x d t 2 = F 1 x r 1 F 2 xd r 2 m d 2 y d t 2 = F 1 y r 1 F 2 y r 2

Tenemos un sistema de dos ecuaciones diferenciales de segundo orden

d 2 x d t 2 = G M 1 ( x 2 + y 2 ) 3/2 x G M 2 ( (xd) 2 + y 2 ) 3/2 (xd) d 2 y d t 2 = G M 1 ( x 2 + y 2 ) 3/2 y G M 2 ( (xd) 2 + y 2 ) 3/2 y

Que se resuelven por procedimientos numéricos con las condiciones iniciales siguientes: En el instante t=0, la partícula se encuentra en la posición (x0, y0) y las componentes de su velocidad inicial son, (v0x, v0y).

Escalas

Antes de resolver el sistema de ecuaciones diferenciales por procedimientos numéricos, es conveniente prepararlas para que el ordenador no maneje números excesivamente grandes o pequeños.

Establecemos un sistema de unidades en el que la longitud se mide en unidades astronómicas, la distancia media entre el Sol y la Tierra. L=una UA=1.496·1011 m y el tiempo en unidades de año, P=un año= 365.26 días=3.156·107 s.

Supongamos el primer cuerpo es el Sol, M1= 1.98·1030 kg y la masa del segundo cuerpo es un múltiplo de la masa del Sol, M2=αM1

En el nuevo sistema de unidades x=Lx', t=P·t', la primera ecuación diferencial se escribe

d 2 x ' dt ' 2 L P 2 = G M 1 (x ' 2 +y ' 2 ) 3/2 x' L L 3 Gα M 1 ( (x'd') 2 +y ' 2 ) 3/2 (x'd') L L 3 d 2 x ' dt ' 2 = G M 1 P 2 L 3 ( x' (x ' 2 +y ' 2 ) 3/2 + α(x'd') ( (x'd') 2 +y ' 2 ) 3/2 )

y de forma similar la segunda ecuación diferencial.

Como L es el semieje mayor de la órbita de la Tierra alrededor del Sol, P es el periodo o tiempo que tarda en dar una vuelta completa y M1 es la masa del Sol. Por la tercera ley de Kepler, el término

G M 1 P 2 L 3 =4 π 2

Volviendo a la notación previa: x e y para la posición y t para el tiempo en el nuevo sistema de unidades. El sistema de ecuaciones diferenciales se escribe

d 2 x d t 2 =4 π 2 ( x ( x 2 + y 2 ) 3/2 + α(xd) ( (xd) 2 + y 2 ) 3/2 ) d 2 y d t 2 =4 π 2 ( y ( x 2 + y 2 ) 3/2 + αy ( (xd) 2 + y 2 ) 3/2 )

Principio de conservación de la energía

La energía total de la partícula es una constante del movimiento.

La energía de la partícula de masa m en el instante inicial t=0 es

E 0 = 1 2 m v 0 2 Gm M 1 r 1 Gm M 2 r 2 E 0 = 1 2 m( v 0x 2 + v 0y 2 ) Gm M 1 x 0 2 + y 0 2 Gm M 2 ( x 0 d) 2 + y 0 2

Cuando E0<0 la partícula permanece confinada en el espacio que rodea a los dos cuerpos. Cuando E0≥0 la partícula escapa al infinito

La energía de la partícula en el instante t es igual a

E= 1 2 m( v x 2 + v y 2 ) Gm M 1 x 2 + y 2 Gm M 2 (xd) 2 + y 2

En el nuevo sistema de unidades establecido, las magnitudes velocidad y posición están relacionadas del siguiente modo:

v=v’·L/P, x=x’·L, y=y’·L, d=d’·L

E m = 1 2 v ' 2 L 2 P 2 G M 1 L ( 1 x ' 2 +y ' 2 + α (x'd') 2 +y ' 2 ) E m = 1 2 v ' 2 L 2 P 2 4 π 2 L 2 P 2 ( 1 x ' 2 +y ' 2 + α (x'd') 2 +y ' 2 ) e ' = E m P 2 L 2 = 1 2 v ' 2 4 π 2 ( 1 x ' 2 +y ' 2 + α (x'd') 2 +y ' 2 )

Volviendo a la notación previa. Definimos una nueva energía e por unidad de masa en este sistema de unidades

e= 1 2 ( v x 2 + v y 2 )4 π 2 ( 1 x 2 + y 2 + α (xd) 2 + y 2 )

Se evalúa el cociente

| e e 0 e 0 |·100

que denominaremos tanto por ciento de error.

Ejemplos

Sea

Ejemplo 1

Ejemplo 2

Ejemplo 3

x0=zeros(1,4);
alpha=1; 
d=2.1;
x0(1)=1;  %abscisa x0
x0(2)=0;  %vx_0, velocidad inicial X
x0(3)=0; %ordenada y0
x0(4)=6.28; %vy_0, velocidad inicial Y
tspan=[0,2];
e0=(x0(2)^2+x0(4)^2)/2-4*pi^2*(1/sqrt(x0(1)^2+x0(3)^2)
+alpha/sqrt((x0(1)-d)^2+x0(3)^2));
fg=@(t,x)[x(2);-4*pi^2*(x(1)/(x(1)^2+x(3)^2)^(3/2)
+alpha*(x(1)-d)/((x(1)-d)^2+x(3)^2)^(3/2)); x(4);
-4*pi^2*(x(3)/(x(1)^2+x(3)^2)^(3/2)
+alpha*x(3)/((x(1)-d)^2+x(3)^2)^(3/2))];
[t,x]=ode45(fg,tspan,x0);
hold on
plot(x(:,1),x(:,3)) 
ef=(x(end,2)^2+x(end,4)^2)/2-4*pi^2*(1/sqrt(x(end,1)^2
+x(end,3)^2)+alpha/sqrt((x(end,1)-d)^2+x(end,3)^2));
fprintf('Variación de energía %1.3f\n', 100*abs((ef-e0)/e0));
xlabel('x')
ylabel('y')
grid on
axis equal
title('Problema de Euler de los tres cuerpos')
Variación de energía 9.111

Referencias

Wild W. J. Euler's three-body problem. Am. J. Phys. 48(4) April 1980, pp. 297-301