In [1]:
import numpy as np
from scipy.integrate import ode, odeint
import matplotlib.pyplot as plt
%matplotlib inline
In [20]:
M = 0.5
t0 = 0 # initial start time
tfinal = 5 # final start time
dt = .005 # time step to solution (dopri5 algorithm uses adaptive)
t_output = np.arange(t0, tfinal, dt) #tiempo en el cual voy a resolver las ecuaciones
y0 = np.array([.25,0.0,1,1]) # initial position and velocity
def ecu(Y, t, M):
x, y, v_x, v_y = Y
d3 = (x**2 + y**2)**1.5
return [v_x, v_y, -M*x/d3, -M*y/d3]
y_result = odeint(ecu, y0, t_output, args=(M,))
x_res, y_res, vx_res, vy_res = y_result.T # extract x and y cols
plt.figure()
plt.plot(x_res, y_res)
Out[20]:
In [3]:
'''
Numerical solution of the one-body problem using dopri5
'''
# Compute RHS of ODE, f(t,Y)
def f(t,Y):
x,y,vx,vy = Y # define individual values
d3 = (x**2+y**2)**1.5 # ||x||^3
M = .5
return np.array([vx,vy,-M*x/d3,-M*y/d3])
# initial parameters
t0 = 0 # initial start time
tfinal = 5 # final start time
dt = .005 # time step to solution (dopri5 algorithm uses adaptive)
y0 = np.array([.25,0.0,0.,.45]) # initial position and velocity
# initiate integrator object
test = ode(f).set_integrator('dopri5',atol=1e-6) # prescribe tolerance for adaptive time step
test.set_initial_value(y0,t0) # set initial time and initial value
fig = plt.figure() # initialize figure
ax = fig.add_subplot(111) # name of the plot
while test.successful() and test.t < tfinal:
# integrate
test.integrate(test.t+dt)
# plot position
ax.plot(test.y[0],test.y[1],'b.',markersize=2)
ax.plot(0,0,'oy',markersize=12)
ax.set_xlim([-.15,.35])
ax.set_ylim([-.12,.15])
ax.set_title('One-body orbit')
In [2]:
def euler_step(u, f, dt):
"""Returns the solution at the next time-step using Euler's method.
Parameters
----------
u : array of float
solution at the previous time-step.
f : function
function to compute the right hand-side of the system of equation.
dt : float
time-increment.
Returns
-------
u_n_plus_1 : array of float
approximate solution at the next time step.
return u + dt * f(u)
"""
return u + dt * f(u)
def f(u):
"""Returns the right-hand side of the rocket flight system of equations.
Parameters
----------
u : array of float
array containing the solution at time n.
Returns
-------
dudt : array of float
array containing the RHS given u.
"""
r1x = u[0]
r1y = u[1]
r2x = u[2]
r2y = u[3]
v1x = u[4]
v1y = u[5]
v2x = u[6]
v2y = u[7]
dist_1_2 = ((r1x - r2x)**2 + (r1y - r2y)**2)**1.5 # ||r1 - r2||^3
return np.array([v1x, v1y, v2x, v2y, G*m2/dist_1_2*(r1x - r2x), G*m2/dist_1_2*(r1y - r2y),\
G*m1/dist_1_2*(r2x - r1x), G*m1/dist_1_2*(r2y - r1y)])
In [8]:
#Parametros del modelo
G = 0.5
m1 = 5.
m2 = 0.
T = 400.0 # final time
dt = 0.1 # time increment
N = int(T/dt) + 1 # number of time-steps
t = np.linspace(0.0, T, N) # time discretization
#Codiciones iniciales
r1x = .25
r1y = 0.
r2x = 0.
r2y = 0.
v1x = 1.
v1y = 1.
v2x = 0.
v2y = 0.
# initialize the array containing the solution for each time-step
u = np.empty((N, 8))
u[0] = np.array([r1x, r1y, r2x, r2y, v1x, v1y, v2x, v2y])# fill 1st element with initial values
# time loop - Euler method
for n in range(N-1):
u[n+1] = euler_step(u[n], f, dt)
In [9]:
plt.plot(u[:,0], u[:,1])
plt.plot(u[:,2], u[:,3])
Out[9]:
In [2]:
# Esta celda es para darle un estilo al notebook. Fue diseñado por Lorena Barba (@LorenaABarba)
from IPython.core.display import HTML
css_file = '/home/franco/Desktop/Proyectos/Pewen/css/IPython_personal_style.css'
HTML(open(css_file, "r").read())
Out[2]:
In [ ]: