\begin{equation*} \ddot{u} + \omega^{2} u = \begin{cases} \frac{F_{0}}{m} & 0 \leqslant t \leqslant t_{1} \\ 0 & t \geqslant 0 \end{cases} \end{equation*}
La solución es
\begin{equation*} u = \begin{cases} A_{1} \sin \bigl( \omega t \bigr) + B_{1} \cos \bigl( \omega t \bigr) + \frac{F_{0}}{k} & 0 \leqslant t \leqslant t_{1} \\ A_{2} \sin \bigl( \omega t \bigr) + B_{2} \cos \bigl( \omega t \bigr) & t \geqslant 0 \end{cases} \end{equation*}
Coeficientes
\begin{align*} A_{1} &= \frac{\dot{u}_{0}}{\omega} \\ B_{1} &= u_{0} - \frac{F_{0}}{k} \\ A_{2} &= \frac{\dot{u}(t_{1})}{\omega} \\ B_{2} &= u(t_{1}) \end{align*}
Usando una traslación de coordenadas en la vibración libre
\begin{equation*} u = \begin{cases} A_{1} \sin \bigl( \omega t \bigr) + B_{1} \cos \bigl( \omega t \bigr) + \frac{F_{0}}{k} & 0 \leqslant t \leqslant t_{1} \\ A_{2} \sin \bigl[ \omega \bigl( t - t_{1} \bigr) \bigr] + B_{2} \cos \bigl[ \omega \bigl( t - t_{1} \bigr) \bigr] & t \geqslant t_{1} \end{cases} \end{equation*}
import numpy as np
import matplotlib.pyplot as plt
w = 42.0 # Tn
m = w/981.0 # Tn/cm/seg^2
k = 10.5 # Tn/cm
u0 = 0.0 # cm
v0 = 0.0 # cm/seg
F0 = 2.0 # Tn
t1 = 1.0 # segundos
omega = np.power(k/m,1.0/2.0)
A1 = v0/omega
B1 = u0 - F0/k
A2 = A1*np.cos(omega*t1) - B1*np.sin(omega*t1)
B2 = A1*np.sin(omega*t1) + B1*np.cos(omega*t1) + F0/k
print 'omega =', omega
print 'A1 =', A1
print 'B1 =', B1
print 'A2 =', A2
print 'B2 =', B2
t = np.linspace(0,3,3000)
posicion = np.where(t<=t1, A1*np.sin(omega*t) + B1*np.cos(omega*t) + F0/k, A2*np.sin(omega*(t-t1)) + B2*np.cos(omega*(t-t1)))
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 = \begin{cases} B_{1} \cos \bigl( \omega t \bigr) + \frac{F_{0}}{k} & 0 \leqslant t \leqslant t_{1} \\ A_{2} \sin \bigl[ \omega \bigl( t - t_{1} \bigr) \bigr] + B_{2} \cos \bigl[ \omega \bigl( t - t_{1} \bigr) \bigr] & t \geqslant t_{1} \end{cases} \end{equation*}
Coeficientes
\begin{align*} B_{1} &= - \frac{F_{0}}{k} \\ A_{2} &= \frac{\dot{u}(t_{1})}{\omega} \\ B_{2} &= u(t_{1}) \end{align*}
import numpy as np
import matplotlib.pyplot as plt
w = 42.0 # Tn
m = w/981.0 # Tn/cm/seg^2
k = 10.5 # Tn/cm
F0 = 2.0 # Tn
t1 = 1.0 # segundos
omega = np.power(k/m,1.0/2.0)
B1 = -F0/k
A2 = -B1*np.sin(omega*t1)
B2 = B1*np.cos(omega*t1) + F0/k
print 'omega =', omega
print 'B1 =', B1
print 'A2 =', A2
print 'B2 =', B2
t = np.linspace(0,3,3000)
posicion = np.where(t<=t1, B1*np.cos(omega*t) + F0/k, A2*np.sin(omega*(t-t1)) + B2*np.cos(omega*(t-t1)))
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 &= \begin{cases} \frac{F_{0}}{m} & 0 \leqslant t \leqslant 1 \\\ 0 & t \geqslant 1 \end{cases} \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 = 42.0/981.0
k = 10.5
omega = np.power(k/m,1.0/2.0)
F0 = 2.0
t1 = 1.0 # segundos
if x <= t1:
return [U[1],-omega**2*U[0] + F0/m]
else:
return [U[1],-omega**2*U[0]]
U0 = [0.0,0.0]
xs = np.linspace(0, 3, 3000)
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()
Periodo
\begin{equation*} T = \frac{2 \pi}{\omega} \end{equation*}
T = (2.0*np.pi)/omega
print 'T =', T
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,3,3000)
posicion = np.where(t<=t1, A1*np.sin(omega*t) + B1*np.cos(omega*t) + F0/k, A2*np.sin(omega*(t-t1)) + B2*np.cos(omega*(t-t1)))
trace = go.Scatter(x = t, y = posicion)
data = [trace]
offline.iplot(data)
Usando la transformada de Hilbert
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import hilbert
t = np.linspace(0,3,3000)
posicion = np.where(t<=t1, A1*np.sin(omega*t) + B1*np.cos(omega*t) + F0/k, A2*np.sin(omega*(t-t1)) + B2*np.cos(omega*(t-t1)))
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,3,3000)
posicion = np.where(t<=t1, A1*np.sin(omega*t) + B1*np.cos(omega*t) + F0/k, A2*np.sin(omega*(t-t1)) + B2*np.cos(omega*(t-t1)))
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]
Desplazamientos
\begin{equation*} u = \begin{cases} \frac{F_{0}}{k} \bigl[ 1 - \cos \bigl( \omega t \bigr) \bigr] & 0 \leqslant t \leqslant t_{1} \\ A_{2} \sin \bigl( \omega t \bigr) + B_{2} \cos \bigl( \omega t \bigr) & t \geqslant 0 \end{cases} \end{equation*}
Valores conocidos
\begin{align*} \omega t &= 2 \pi \frac{t}{T} \\ \sqrt{A_{2}^{2} + B_{2}^{2}} &= \frac{F_{0}}{k} \sqrt{\bigl[ \sin \bigl( \omega t \bigr) \bigr]^{2} + \bigl[ 1 - \cos \bigl( \omega t \bigr) \bigr]^{2}} = \frac{F_{0}}{k} \sqrt{2 \bigl[ 1 - \cos \bigl( \omega t \bigr) \bigr]} = \frac{F_{0}}{k} \big | 2 \sin \bigl( \frac{\omega t}{2} \bigr) \big | \end{align*}
Reemplazando
\begin{equation*} u = \begin{cases} \frac{F_{0}}{k} \bigl[ 1 - \cos \bigl( 2 \pi \frac{t}{T} \bigr) \bigr] & 0 \leqslant t \leqslant t_{1} & \mbox{respuesta forzada} \\ \frac{2 F_{0}}{k} & t \geqslant t_{1} & \mbox{respuesta forzada} \\ \frac{2 F_{0}}{k} \big | \sin \bigl( \pi \frac{t}{T} \bigr) \big | & t \geqslant 0 & \mbox{respuesta libre} \end{cases} \end{equation*}
Reescribiendo
\begin{equation*} u = \begin{cases} \frac{F_{0}}{k} \bigl[ 1 - \cos \bigl( 2 \pi \frac{t}{t_{1}} \frac{t_{1}}{T} \bigr) \bigr] & 0 \leqslant \frac{t_{1}}{T} \leqslant \frac{1}{2} & \mbox{respuesta forzada} \\ \frac{2 F_{0}}{k} & \frac{t_{1}}{T} \geqslant \frac{1}{2} & \mbox{respuesta forzada} \\ \frac{2 F_{0}}{k} \big | \sin \bigl( \pi \frac{t}{t_{1}} \frac{t_{1}}{T} \bigr) \big | & \frac{t_{1}}{T} \geqslant 0 & \mbox{respuesta libre} \end{cases} \end{equation*}
Para simplificar $\frac{t}{t_{1}} = 1$
t = np.linspace(0,3,3000)
u0 = 2*np.abs(np.sin(np.pi*t))
u1 = np.where(t<=0.5, 1 - np.cos(2*np.pi*t), 2)
plt.figure(figsize=(11,8.5))
plt.plot(t,u0, '--', label='respuesta libre')
plt.plot(t,u1, label='respuesta forzada')
plt.legend(loc='upper right')
plt.xlabel(r'$\frac{t_{1}}{T}$')
plt.ylabel(r'$R_{d}$')
plt.grid(True)
plt.show()
El espectro de respuesta es la envolvente
\begin{equation*} u_{\max} = \begin{cases} \frac{2 F_{0}}{k} \big | \sin \bigl( \pi \frac{t_{1}}{T} \bigr) \big | & 0 \leqslant \frac{t_{1}}{T} \leqslant \frac{1}{2} & \mbox{respuesta libre} \\ \frac{2 F_{0}}{k} & \frac{t_{1}}{T} \geqslant \frac{1}{2} & \mbox{respuesta forzada} \end{cases} \end{equation*}
t = np.linspace(0,3,3000)
u = np.where(t<=0.5, 2*np.abs(np.sin(np.pi*t)), 2)
plt.figure(figsize=(11,8.5))
plt.plot(t,u)
plt.xlabel(r'$\frac{t_{1}}{T}$')
plt.ylabel(r'$R_{d}$')
plt.grid(True)
plt.show()
Factor de amplificación máximo
valor_pico = np.where(u>= max(u))[0]
print 't1/T =', 1/(2.0*t[valor_pico[0]])
print 'Rd =', u[valor_pico[0]]