El péndulo doble es uno de los primeros modelos que describen la teoría del caos en los sistemas dinámicos. El péndulo doble es la unión de dos péndulos simples como se puede observar en la siguiente imagen respectivamente.
 |
| DCL del Péndulo Doble |
Usando la segunda ley de Newton obtenemos las siguientes ecuaciones de movimiento del péndulo doble respectivamente.
Para la masa ${\large M_{1}}$ tenemos las siguientes ecuaciones de movimiento
$${\large \begin{align*} M_{1}\frac{d^2 x_{1}}{dt^2} &= -T_{1}\sin{\theta_{1}} + T_{2}\sin{\theta_{2}} \\ M_{1}\frac{d^2 y_{1}}{dt^2} &= T_{1}\cos{\theta_{1}} - T_{2}\cos{\theta_{2}} - M_{1}g \end{align*} }$$
Para la masa ${\large M_{2}}$ tenemos las siguientes ecuaciones de movimiento
$${\large \begin{align*} M_{2}\frac{d^2 x_{2}}{dt^2} &= -T_{2}\sin{\theta_{2}} \\ M_{2}\frac{d^2 y_{2}}{dt^2} &= T_{2}\cos{\theta_{2}} - M_{2}g \end{align*} }$$
Para resolver las ecuaciones de movimiento tenemos que encontrar los desplazamientos angulares ${\large \theta_{1}}$ y ${\large \theta_{2}}$. Por consiguiente se obtienen las siguientes ecuaciones diferenciales respectivas.
$${\large \begin{align*} \ddot{\theta_{1}} &= \frac{-g(2M_{1} + M_{2})\sin{(\theta_{1})} -gM_{2}\sin{(\theta_{1} - 2\theta_{2})} - 2L_{2}M_{2}\dot{\theta_{2}}^2\sin{(\theta_{1} -\theta_{2})} - \dot{\theta_{1}}^2 L_{1}\sin{(2(\theta_{1} - \theta_{2}))}}{L_{1}(2M_{1}+M_{2} -M_{2}\cos{2(\theta_{1} - \theta_{2})})} \\ \ddot{\theta_{2}} &= \frac{2L_{1}(M_{1} + M_{2})\dot{\theta_{1}}^2 \sin{(\theta_{1} - \theta_{2})} + 2g(M_{1} + M_{2})\cos{(\theta_{1})}\sin{(\theta_{1} - \theta_{2})} + L_{2}M_{2}\dot{\theta_{2}}^2 \sin{(2(\theta_{1} - \theta_{2}))}}{L_{2}(2M_{1}+M_{2} -M_{2}\cos{(2(\theta_{1} - \theta_{2})})}\end{align*}}$$
Como las ecuaciones diferenciales anteriores son difíciles de resolver de forma analítica entonces se procede a resolver de forma numérica.
from numpy import sin, cos
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integrate
from IPython.display import HTML
import matplotlib.animation as animation
from collections import deque
def pendouble(L1, L2, M1, M2, th1, th2, w1, w2, T, N):
""" T = es el tiempo de simulación
[M1, M2] son las masas del péndulo
[L1, L2] son las longitudes de las barras
[th1, th2] son los ángulos iniciales del péndulo
[w1, w2] son las velocidades iniciales del péndulo
N es el número de pasos de la curva """
def derivs(state, t):
dydx = np.zeros_like(state)
dydx[0] = state[1]
delta = state[0] - state[2]
den1 = (2*M1+M2) * L1 - M2 * L1 * cos(2*delta)
dydx[1] = (- 9.8*(2*M1+M2)*sin(state[0]) - M2*9.8*sin(delta - state[2])
- 2*sin(delta)*M2*L2*state[3]*state[3]
- state[1]*state[1]*L1*sin(2*delta) ) / (den1)
dydx[2] = state[3]
den2 = (L2/L1) * den1
dydx[3] = ( 2*(M1+M2)*L1*state[1]*state[1]*sin(delta) +
2*9.8*(M1+M2)*cos(state[0])*sin(delta) +
state[3]*state[3]*L2*M2*sin(2*delta) ) / (den2)
return dydx
# crear una matriz de tiempo
dt = 0.02
t = np.arange(0, T, dt)
# estado inicial
state = np.radians([th1, w1, th2, w2])
# integrando la EDO usando scipy.integrate
y = integrate.odeint(derivs, state, t)
x1 = L1*sin(y[:, 0])
y1 = -L1*cos(y[:, 0])
x2 = L2*sin(y[:, 2]) + x1
y2 = -L2*cos(y[:, 2]) + y1
fig = plt.figure(figsize=(16, 10))
ax = fig.add_subplot(autoscale_on=False, xlim=(-L1-L2, L1+L2), ylim=(-L1-L2, 1.0))
ax.set_aspect('equal')
ax.text(-2.0, 1.100, "Simulación del Pendulo Doble para theta1 = {}$^\circ$, theta2 {}$^\circ$".format(th1, th2)
, size = 20)
ax.grid()
O = ax.plot(0,0, color='r', marker='o', markeredgecolor='k', lw = 6, ms = 20)
o1, = ax.plot([],[], color='r', marker='o', markeredgecolor='k', lw = 6, ms = 20)
o2, = ax.plot([],[], color='r', marker='o', markeredgecolor='k', lw = 6, ms = 20)
line, = ax.plot([], [], 'o-', lw=4, color='b')
trace, = ax.plot([], [], ',-', lw=2, color='k')
time_template = 'tiempo = %.1fs'
time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes, size= 20)
history_x, history_y = deque(maxlen=N), deque(maxlen=N)
def animate(i):
thisx = [0, x1[i], x2[i]]
thisy = [0, y1[i], y2[i]]
p1x = [x1[i], x1[i]]
p1y = [y1[i], y1[i]]
p2x = [x2[i], x2[i]]
p2y = [y2[i], y2[i]]
if i == 0:
history_x.clear()
history_y.clear()
history_x.appendleft(thisx[2])
history_y.appendleft(thisy[2])
line.set_data(thisx, thisy)
o1.set_data(p1x, p1y)
o2.set_data(p2x, p2y)
trace.set_data(history_x, history_y)
time_text.set_text(time_template % (i*dt))
return line, trace, time_text
ani = animation.FuncAnimation(fig, animate, len(y), interval=dt*1000, blit=True)
return HTML(ani.to_html5_video())
Aplicación del código fuente
Aplicando el codigo fuente del péndulo doble obtenemos una animación para el pendulo doble con sus respectivas condiciones iniciales
 |
| ${\large pendouble(1, 1, 1, 1, 120, 60, 0, 0, 20, 1000)}$ |
Como se menciono anteriormente, el péndulo doble es uno de los primeros modelos que describen la teoría del caos. Por ejemplo para las siguientes condiciones iniciales se puede observar este comportamiento caótico.
${\large L_{1} = L_{2} = 1 \hspace{0.2cm} m }$
${\large w_{1} = w_{2} = 0 \hspace{0.2cm} rad/s }$
${\large \theta_{1} = 1\hspace{0.2cm} rad }$
${\large M_{1} = 1.5 \hspace{0.2cm} Kg }$
${\large M_{2} = 2.0 \hspace{0.2cm} Kg }$
 |
| Diferentes $\theta_{2}$ para cada simulación |
Los diferentes $\theta_{2}$ para los péndulos dobles para sus correspondientes colores en la animación anterior.
${\large \theta_{2} = 1.998 \hspace{0.2cm} rad }$
${\large \theta_{2} = 2.0 \hspace{0.2cm} rad }$
${\large \theta_{2} = 2.002 \hspace{0.2cm} rad }$
Comentarios
Publicar un comentario