Se dice que un sistema cualquiera, mecánico, eléctrico, neumático, etc., es un oscilador armónico si, cuando se deja en libertad fuera de su posición de equilibrio, vuelve hacia ella describiendo oscilaciones sinusoidales, o sinusoidales amortiguadas en torno a dicha posición estable.
En realidad esto es el estudio de oscilaciones.
In [1]:
from IPython.display import YouTubeVideo
YouTubeVideo('k5yTVHr6V14')
Out[1]:
Los sistemas mas sencillos a estudiar en oscilaciones son el sistema masa-resorte y el péndulo simple.
La solución a este sistema masa-resorte se explica en términos de la segunda ley de Newton. Para este caso, si la masa permanece constante y solo consideramos la dirección en $x$. Entonces,
\begin{equation}
F = m \frac{d^2x}{dt^2}.
\end{equation}
¿Cuál es la fuerza? Ley de Hooke! \begin{equation} F = -k x, \quad k > 0. \end{equation}
Vemos que la fuerza se opone al desplazamiento y su intensidad es proporcional al mismo. Y $k$ es la constante elástica o recuperadora del resorte.
Entonces, un modelo del sistema masa-resorte está descrito por la siguiente ecuación diferencial:
cuya solución se escribe como \begin{equation} x(t) = A \cos(\omega_{o} t) + B \sin(\omega_{o} t) \end{equation} Y su primera derivada (velocidad) sería \begin{equation} \frac{dx(t)}{dt} = \omega_{0}[- A \sin(\omega_{0} t) + B\cos(\omega_{0}t)] \end{equation}
Ver en el tablero que significa solución de la ecuación diferencial.
Esta instrucción es para que las gráficas aparezcan dentro de este entorno.
In [2]:
%matplotlib inline
Esta es la librería con todas las instrucciones para realizar gráficos.
In [3]:
import matplotlib.pyplot as plt
In [4]:
import matplotlib as mpl
label_size = 14
mpl.rcParams['xtick.labelsize'] = label_size
mpl.rcParams['ytick.labelsize'] = label_size
Y esta es la librería con todas las funciones matemáticas necesarias.
In [6]:
import numpy as np
In [19]:
# Definición de funciones a graficar
A, B, w0 = .5, .1, .5 # Parámetros
t = np.linspace(0, 50, 100) # Creamos vector de tiempo de 0 a 50 con 100 puntos
x = A*np.cos(w0*t)+B*np.sin(w0*t) # Función de posición
dx = w0*(-A*np.sin(w0*t)+B*np.cos(w0*t)) # Función de velocidad
# Gráfico
plt.figure(figsize = (7, 4)) # Ventana de gráfica con tamaño
plt.plot(t, x, '-', lw = 1, ms = 4,
label = '$x(t)$') # Explicación
plt.plot(t, dx, 'ro-', lw = 1, ms = 4,
label = r'$\dot{x(t)}$')
plt.xlabel('$t$', fontsize = 20) # Etiqueta eje x
plt.show()
In [22]:
# Colores, etiquetas y otros formatos
plt.figure(figsize = (7, 4))
plt.scatter(t, x, lw = 0, c = 'red',
label = '$x(t)$') # Gráfica con puntos
plt.plot(t, x, 'r-', lw = 1) # Grafica normal
plt.scatter(t, dx, lw = 0, c = 'b',
label = r'$\frac{dx}{dt}$') # Con la r, los backslash se tratan como un literal, no como un escape
plt.plot(t, dx, 'b-', lw = 1)
plt.xlabel('$t$', fontsize = 20)
plt.legend(loc = 'best') # Leyenda con las etiquetas de las gráficas
plt.show()
Y si consideramos un conjunto de frecuencias de oscilación, entonces
In [20]:
frecuencias = np.array([.1, .2 , .5, .6]) # Vector de diferentes frecuencias
plt.figure(figsize = (7, 4)) # Ventana de gráfica con tamaño
# Graficamos para cada frecuencia
for w0 in frecuencias:
x = A*np.cos(w0*t)+B*np.sin(w0*t)
plt.plot(t, x, '*-')
plt.xlabel('$t$', fontsize = 16) # Etiqueta eje x
plt.ylabel('$x(t)$', fontsize = 16) # Etiqueta eje y
plt.title('Oscilaciones', fontsize = 16) # Título de la gráfica
plt.show()
Estos colores, son el default de matplotlib, sin embargo existe otra librería dedicada, entre otras cosas, a la presentación de gráficos.
In [26]:
import seaborn as sns
sns.set(style='ticks', palette='Set2')
In [31]:
frecuencias = np.array([.1, .2 , .5, .6])
plt.figure(figsize = (7, 4))
for w0 in frecuencias:
x = A*np.cos(w0*t)+B*np.sin(w0*t)
plt.plot(t, x, 'o-',
label = '$\omega_0 = %s$'%w0) # Etiqueta cada gráfica con frecuencia correspondiente (conversion float a string)
plt.xlabel('$t$', fontsize = 16)
plt.ylabel('$x(t)$', fontsize = 16)
plt.title('Oscilaciones', fontsize = 16)
plt.legend(loc='center left', bbox_to_anchor=(1.05, 0.5), prop={'size': 14})
plt.show()
Si queremos tener manipular un poco mas las cosas, hacemos uso de lo siguiente:
In [33]:
from ipywidgets import *
In [35]:
def masa_resorte(t = 0):
A, B, w0 = .5, .1, .5 # Parámetros
x = A*np.cos(w0*t)+B*np.sin(w0*t) # Función de posición
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.plot(x, [0], 'ko', ms = 10)
ax.set_xlim(xmin = -0.6, xmax = .6)
ax.axvline(x=0, color = 'r')
ax.axhline(y=0, color = 'grey', lw = 1)
fig.canvas.draw()
interact(masa_resorte, t = (0, 50,.01));
La opción de arriba generalmente será lenta, así que lo recomendable es usar interact_manual.
In [38]:
def masa_resorte(t = 0):
A, B, w0 = .5, .1, .5 # Parámetros
x = A*np.cos(w0*t)+B*np.sin(w0*t) # Función de posición
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.plot(x, [0], 'ko', ms = 10)
ax.set_xlim(xmin = -0.6, xmax = .6)
ax.axvline(x=0, color = 'r')
ax.axhline(y=0, color = 'grey', lw = 1)
fig.canvas.draw()
interact_manual(masa_resorte, t = (0, 50,.01));
Ahora, si fijamos nuestra atención al movimiento de un péndulo simple (oscilaciones pequeñas), la ecuación diferencial a resolver tiene la misma forma:
\begin{equation} \frac{d^2 \theta}{dt^2} + \omega_{0}^{2}\, \theta = 0, \quad\mbox{donde}\quad \omega_{0}^2 = \frac{g}{l}. \end{equation}La diferencia más evidente es como hemos definido a $\omega_{0}$. Esto quiere decir que,
\begin{equation} \theta(t) = A\cos(\omega_{0} t) + B\sin(\omega_{0}t) \end{equation}Si graficamos la ecuación de arriba vamos a encontrar un comportamiento muy similar al ya discutido anteriormente. Es por ello que ahora veremos el movimiento en el plano $xy$. Es decir,
\begin{align} x &= l \sin(\theta), \quad y = l \cos(\theta) \end{align}
In [39]:
# Podemos definir una función que nos entregue theta dados los parámetros y el tiempo
def theta_t(a, b, g, l, t):
omega_0 = np.sqrt(g/l)
return a * np.cos(omega_0 * t) + b * np.sin(omega_0 * t)
In [41]:
# Hacemos un gráfico interactivo del péndulo
def pendulo_simple(t = 0):
fig = plt.figure(figsize = (5,5))
ax = fig.add_subplot(1, 1, 1)
x = 2 * np.sin(theta_t(.4, .6, 9.8, 2, t))
y = - 2 * np.cos(theta_t(.4, .6, 9.8, 2, t))
ax.plot(x, y, 'ko', ms = 10)
ax.plot([0], [0], 'rD')
ax.plot([0, x ], [0, y], 'k-', lw = 1)
ax.set_xlim(xmin = -2.2, xmax = 2.2)
ax.set_ylim(ymin = -2.2, ymax = .2)
fig.canvas.draw()
interact_manual(pendulo_simple, t = (0, 10,.01));
Realmente lo que se tiene que resolver es,
\begin{equation} \theta(t) = \theta(0) \cos(\omega_{0} t) + \frac{\dot{\theta}(0)}{\omega_{0}} \sin(\omega_{0} t) \end{equation}Actividad. Modificar el programa anterior para incorporar las condiciones iniciales.
In [120]:
# Solución:
def theta_t():
return
In [ ]:
def pendulo_simple(t = 0):
fig = plt.figure(figsize = (5,5))
ax = fig.add_subplot(1, 1, 1)
x = 2 * np.sin(theta_t( , t))
y = - 2 * np.cos(theta_t(, t))
ax.plot(x, y, 'ko', ms = 10)
ax.plot([0], [0], 'rD')
ax.plot([0, x ], [0, y], 'k-', lw = 1)
ax.set_xlim(xmin = -2.2, xmax = 2.2)
ax.set_ylim(ymin = -2.2, ymax = .2)
fig.canvas.draw()
interact_manual(pendulo_simple, t = (0, 10,.01));
La posición y velocidad para el sistema masa-resorte se escriben como:
In [43]:
k = 3 #constante elástica [N]/[m]
m = 1 # [kg]
omega_0 = np.sqrt(k/m)
x_0 = .5
dx_0 = .1
In [58]:
t = np.linspace(0, 50, 300)
In [59]:
x_t = x_0 *np.cos(omega_0 *t) + (dx_0/omega_0) * np.sin(omega_0 *t)
dx_t = -omega_0 * x_0 * np.sin(omega_0 * t) + dx_0 * np.cos(omega_0 * t)
In [48]:
plt.figure(figsize = (7, 4))
plt.plot(t, x_t, label = '$x(t)$', lw = 1)
plt.plot(t, dx_t, label = '$\dot{x}(t)$', lw = 1)
#plt.plot(t, dx_t/omega_0, label = '$\dot{x}(t)$', lw = 1) # Mostrar que al escalar, la amplitud queda igual
plt.legend(loc='center left', bbox_to_anchor=(1.01, 0.5), prop={'size': 14})
plt.xlabel('$t$', fontsize = 18)
plt.show()
In [60]:
plt.figure(figsize = (5, 5))
plt.plot(x_t, dx_t/omega_0, 'ro', ms = 2)
plt.xlabel('$x(t)$', fontsize = 18)
plt.ylabel('$\dot{x}(t)/\omega_0$', fontsize = 18)
plt.show()
In [61]:
plt.figure(figsize = (5, 5))
plt.scatter(x_t, dx_t/omega_0, cmap = 'viridis', c = dx_t, s = 8, lw = 0)
plt.xlabel('$x(t)$', fontsize = 18)
plt.ylabel('$\dot{x}(t)/\omega_0$', fontsize = 18)
plt.show()
In [53]:
k = 3 #constante elástica [N]/[m]
m = 1 # [kg]
omega_0 = np.sqrt(k/m)
In [79]:
t = np.linspace(0, 50, 50)
In [63]:
x_0s = np.array([.7, .5, .25, .1])
dx_0s = np.array([.2, .1, .05, .01])
cmaps = np.array(['viridis', 'inferno', 'magma', 'plasma'])
In [65]:
plt.figure(figsize = (6, 6))
for indx, x_0 in enumerate(x_0s):
x_t = x_0 *np.cos(omega_0 *t) + (dx_0s[indx]/omega_0) * np.sin(omega_0 *t)
dx_t = -omega_0 * x_0 * np.sin(omega_0 * t) + dx_0s[indx] * np.cos(omega_0 * t)
plt.scatter(x_t, dx_t/omega_0, cmap = cmaps[indx],
c = dx_t, s = 10,
lw = 0)
plt.xlabel('$x(t)$', fontsize = 18)
plt.ylabel('$\dot{x}(t)/\omega_0$', fontsize = 18)
#plt.legend(loc='center left', bbox_to_anchor=(1.05, 0.5))
Trayectorias del oscilador armónico simple en el espacio fase $(x,\, \dot{x}\,/\omega_0)$ para diferentes valores de la energía.