Autor: Eduardo Vélez Moreno
Pendulo Doble
Librerías
In [1]:
#Librerías
import numpy as np
from numpy import sin, cos
from [Link] import odeint
import [Link] as plt
Parametros del sistema
In [2]:
#Características del sistema
g = 9.8 # Aceleración de la gravedad
l1 = 1.0 # Longitud del primer pendulo
l2 = 1.0 # Longitud del segundo pendulo
m1 = 1.0 # Masa del pendulo 1 in kg
m2 = 1.0 # Masa del pendulo 2 in kg
Condiciones iniciales para pequeñas oscilaciones
In [3]:
#Condiciones iniciales para pequeñas oscilaciones
x1 = 10.0 # Posicion angular del pendulo 1
x2 = 0.0 # Velocidad angular del pendulo 1
x3 = -10.0 # Posicion angulor del pendulo 2
x4 = 0.0 # Velocidad angular del pendulo 2
Condiciones iniciales para grandes oscilaciones
In [4]:
#Condiciones iniciales para oscilaciones grandes
x5 = 120.0 # Posicion angular del pendulo 1
x6 = 0.0 # Velocidad angular del pendulo 1
x7 = -20 # Posicion angulor del pendulo 2
x8 = 0.0 # Velocidad angular del pendulo 2
Vector de estado
In [5]:
#Vector de estados
estado1 = [Link]([x1, x2, x3, x4]) #Pequelos oscilaciones
estado2 = [Link]([x5, x6, x7, x8]) #Oscilaciones grandes
Vector de tiempo
In [6]:
#Tiempo de simulacion
T=0.05 # Periodo de tiempo
t=[Link](0.0,30,T) # Vector de tiempo
Funcion que describe el Modelo no lineal
q
In [7]:
def nolineal(estado, t):
# Vector de ODES:
# Crea una vector del tamaño de estado.
# Servira para guardar las ecuaciones diferenciales de primer orden
dydx = np.zeros_like(estado)
# ODE 1
# Velocidad angular del primer pendulo : derivada de x_1 = x_2
dydx[0] = estado[1]
# ODE 2
# Aceleracion angular del primer pendulo : derivada de x_2
den1=(m1+m2)-m2*(cos(estado[0]-estado[2]))**2 #Denominador comun
dydx[1] = -(m2*((estado[1])**2)*sin(estado[0]-estado[2])*cos(estado[0]-estado[2]
+(m2*g*sin(estado[2])*cos(estado[0]-estado[2]))/(den1*l1) \
-(m2*l2*((estado[3])**2)*sin(estado[0]-estado[2]))/(den1*l1) \
-((m1+m2)*g*sin(estado[0]))/(den1*l1)
# ODE 3
# Velocidad angular del segundo pendulo : derivada de x_3 = x_4
dydx[2] = estado[3]
# ODE 4
# Aceleracion angular del segundo pendulo : derivada de x_4
den2 = l2*den1
dydx[3] = +((m1+m2)*l1*(estado[1]**2)*sin(estado[0]-estado[2]))/den2 \
+(m2*(estado[3]**2)*sin(estado[0]-estado[2])*cos(estado[0]-estado[2]))
+((m1+m2)*g*sin(estado[0])*cos(estado[0]-estado[2]))/den2 \
-((m1+m2)*g*sin(estado[2]))/den2
return dydx
Funcion que describe el modelo lineal
In [8]:
def lineal(estado, t):
# Vector de ODES:
# Crea una vector del tamaño de estado.
# Servira para guardar las ecuaciones diferenciales de primer orden
dydx2 = np.zeros_like(estado)
# ODE 1
# Velocidad angular del primer pendulo : derivada de x_1 = x_2
dydx2[0] = estado[1]
# ODE 2
# Aceleracion angular del primer pendulo : derivada de x_2
den1=m1*l1 #Denominador comun
dydx2[1] = (-(m1+m2)*g*estado[0]+m2*g*estado[2])/den1
# ODE 3
# Velocidad angular del segundo pendulo : derivada de x_3 = x_4
dydx2[2] = estado[3]
# ODE 4
# Aceleracion angular del segundo pendulo : derivada de x_4
den2 = m1*l2
dydx2[3] = ((m1+m2)*g*estado[0]-(m1+m2)*g*estado[2])/den2
return dydx2
Creacion de graficos para pequeñas oscilaciones
In [9]:
#Posicion y velocidad del pendulo 1
y1 = odeint(nolineal, estado1, t)
y3 = odeint(lineal, estado1, t)
fig, axs = [Link](4, 2,figsize=(15,25))
#Grafica 1 'Posicion y Velocidad para la masa 1: Modelo no lineal'
axs[0, 1].set_title('Posicion y Velocidad para la masa 1: Modelo no lineal')
axs[0, 1].plot(t,y1[:,0])
axs[0, 1].plot(t,y1[:,1])
#Grafica 2 'Posicion y Velocidad para la masa 1: Modelo lineal'
axs[0, 0].set_title('Posicion y Velocidad para la masa 1: Modelo lineal')
axs[0, 0].plot(t,y3[:,0])
axs[0, 0].plot(t,y3[:,1])
#Grafica 3 'Plano de fase: Posicion vs Velocidad, masa 1: Modelo no lineal'
axs[1, 1].set_title('Plano de fase: Posicion vs Velocidad, masa 1: Modelo no lineal'
axs[1, 1].plot(y1[:,0],y1[:,1],color=[0.1,0.1,.1])
#Grafica 4 'Plano de fase, Posicion vs Velocidad, masa 1: Modelo lineal'
axs[1, 0].set_title('Plano de fase, Posicion vs Velocidad, masa 1: Modelo lineal')
axs[1, 0].plot(y3[:,0],y3[:,1],color=[0.1,0.1,.1])
#Grafica 5 'Posicion y Velocidad para la masa 2: Modelo no lineal'
axs[2, 1].set_title('Posicion y Velocidad para la masa 2: Modelo no lineal')
axs[2, 1].plot(t,y1[:,2])
axs[2, 1].plot(t,y1[:,3])
#Grafica 6 'Posicion y Velocidad para la masa 2: Modelo lineal'
axs[2, 0].set_title('Posicion y Velocidad para la masa 2: Modelo lineal')
axs[2, 0].plot(t,y3[:,2])
axs[2, 0].plot(t,y3[:,3])
#Grafica 7 'Plano de fase de la masa 2: Modelo no lineal'
axs[3, 1].set_title('Plano de fase: Posicion vs Velocidad, masa 2: Modelo no lineal'
axs[3, 1].plot(y1[:,2],y1[:,3],color=[0.1,0.1,.1])
#Grafica 8 'Plano de fase de la masa 2: Modelo lineal'
axs[3, 0].set_title('Plano de fase: Posicion vs Velocidad, masa 2: Modelo lineal')
axs[3, 0].plot(y3[:,2],y3[:,3],color=[0.1,0.1,.1])
[Link]()
Creación de graficos para grandes oscilaciones
In [10]:
#Posicion y velocidad del pendulo 1
y2 = odeint(nolineal, estado2, t)
y4 = odeint(lineal, estado2, t)
fig, axs = [Link](4, 2,figsize=(15,25))
#Grafica 1 'Posicion y Velocidad para la masa 1: Modelo no lineal'
axs[0, 1].set_title('Posicion y Velocidad para la masa 1: Modelo no lineal')
axs[0, 1].plot(t,y2[:,0])
axs[0, 1].plot(t,y2[:,1])
#Grafica 2 'Posicion y Velocidad para la masa 1: Modelo lineal'
axs[0, 0].set_title('Posicion y Velocidad para la masa 1: Modelo lineal')
axs[0, 0].plot(t,y4[:,0])
axs[0, 0].plot(t,y4[:,1])
#Grafica 3 'Plano de fase: Posicion vs Velocidad, masa 1: Modelo no lineal'
axs[1, 1].set_title('Plano de fase: Posicion vs Velocidad, masa 1: Modelo no lineal'
axs[1, 1].plot(y2[:,0],y2[:,1],color=[0.1,0.1,.1])
#Grafica 4 'Plano de fase, Posicion vs Velocidad, masa 1: Modelo lineal'
axs[1, 0].set_title('Plano de fase, Posicion vs Velocidad, masa 1: Modelo lineal')
axs[1, 0].plot(y4[:,0],y4[:,1],color=[0.1,0.1,.1])
#Grafica 5 'Posicion y Velocidad para la masa 2: Modelo no lineal'
axs[2, 1].set_title('Posicion y Velocidad para la masa 2: Modelo no lineal')
axs[2, 1].plot(t,y2[:,2])
axs[2, 1].plot(t,y2[:,3])
#Grafica 6 'Posicion y Velocidad para la masa 2: Modelo lineal'
axs[2, 0].set_title('Posicion y Velocidad para la masa 2: Modelo lineal')
axs[2, 0].plot(t,y4[:,2])
axs[2, 0].plot(t,y4[:,3])
#Grafica 7 'Plano de fase de la masa 2: Modelo no lineal'
axs[3, 1].set_title('Plano de fase: Posicion vs Velocidad, masa 2: Modelo no lineal'
axs[3, 1].plot(y2[:,2],y2[:,3],color=[0.1,0.1,.1])
#Grafica 8 'Plano de fase de la masa 2: Modelo lineal'
axs[3, 0].set_title('Plano de fase: Posicion vs Velocidad, masa 2: Modelo lineal')
axs[3, 0].plot(y4[:,2],y4[:,3],color=[0.1,0.1,.1])
#[Link]()
[<[Link].Line2D at 0x2b9157f6d00>]
Out[10]:
Creacion de las animaciones
In [11]: import [Link] as animation
from [Link] import PillowWriter
%matplotlib
#Descomentar para oscilaciones grandes
#y1=y2;
#y3=y4;
#Modelo no lineal pequeñas oscilaciones
xa1 = l1*sin(y1[:, 0])
ya1 = -l1*cos(y1[:, 0])
xa2 = l2*sin(y1[:, 2]) + xa1
ya2 = -l2*cos(y1[:, 2]) + ya1
#Modelo lineal pequeñas oscilaciones
xa5 = l1*sin(y3[:, 0])
ya5 = -l1*cos(y3[:, 0])
xa6 = l2*sin(y3[:, 2]) + xa1
ya6 = -l2*cos(y3[:, 2]) + ya1
fig = [Link]()
[Link]('Pendulo doble: No lineal vs Lineal')
gs = fig.add_gridspec(3, 2, hspace=0, wspace=0)
[Link](sharex='col', sharey='row')
ax = fig.add_subplot(321, autoscale_on=True, xlim=(-3, 3), ylim=(-3, 3))
[Link]('off')
ax.set_aspect('equal')
ax2 = fig.add_subplot(322, autoscale_on=True, xlim=(-3, 3), ylim=(-3, 3))
[Link]('off')
ax2.set_aspect('equal')
ax3 = fig.add_subplot(323, autoscale_on=True, xlim=(0, 30), ylim=(-10, 10))
[Link]('off')
ax4 = fig.add_subplot(324, autoscale_on=True, xlim=(0, 30), ylim=(-10, 10))
[Link]('off')
ax5 = fig.add_subplot(325, autoscale_on=True, xlim=(0, 30), ylim=(-10, 10))
[Link]('off')
ax6 = fig.add_subplot(326, autoscale_on=True, xlim=(0, 30), ylim=(-10, 10))
[Link]('off')
#[Link]()
#Pendulo pequeñas oscilaciones modelo no lineal
line, = [Link]([], [], 'o-', lw=1)
path1line, = [Link]([], [], 'green', lw=.5)
path2line, = [Link]([], [], 'black', lw=.5)
#Posicion masa1
m1pos, = [Link]([], [], 'orange', lw=1)
#Velocidad masa1
m1vel, = [Link]([], [],'r', lw=1)
#Posicion masa2
m2pos, = [Link]([], [], 'orange', lw=1)
#Posicion masa2
m2vel, = [Link]([], [], 'r', lw=1)
#Pendulo pequeñas oscilaciones modelo lineal
line2, = [Link]([], [], 'o-', lw=1)
path1line2, = [Link]([], [], 'green', lw=.5)
path2line2, = [Link]([], [], 'black', lw=.5)
#Posicion masa2
m2pos1, = [Link]([], [], 'orange', lw=1)
#Velocidad masa1
m2vel1, = [Link]([], [],'r', lw=1)
#Posicion masa2
m2pos2, = [Link]([], [], 'orange', lw=1)
#Posicion masa2
m2vel2, = [Link]([], [], 'r', lw=1)
def init():
line.set_data([], [])
line2.set_data([], [])
return [line,line2]
def animate(i):
#Pendulo masa 1
thisx = [0, xa1[i], xa2[i]]
thisy = [0, ya1[i], ya2[i]]
line.set_data(thisx, thisy)
path1line.set_data(xa1[:i], ya1[:i])
path2line.set_data(xa2[:i], ya2[:i])
#Trayectoria
thisx1 = [0, xa5[i], xa6[i]]
thisy1 = [0, ya5[i], ya6[i]]
line2.set_data(thisx1, thisy1)
path1line2.set_data(xa5[:i], ya5[:i])
path2line2.set_data(xa6[:i], ya6[:i])
m1pos.set_data(t[:i], y1[:,0][:i])
m1vel.set_data(t[:i], y1[:,1][:i])
m2pos1.set_data(t[:i], y3[:,0][:i])
m2vel1.set_data(t[:i], y3[:,1][:i])
m2pos.set_data(t[:i], y1[:,2][:i])
m2vel.set_data(t[:i], y1[:,3][:i])
m2pos2.set_data(t[:i], y3[:,2][:i])
m2vel2.set_data(t[:i], y3[:,3][:i])
return [line,line2,m1pos,m1vel,m2pos,m2vel,path1line,path2line,path1line2,path2l
ani = [Link](fig, animate, [Link](1, len(y1)),
interval=25, blit=True, init_func=init)
#[Link]("[Link]", writer=PillowWriter(fps=24))
[Link]()
Using matplotlib backend: TkAgg
C:\Program Files\Python39\lib\site-packages\ipykernel\[Link]: RuntimeWarn
ing: coroutine 'Kernel.do_one_iteration' was never awaited
[Link]()
RuntimeWarning: Enable tracemalloc to get the object allocation traceback
