IMT1001-2023: Introduccion a la Ingenieria Matematica
Modulo Analisis Numerico
Prof. Manuel A. Sánchez
In [247]:
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import pandas as pd
from scipy.optimize import newton
from IPython.display import display, HTML
display(HTML("""<style>.output {display: flex;align-items: center;text-align: center;}</style>"""))
Metodo de Euler explicito¶
In [248]:
def Forward_Euler(odefun, t_span, y0, h=None, return_trajectory=False):
'''
Metodo de Euler explicito
Resuelve: dy/dy = odefun(t,y), y(t0) = y0
Input : odefun, t_span = (tiemp inicial, tiempo final)
y0 condicion inicial, h paso en t,
return_trajectory returnar o no toda la trajectoria de la solucion
Output : (t,y)
'''
t0, tf = t_span
t = t0; y = y0
if h is None: h = (tf-t0)/100
yh = []; tn = []
if return_trajectory: yh.append(y0); tn.append(0.0)
while t<tf:
y = ForwardEuler_singlestep(odefun,h,t,y)
t+=h
if return_trajectory: yh.append(y); tn.append(t)
if return_trajectory:
return (tn,np.array(yh))
else:
return t,y
def ForwardEuler_singlestep(fun,dt,t0,y0):
yout = y0 + dt*fun(t0, y0)
return yout
Ejemplo 1:¶
Considere el problema de Cauchy: $$ \frac{dy(t)}{dt} = \cos(2 y(t)), \quad t\in[0,1],\quad \mbox{condicion inicial}\,\,y(0) = 0 $$
In [7]:
fun_ex1 = lambda t, y: np.cos(2*y)
funprime_ex1 = lambda t, y: -2.0*np.sin(2*y)
y_ex1 = lambda t: 0.5*np.arcsin((np.exp(4*t)-1)/(np.exp(4*t)+1))
t_ex1, yh_ex1 = Forward_Euler(odefun=fun_ex1,
t_span=(0,1),
y0=0,
h=2**(-4),
return_trajectory=True)
_,_ = plot_trayectory(t_ex1, yh_ex1, exact=y_ex1)
plt.show()
Ejemplo 2: estabilidad¶
Considere el problema $$ y' = \lambda y,\quad t\in(0,\infty),\quad\text{y} \quad y(0) = y_0 $$
In [63]:
lamb = -1.0; fun_ex2 = lambda t,y: lamb*y
y0 = np.array([1.0])
y_ex2 = lambda t: y0*np.exp(lamb*t)
t_ex2, yh_ex2 = Forward_Euler(odefun=fun_ex2,
t_span=(0,30),
y0=y0,
h=1.99,#h=2.01,1.99
return_trajectory=True)
_,_ = plot_trayectory(t_ex2, yh_ex2, exact=y_ex2)
plt.show()
In [65]:
lamb = -1.0; fun_ex2 = lambda t,y: lamb*y
y0 = np.array([1.0])
y_ex2 = lambda t: y0*np.exp(lamb*t)
t_ex2i, yh_ex2i = Backward_Euler(odefun=fun_ex2,
t_span=(0,30),
y0=y0,
h=3.0,#h=2.01,1.99
return_trajectory=True)
_,_ = plot_trayectory(t_ex2i, yh_ex2i, exact=y_ex2)
plt.show()
Ejemplo 3: estabilidad¶
$$ y' = \arctan(3y) - 3y + t, \quad t>0, \qquad y(0) = 1 $$Como $\partial f/ \partial y = 3/(1+9y^{2}) -3<0$, $\lambda_{\max} = \max|f_{y}| = 3$. As\'i para Euler explicito podemos escoger $h<2/3$
In [68]:
fun_ex3 = lambda t,y: 1.0-y**2# np.arctan(3*y)-3*y+t
funprime_ex2 = lambda t, y: -2*y
y0 = np.array([(np.exp(1.0)-1)/(np.exp(1.0)+1)])
y_ex3 = lambda t: (np.exp(2*t+1)-1.0)/(np.exp(2*t+1)+1)
t_ex3, yh_ex3 = Forward_Euler(odefun=fun_ex3,
t_span=(0,300),
y0=y0,
#h=30.0/14.0,
#h=30.0/16.0,
h=0.5,
return_trajectory=True)
In [69]:
_,_ = plot_trayectory(t_ex3, yh_ex3, exact=y_ex3)
plt.show()
Metodo de Euler implicito¶
In [257]:
def Backward_Euler(odefun, t_span, y0, odefunprime=None, h=None, return_trajectory=False):
'''
Metodo de Euler implcito
Resuelve: dy/dy = odefun(t,y), y(t0) = y0
Input : odefun, t_span = (tiemp inicial, tiempo final)
y0 condicion inicial, h paso en t,
return_trajectory returnar o no toda la trajectoria de la solucion
Output : (t,y)
'''
t0, tf = t_span
t = t0; y = y0
if h is None: h = (tf-t0)/100
yh = []; tn = []
if return_trajectory: yh.append(y0); tn.append(0.0)
while t<tf:
y = BackwardEuler_singlestep(odefun,h,t,y, dfun=odefunprime)
t+=h
if return_trajectory: yh.append(y); tn.append(t)
if return_trajectory:
return (tn,np.array(yh))
else:
return (t,np.array(y))
def BackwardEuler_singlestep(fun,dt,t0,y0,dfun=None):
def funnonlinear(y):
return y - dt*fun(t0,y) - y0
if dfun is not None:
def funprimenonlinear(y):
return 1 - dt*dfun(t0,y)*fun(t0,y)
yout = newton(funnonlinear, y0, fprime=funprimenonlinear)
else:
yout = newton(funnonlinear, y0)
return yout
In [32]:
t_ex12, yh_ex12 = Backward_Euler(odefun=fun_ex1,
t_span=(0,1),
y0=0,
h=2**(-4),
odefunprime=funprime_ex1,
return_trajectory=True)
_,_ = plot_trayectory(t_ex12, yh_ex12,exact=y_ex1)
plt.show()
Metodo de Heun o Euler modificado¶
In [78]:
def Heun(odefun, t_span, y0, h=None, return_trajectory=False):
'''
MEtodo de Heum o De Euler mejorado o Euler modificado
Resuelve: dy/dy = odefun(t,y), y(t0) = y0
Input : odefun, t_span = (tiemp inicial, tiempo final)
y0 condicion inicial, h paso en t,
return_trajectory returnar o no toda la trajectoria de la solucion
Output : (t,y)
'''
t0, tf = t_span
t = t0; y = y0
if h is None: h = (tf-t0)/100
yh = []; tn = []
if return_trajectory: yh.append(y0); tn.append(0.0)
while t<tf:
y = Heun_singlestep(odefun,h,t,y)
t+=h
if return_trajectory: yh.append(y); tn.append(t)
if return_trajectory:
return (tn,yh)
else:
return (t,y)
def Heun_singlestep(fun,dt,t0,y0):
f0 = fun(t0,y0)
z = y0 + dt * f0
yout = y0 + 0.5*dt*(f0 + fun(t0+dt, z))
return yout
Estudio de convergencia de metodos de Euler¶
In [55]:
t0 = 0.0; tf = 1.0
# number of subintervals
NN = [2,4,8,16,32,64,128, 256, 512, 1024]
h = [(tf-t0)/N for N in NN]
Error_fe = []
Error_be = []
Error_me = []
print(y_ex1(tf))
for hj in h:
# Solution forward Euler
t_fe, yh_fe = Forward_Euler(odefun=fun_ex1,
t_span=(t0,tf),
y0=0,
h=hj)
# Solution backward Euler
t_be, yh_be = Backward_Euler(odefun=fun_ex1,
t_span=(t0,tf),
y0=0,
h=hj)
# Solution Modified Euler
t_me, yh_me = Heun(odefun=fun_ex1,
t_span=(t0,tf),
y0=0,
h=hj)
# Compute Error
Error_fe.append(abs(y_ex1(tf)-yh_fe))
Error_be.append(abs(y_ex1(tf)-yh_be))
Error_me.append(abs(y_ex1(tf)-yh_me))
# rates of convergence
r_fe = np.zeros(len(h))
r_be = np.zeros(len(h))
r_me = np.zeros(len(h))
for j in range(len(h)-1):
r_fe[j+1] = np.log(Error_fe[j+1]/Error_fe[j])/np.log(h[j+1]/h[j])
r_be[j+1] = np.log(Error_be[j+1]/Error_be[j])/np.log(h[j+1]/h[j])
r_me[j+1] = np.log(Error_me[j+1]/Error_me[j])/np.log(h[j+1]/h[j])
d = {'h': h, 'F. Euler Error': Error_fe, 'p FE': r_fe,
'B. Euler Error': Error_be, 'p BE': r_be,
'M. Euler Error': Error_me, 'p ME': r_me, }
df = pd.DataFrame(data=d)
In [56]:
df.style.format({'h':'{:.2e}',
'F. Euler Error':'{:.2e}', 'p FE':'{:.2f}',
'B. Euler Error':'{:.2e}', 'p BE':'{:.2f}',
'M. Euler Error':'{:.2e}', 'p ME':'{:.2f}'}
).set_caption("Convergence Table for first order Euler methods")
Out[56]:
In [58]:
fig, ax = plt.subplots(1,1,figsize=(8,6))
ax.loglog(d['h'], d['F. Euler Error'], ':o', label='Forward Euler')
ax.loglog(d['h'], d['B. Euler Error'], ':*', label='Backward Euler')
ax.loglog(d['h'], d['M. Euler Error'], ':s', label='Modified Euler')
ax.legend(); plt.show()
In [19]:
t_cn, yh_cn = Crank_Nicolson(odefun=fun_ex1,
t_span=(0,1),
y0=np.array([0.0]),
h=2**(-4),
odefunprime=funprime_ex1,
return_trajectory=True)
_,_ = plot_trayectory(t_cn, yh_cn,exact=y_ex1)
plt.show()
In [21]:
t_rk, yh_rk = RK3(odefun=fun_ex1,
t_span=(0,1),
y0=np.array([0.0]),
h=2**(-4),
return_trajectory=True)
_,_ = plot_trayectory(t_rk, yh_rk, exact=y_ex1)
plt.show()
Test de convergencia de metodos Runge Kutta explicitos¶
In [52]:
df.style.format({'h':'{:.2e}',
'RK 2 Error':'{:.2e}', 'p RK2':'{:.2f}',
'RK 3 Error':'{:.2e}', 'p RK3':'{:.2f}',
'RK 4 Error':'{:.2e}', 'p RK4':'{:.2f}'}
).set_caption("Tabla de convergencia metodos de Runge-Kutta explicitos")
Out[52]:
In [53]:
fig, ax = plt.subplots(1,1,figsize=(8,6))
ax.loglog(d['h'], d['RK 2 Error'], ':o', label='RK2')
ax.loglog(d['h'], d['RK 3 Error'], ':*', label='RK3')
ax.loglog(d['h'], d['RK 4 Error'], ':s', label='RK4')
ax.legend()
plt.show()
In [26]:
def lorenz(t,y):
sigma = 10.0
beta = 8.0/3.0
rho = 28.0
f = [sigma*(y[1]-y[0]),
y[0]*(rho-y[2])-y[1],
y[0]*y[1] - beta*y[2]]
return np.array(f)
In [27]:
t_rk, yh_rk = RK4(odefun=lorenz,
t_span=(0,20),
y0=np.array([-8.0,8.0,27.0]),
h=0.01,
return_trajectory=True)
ax = plt.figure(figsize=(8,6)).add_subplot(projection='3d')
ax.plot(yh_rk[:,0], yh_rk[:,1], yh_rk[:,2], 'C04')
# plt.axis('off')
plt.show()
Ejemplo: Modelo SIR para epidemias¶
In [81]:
# !pip install graphviz
In [80]:
# import graphviz
# from graphviz import Digraph
# g = Digraph('G', engine="neato", filename='ex.gv',format='pdf')
# g.attr(size='7')
# g.node('S(t)',pos='1,1!', style='filled', fillcolor='#4285F4')
# g.node('I(t)',pos='3,1!', style='filled', fillcolor='#DB4437')
# g.node('R(t)',pos='5,1!', style='filled', fillcolor='#0F9D58')
# # yellow = #F4B400
# g.edge('S(t)','I(t)', label = '<β>')
# g.edge('I(t)','R(t)', label = '<γ>')
# g.render()
# g
In [249]:
def SIR_ode(t, y):
# SIR Model
# dS/dt = f1; dS/dt = f2; dS/dt = f3
# variables
S = y[0]; I = y[1]; R = y[2]
# parameters
R_0 = 1.64
gamma = 1.0/11.0
beta = R_0*gamma/S0
# right-hand side
f1 = -beta*S*I
f2 = beta*S*I -gamma*I
f3 = gamma*I
return np.array([f1, f2, f3])
# Initial conditions
S0 = 254.0; I0 = 7.0; R0 = 0.0
In [250]:
tE, yhE = Forward_Euler(odefun=SIR_ode,
t_span=(0,100),
y0=np.array([S0, I0, R0]),
h=0.01,
return_trajectory=True)
t, yh = RK4(odefun=SIR_ode,
t_span=(0,100),
y0=np.array([S0, I0, R0]),
h=0.01,
return_trajectory=True)
In [251]:
plt.figure(figsize=(8,6))
plt.plot(tE, yhE[:,0], 'b', label = r'$S(t)$ FE')
plt.plot(tE, yhE[:,1], 'r', label = r'$I(t)$ FE')
plt.plot(tE, yhE[:,2], 'g', label = r'$R(t)$ FE')
plt.legend(loc='best')
plt.xlabel('t')
plt.grid()
plt.show()
In [252]:
plt.figure(figsize=(8,6))
plt.plot(t, yh[:, 0], 'b', label=r'$S(t)$ RK4')
plt.plot(t, yh[:, 1], 'r', label=r'$I(t)$ RK4')
plt.plot(t, yh[:, 2], 'g', label=r'$R(t)$ RK4')
plt.legend(loc='best')
plt.xlabel('t')
plt.grid()
plt.show()
Sistemas Hamiltonianos¶
Ejemplo: Pendulo simple¶
\begin{eqnarray} \frac{d p}{dt} & = &- (\frac{g}{ L}) \sin(q) \\ \frac{d q}{dt} & = & p \end{eqnarray}Hamiltonian function: $H(p,q) = \frac{1}{2} p^{2} - (\frac{g}{L})\cos (q)$. Note that \begin{eqnarray} \frac{\partial H}{\partial q} & = & (\frac{g}{L})\sin(q) \\ \frac{\partial H}{\partial p} & = & p \end{eqnarray}
Symplectic Euler method: \begin{eqnarray} p_h^{n+1} & = & p_h^{n} + h(- (\frac{g}{L})\sin(q_h^{n})) \\ q_h^{n+1} & = & q_h^{n} + h( p_h^{n+1}) \end{eqnarray}
In [348]:
def Symplectic_Euler(odefun, t_span, y0, odefunprime=None, h=None, return_trajectory=False):
'''
Metodo de Euler simplectico
Resuelve: dy/dy = odefun(t,y), y(t0) = y0
Input : odefun, t_span = (tiemp inicial, tiempo final)
y0 condicion inicial, h paso en t,
return_trajectory returnar o no toda la trajectoria de la solucion
Output : (t,y)
'''
t0, tf = t_span
t = t0; y = y0
if h is None: h = (tf-t0)/100
yh = []; tn = []
if return_trajectory: yh.append(y0); tn.append(0.0)
while t<tf:
y = SymplecticEuler_singlestep(odefun,h,t,y, separable=True, dfun=odefunprime)
t+=h
if return_trajectory: yh.append(y); tn.append(t)
if return_trajectory:
return (tn,np.array(yh))
else:
return (t,y)
def SymplecticEuler_singlestep(fun_Ham,dt,t0,y0,separable=True,dfun=None):
Hp, Hq = fun_Ham
if len(y0)%2 == 0:
n = int(len(y0)/2)
p0 = y0[:n]
q0 = y0[n:]
if separable: # the method is explicit!
p = p0 - dt*Hq(t0, p0, q0)
q = q0 + dt*Hp(t0, p ,q0)
else:
def funnonlinear(p):
return p + dt*Hq(t0,p,q0) - p0
p = newton(funnonlinear, p0)
q = q0 + dt*Hp(t0+dt, p ,q0)
return np.concatenate((p,q))
else:
raise NameError('dimension is not even')
In [349]:
def Math_pendulum(t, y):
# parameters
global g, L
g = 9.8; L = 1
p,q = y
f1 = -(g/ L)*np.sin(q)
f2 = p
return np.array([f1,f2])
def Math_pendulum_Hq(t, p, q):
return (g/ L)*np.sin(q)
def Math_pendulum_Hp(t, p, q):
return p
def Energy_pendulum(p,q):
return 0.5*p**2-(g/ L)*np.cos(q)
# Condiciones iniciales
p0 = 0.0; q0 = np.pi/3
In [350]:
h = 0.01; T = 10
# Forward Euler solution
t_fe, y_fe = Forward_Euler( odefun=Math_pendulum,
t_span=(0,T),
y0=np.array([p0,q0]),
h=h,
return_trajectory=True)
p_fe = y_fe[:,0]; q_fe = y_fe[:,1]
# # Backward Euler solution
t_be, y_be = Backward_Euler( odefun=Math_pendulum,
t_span=(0,T),
y0=np.array([p0,q0]),
h=h,
return_trajectory=True)
p_be = y_be[:,0]; q_be = y_be[:,1]
# # Symplectic Euler solution
t_se, y_se = Symplectic_Euler( odefun=(Math_pendulum_Hp, Math_pendulum_Hq),
t_span=(0,T),
y0=np.array([p0,q0]),
h=h,
return_trajectory=True)
p_se = y_se[:,0]; q_se = y_se[:,1]
In [351]:
plt.figure(figsize=(8,6))
plt.plot(q_fe[:500],p_fe[:500],'oC0')
plt.plot(q_be[:500],p_be[:500],'oC1')
plt.plot(q_se[:500],p_se[:500],'oC2')
plt.legend(['explicit Euler','implicit Euler', 'symplectic Euler'])
plt.show()
In [265]:
H_fe = Energy_pendulum(p_fe, q_fe)
H_be = Energy_pendulum(p_be, q_be)
H_se = Energy_pendulum(p_se, q_se)
plt.figure(figsize=(8,6))
plt.plot(t_fe, H_fe,'--C0')
plt.plot(t_be, H_be,'--C1')
plt.plot(t_se, H_se,'--C2')
plt.legend(['explicit Euler','implicit Euler', 'symplectic Euler'])
plt.show()
Simulacion del pendulo¶
In [267]:
#ani = animation.FuncAnimation(fig, animate, frames=10, interval=interval, blit=True, init_func=init);
ani = animation.FuncAnimation(fig, animate, frames=100, interval=interval, blit=False, init_func=init);
# equivalent to rcParams['animation.html'] = 'html5'
HTML(ani.to_jshtml())
Out[267]:
The three-body problem¶
Notation and assumptions:
- $M_s$: mass of the Sun
- $M_e$: mass of the Earth
- $M_m$: mass of Mars
- $M_s \approx 330000 M_e$, $M_m \approx 0.1 M_e$
- Center of gravity of the three bodies approximately coincides with the center of the Sun (Sun remain still in the model)
- The three objects remain in the same plane described by their initial positions (2D problem)
- $\mathbf{x}_e = (x_e, y_e)$ denote the earth position.
- $\mathbf{x}_m = (x_m, y_m)$ denote the earth position.
Universal gravitational law
\begin{eqnarray} M_e\frac{d^{2} \mathbf{x}_e}{dt^{2}} & = & G M_m M_e \frac{\mathbf{x}_m-\mathbf{x}_e}{\|\mathbf{x}_m-\mathbf{x}_e\|^{3}} -G M_e M_s\frac{\mathbf{x}_e}{\|\mathbf{x}_e\|^{3}} \\ M_m\frac{d^{2} \mathbf{x}_m}{dt^{2}} & = & G M_e M_m \frac{\mathbf{x}_e-\mathbf{x}_m}{\|\mathbf{x}_e-\mathbf{x}_m\|^{3}} -G M_m M_s\frac{\mathbf{x}_m}{\|\mathbf{x}_m\|^{3}} \end{eqnarray}non-dimensional
\begin{eqnarray} \frac{d^{2} \mathbf{x}_e}{dt^{2}} & = & 4 \pi^2\left( \frac{M_m}{M_s} \frac{\mathbf{x}_m-\mathbf{x}_e}{\|\mathbf{x}_m-\mathbf{x}_e\|^{3}} -\frac{\mathbf{x}_e}{\|\mathbf{x}_e\|^{3}} \right)\\ \frac{d^{2} \mathbf{x}_m}{dt^{2}} & = & 4 \pi^2\left( \frac{M_e}{M_s} \frac{\mathbf{x}_e-\mathbf{x}_m}{\|\mathbf{x}_e-\mathbf{x}_m\|^{3}} -\frac{\mathbf{x}_m}{\|\mathbf{x}_m\|^{3}} \right) \end{eqnarray}We introduce the variables \begin{eqnarray} \mathbf{q} &=& \left(q_1, q_2, q_3, q_4\right) = \left(x_e, y_e, x_m, y_m\right) \\ \mathbf{p} &=& \left(p_1, p_2, p_3, p_4\right) = \left(M_e\frac{d x_e}{dt}, M_e\frac{d y_e}{dt}, M_m\frac{d x_m}{dt}, M_m \frac{d y_m}{dt} \right) \end{eqnarray}
Hamiltonian system: \begin{eqnarray} \frac{d q_1}{dt} & = & \frac{p_1}{M_e}\\ \frac{d q_2}{dt} & = & \frac{p_2}{M_e}\\ \frac{d q_3}{dt} & = & \frac{p_3}{M_m}\\ \frac{d q_4}{dt} & = & \frac{p_4}{M_m}\\ \frac{1}{M_e}\frac{d p_1}{dt} & = & 4 \pi^2\left( \frac{M_m}{M_s} \frac{q_3-q_1}{\|(q_3,q_4)-(q_1, q_2)\|^{3}} -\frac{q_1}{\|(q_1, q_2)\|^{3}} \right) \\ \frac{1}{M_e}\frac{d p_3}{dt} & = & 4 \pi^2\left( \frac{M_m}{M_s} \frac{q_4-q_2}{\|(q_3,q_4)-(q_1, q_2)\|^{3}} -\frac{q_2}{\|(q_1, q_2)\|^{3}} \right)\\ \frac{1}{M_m}\frac{d p_3}{dt} & = & 4 \pi^2\left( \frac{M_e}{M_s} \frac{q_1-q_3}{\|(q_1,q_2)-(q_3, q_4)\|^{3}} -\frac{q_3}{\|(q_3, q_4)\|^{3}} \right) \\ \frac{1}{M_m}\frac{d p_4}{dt} & = & 4 \pi^2\left( \frac{M_e}{M_s} \frac{q_2-q_4}{\|(q_1,q_2)-(q_3, q_4)\|^{3}} -\frac{q_4}{\|(q_3, q_4)\|^{3}} \right) \end{eqnarray}
Hamiltonian: $ \mathbf{p} = (p_1,p_2,p_3,p_4), \mathbf{q} = (q_1,q_2,q_3,q_4)$
\begin{equation} H(\mathbf{p}, \mathbf{q}) = \frac{1}{2}\left(\frac{p_1^{2}}{M_e}+\frac{p_2^{2}}{M_e}+\frac{p_3^{2}}{M_m}+\frac{p_4^{2}}{M_m} \right) - \frac{4 \pi^2}{M_s}\left( M_{e} M_{m} \frac{1}{\|(q_1,q_2)-(q_3, q_4)\|} + M_{e} M_{s} \frac{1}{\|(q_1, q_2)\|} + M_{m} M_{s} \frac{1}{\|(q_3, q_4)\|} \right) \end{equation}In [329]:
global M_e, M_m, M_s
M_e = 1.0
M_m = 0.1
M_s = 330000.0
def threebody_problem_Hp(t,p,q):
f1 = p[0]/M_e
f2 = p[1]/M_e
f3 = p[2]/M_m
f4 = p[3]/M_m
return np.array([f1,f2,f3,f4])
def threebody_problem_Hq(t,p,q):
f1 = -M_e*4*np.pi**2*( M_m/M_s * (q[2]-q[0])/((q[2]-q[0])**2+(q[3]-q[1])**2 )**(3/2) - q[0]/(q[0]**2+q[1]**2)**(3/2))
f2 = -M_e*4*np.pi**2*( M_m/M_s * (q[3]-q[1])/((q[2]-q[0])**2+(q[3]-q[1])**2 )**(3/2) - q[1]/(q[0]**2+q[1]**2)**(3/2))
f3 = -M_m*4*np.pi**2*( M_e/M_s * (q[2]-q[0])/((q[2]-q[0])**2+(q[3]-q[1])**2 )**(3/2) - q[2]/(q[2]**2+q[3]**2)**(3/2))
f4 = -M_m*4*np.pi**2*( M_e/M_s * (q[3]-q[1])/((q[2]-q[0])**2+(q[3]-q[1])**2 )**(3/2) - q[3]/(q[2]**2+q[3]**2)**(3/2))
return np.array([f1,f2,f3,f4])
def threebody_problem(t,y):
if len(y)%2 == 0:
n = int(len(y)/2)
p = y[:n]
q = y[n:]
f2 = threebody_problem_Hp(t,p,q)
f1 = -threebody_problem_Hq(t,p,q)
return np.concatenate((f1,f2))
else:
raise NameError('dimension is not even')
In [355]:
# Condicion inicial
p0 = np.array([0.0,-5.1, 0.0, -0.46])
q0 = np.array([1.0,0.0,1.52,0.0])
h = 10/20000; T = 10
# Forward Euler solution
t_fe, y_fe = Forward_Euler( odefun=threebody_problem,
t_span=(0,T),
y0=np.concatenate((p0,q0)),
h=h,
return_trajectory=True)
p_fe = y_fe[:,:4]; q_fe = y_fe[:,4:]
# # Symplectic Euler soluton
t_se, y_se = Symplectic_Euler( odefun=(threebody_problem_Hp, threebody_problem_Hq),
t_span=(0,T),
y0=np.concatenate((p0,q0)),
h=h,
return_trajectory=True)
p_se = y_se[:,:4]; q_se = y_se[:,4:]
In [356]:
%matplotlib inline
earth_pos_x_fe = q_fe[:,0]
earth_pos_y_fe = q_fe[:,1]
mars_pos_x_fe = q_fe[:,2]
mars_pos_y_fe = q_fe[:,3]
earth_pos_x = q_se[:,0]
earth_pos_y = q_se[:,1]
mars_pos_x = q_se[:,2]
mars_pos_y = q_se[:,3]
fig1, ax1 = plt.subplots(1,2, figsize=(12,6));
fig1.suptitle('Comparison of orbits. Forward euler vs Symplectic Euler')
ax1[0].plot(0,0, 'oC1', markersize = 30)
ax1[0].plot(earth_pos_x_fe,earth_pos_y_fe,'oC0')
ax1[0].plot(mars_pos_x_fe,mars_pos_y_fe,'oC3')
ax1[1].plot(0,0, 'oC1', markersize = 30)
ax1[1].plot(earth_pos_x,earth_pos_y,'oC0')
ax1[1].plot(mars_pos_x,mars_pos_y,'oC3')
plt.show()
In [375]:
matplotlib.rcParams['animation.embed_limit'] = 4
ani = animation.FuncAnimation(fig3, animate0, frames=500, interval=interval, blit=True, init_func=init0);
# equivalent to rcParams['animation.html'] = 'html5'
HTML(ani.to_jshtml())
Out[375]:
In [ ]: