El Problema De Los Tres Cuerpos Restringido O De Euler. Leonardo L´ opez opez Hern´ andez. andez. F´ısica ısi ca Computac Comp utacion ional. al. Laboratorio De Ense˜ nanza nanza En C´ omputo omp uto En F´ısica. ısi ca.
Abstract Se estudia el problema de los tres cuerpos restringido o de Euler. Se plantean las ecuaciones del sistema utilizando utilizando la formulaci´ formulaci´ on on de Hamilton-La Hamilton-Lagrang grange, e, dichas dichas ecuaciones ecuaciones se resuelven resuelven num´ num´ericamen ericamente te por el m´ etodo etodo de Runge-Kutta de orden cuatro, encontrando con ello las trayectorias trayectorias que describe una pari´ pari´ıcula localizada localizada en alg´ un punto del espacio en presencia de dos masas relativamente grandes en comparaci´on con un la part´ıcula ıcul a anteri ant erior. or.
1.
Intr Introdu oducc cci´ i´ on on
2. 2.1. 2.1.
Uno los problemas que m´ as as ha inquietado inquiet ado a f´ısicos y matem´ aticos es el problema de los tres cuerpos. aticos Tiene su origen en la ley de gravitaci´ on on universal y la segunda ley de Newton, la primera dicta que dos cuerpos espaciados interact´ uan mediante una fuerza uan proporc proporcion ional al al product productoo de sus masas e inve inversa rsa-mente mente proporc proporcina inall al cuadra cuadrado do de sus distanc distancias ias,, la segunda estipula que, dada una fuerza sobre un cuerpo, cuerp o, ´esta esta produce produ ce acelera a celeraci´ ci´ on on sobre el mismo.
Marco arco Te´ orico orico El Probl Problema ema De Eu Euler ler..
Cons Consid id´´eres e resee el prob proble lema ma de los tres tres cuer cuerpos pos rest restri ring ngid idoo o de Eule Euler, r, esto esto es, es, dos dos ma masa sass fijas fijas espaciadas una distancia d, por comodidad colocadas sobre el eje x y con el origen en una de ellas; adem´ as, una tercera masa peque˜ as, na na en una posici´ on on arbitraria como en la figura.
Dado esto, el problema de los tres cuerpos consiste en dete determ rmin inar ar las posic posicion iones es espa espaci ciale aless de tres tres cuerpos en cualquier instante dadas sus posiciones y veloc velocid idad ades es inic inicia iale les, s, cons consid ider eran ando do que que entr entree tales tales cuerpos cuerpos solo solo hay hay inter interacc acci´ i´ on on gravit gravitaci aciona onall y adem´ as, as, que la aceleraci´on on producida por la fuerza de interacci´on on es igual a la derivada de la velocidad. La solu soluci ci´ on o´ n a este este prob proble lema ma es muy comp compli lica ca-da, da, sin sin em emba barg rgo, o, exis existe ten n plan plante team amie ien ntos tos m´ as simples, por ejemplo, el problema de Euler, el cual se basa en colocar dos masas relativamente grandes, fijas fijas en una una posi posici ci´ on o´ n y otra otra ma masa sa de ma magn gnit itud ud peque˜ na n a comp compar arad adaa a las ante anteri rior ores es con con liber liberta tad d de mo mov verse erse,, el prob problem lemaa cons consis iste te enton entonce ces, s, en determinar las posiciones y velocidades a cualquier tiempo de la part´ part´ıcula que se mueve mueve “libremente”.
Consid´ Con sid´erese eres e adem´ adem as a´s que la unica u ´nica interacci´ on on entre las tres masas es del tipo gravitatorio. Las masas on, on, as´ı que, m1 y m2 se suponen fijas en su posici´ la unic u ´ nicaa ma masa sa que que se pued puedee mo mov ver debi debido do a la presencia de las otras dos es m.
La Lass cons consid ider erac acio ione ness ante anteri rior ores es redu reduce cen n el sist sisteema de 6 ecuaciones diferenciales de segundo orden no linea lineale less a sola solame ment ntee 2, lo cual cual simp simplifi lifica ca los los c´alculos. alculos.
Es nece necesa sari rioo enco encont ntra rarr las las ecua ecuaci cion ones es de mo movi vi-miento miento en coordenadas coordenadas cartesianas cartesianas que caracteriza caracterizan n el sistema, utilizando la formulaci´ on on de HamiltonLagrange [? [?].
Figura Figura 1: Configuraci´ Configuraci´ on on de las masas.
1
Para ello n´ otese que el Lagrangiano del sistema Realizando los c´ alculos se obtiene que: esta dado por: P x = mx˙ L = T − V
(1)
(10)
P y = my˙
(11)
Mientras que el Hamiltoniano del sistema ser en coordenas cartesianas: Introduciendo lo anterior en la ecuaci´ o n (2) se escribe finalmente el Hamiltoniano del sistema como: ˙ x + yP ˙ y−L (2) H = xP P y2 P x2 Gmm1 H = + − 2m 2m x2 + y 2
La energ´ıa cin´etica del sistema es solamente la energ´ıa de la masa m, la cual esta dada por: T =
1 m(x˙ 2 + y˙ 2 ) 2
−
(xGmm d) + y −
2 2
2
(12)
(3)
N´o tese que ´este Hamiltoniano tiene la forma La energ´ıa potencial del sistema es solamente la de H = T + V , por lo cual el sistema es conservativo. energ´ıa potencial de la masa m, ´esta a su vez es la o n del Lagrangiano (eq.7), suma de la energ´ıa potencial debido a la fuerza sobre Tomando la expresi´ m por efecto de la masa m1 y la masa m2 , consid´ ere- se encuentran las ecuaciones de “Euler-Lagrange” on se primero, la contribuci´ on de m y m1 . De la figura para x y para y, utilizando la expresi´ puede verse que la energ´ıa potencial asociada al sisd ∂L ∂L ( ) (13) − tema m, m1 es: dt ∂ q ˙1 ∂q i V 1 = −
Gmm x +y 1
(4) Haciendo c´ alculos se tienen las siguientes dos ecuaciones diferenciales de orden dos, no lineales y acoDonde G es la constante de gravitaci´on universal, pladas: 3 cuyo valor es G = 6,693x10 11 kgm s2 ; Ahora, la cond2 x 1 2 x − ((x d)Gm = − (x2 Gm 2 +y 2 )3/2 (x − d) tribuci´ on al potencial del sistema m, m2 puede endt2 +y 2 )3/2 d2 y 1 2 contrarse como sigue: la masa m, esta a una distan= − (x2Gm y − ((x d)Gm 2 +y2 )3/2 y dt2 +y 2 )3/2 cia radial dada por x22 + y22 , donde x2 = x − d y (14) y2 = y, por tanto se escribe la energ´ıa potencial de este subsistema como: ´ Estas ecuaciones se resolver´an utilizando el m´etodo de Runge − Kutta de orden cuatro. Gmm2 (5) V 2 = − (x − d)2 + y 2 2
2
−
∗
−
−
La energ´ıa potencial total est´ a dada por la suma de 2.2. V 1 y V 2 , de tal forma que: V = −
Gmm x +y 1
−
(xGmm d) + y 2 2
M´ etodo de Runge-Kutta de orden cuatro. [?][?]
(6)
Uno de los m´etodos num´ ericos m´ as utilizados para la resoluci´ on de ecuaciones diferenciales es el Por tanto el Lagrangiano del sistema se escribe como: m´etodo de Runge-Kutta de orden cuatro. 2
2
−
2
Este m´ etodo consiste en aproximar la soluci´ on (7) y (t) de una ecuaci´on diferencial dada a un problema 2 − de la forma dy dt = f (t, y), para t en el intervalo [a, b], on inicial y (a) = α. Los momentos generalizados en x y en y est´ an dados sujeto a la condici´ por: En forma expl´ıcita este m´ etodo consiste en ha∂L P x = (8) cer la divisi´on del intervalo temporal y para cada ∂ x˙ punto obtenido hacer las iteraciones de acuerdo a lo ∂L (9) siguiente: P y = ˙ 1 Gmm1 + L = m(x˙ 2 + y˙ 2 )+ 2 x2 + y2
(xGmm d) + y 2 2
∂ y
2
ω = α K = hf (t , ω ) K = hf (t + , ω + K ) K = hf (t + , ω + K ) K = hf (t , ω + K ) ω = ω + (K + 2K + 2K + K )
ω (a) = α ω (a) = α ...... = ...... . . . . . . = . . . . . . . . . . . . = . . . . . . ω (a) = α
0
1
i
2
i
3
i
4
i+1
i
h 2 h 2
i+1 1 i 6
i i
i
1 2 1 2
1
(15)
2
3
1
2
3
4
dup dt
1
1
2
p
1
1
2
p
1
1
2
2
(18)
p
= hf i (t j , ω1,j , ω2,j , . . . , ω p,j ) h 1 1 2,i = hf i (t j + 2 , ω1,j + 2 K 1,1 , . . . , ω p,j + 2 K 1,p ) h 1 1 3,i = hf i (t j + 2 , ω1,j + 2 K 2,1 , . . . , ω p,j + 2 K 2,p ) 4,i = hf i (t j + h, ω1,j + K 3,1 , . . . , ω p,j + K 3,m ) ωi,j+1 = ωi,j + 16 (K 1,i + 2K 2,i + 2K 3,i + K 4,i ) (19) 1,j
Un sistema de ecuaciones de orden p, sujetas a p condiciones iniciales puede expresarse como:
du1 dt du2 dt
2,0
K K K K
M´ etodo de Runge-Kutta de orden cuatro para un sistemas de ecuaciones. [?]
= f (t, u , u , . . . , u ) = f (t, u , u , . . . , u ) .... = .... . . . . = . . . . . . . . = . . . . = f (t, u , u , . . . , u )
1
p,0
Para cada i = 0, 1, . . . , n − 1. La notaci´ on de K i se introduce para prescindir de las anidaciones sucesivas en la segunda variable de f (t, y).
2.2.1.
1,0
Para i = 1, 2, . . . , p.
3.
Desarrollo
(16) Para trabajar de forma m´as c´ omoda es necesario hacer unos ajustes a las ecuaciones a resolver (eq. 13); para ello consid´erese establecer un sistema en el cual se eviten n´ umeros muy grandes o demasiado pequeos.
p
u (a) = α u (a) = α ...... = ...... . . . . . . = . . . . . . . . . . . . = . . . . . . u (a) = α
T´ o mese tal sistema como el sistema en donde las distancias se midan en unidades astron´ omicas y 2 2 el tiempo en unidades de aos. Adem’as, es necesasio (17) llevar a cabo un reescalamiento para evitar n´ umeros demasiado pequeos como G, para ello consi´edese un cambio de variable en el lagrangiano, de tal t forma que x =Lx , y =y y t = P , donde L es el semieje p p mayor de la o´rbita de m2 a m1 (o viceversa), y P es el periodo que tarda en completar una Para t ∈ [a, b]. Dado un entero n, la partici´ on ´orbita; sup´ ongase adem´ as, que se cumple la relaci´on del intervalo [a, b] esta dada por t j = a + jh , para m2 = αm1 , entonces se tiene que: j = 1, 2, . . . , n; mientras que las aproximaciones para ui (t j ) ser´an dadas por ωi,j . El conjunto de ecuaciones (13), se resolvieron uti1
1
lizando el m´etodo de Runge-Kutta,con p = 2, haExtendiendo el m´etodo a p ecuaciones se tendr´a en- ciendo G igual a 1, con el argumento de que los tonces que: resultados ser´an semejantes, solo que mas grandes, utilizando un sencillo cambio de notaci´ on, u1 = x, u2 = y , u3 = P x y u4 = P y , con lo cual se tiene el 3
siguiente sistema de ecuaciones:
3
= um3 = um4 1 u1 = − (u2mm +u2 )( 3/2)
3
=
u˙ u˙ u˙ u˙
1 2
mm2 (u1 −d) − ((u1 −d)2 +u22 )( 3/2) 1 2 mm1 u2 mm2 u2 − 2 − (u1 +u22 )( 3/2) ((u1 −d)2 +u22 )( 3/2)
x21=h*f(t+h/2,w1+x11/2,w2+x12/2,w3+x13/2,w4+x14/2) x22=h*g(t+h/2,w1+x11/2,w2+x12/2,w3+x13/2,w4+x14/2) x23=h*u(t+h/2,w1+x11/2,w2+x12/2,w3+x13/2,w4+x14/2) x24=h*v(t+h/2,w1+x11/2,w2+x12/2,w3+x13/2,w4+x14/2) x31=h*f(t+h/2,w1+x21/2,w2+x22/2,w3+x23/2,w4+x24/2) x32=h*g(t+h/2,w1+x21/2,w2+x22/2,w3+x23/2,w4+x24/2) x33=h*u(t+h/2,w1+x21/2,w2+x22/2,w3+x23/2,w4+x24/2) x34=h*v(t+h/2,w1+x21/2,w2+x22/2,w3+x23/2,w4+x24/2)
(20)
x41=h*f(t+h,w1+x31,w2+x32,w3+x33,w4+x34) x42=h*g(t+h,w1+x31,w2+x32,w3+x33,w4+x34) x43=h*u(t+h,w1+x31,w2+x32,w3+x33,w4+x34) x44=h*v(t+h,w1+x31,w2+x32,w3+x33,w4+x34)
Se implemento el programa Euler.f90 en lenguaje fortran90, el c´ odigo es el siguiente:
w1=w1+(x11+2*x21+2*x31+x41)/6 w2=w2+(x12+2*x22+2*x32+x42)/6 w3=w3+(x13+2*x23+2*x33+x43)/6 w4=w4+(x14+2*x24+2*x34+x44)/6
program Euler !este programa resuelve el sistema de cuatro ecuaciones diferenciales !del problema de los tres cuerpos de Euler, no lineales y acopladas !sujeto a las condiciones iniciales y10 y y20 para las posiciones, !y30 y y04 para las velocidades.
!salida en pantalla de los resultados print*, ’ La solucion numerica para x, y, xprima y yprima, es: ’ print*, w1, w2, w3, w4 write(1,*) t, w1 write(2,*) t, w2 write(3,*) t, w3 write(4,*) t, w4 write(5,*) t, w1, w2 write(6,*) w1, w2
implicit none !definicion de variables reales real (8) :: x11, x12, x13, x14, x21, x22, x23, x24, x31, x32, x33, x34 real (8) :: x41, x42, x43, x44, f, g, u, v, t, a, b, h, d, Px, Py real (8) :: w1, w2, w3, w4, y10, y20, y30, y40, m1, m2, alpha real, parameter :: pi=3.14159 integer :: i, n
end do print*
!descripcion del programa print*, ’————————————————————’ print*, ’ El Problema De Los Tres Cuerpos De Euler ’ print*, ’————————————————————’ print* print*, ’Este programa resuelve el problema de los’ print*, ’tres cuerpos de Euler, dada una de las’ print*, ’masas como multiplo de la otra m2=alpha*m1’ print*, ’dadas las posiciones y velocidades iniciales’ print*, ’de la masa peque˜ na; utilizando el metodo’ print*, ’de Runge-Kutta.’
!salida de datos en graficas de gnuplot call system (’gnuplot trescuerpos1 .gp ) callsystem( gnuplottrescuerpos 2 .gp ) callsystem( gnuplottrescuerpos 3 .gp ) callsystem( gnuplottrescuerpos 4 .gp ) callsystem( gnuplottrescuerpos 5 .gp ) callsystem( gnuplottrescuerpos 6 .gp )
print ∗ end program Euler
!introduccion de los valores a usar print*, ’ Da los valores de las masas m1 y m2 ’ read*, m1, m2 print*, ’ Da el intervalo de tiempo ’ read*, a, b print*, ’ Da el numero de divisiones de tu intervalo ’ read*, n print*, ’ Da las posiciones iniciales de la masa m ’ read*, y10, y20 print*, ’Da las velocidades iniciales de la masa m’ read*, y30, y40 print*, ’ Da la separacion entre m1 y m2 ’ read*, d
real(8)functionf (t , x , y , P x , P y) real(8) :: t , x , y , P x , P y f = P x end function f real(8)functiong(t , x , y , P x , P y) real(8) :: t , x , y , P x , P y g = Py end function g real(8)functionu(t , x , y , P x , P y) real(8) :: t , x , y , P x , P y u=-((4*pi*pi*x/(x**2+y**2)**(3/2))-(4*pi*pi*alpha*(x-d)/(((xd)**2+y**2)**(3/2)))) end function u
!apertura de ficheros open(unit=1, file=’x.dat’, status=’replace’) open(unit=2, file=’xprima.dat’, status=’replace’) open(unit=3, file=’y.dat’, status=’replace’) open(unit=4, file=’yprima.dat’, status=’replace’) open(unit=5, file=’txy.dat’, status=’replace’) open(unit=6, file=’xy.dat’, status=’replace’)
real(8)functionv(t , x , y , P x , P y) real(8) :: t , x , y , P x , P y v=-(4*pi*pi*y/((x**2+y**2)**(3/2)))-(4*pi*pi*alpha*y/(((xd)**2+y**2)**(3/2))) end function v
!definicion de algunos valores alpha=m1/m2 h=(b-a)/n t=a
Las ecuaciones
w1=y10 w2=y20 w3=y30 w4=y40
4.
Resultados.
!inicio del ciclo del algoritmo de Runge-Kutta do i= 1, n
5.
Conclusiones.
t=a+(i-1)*h
Referencias
x11=h*f(t,w1,w2,w3,w4) x12=h*g(t,w1,w2,w3,w4) x13=h*u(t,w1,w2,w3,w4) x14=h*v(t,w1,w2,w3,w4)
[1] Classical Dynamics Of Particles And Systems, Stephen T. Thornton, Jerry B. Marion, Fith 4
Edition, Thomson Brooks/Cole 2004. [2] An´ alisis Num´ erico, Richard L. Burden, J. Douglas Faires, 7ma. Edici´ on, Editorial Thomson Learning, M´exico D.F. 2004. [3] An Introduction To Computational Physics, Tao Pang, Second Edition, Editorial Cambrige, New York 2006.
5
