Simulaciones en el problema restringido circular

plano de los tres cuerpos

Fernando Chamizo

6 de junio de 2017

Resumen
Después de algunas consideraciones teóricas básicas sobre los puntos de Lagrange y las
curvas de velocidad cero, se estudia experimentalmente el comportamiento de algunas de
las órbitas que aparecen en el problema restringido circular plano de los tres cuerpos. El
código empleado para generar las figuras y los datos está disponible en:
http://www.uam.es/fernando.chamizo/physics/files/prog_calc_astr.tar.gz

Índice
1. El modelo teórico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

2. Curvas de velocidad cero . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

3. Esquema del código para el cálculo de las órbitas . . . . . . . . . . . . . . . . 7

4. Confinamiento alrededor de las masas . . . . . . . . . . . . . . . . . . . . . . . 10

5. Algunas órbitas periódicas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

6. Animaciones en el sistema inercial centro de masas . . . . . . . . . . . . . . 16

Referencias . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

1
1. El modelo teórico
Dos masas m1 y m2 bajo la acción de la atracción gravitatoria admiten órbitas elı́pticas con
foco común en el centro de masas, que supondremos el origen, y ejes mayores contenidos en una
misma recta. En el caso especial de excentricidad nula tenemos órbitas circulares recorridas
con velocidad angular constante ω0 dada por la tercera ley de Kepler
m1 + m2
(1) ω02 = G 3
r12

con r12 la distancia entre ambas masas. Esta situación se ajusta aproximadamente a los movi-
mientos de los planetas con la excepción de Mercurio ya que no presentan grandes excentrici-
dades:
' ♀ ♁ ♂ X Y Z [
e 0.2056 0.0068 0.0167 0.0934 0.0484 0.0542 0.0472 0.0086
El problema restringido circular de los tres cuerpos consiste en introducir una tercera ma-
sa m3 cuya influencia gravitatoria sobre m1 y m2 es despreciable. Por ejemplo, podrı́a corres-
ponder a una sonda espacial suficientemente cerca de Marte para que el sistema Sol-Marte sea
predominante.

Si r(t) describe el movimiento de m3 , las ecuaciones de Newton exigen

Gm1 r13 Gm2 r23
(2) r̈ = − 3 − 3 .
r13 r23
P
Escribamos r = xi ei donde R = {e1 , e2 , e3 } es un sistema de referencia que gira con
velocidad angular ω. Empleando ėi = ω × ei se deduce
X

(3) r̈ = ẍi ei + 2ẋi ω × ei + ẋi ω × (ω × ei ) .

El segundo sumando es la aceleración de Coriolis y el tercero la aceleración centrı́fuga. Digamos

que el plano orbital de m1 y m2 es el plano XY , entonces eligiendo e3 = b z y ω = ω0 e3 , en el
sistema de referencia R estas masas permanecerán estáticas: sus coordenadas serán constantes.
Escojamos ahora un sistema de unidades de masas y longitudes en el que

(4) G(m1 + m2 ) = 1 y r12 = 1

entonces (1) asegura ω0 = 1. Si denominamos (−µ, 0) y (1 − µ, 0) a las posiciones estáticas de

m1 y m2 , se cumple (recordemos que el centro de masas se habı́a supuesto el origen)
m2 m1
(5) µ= y 1−µ= .
m1 + m2 m1 + m2

2
Debido a la simetrı́a, no hay restricción en suponer 0 ≤ µ ≤ 1/2. Con esta notación, las
coordenadas en R del segundo miembro de (2) son
1 − µ µ
q q
(6) ∇ + con r1 = (x1 + µ)2 + x22 + x23 y r2 = (x1 − 1 + µ)2 + x22 + x23 .
r1 r2
Separando las coordenadas del primer miembro de (2) con (3), se deducen las ecuaciones de
movimiento en R
∂U ∂U ∂U 1 1−µ µ
(7) ẍ1 = 2ẋ2 + , ẍ2 = −2ẋ1 + , ẍ3 = con U = (x21 + x22 ) + + .
∂x1 ∂x2 ∂x3 2 r1 r2

Los puntos de Lagrange son aquellos en los que la masa m3 puede permanecer con velocidad
relativa nula con respecto a m1 y m2 . En el sistema no rotado, verı́amos que las tres masas
avanzan con la misma velocidad angular en sus órbitas circulares. De acuerdo con (7) los
puntos de Lagrange vienen determinados por la ecuación ∇U = ~0. La tercera coordenada lleva
a x3 = 0, como es lógico pensando en la dirección de las fuerzas. En consonancia con ello, en
este trabajo nos centraremos sobre todo en el problema plano en que m3 se mueve en el plano
orbital de m1 y m2 identificado con XY . Consecuentemente en lo sucesivo escribiremos x e y
en vez de x1 y x2 . De este modo el potencial efectivo U es
( p
1 2 2 1−µ µ r1 = (x + µ)2 + y 2 ,
(8) U = (x + y ) + + con p
2 r1 r2 r2 = (x − 1 + µ)2 + y 2 .

Reescribiendo esta fórmula como

1
1
1
(9) U = (1 − µ) r1−1 + r12 + µ r2−1 + r22 − µ(1 − µ),
2 2 2
las ecuaciones para los puntos de Lagrange son

(1 − µ) 1 − r−3 (x + µ) + µ 1 − r−3 (x − 1 + µ) = 0,



1 2
(10)
(1 − µ) 1 − r−3
y + µ 1 − r−3
y = 0.
1 2

Las soluciones que se ven a simple vista son las que son menos naturales desde el punto de
vista fı́sico: r1 = r2 = 1 que dan lugar a los puntos
√ √
1 3
1 3

(11) L4 = − µ, y L5 = − µ, − .
2 2 2 2
Geométricamente corresponden a que las tres masas formen un triángulo equilátero. El resto
de las soluciones de (10), tienen y = 0, es decir, corresponden al caso natural fı́sicamente en
que las tres masas están alineadas y las fuerzas sobre m3 se compensan con la fuerza centrı́fuga.

3
La numeración habitual es llamar L1 a la solución entre m1 y m2 , L2 a la que está a la derecha
de m2 y L3 el que está a la izquierda de m1 . Son soluciones de ecuaciones de quinto grado
[Cor98] que en general no admiten soluciones explı́citas. Una situación natural se da cuando
m1 es mucho mayor que m2 (por ejemplo en el citado ejemplo del Sistema Solar), lo que se
traduce por (5) en que µ es pequeño. En esa situación, se tienen las aproximaciones

L1 ≈ 1 − (µ/3)1/3 , 0 , L2 ≈ 1 + (µ/3)1/3 , 0 y L3 ≈ − 1 − 5µ/12, 0 .

(12)

Es fácil dar explicaciones sencillas de estas aproximaciones sin pasar al sistema de referencia R.
Por ejemplo, para L1 , una partı́cula en una órbita circular de radio r con frecuencia angular
ω0 sufre con nuestras unidades una fuerza centrı́fuga ω02 r = r y una atracción gravitatoria por
m1 y m2 dada por (1 − µ)/(r − µ)2 y µ/(r − 1 + µ)2 . El balance de fuerzas para r = 1 −  bajo
el ansatz µ = O(2 ) es
1−µ µ
(13) = +r que implica 1 + 2 + O(2 ) = µ−2 + 1 − .
(r − µ)2 (r − 1 + µ)2

Por tanto  = (µ/3)1/3 + O(µ2/3 ), que es coherente con el ansatz.

La siguiente función en SageMath1 genera la lista de los puntos de Lagrange “exactos”
resolviendo (10) con y = 0 y usando las expresiones (11):

def l l a g (mu) :
l a g = [ 0 , 0 , 0 , ((1 −2∗mu) / 2 , s q r t ( 3 ) / 2 ) , ((1 −2∗mu) / 2 , −s q r t ( 3 ) / 2 ) ]

f u = (1−mu) ∗ ( x−1+mu) ˆ 2 ∗ ( ( x+mu) ˆ3−1)+mu∗ ( x+mu) ˆ 2 ∗ ( ( x−1+mu) ˆ3+1)

l a g [ 0 ] = ( f i n d r o o t ( fu , −mu,1 −mu) , 0 )

f u = (1−mu) ∗ ( x−1+mu) ˆ 2 ∗ ( ( x+mu) ˆ3−1)+mu∗ ( x+mu) ˆ 2 ∗ ( ( x−1+mu) ˆ3−1)

l a g [ 1 ] = ( f i n d r o o t ( fu , 1−mu, 2 ) , 0 )

f u = (1−mu) ∗ ( x−1+mu) ˆ 2 ∗ ( ( x+mu) ˆ3+1)+mu∗ ( x+mu) ˆ 2 ∗ ( ( x−1+mu) ˆ3+1)

l a g [ 2 ] = ( f i n d r o o t ( fu , −2, −mu) , 0 )

return l a g

Parte del código citado después hará referencia a esta función.

2. Curvas de velocidad cero

A partir de (7) se obtiene
∂U ∂U ∂U dU
(14) ẋ1 ẍ1 + ẋ2 ẍ2 + ẋ3 ẍ3 = ẋ1 + ẋ2 + ẋ3 = .
∂x1 ∂x2 ∂x3 dt
1
Es un sistema algebraico computacional que combina varios paquetes contrastados escritos en diferentes
lenguajes a través de una sintaxis basada en Python.

4
Lo cual implica que la llamada constante de Jacobi

(15) cJ = 2U − (ẋ21 + ẋ22 + ẋ23 )

es una constante de movimiento.

Es posible traducir esta ley de conservación en una formulación Hamiltoniana [Ver09]. Ası́ en
el caso plano las ecuaciones de movimiento pueden derivarse de las ecuaciones de Hamilton
con
1 2 1−µ µ
px + p2y + ypx − xpy −


(16) H(x, y, px , py ) = −
2 r1 r2
donde px = ẋ − y, py = ẏ + x, son los momentos conjugados a x e y. La conservación de cJ con
x3 = 0 equivale a la conservación de la “energı́a” H.
Según la definición de la constante de Jacobi, las curvas implı́citas cJ = 2U están formadas
por todos los puntos en los que ẋ = ẏ = 0. Se dice que son curvas de velocidad cero. Hay que
aclarar que esta velocidad es relativa al sistema de referencia R, por tanto vistas “desde fuera”
corresponden a puntos en los que m3 se mueve con velocidad angular igual a la de m1 y m2 .
Estas curvas no corresponden a órbitas de m3 porque no satisfacen (7), lo cual es consecuencia
indirecta de que solo hay cinco puntos de Lagrange. La importancia de las curvas de velocidad
cero radica en que establecen la frontera entre las zonas cJ > 2U que son zonas prohibidas al
movimiento por corresponder a velocidades imaginarias, y las zonas permitidas cJ > 2U que
podrı́a alcanzar m3 en una trayectoria con condiciones iniciales correspondientes a cJ .
Las curvas de velocidad cero correspondientes a las diferentes constantes de Jacobi no son
entonces otra cosa que las curvas de nivel de U . Por otro lado, según vimos tras (7), los puntos
de Lagrange son puntos crı́ticos para U . Es fácil ver que en L4 y L5 se alcanza un mı́nimo ya
que son puntos crı́ticos y en ellos

∂2U 3 ∂ 2 U ∂ 2 U  ∂ 2 U 2 3 9 27 27
(17) 2
= > 0, 2 2
− = · − (1 − 2µ)2 = µ(1 − µ) > 0.
∂x 4 ∂x ∂y ∂x∂y 4 4 16 4
Por tanto las curvas de nivel que las contienen se reducen a un punto. Por otra parte, L1 , L2 y
L3 son puntos de silla de U [MD99, §3.7.1] y entonces existen curvas de velocidad cero crı́ticas
que pasan por ellos. La siguiente función en SageMath dada una lista de ı́ndices L dibuja estas
curvas crı́ticas donde i − 1 corresponde a Li .
def z e r o v e l c ( a , L) :
lag = l l a g (a)
f = x ˆ2 + y ˆ2 +2∗((1−a ) / s q r t ( ( x+a ) ˆ2+y ˆ 2 )+a / s q r t ( ( x−1+a ) ˆ2+y ˆ 2 ) )
P = p o i n t ( [ ( 0 , 0 ) ] , s i z e =0)
for k in L :
c = f ( x=l a g [ k ] [ 0 ] , y=l a g [ k ] [ 1 ] )
P += i m p l i c i t p l o t ( f −c , ( x , − 2 , 2) , ( y , − 2 , 2) , p l o t p o i n t s =400)
return P

5
Las siguientes figuras se obtienen invocando a zero_vel_c(0.1,L) con L=[0], L=[1], L=[2]
y L=[0,1,2] (esto es, la superposición de las figuras anteriores).

Para interpretar mejor los resultados se dibujan en rojo los cinco puntos de Lagrange y en
verde las dos masas m1 y m2 . Además se indican en azul en tamaño pequeño el origen y el
(1, 0) para que sirvan de referencia.
Es interesante ver la evolución de estas curvas cuando se disminuye el valor de µ, por
ejemplo cuando cambiamos µ = 0.1 por µ = 0.005 son

Para dibujar familias de curvas de velocidad cero (no crı́ticas) se emplea la siguiente función
que ajusta autmáticamente el rango aprovechando que en L4 se alcanza un mı́nimo y buscando
un máximo en los bordes del marco de la figura.
def z e r o v e l (mu, l =1.6) :
P = p o i n t ( [ ( 0 , 0 ) ] , s i z e =0)
f = x ˆ2 + y ˆ2 +2∗((1−mu) / s q r t ( ( x+mu) ˆ2+y ˆ 2 )+mu/ s q r t ( ( x−1+mu) ˆ2+y ˆ 2 ) )
cmin = f ( x=(1−2∗mu) / 2 , y=s q r t ( 3 ) / 2 )
cmax = min ( f ( x=−l , y=0) , f ( x=0 , y=−l ) )
N = 12
h = ( cmax−cmin ) /N
for c in s r a n g e ( cmin+h / 2 , cmax , h ) :
P += i m p l i c i t p l o t ( f −c , ( x,− l , l ) , ( y,− l , l ) , p l o t p o i n t s =400)
return P

El número de curvas viene dado por N y está arbitrariamente fijado a 12. el parámetro opcional
l especifica el tamaño del cuadrado dentro del que se dibujan las curvas.

6
 Aquı́ se muestran algunos ejemplos para diferentes valores del único argumento obligato-
rio µ.

µ = 0.5 µ = 0.1 µ = 0.05 µ = 0.01

La última figura muestra un hecho fı́sicamente claro: cuando la masa m2 es despreciable y

no estamos en sus cercanı́as, el sistema tiende a tener simetrı́a radial y las curvas de velocidad
cero son casi circulares.
En [Dan88, p. 259] (ver también [MD99, p. 81]) se puede ver una representación tridimen-
sional de las curvas de velocidad cero como curvas de nivel de 2U .

3. Esquema del código para el cálculo de las órbitas

La existencia de la formulación Hamiltoniana permitirı́a en principio el uso de integradores
simplécticos pero como el Hamiltoniano (16) no es separable podrı́an ser costosos numérica-
mente por ser implı́citos (parece que hay alternativas explı́citas aumentando artificialmente
la dimensión y componiendo soluciones [Tao16]). En aras de la simplicidad, usamos el clásico
esquema explı́cito Runge-Kutta RK4 que, dicho sea de paso, hemos visto que se emplea en algún
artı́culo de investigación reciente [Sic10]. La implementación de esta parte es en C.
Las coordenadas en el espacio de fases de velocidades (posiciones y velocidades) están
representadas por la estructura ph_sp que responde a

typedef struct ph sp {
double t x;
double t y;
double t vx ;
double t vy ;
} ph sp ;

Para entender el esquema del código C, examinemos la rutina main:

1 int main ( void ) {
2 int i ;
3 FILE ∗ f p ;
4
5 ph sp ∗ c o o r d r k 4 = ( ph sp ∗ ) c a l l o c ( N STEPS+1, sizeof ( ph sp ) ) ;

7
 6
7 runge kutta4 ( coord rk4 ) ;
8
9 f p = f o p e n ( " r3b . dat " , "w" ) ;
10 if (NULL == f p ) {
11 f p r i n t f ( s t d e r r , " File failed \n" ) ;
12 return FAILURE ;
13 }
14 f p r i n t f ( fp , " %.15e\n %.15e\n %d\n" , MU, CJ , N STEPS+1) ;
15 for ( i =0; i <N STEPS+1; i ++){
16 f p r i n t f ( fp , " %e , %e\n" , ( double ) c o o r d r k 4 [ i ] . x ,
,→ ( double ) c o o r d r k 4 [ i ] . y ) ;
17 }
18 f c l o s e ( fp ) ;
19
20 f r e e ( coord rk4 ) ;
21
22 return SUCCESS ;
23 }

La lı́nea 5 reserva el espacio en memoria para almacenar los datos de la órbita y la 7 llama al
método numérico de Runge-Kutta. El resultado quedará acumulado en la estructura de datos
coord_rk4. Las lı́neas 9–18 extraen las coordenadas x e y y las escribe en el fichero de datos
r3b.dat que se encabeza con los valores de µ, la constante de Jacobi cJ y el número de pasos
más uno. Los dos primeros los utilizará el código que dibuja los datos y el último evitará errores
de lectura. Finalmente, en la lı́nea 20 se libera la memoria.

La rutina runge_kutta4 incorpora el método Runge-Kutta 4 habitual para resolver una

ecuación diferencial y 0 (x) = f (x, y) [SB02, (7.2.1.14)]. La única dificultad en la implementación
es que hay tratar el problema vectorial. En relación con ello, se añade una función que halla
combinaciones lineales de la forma ~v + λw ~ pero nada de esto es suficientemente reseñable para
copiar aquı́ el código concreto. La función f está aquı́ determinada por la dinámica (7), que
en el espacio de fases de velocidades pasa a ser de primer orden. La función que evalúa f se
denomina dynamics y responde al código

ph sp dynamics ( ph sp ph ) {
ph sp outph ;
d o u b l e t r 1 c = ph . x+MU;
d o u b l e t r 2 c = ph . x−1.0+MU;
r 1 c = s q r t ( r 1 c ∗ r 1 c+ph . y∗ph . y ) ;
r 1 c ∗= r 1 c ∗ r 1 c ;
r 2 c = s q r t ( r 2 c ∗ r 2 c+ph . y∗ph . y ) ;
r 2 c ∗= r 2 c ∗ r 2 c ;
/* velocity */
outph . x = ph . vx ;
outph . y = ph . vy ;
/* Force */
outph . vx = −(1.0−MU) ∗ ( ph . x+MU) / r1c−MU∗ ( ph . x−1.0+MU) / r 2 c ;
outph . vy = −(1.0−MU) ∗ph . y/ r1c−MU∗ph . y/ r 2 c ;

8
 /* Efective force ( Coriolis correction ) */
outph . vx += 2 . 0 ∗ ph . vy + ph . x ;
outph . vy += −2.0∗ph . vx + ph . y ;
return outph ;
}

Los valores que determinan la órbita aparecen como parámetros a través de la directiva
#define dentro de un fichero de cabecera .h que, como veremos, es generado automáticamente
por un programa

SageMath.

Los parámetros son µ y las posiciones y velocidades iniciales
x(0), y(0) y ẋ(0), ẏ(0) (indicados a la izquierda para un ejemplo de órbita horseshoe que
veremos más adelante)

# define MU 0 . 0 0 0 9 5 3 8 7 5
# define X0 −0.97668
# define Y0 0 . 0 # define N ORB 3 0 . 0
# define VX0 0 . 0 # define DT 0 . 0 1
# define VY0 −0.06118

Por otra parte, hay dos parámetros relativos al método numérico (columna de la derecha),
N_ORB indica el tiempo en unidades de órbitas de las masas m1 y m2 , mientras que DT indica
la discretización en el tiempo. A partir de ellos se define el número de pasos con
# define N STEPS ( ( int ) ( 2 . 0 ∗ M PI∗N ORB/DT) )

que es conveniente utilizar directamente en la rutina runge_kutta4.

Para mayor comodidad, se interactúa con el código C anterior a través de una función,
orbit, definida en SageMath (de hecho en la parte no gráfica es en su mayorı́a Python puro)
que tiene tres argumentos, el último de ellos opcional. El primero es una lista con los parámetros
orbitales y el segundo con los del método numérico, ordenados como se ha indicado antes. Ası́,
la órbita puesta como ejemplo se invoca con
Lorb = [ 0 . 0 0 0 9 5 3 8 7 5 , −0.97668 , 0 . 0 , 0 . 0 , − 0 . 0 6 1 1 8 ]
Lnum = [ 3 0 . 0 , 0 . 0 1 ]
o r b i t ( Lorb , Lnum )

El tercer parámetro opcional es una lista de la forma [xmin,xmax,ymin,ymax,filen] que

indica las dimensiones de la ventana que se va dibujar y el nombre del fichero en el que se
guardará la figura. La definición de orbit es
def o r b i t ( Lor , Lnu , Lpl = [ − 1 . 6 , 1 . 6 , − 1 . 6 , 1 . 6 , ’./ r3b . eps ’ ] ) :
d a t a t o f i l e ( Lor , Lnu )
r u n a n d p l o t ( Lpl )

donde data_to_file crea el fichero de cabecera .h para que lo emplee el programa en C

y run_and_plot llama a dicho programa y lee y representa los resultados almacenados en

9
el fichero r3b.dat. Además de los datos se representa también la curva de velocidad cero
correspondiente a la órbita para que sean claras las zonas prohibidas.
Una pequeña nota técnica es que para evitar imágenes muy “pesadas” que no añaden
información, en run_and_plot solo se aprovechan los datos de r3b.dat que podrı́an dar lugar
a pı́xeles distintos con resolución razonable. Esto se hace con el if de
for k in r a n g e ( n s t e p s ) :
l i n e=f . r e a d l i n e ( )
xy = l i n e . s p l i t ( ’,’ )
xy = ( f l o a t ( xy [ 0 ] ) , f l o a t ( xy [ 1 ] ) )
if abs ( l x y [0] − xy [ 0 ] ) / r e s x + abs ( l x y [1] − xy [ 1 ] ) / r e s y > 1 :
l x y = xy
Ldat . append ( xy )

Aquı́ res_x y res_y se calculan automáticamente para que haya a lo más 800 pı́xeles por
coordenada.

4. Confinamiento alrededor de las masas

Como ya hemos mencionado, en L1 = (x1 , 0) el potencial U alcanza un punto de silla y la
presencia a izquierda y derecha, en (−µ, 0) y (1 − µ, 0), de las singularidades correspondientes
a las masas m1 y m2 causa que la curva de velocidad cero crı́tica que pasa por L1 tenga dos
componentes conexas, una exterior con forma redondeada y otra interior con forma de ocho
tumbado (aquı́ [Dan88, Fig. 8.6] puede dar cierta intuición de por qué ocurre esto). Según µ
varı́a entre 1/2 y 0, la coordenada x1 se traslada hacia la derecha y consecuentemente el lazo
derecho de la componente interior se va reduciendo.

µ = 0.4 µ = 0.3 µ = 0.1 µ = 0.01

El valor de la constante de Jacobi para estas curvas crı́ticas es
2(1 − µ) 2µ
(18) c1 (µ) = x21 + + .
x1 + µ 1 − µ − x1

Cuando µ → 0+ se tiene x1 → 1− y por tanto c1 (0+ ) = 3+ . En el otro extremo, si µ → 0.5−

entonces x1 → 0+ y c1 (0.5− ) = 4− . Para comprobar que c1 (µ) crece con µ, llamemos F (µ, x1 )

10


al segundo miembro de (18). Tenemos que ver que la derivada de c1 (µ) = F µ, x1 (µ) no se
anula para 0 < µ < 1/2 y ası́ c1 será monótona. Utilizando que, por definición, en L1 se cumple
∂F/∂x1 = 0, la anulación de c01 (µ) implicarı́a
∂F ∂F dx1 ∂F 2(1 + x1 ) 2(1 − x1 )
(19) c01 (µ) = + = =− 2
+ = 0,
∂µ ∂x1 dµ ∂µ (x1 + µ) (1 − µ − x1 )2
y de aquı́
1 1 2T 2T
(20) µ= + T− donde se ha escrito x1 = con 0 ≤ T ≤ 1.
2 2 1 + T2 1 + T2
Sustituyendo en la primera ecuación de (10) con x = x1 e y = 0 y extrayendo un factor
T (T − 1), se llega una ecuación de tercer grado que no tiene solución en T ∈ (0, 1).

En [MD99, (3.94)] se menciona la aproximación asintótica c1 (µ) ∼ 3 + (9µ)2/3 − 10µ/3 para

µ → 0+ pero los siguientes datos muestran que en realidad es un muy buen sustituto de la
fórmula no explı́cita (18) en todo el rango 0 ≤ µ ≤ 1/2:

c1 (0.50) = 4.00000 c1 (0.20) = 3.80465

∼ 4.05901 ∼ 3.81306
c1 (0.40) = 3.98091 c1 (0.10) = 3.59695
∼ 4.01559 ∼ 3.59884
c1 (0.30) = 3.92015 c1 (0.01) = 3.16764
∼ 3.93899 ∼ 3.16750

Si consideramos una curva de velocidad cero para un valor de la constante de Jacobi cJ

mayor que c1 , el ocho se romperá en dos óvalos alrededor en m1 y m2 . Si por el contrario, cJ es
algo menor que c1 entonces la “cintura” del ocho se ensanchará pero en este caso hay un lı́mite:
cuando cJ llegue al valor de c2 correspondiente a la curva de velocidad cero que contiene a L2 ,
entonces ya dejará de haber dos componentes conexas.

c1 < cJ c2 < cJ < c1 cJ < c2

Las figuras corresponden a µ = 0.2 que implica c1 = 3.8047 y c2 = 3.5524. Los valores elegidos
de cJ son respectivamente 3.9, 3.7 y 3.5.

11
 En los dos primeros casos tenemos órbitas confinadas alrededor de m1 y m2 porque, recor-
demos, las curvas de velocidad cero establecen fronteras para las regiones permitidas.

Ya una situación simétrica da una idea de la complejidad de los movimientos en el problema

de los tres cuerpos, incluso en la versión restringida que manejamos aquı́. Si tomamos µ = 1/2
con (x0 , y0 ) = (0, 0) y (vx0 , vy0 ) = (0.5, 0) y esperamos cuatro órbitas de m1 y m2 entonces la
trayectoria de m3 es

Lorb = [ 0 . 5 , 0 . 0 , 0 . 0 , 0 . 5 , 0 . 0 ]
Lnum = [ 5 . 0 , 0 . 0 0 0 1 ]

c1 = 4, cJ = 3.75, c2 = 3.4568

Si nos fijamos en la salida desde el punto inicial, que es el origen, observamos un desvı́o en
la trayectoria hacia abajo a pesar de la velocidad inicial es horizontal. La explicación es que la
simetrı́a y ↔ −y es solo aparente pues en el sistema inercial centro de masas m1 y m2 están
rotando. La desviación se debe al término de Coriolis en (3).

Si acercamos cJ a c1 , la región permitida tenderá a estrangularse y la comunicación entre m1

y m2 será más inusual. Eso es lo que se trata de ilustrar con la primera de las siguientes figuras,
en la que c1 − cJ = 0.021975. Nótese el retroceso inicial debido a que vx0 no es suficientemente
grande.

(x0 , y0 ) = (0, 0) (vx0 , vy0 ) = (0.6, 0.12) (x0 , y0 ) = (0.1, 0) (vx0 , vy0 ) = (0.85, 0.57)
µ = 0.4, cJ = 3.958933 µ = 0.29, cJ = 3.558268
Si, por el contrario, cJ es ligeramente menor que c2 , se abrirá un hueco hacia la región exterior
no acotada y m3 podrı́a escapar de los alrededores de m1 y m2 , como muestra la segunda figura
en la que c2 − cJ = 0.003822.

12
 En el caso c1 < cJ el movimiento está confinado al óvalo que contiene a una de las masas
siempre que la posición inicial esté en él. Tomemos para todos los ejemplos siguientes µ = 0.2
y cJ = 4.098992. En este caso los dos óvalos tiene forma casi circular debido a que cJ y c1 no
están próximos. Por (15), la constante de Jacobi depende del módulo de la velocidad pero no
de su dirección, sin embargo las órbitas, por supuesto, pueden ser bien distintas. Partiendo del
mismo punto (x0 , y0 ) = (0.1, 0.08), para el tiempo de tres órbitas de m1 y m2 se obtiene, con
las velocidades iniciales indicadas,

(vx0 , vy0 ) = (0.691587, 1.07708) (vx0 , vy0 ) = (−0.532668, 1.16390)

Para conseguir que m3 esté confinada alrededor de m2 podemos tomar un (x0 , y0 ) cercano
a (1 − µ, 0), la singularidad de U , y después compensar con una velocidad inicial adecuada
en (15) para que cJ sea 4.098992 como antes. Por ejemplo, en la figura de la derecha se ha
tomado (x0 , y0 ) = (0.9, 0) que habrı́a dado la constante 6.264545 con velocidad inicial nula y
que con vx20 + vy20 = 2.16557 se convierte en la cJ de más arriba.

(vx0 , vy0 ) = (1.12331, 0.613652) (vx0 , vy0 ) = (1.29144, 0.705515)

La figura de la izquierda es un ejemplo más correspondiente a (x0 , y0 ) = (0.1, 0.08). Nótese

que los puntos cuspidales están en el borde del óvalo, esto se debe a que en ellos ẋ = ẏ = 0 y
entonces tienen que pertenecer a la curva de velocidad cero.

13
5. Algunas órbitas periódicas
Resulta un poco chocante que en un sistema caótico como el de los tres cuerpos se puedan
encontrar soluciones periódicas (y lo es mucho más que en 2000 se hallara una sencilla para
el problema sin restricciones [CM00]) pero desde el punto de vista numérico no lo es tanto
pensando en métodos de disparo, esto es, intentando por aproximaciones sucesivas lanzar m3
de modo que la posición y velocidades finales se aproximen a las iniciales. Habitualmente se
emplea la simetrı́a (t, x, y, ẋ, ẏ) ↔ (−t, x, −y, −ẋ, ẏ) del problema, por ello las condiciones
iniciales son a menudo de la forma y0 = vx0 = 0 [BM05].
En [Per17] se prueba una especie de ley de Kepler para las órbitas periódicas y entre los
ejemplos que se mencionan allı́ representamos los dos siguientes:

x0 = 0.732035 x0 = 1.1378
y0 = 0.713203 y0 = 0
vx0 = 0.112027 vx0 = 0
vy0 = −0.114985 vy0 = −0.174265

cJ = 2.974994 cJ = 3.022873

Los valores de µ son 9.53875 · 10−4 en la primera figura, correspondiente al sistema Sol-Júpiter,
y 7.80369 · 10−5 en la segunda, correspondiente al sistema Júpiter-Ganı́medes.

Los siguientes ejemplos, para µ = 0.002 y µ = 0.2, respectivamente, los encontramos al

azar experimentando.

(x0 , y0 ) = (−0.9, 0) (x0 , y0 ) = (0.1, 0.08)

(vx0 , vy0 ) = (0, −0.058) (vx0 , vy0 ) = (1.0, 0.8)

La curiosa estrella de siete puntas apareció mientras se estudiaban las órbitas confinadas.

14
 Entre las órbitas menos intuitivas en la variante del problema de los tres cuerpos que nos
ocupa, están las llamadas órbitas de herradura (horseshoe orbits). Nos centramos aquı́ en el
caso periódico. Una definición informal (cf. [BM05]) es que son aquellas con forma de “C”
encerrando a L3 , L4 y L5 pero no a L1 ni a L2 . Este fenómeno aparece para valores de µ
pequeños y con cJ y c1 cercanos a 3 (recordemos que c1 ≥ 3 es el valor de 2U en L1 ).
Dos ejemplos, con µ = 9.53875 · 10−4 correspondiente al sistema Sol-Júpiter, que aparecen
en [Tay81] y están representadas en [MD99, Fig. 3.17], son:

x0 = −0.97668 x0 = −1.02745
y0 = 0 y0 = 0
vx 0 = 0 vx0 = 0
vy0 = −0.06118 vy0 = 0.04032

cJ = 2.99892672 cJ = 3.00148629

En el primer caso las curvas de velocidad cero se reducen prácticamente a L4 y L5 , mientras

que en el segundo, está próxima a la curva crı́tica que pasa por L3 , lo cual impide que haya
oscilaciones como las del primer caso porque la curva de velocidad cero delimita una zona
prohibida.
Estas curiosas órbitas no son una especulación teórica, de hecho el asteroide 2002 AA29
sigue una órbita de herradura en el sistema Sol-Tierra. Lo que verı́a un observador inercial es
que m3 sigue órbitas más o menos circulares pero que la cabo de un número suficientemente
grande de órbitas, su distancia angular a m2 ha variado ostensiblemente hasta poder estar
bastante cerca sin ser capturado. En [MD81] se introducen fórmulas aproximadas cuando µ es
muy pequeño para esta distancia mı́nima y para el grosos de la “herradura”.
El trabajo [BM05] está dedicado a producir familias de órbitas de herradura. De entre
las 27 que se muestran para µ = 10−4 , destacamos las tres siguientes:

(x0 , y0 ) = (−0.93449350, 0) (x0 , y0 ) = (−1.05221706, 0) (x0 , y0 ) = (−1.06764576, 0)

(vx0 , vy0 ) = (0, −0.14032108) (vx0 , vy0 ) = (0, 0.07736088) (vx0 , vy0 ) = (0, 0.11629385)

15
Los valores de la constante de Jacobi cJ y de c1 y c2 (el valor de 2U en L2 ) están recogidos en
la tabla:
c1 3.00898924 c1 3.00898924 c1 3.00898924
cJ 2.99390329 cJ 3.00201256 cJ 2.99970858
c2 3.00885590 c2 3.00885590 c2 3.00885590
Como habı́amos anticipado, son valores próximos entre sı́ y a tres.

6. Animaciones en el sistema inercial centro de masas

Las imágenes que hemos mostrado hasta ahora han sido siempre en el sistema de referencia
rotatorio. En esta sección se explica cómo crear fotogramas en el sistema de referencia estático
e integrarlos en un fichero GIF animado que permite una visualización cómoda de diversas
situaciones relativas al caso particular del problema de los tres cuerpos que estamos estudiando.
La función que lleva a cabo el proceso es make_animation, definida como
def make animation ( Lorb , Lnum , Lplo , n f r , d e l a =100 , f l a g = F a l s e ) :
Ldat = r e t r i e v e d a t a f r a m e s ( Lorb , Lnum , n f r )
make frames ( Lorb [ 0 ] , Ldat , Lplo , f l a g )
f i l e n = o s . path . s p l i t e x t ( Lplo [ − 1 ] ) [ 0 ]
s u b p r o c e s s . c a l l ( ’./ animate ’+s t r ( d e l a )+’ ’+f i l e n , s h e l l=True )

Los tres primeros argumentos son los datos orbitales, del algoritmo numérico y de la repre-
sentación que ya han aparecido con anterioridad. Hay un argumento obligatorio más, nfr que
es el número de fotogramas. Finalmente, están los argumentos dela y f_lag, el primero es el
tiempo en milisegundos entre fotogramas y el segundo tiene la información de si se representan
los ejes y los puntos de Lagrange dependiendo respectivamente de si los bits 1 y 0 están activos
o no. Es decir, f_lag=3 pinta ambos y f_lag=2 pinta los ejes pero no los puntos de Lagrange.
Con retrieve_data_frames se llama al programa C que hace los cálculos con el método
de Runge-Kutta y del fichero de datos resultante se eligen nfr de ellos. La parte relevante del
código es:
k = 0
# retrieve only nfr values
for n f in r a n g e ( n f r ) :
ind = f l o o r ( nf ∗ nsteps / n f r )
while ( k!= i n d ) :
line = f . readline ()
k +=1
line = f . readline ()

La fórmula para ind asegura que las nfr posiciones que indican los datos recuperados estén
uniformemente distribuidos en el tiempo (aquı́ nsteps es el número de lı́neas del fichero de
datos). La salida de retrieve_data_frames es una lista de nfr pares dados por tiempos y

16
posiciones en el sistema no rotado. Esta lista la utiliza make_frames creando nfr imágenes
temporales temp_frame*.png con un punto azul en la posición de m3 tomada de Ldat y dos
puntos verdes en las posiciones de m1 y m2 . Además, make_frames crea en el fichero indicado
en Lplo una superposición de todos los fotogramas. Finalmente, se llama al script de Shell
animate que a través de la suite de código abierto ImageMagick compone el GIF animado con
el nombre tomado del último elemento de L_plo, salvo que la extensión se modifica a .gif. El
contenido del script es
# !/ bin / bash

c o n v e r t −d e l a y $1 −l o o p 1 temp frame ∗ . png −d e c o n s t r u c t −l a y e r s o p t i m i z e $2 . g i f

rm temp frame ∗ . png
du −h $2 . g i f

Los modificadores -deconstruct -layers optimize son importantes para crear ficheros de
tamaño reducido.

Como ejemplo, ejecutando

Lorb = [ 0 . 2 9 , 0 . 5 4 2 7 0 0 8 , −0.2574309 , −0.1441290 , − 0.551 2641]
Lnum = [ 0 . 2 9 5 5 , 0 . 0 0 0 1 ]
Lplo = [ − 0 . 5 7 , 0 . 7 , − 0 . 5 5 , 0 . 7 , ’ ../../ images / frames . eps ’ ]
make animation ( Lorb , Lnum , Lplo , 1 0 0 , 2 0 , 3 )

tendremos la animación de un tirón gravitatorio (gravitational slingshot) por parte de m2 que

lanza m3 hacia m1 . En http://www.uam.es/fernando.chamizo/oscuro/images/gravsl.gif
está el fichero resultante. Para mostrar algunas imágenes aquı́, podemos modificar animate de
modo que no borre los fotogramas temporales. Para que el tirón gravitatorio sea los suficien-
temente fuerte como para cambiar totalmente la trayectoria de m3 , tiene que pasar rozando a
m3 y eso es lo que vemos en el último de los siguientes cuatro fotogramas, donde m3 (azul),
casi coincide con m2 (verde).

Después m3 se dirige hacia m1 y la orbitará cierto número de veces. Mostrando tan pocos
fotogramas es difı́cil hacerse una idea del fenómeno, sobre todo teniendo en cuenta la gran

17
velocidad que adquiere m3 cuando está cerca de m2 , por ello resulta ilustrativa la última
imagen que incluye las posiciones de las tres masas a lo largo de los cien fotogramas.

0.6

0.4

0.2

-0.4 -0.2 0.2 0.4 0.6

-0.2

-0.4

La separación entre los puntos azules en esta última imagen, da una idea de la velocidad de m2
que hemos mencionado. Según se acerca al punto del tirón los puntos se van espaciando, tanto
que 100 fotogramas se vuelven pocos para ubicarlo bien. Después m3 sigue bajo la influencia
gravitatoria de m2 que se va alejando por la parte superior, lo que causa en cambio en la
curvatura de la trayectoria. En los momentos a medio camino entre la influencia de m1 y m2
la velocidad disminuye y después los puntos vuelven a espaciarse cuando m3 entra claramente
bajo la acción de m1 .

También hemos hecho un script para realizar vı́deos con el software libre ffmpeg a partir
de los fotogramas generados por make_animation.
# !/ bin / bash

# Comment ’rm temp_frame *. png ’ in animate and

# run ’ make_animation ( Lorb , Lnum , Lplo , 300 , 25 , 3) ’
rm output1 . mp4
ffmpeg −f r a m e r a t e 24 −s t a r t n u m b e r 0 − i temp frame %03d . png output1 . mp4
rm t e m p p l a y l i s t . t x t
for i i n { 1 . . 2 } ; do printf " file ’ %s ’\ n" output1 . mp4 >> t e m p p l a y l i s t . t x t ; done
rm output . mp4
ffmpeg −f c o n c a t − i t e m p p l a y l i s t . t x t −c copy output . mp4
rm output1 . mp4
echo ’ ’
du −h output . mp4

Un ejemplo correspondiente a una de las órbitas periódicas (la de forma de triángulo) de

la sección anterior se puede descargar en http://www.uam.es/fernando.chamizo/oscuro/
images/output.mp4. El 2 en el bucle for indica el número de periodos. A pesar de que
µ = 0.002 es pequeño, con un poco de atención se observan pequeños desplazamientos de
la masa central.

18
Referencias
[BM05] E. Barrabés and S. Mikkola. Families of periodic horseshoe orbits in the restricted
three-body problem. Astron. Astrophys., 432:1115–1129, mar 2005.

[CM00] A. Chenciner and R. Montgomery. A remarkable periodic solution of the three-body

problem in the case of equal masses. Ann. of Math. (2), 152(3):881–901, 2000.

[Cor98] N. J. Cornish. The Lagrange points. 1998. WMAP Education and Outreach
https://map.gsfc.nasa.gov/ContentMedia/lagrange.pdf.

[Dan88] J. M. A. Danby. Fundamentals of celestial mechanics. Willmann-Bell, Inc., Richmond,

VA, second edition, 1988.

[MD81] C. D. Murray and S. F. Dermott. The dynamics of tadpole and horseshoe orbits. I -
Theory. Icarus, 48:1–22, October 1981.

[MD99] C. D. Murray and S. F. Dermott. Solar system dynamics. Cambridge University

Press, Cambridge, 1999.

[Per17] O. Perdomo. On the period of the periodic orbits of the restricted three body problem.
Celestial Mech. Dynam. Astronom., doi 10.1007/s10569-017-9766-8, 2017.

[SB02] J. Stoer and R. Bulirsch. Introduction to numerical analysis, volume 12 of Texts in

Applied Mathematics. Springer-Verlag, New York, third edition, 2002. Translated
from the German by R. Bartels, W. Gautschi and C. Witzgall.

[Sic10] B. Sicardy. Stability of the triangular Lagrange points beyond Gascheau’s value.
Celestial Mech. Dynam. Astronom., 107(1-2):145–155, 2010.

[Tao16] M. Tao. Explicit symplectic approximation of nonseparable Hamiltonians: Algorithm

and long time performance. Phys. Rev. E (3), 94:043303, 2016.

[Tay81] D. B. Taylor. Horseshoe periodic orbits in the restricted problem of three bodies for
a sun-Jupiter mass ratio. Astron. Astrophys., 103:288–294, nov 1981.

[Ver09] C. Verzijl. Notes on the r3bp project. https://tinyurl.com/lhezete, 2009.

19