Vibración forzada sin amortiguamiento¶
\begin{equation*} m \ddot{u} + k u = F_{0} \sin \bigl( \bar{\omega} t \bigr) \end{equation*}
Reescribiendo la ecuación
\begin{equation} \ddot{u} + \omega^{2} u = \frac{F_{0}}{m} \sin \bigl( \bar{\omega} t \bigr) \end{equation}
La solución es
\begin{equation*} u = A_{1} \sin \bigl( \omega t \bigr) + A_{2} \cos \bigl( \omega t \bigr) + G \sin \bigl( \bar{\omega} t \bigr) \end{equation*}
Coeficientes
\begin{align*} \gamma &= \frac{\bar{\omega}}{\omega} \\ A_{1} &= \frac{\dot{u}_{0}}{\omega} - \frac{F_{0}}{k} \frac{\gamma}{1 - \gamma^{2}}\\ A_{2} &= u_{0} \\ G &= \frac{F_{0}}{k} \frac{1}{1 - \gamma^{2}} \end{align*}
Ejemplo en clase¶
Usando
\begin{equation*} u = A \sin \bigl( \omega t \bigr) + B \cos \bigl( \omega t \bigr) + G \sin \bigl( \bar{\omega} t \bigr) \end{equation*}
import numpy as np
import matplotlib.pyplot as plt
w = 25.0 # Tn
m = w/981.0 # Tn/cm/seg^2
k = 7.0 # Tn/cm
u0 = 0.0 # cm
v0 = 0.0 # cm/seg
beta = 0.07
F0 = 2.0
omega_bar = 18.1
omega = np.power(k/m,1.0/2.0)
gamma = omega_bar/omega
G = (F0/k)*(1.0)/(1.0 - gamma**2)
A = (v0/omega) - (F0/k)*(gamma)/(1.0 - gamma**2)
B = u0
print 'omega =', omega
print 'gamma =',gamma
print 'G =', G
print 'A =', A
print 'B =', B
t = np.linspace(0,15,15000)
posicion = A*np.sin(omega*t) + B*np.cos(omega*t) + G*np.sin(omega_bar*t)
plt.figure(figsize=(19,8.5))
plt.plot(t,posicion)
plt.xlabel(r'$t$ seg.')
plt.ylabel(r'$u$ cm.')
plt.grid(True)
plt.show()
Cuando $u_{0} = 0$ y $\dot{u}_{0} = 0$
\begin{equation*} u = A_{1} \sin \bigl( \omega t \bigr) + G \sin \bigl( \bar{\omega} t \bigr) \end{equation*}
Coeficientes
\begin{align*} \gamma &= \frac{\bar{\omega}}{\omega} \\ A_{1} &= - \frac{F_{0}}{k} \frac{\gamma}{1 - \gamma^{2}}\\ G &= \frac{F_{0}}{k} \frac{1}{1 - \gamma^{2}} \end{align*}
import numpy as np
import matplotlib.pyplot as plt
w = 25.0 # Tn
m = w/981.0 # Tn/cm/seg^2
k = 7.0 # Tn/cm
beta = 0.07
F0 = 2.0
omega_bar = 18.1
omega = np.power(k/m,1.0/2.0)
gamma = omega_bar/omega
G = (F0/k)*(1.0)/(1.0 - gamma**2)
A = -(F0/k)*(gamma)/(1.0 - gamma**2)
print 'omega =', omega
print 'gamma =',gamma
print 'G =', G
print 'A =', A
t = np.linspace(0,15,15000)
posicion = A*np.sin(omega*t) + G*np.sin(omega_bar*t)
plt.figure(figsize=(19,8.5))
plt.plot(t,posicion)
plt.xlabel(r'$t$ seg.')
plt.ylabel(r'$u$ cm.')
plt.grid(True)
plt.show()
Usando odeint
\begin{align*} u^{\prime} &= z \\ z^{\prime} + \omega^{2} u &= \frac{F_{0}}{m} \sin \bigl( \bar{\omega} t \bigr) \end{align*}
sujeto a $u(0) = 0$ y $u^{\prime}(0) = 0$
from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt
def dU_dx(U, x):
# Here U is a vector such that y=U[0] and z=U[1]. This function should return [z', u']
m = 25.0/981.0
k = 7.0
omega = np.power(k/m,1.0/2.0)
F0 = 2.0
omega_bar = 18.1
return [U[1],-omega**2*U[0] + (F0/m)*np.sin(omega_bar*x)]
U0 = [0.0,0.0]
#U0 = [6.5,32.0]
xs = np.linspace(0, 15, 15000)
Us = odeint(dU_dx, U0, xs)
ys = Us[:,0]
plt.figure(figsize=(19,8.5))
plt.plot(xs, ys)
plt.xlabel(r'$t$ seg.')
plt.ylabel(r'$u$ cm.')
plt.grid(True)
plt.show()
Valor máximo aproximado
valor_pico = np.where(posicion >= max(posicion))[0]
print 't pico =', t[valor_pico[0]]
print 'valor pico =', posicion[valor_pico[0]]
Gráfico interactivo
import plotly.offline as offline
import plotly.graph_objs as go
import numpy as np
offline.init_notebook_mode()
t = np.linspace(0,15,15000)
posicion = A*np.sin(omega*t) + B*np.cos(omega*t) + G*np.sin(omega_bar*t)
trace = go.Scatter(x = t, y = posicion)
data = [trace]
offline.iplot(data)
Desplazamientos máximos¶
Usando la transformada de Hilbert
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import hilbert
t = np.linspace(0,15,15000)
posicion = A*np.sin(omega*t) + B*np.cos(omega*t) + G*np.sin(omega_bar*t)
analytic_signal = hilbert(posicion)
amplitude_envelope = np.abs(analytic_signal)
instantaneous_phase = np.unwrap(np.angle(analytic_signal))
instantaneous_frequency = (np.diff(instantaneous_phase)/(2.0*np.pi))
fig = plt.figure(figsize=(19,8.5))
ax0 = fig.add_subplot(211)
ax0.plot(t, posicion, label='posicion')
ax0.plot(t, amplitude_envelope, label='envolvente')
ax0.set_xlabel("tiempo en segundos")
ax0.legend()
ax1 = fig.add_subplot(212)
ax1.plot(t[1:], instantaneous_frequency)
ax1.set_xlabel("tiempo en segundos")
plt.show()
Usando argrelextrema
import numpy as np
from scipy.signal import argrelextrema
import matplotlib.pyplot as plt
t = np.linspace(0,15,15000)
posicion = A*np.sin(omega*t) + B*np.cos(omega*t) + G*np.sin(omega_bar*t)
maximo = argrelextrema(posicion, np.greater)[0] #array of indexes of the locals maxima
x = [t[maximo] for i in maximo][0]
y = [posicion[maximo] for i in maximo][0]
plt.figure(figsize=(19,8.5))
plt.plot(t,posicion)
plt.plot(x, y, 'ro')
plt.xlabel(r'$t$ seg.')
plt.ylabel(r'$u$ cm.')
plt.show()
Desplazamiento máximo
for i in range(np.size(x)):
print i, '\t', round(x[i],3), '\t', round(y[i],3)
Periodo aproximado
for i in range(np.size(x) - 1):
print i+1, '\t', x[i+1]-x[i]
Resonancia¶
\begin{equation*} u = - \frac{F_{0}}{k} \frac{\gamma}{1 - \gamma^{2}} \sin \bigl( \omega t \bigr) + \frac{F_{0}}{k} \frac{1}{1 - \gamma^{2}} \sin \bigl( \bar{\omega} t \bigr) \end{equation*}
Se produce cuando $\gamma = 1$
\begin{equation*} \lim_{\gamma \to 1} u = - \frac{F_{0}}{k} \frac{\gamma}{1 - \gamma^{2}} \sin \bigl( \omega t \bigr) + \frac{F_{0}}{k} \frac{1}{1 - \gamma^{2}} \sin \bigl( \gamma \omega t \bigr) = - \frac{F_{0}}{2 k} \omega t \cos \bigl( \omega t \bigr) + \frac{F_{0}}{2 k} \sin \bigl( \omega t \bigr) \end{equation*}
t = np.linspace(0,15,15000)
posicion = (-F0/k)*omega*t*np.cos(omega*t) + (F0/k)*np.sin(omega*t)
plt.figure(figsize=(19,8.5))
plt.plot(t,posicion)
plt.xlabel(r'$t$ seg.')
plt.ylabel(r'$u$ cm.')
plt.grid(True)
plt.show()
