In [65]:
import numpy as np
from math import sin
#Para resolver la ecución diff
from scipy.integrate import odeint
#Para los plot
import matplotlib.pyplot as plt
from matplotlib import animation
#Para mostrar la animación inline
from JSAnimation import IPython_display
#Magic para hacer los plot inline
%matplotlib inline
In [52]:
#Constantes del problema
g = 9.8 #m/s^2
l = 1 #m
g_l = g/l
y0 = [4, 0] #Condiciones iniciales. [theta] = rad, omega(0) = 0
t_output = np.arange(0, 3, 0.01) #tiempo
def ecu(x_vec, t):
theta, omega = x_vec # tuple unpacking
return [omega, -g_l*sin(theta)]
y_result = odeint(ecu, y0, t_output)
tetha_result, omega_result = y_result.T # extract x and y cols
In [58]:
x = np.sin(tetha_result)*l
y = -np.cos(tetha_result)*l
x2 = np.arange(-1, 1, 0.0001)
y2 = np.sqrt(1 - x2**2)
y3 = -y2
plt.figure(figsize = (8, 8))
plt.plot(x, y, linewidth = 2, label = 'Mov. pendulo')
plt.plot(x2, y2, color = 'red', ls = '--', linewidth = 3, label = r'$1 = x^{2} + y^{2}$')
plt.plot(x2, y3, color = 'red', ls = '--', linewidth = 3)
plt.grid(True)
plt.ylim(ymin = -1.2, ymax = 1.2)
plt.xlim(xmin = -1.2, xmax = 1.2)
plt.title(r'$Pendulo\; Simple$', fontsize=20)
plt.xlabel(r'$Posicion\; Horizonta$l', fontsize=20)
plt.ylabel(r'$Posicion\; Vertical$', fontsize=20)
plt.legend(loc = "best", fontsize=20)
plt.show()
In [68]:
# Primero seteamos la figure, los ejes y los elementos del plot que se van a animar
fig = plt.figure(figsize = (6, 6))
ax = plt.axes(xlim=(-1.2, 1.2), ylim=(-1.2, 1.2))
ax.grid()
line, = ax.plot([], [], 'o-', lw=2)
# Background de cada frame
def init():
line.set_data([], [])
return line,
# Función que hace la animación
def animate(i):
thisx = [0, x[i]]
thisy = [0, y[i]]
line.set_data(thisx, thisy)
return line,
animation.FuncAnimation(fig, animate, init_func=init, frames=len(x), interval=20, blit=True)
#anim.save('pendulo_simple.mp4', fps=30, extra_args=['-vcodec', 'libx264'])
Out[68]:
In [60]:
h = l*(1 - np.cos(tetha_result))
v = omega_result*l
e_k = 1/2.*v**2
e_p = g*h
e_total = e_k + e_p
plt.figure(figsize = (10,8))
plt.plot(t_output, e_k, linewidth = 2, label = r'$Energ\acute{\imath}a\; Cinetica$')
plt.plot(t_output, e_p, linewidth = 2, label = r'$Energ\acute{\imath}a\; Potencial$')
plt.plot(t_output, e_total, linewidth = 2)
plt.grid(True)
plt.legend(loc = "lower right", fontsize=20)
plt.title(r'$Energ\acute{\imath}a\; del\; Sistema$', fontsize=20)
plt.xlabel('$Tiempo$', fontsize=20)
plt.ylabel(r'$Energ \acute{i} a $', fontsize=20)
plt.show()
Observar que el movimiento del péndulo es un circulo perfecto, coincidiendo con la teoría. Ademas, la energía total del sistema se conserva.
PNER IMAGEN DEL PENDULO
In [6]:
N = 1000 # number of steps to take
y = np.zeros ([4])
L_o = 1.0
L = 1.0
v_o = 0.0
theta_o = 0.3
omega_o = 0.0
# set initial state
y[0] = L
y[1] = v_o
y[2] = theta_o
y[3] = omega_o
time = np.linspace(0 , 25 , N)
k = 3.5 # spring constant, in N/m
m = 0.2 # mass, in kg
gravity = 9.8 # g, in m/sˆ2
def spring_pendulum(y, time):
g0 = y[1]
g1 = (L_o + y[0]) * y[3] * y[3] - k/m*y[0] + gravity * np.cos(y[2])
g2 = y[3]
g3 = -(gravity * np.sin(y[2]) + 2.0 * y[1] * y[3])/(L_o + y[0])
return np.array([g0,g1,g2,g3])
# Now we do the calculations.
answer = odeint(spring_pendulum,y,time)
# Now graph the results.
# rather than graph in terms oft, I’m going
# to graph the track the mass takes in 2D.
# This will require that I change L, theta data
# to x, y data .
xdata = (L_o + answer[:,0]) * np.sin(answer[:, 2])
ydata = -(L_o + answer[:,0]) * np.cos(answer[:,2])
plt.figure(figsize=(10,10))
plt.plot(xdata, ydata, 'r-')
plt.xlabel('Horizontal position')
plt.ylabel('Vertical position')
Out[6]:
In [1]:
# Execute this cell to load the notebook's style sheet, then ignore it
from IPython.core.display import HTML
css_file = '/home/franco/Desktop/Congresos_y_Cursos/Numerical-mooc/numericalmoocstyle.css'
HTML(open(css_file, "r").read())
Out[1]:
In [ ]: