Ejercicio: 2Eva2026PAOI_T2 EDO Problema de dos cuerpos
literal a
Simplificar las expresiones:
G = 4.982 x 10-19
m1 = 5.9722 x 1024
m2 = 7.3490 x 1022
Se empieza reemplazando las constantes para simplificar la ecuación..
\frac{d^2 x}{dt^2}=\frac{4.982 ( 10^{-19}) 5.9722 (10^{24})}{\left( \sqrt{x^2+y^2}\right)^3}\left(-x\right) \frac{d^2 y}{dt^2}=\frac{4.982 (10^{-19}) 5.9722 x (10^{24})}{\left( \sqrt{x^2+y^2}\right)^3}\left(-y\right)quedando de la siguiente forma:
\frac{d^2 x}{dt^2}=-2.9754(10^6)\frac{x}{\left( \sqrt{x^2+y^2}\right)^3} \frac{d^2 y}{dt^2}=-2.9754(10^6) \frac{y}{\left( \sqrt{x^2+y^2}\right)^3}
Condiciones iniciales cuando t0 =0
x0 = 384.4
y0 = 0
vx0 = 0
vy0 = 88.128
Genera la tabla a desarrollar:
| ti | xi | yi | vxi | vyi |
|---|---|---|---|---|
| 0 | 384.4 | 0 | 0 | 88.128 |
literal b
El sistema de ecuaciones es de 2da derivada, por lo que se aplica el método de Runge-Kutta de 2do orden a cada expresión.
f_x(t,x,y,vx,vy) = v_x
f_y(t,x,y,vx,vy)= v_y
g_x(t,x,y,vx,vy)=-2.9754(10^6)\frac{x}{\left( \sqrt{x^2+y^2}\right)^3}
g_y(t,x,y,vx,vy) = -2.9754(10^6) \frac{y}{\left( \sqrt{x^2+y^2}\right)^3}
Los pasos del método se adaptan al ejercicio:
K1_{x} = h (vx)
K1_{y} = h (vy)
K1_{vx} = h \left( -2.9754(10^6)\frac{x}{\left( \sqrt{x^2+y^2}\right)^3} \right)
K1_{vy} = h \left( -2.9754(10^6) \frac{y}{\left( \sqrt{x^2+y^2}\right)^3} \right)
K2_{x} = h (vx+K1_{vx})
K2_{y} = h (vy+K1_{vy})
K2_{vx} = h \left( -2.9754(10^6)\frac{x+K1_x}{\left( \sqrt{(x+K1_x)^2+(y+K1_y)^2}\right)^3} \right)
K2_{vy} = h \left( -2.9754(10^6) \frac{y+K1_y}{\left( \sqrt{(x+K1_x)^2+(y+K1_y)^2}\right)^3} \right)
x_{i+1}=x_i+\frac{K1_{x}+K2_{x}}{2}
y_{i+1}=y_i+\frac{K1_{y}+K2_{y}}{2}
vx_{i+1}=vx_i+\frac{K1_{vx}+K2_{vx}}{2}
vy_{i+1}=vy_i+\frac{K1_{vy}+K2_{vy}}{2}
literal c
| ti | xi | yi | vxi | vyi |
|---|---|---|---|---|
| 0 | 384.4 | 0 | 0 | 88.128 |
itera =0
K1_{x} = 0.1 (0) = 0
K1_{y} = 0.1 (88.128) =8.8128
K1_{vx} = 0.1 \left( -2.9754(10^6)\frac{384.4}{\left( \sqrt{384.4^2+0^2}\right)^3} \right) =-2.0136
K1_{vy} = 0.1 \left( -2.9754(10^6) \frac{0}{\left( \sqrt{384.4^2+0^2}\right)^3} \right) = 0
K2_{x} = 0.1 (-2.0136) = -0.2014
K2_{y} = 0.1 (88.128+0) = 8.8128
K2_{vx} = 0.1 \left( -2.9754(10^6)\frac{384.4+0}{\left( \sqrt{(384.4+0)^2+(0+8.8128)^2}\right)^3} \right) = -2.0121
K2_{vy} = 0.1 \left( -2.9754(10^6) \frac{0+8.8128}{\left( \sqrt{(384.4+0)^2+(0+8.8128)^2}\right)^3} \right) =-0.0461
actualiza las variables de la iteración:
x_{1}=384.4+\frac{0-0.2014}{2} =384.30
y_{1}=0+\frac{8.8128+8.8128}{2} = 8.8128
vx_{1}=0+\frac{-2.0136-2.0121}{2} = -2.0128
vy_{1}=88.128+\frac{0+-0.0461}{2} = 88.105
| ti | xi | yi | vxi | vyi |
|---|---|---|---|---|
| 0 | 384.4 | 0 | 0 | 88.128 |
| 0.1 | 384.3 | 8.8128 | -2.0128 | 88.105 |
itera=1
K1_{x} = 0.1 (-2.0128) = -0.2013
K1_{y} = 0.1 (88.105) = 8.8105
K1_{vx} = 0.1 \left( -2.9754(10^6)\frac{384.3}{\left( \sqrt{384.3^2+8.8128^2}\right)^3} \right) =-2.0131
K1_{vy} = 0.1 \left( -2.9754(10^6) \frac{8.8128}{\left( \sqrt{384.3^2+8.8128^2}\right)^3} \right) = -0.0462
K2_{x} = 0.1 (-2.0128-2.0131) = -0.4026
K2_{y} = 0.1 (88.105-0.0462) = 8.8059
K2_{vx} = 0.1 \left( -2.9754(10^6)\frac{384.3-0.2013}{\left( \sqrt{(384.3-0.2013)^2+(8.8128+8.8105)^2}\right)^3} \right) = -2.0105
K2_{vy} = 0.1 \left( -2.9754(10^6) \frac{8.8128+8.8105}{\left( \sqrt{(384.3-0.2013)^2+(8.8128+8.8105)^2}\right)^3} \right) = -0.0922
x_{2}=384.3+\frac{-0.2013-0.4026}{2} = 384.0
y_{2}=8.8128+\frac{8.8105+8.8059}{2}= 17.621
vx_{2}=-2.0128+\frac{-2.0131-2.0105}{2} = -4.0246
vy_{2}=88.105+\frac{-0.0462-0.0922}{2}= 88.036
| ti | xi | yi | vxi | vyi |
|---|---|---|---|---|
| 0 | 384.4 | 0 | 0 | 88.128 |
| 0.1 | 384.3 | 8.8128 | -2.0128 | 88.105 |
| 0.2 | 384.0 | 17.621 | -4.0246 | 88.036 |
literal d
Algoritmo en Python
# 2Eva2026PAOI_T2 EDO Problema de dos cuerpos
# Trayectoria dos cuerpos en espacio
import numpy as np
G = 6.6740e-11*((60*60*24)**2/(1000000**3)) # constante gravitacion
m1 = 5.9722e24 # masa tierra
m2 = 7.3490e22 # masa luna Cuerpo2
r = lambda x,y: np.sqrt(x**2 + y**2)
# ecuacion cuerpo 2
fx = lambda t,x,y,vx,vy: vx
gx = lambda t,x,y,vx,vy: G*m1*(0-x)/(r(x,y))**3
fy = lambda t,x,y,vx,vy: vy
gy = lambda t,x,y,vx,vy: G*m1*(0-y)/(r(x,y))**3
# condiciones iniciales
x0 = 3.84400e8/1000000 # coordenadas Cuerpo2
y0 = 0
vx0 = 0 # velocidad Cuerpo2
vy0 = 1.02e3*(60*60*24)/1000000
t0 = 0
h = 0.1
muestras = 250+1
# Algoritmo como función
def rungekutta2_fg(fx,fy,gx,gy,t0,x0,y0,vx0,vy0,h,muestras):
''' solucion a EDO d2y/dx2 con Runge-Kutta 2do Orden,
adaptado al ejercicio
'''
tamano = muestras + 1
tabla = np.zeros(shape=(tamano,5+8),dtype=float)
tabla[0,:5] = [t0,x0,y0,vx0,vy0]
ti = t0 # valores iniciales
xi = x0
yi = y0
vxi = vx0
vyi = vy0
for i in range(1,tamano,1):
K1x = h * fx(ti,xi,yi,vxi,vyi)
K1y = h * fy(ti,xi,yi,vxi,vyi)
K1vx = h * gx(ti,xi,yi,vxi,vyi)
K1vy = h * gy(ti,xi,yi,vxi,vyi)
K2x = h * fx(ti+h,xi+K1x,yi+K1y,vxi+K1vx,vyi+K1vy)
K2y = h * fy(ti+h,xi+K1x,yi+K1y,vxi+K1vx,vyi+K1vy)
K2vx = h * gx(ti+h,xi+K1x,yi+K1y,vxi+K1vx,vyi+K1vy)
K2vy = h * gy(ti+h,xi+K1x,yi+K1y,vxi+K1vx,vyi+K1vy)
xi = xi + (K1x+K2x)/2
yi = yi + (K1y+K2y)/2
vxi = vxi + (K1vx+K2vx)/2
vyi = vyi + (K1vy+K2vy)/2
ti = ti + h
tabla[i] = [ti,xi,yi,vxi,vyi,
K1x, K1y, K1vx, K1vy,
K2x, K2y, K2vx, K2vy ]
return(tabla)
# PROCEDIMIENTO
tabla = rungekutta2_fg(fx,fy,gx,gy,
t0,x0,y0,vx0,vy0,h,muestras)
n = len(tabla)
# SALIDA
print('Trayectoria de dos cuerpos')
np.set_printoptions(precision=4)
print('EDO f,g con Runge-Kutta 2 Orden')
print('i ','[ ti, xi, yi, vxi, vyi',']')
print(' [ K1x, K1y, K1vx, K1vy ]')
print(' [ K2x, K2y, K2vx, K2vy ]')
tamano = 5
for i in range(0,tamano,1):
txt = ' '
if i>=10:
txt = ' '
print(str(i),tabla[i,:5])
print(txt,tabla[i,5:9])
print(txt,tabla[i,9:])
# GRAFICA
import matplotlib.pyplot as plt
ti = tabla[:,0]
xi = tabla[:,1]
yi = tabla[:,2]
plt.plot(xi,yi,color='blue', label='Trayectoria')
plt.plot(0,0,'o',color='blue',label='Cuerpo1')
plt.plot(xi[0],yi[0],'*',color='red',label='Cuerpo2[0]')
plt.plot(xi[-1],yi[-1],'o',color='red',label='Cuerpo2[n]')
# entorno de gráfica
plt.axhline(0,color='gray',linestyle='dashed')
plt.axvline(0,color='gray',linestyle='dashed')
plt.xlabel('x [1000Km]')
plt.ylabel('y [1000Km]')
plt.title('Trayectoria 2 Cuerpos. h='+str(h)+', muestras='+str(muestras)+', vy0:'+str(vy0))
plt.legend(loc='lower left')
plt.tight_layout()
plt.show()
resultados.txt
Trayectoria de dos cuerpos
EDO f,g con Runge-Kutta 2 Orden
i [ ti, xi, yi, vxi, vyi ]
[ K1x, K1y, K1vx, K1vy ]
[ K2x, K2y, K2vx, K2vy ]
0 [ 0. 384.4 0. 0. 88.128]
[0. 0. 0. 0.]
[0. 0. 0. 0.]
1 [ 1.0000e-01 3.8430e+02 8.8128e+00 -2.0128e+00 8.8105e+01]
[ 0. 8.8128 -2.0136 0. ]
[-0.2014 8.8128 -2.0121 -0.0461]
2 [ 2.0000e-01 3.8400e+02 1.7621e+01 -4.0246e+00 8.8036e+01]
[-0.2013 8.8105 -2.0131 -0.0462]
[-0.4026 8.8059 -2.0105 -0.0922]
3 [ 3.0000e-01 3.8349e+02 2.6420e+01 -6.0343e+00 8.7920e+01]
[-0.4025 8.8036 -2.0115 -0.0923]
[-0.6036 8.7943 -2.0078 -0.1383]
4 [ 0.4 382.7905 35.2051 -8.0407 87.7591]
[-0.6034 8.792 -2.0088 -0.1384]
[-0.8043 8.7782 -2.0041 -0.1843]
literal e
El resultado con el algoritmo se aproxima a la orbita de un cuerpo mas pequeño orbitando sobre uno mucho mas grande. Semejante a lo sugerido en el enunciado sobre el movimiento de la luna alrededor de la tierra. La simplificación no considera el movimiento de la tierra que se observa por ejemplo en las mareas, que es un modelo mas complejo.
literal f
Solo para comprobar que el algoritmo presentado considera la condición de iniciar con una velocidad tangencial menor, hace que la "luna" pierda su órbita y comience a alejarse de la tierra.
El resultado gráfico con el algoritmo es:
