### ¿Cómo funciona la suspensión de un auto?
> Una primer aproximación al modelo de la suspensión de un automovil es considerar el oscilador armónico amortiguado. El cual se representa a través de la siguiente ecuación diferencial. 

\begin{equation}
m\ddot{x} + k x + B \dot{x} = 0
\end{equation}
donde $k$ es la constante del muelle, y $B$ la constante de amortiguación. 

Referencia: 
   - https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html
   - https://docs.scipy.org/doc/scipy/reference/index.html

___
<div>




 </div>

Esta es una ecuación diferencial ordinaria (EDO). En python existe una función llamada _odeint_ del paquete _integrate_ de la libreria _scipy_, que permite integrar sistemas del tipo 
\begin{equation}
\frac{dy}{dt} = f(x,y)
\end{equation}
con condiciones iniciales $y(0) = y_{0}$. Ahora bien, si nos fijamos bien, la ecuación diferencial que tenemos es de segundo orden. No hay problema. La podemos simplificar como un sistema de ecuaciones de primer orden como sigue:

\begin{align}
\dot{x} & = y \\
\dot{y} & = -\frac{k}{m} x - \frac{B}{m} y  
\end{align}


```python
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint 
%matplotlib inline
```


```python
k = 3.0 # Constante del muelle
m = 1.0   # Masa
B = .5 # Constante de amortiguación

def armonico(variables, t):
    x, y = variables
    return [y, -k * x / m - B / m * y]

inicial = [0.6, 0.4] # Vector de posición inicial y velocidad inicial
# condiciones iniciales x(t=0)=0.6 [m] y(t=0)=0.4 [m/s]

#tiempo = np.linspace(0, 15)  # Dominio temporal de 0 a 15 
tiempo = np.arange(0, 20, .01)

resultado = odeint(armonico, inicial, tiempo) 
#El sistema se resuelve  con 
#odeint(sistema, condiciones iniciales, rango donde graficaremos)
```


```python
xx, yy = resultado.T  # extraer posición y velocidad. 
```


```python
import matplotlib as mpl
label_size = 14
mpl.rcParams['xtick.labelsize'] = label_size 
mpl.rcParams['ytick.labelsize'] = label_size 
```


```python
plt.plot(tiempo, xx, c = 'r', label="Posicion") 
plt.plot(tiempo, yy, c = 'k', label="Velocidad")
plt.legend(loc = 'best', prop={'size': 14})    
plt.xlabel('tiempo', fontsize = 14) 
plt.show()
```


```python
omega0 = k/m
```


```python
plt.figure(figsize = (6,6))
plt.scatter(xx, yy/omega0, lw = 0, s = 3, cmap = 'viridis', c = xx)
plt.show()
```


```python
from ipywidgets import *
```


```python
def amortiguado(t = 0):
    fig = plt.figure()
    ax = fig.add_subplot(1, 1, 1)
    ax.plot(tiempo, xx,  'k-')
    ax.plot(tiempo[t], xx[t],  'ro')
    ax.text(4, .6, 'tiempo =    %s'%tiempo[t])
    ax.text(4, .5, 'posición =  %s'%xx[t])
    fig.canvas.draw()
t_f = len(tiempo) 
interact_manual(amortiguado, t = (0, t_f, 1));
```

### Casos

Teníamos
\begin{equation}
m\ddot{x} + k x + B \dot{x} = 0
\end{equation}

si recordamos que $\omega_0 ^2 = k/m$ y definimos $B/m\equiv 2\Gamma$, tendremos

\begin{equation}
\ddot{x} + 2\Gamma \dot{x}+ \omega_0^2 x = 0
\end{equation}

### Amortiguado

Si $\omega_0^2 > \Gamma^2$ se tiene movimiento oscilatorio amortiguado. 


```python
omega0 = k/m
Gamma = B/(2*m)
```


```python
omega0**2, Gamma**2
```




    (9.0, 0.0625)




```python
omega0**2 > Gamma**2
```




    True



Entonces, el primer caso que ya habíamos presentado corresponde a movimiento amortiguado. 


```python
plt.plot(tiempo, xx, c = 'r', label="Posicion") 
plt.plot(tiempo, yy, c = 'k', label="Velocidad")
plt.legend(loc = 'best', prop={'size': 14})    
plt.xlabel('tiempo', fontsize = 14) 
plt.show()
```

### Sobreamortiguado

Si $\omega_0^2 < \Gamma^2$ se tiene movimiento oscilatorio amortiguado. 


```python
k = .1 # Constante del muelle
m = 1.0   # Masa
B = .5 # Constante de amortiguación
```


```python
omega0 = k/m
Gamma = B/(2*m)
```


```python
omega0**2, Gamma**2
```




    (0.010000000000000002, 0.0625)




```python
omega0**2 < Gamma**2
```




    True




```python
inicial = [0.6, 0.4] 
tiempo = np.arange(0, 20, .01)
resultado = odeint(armonico, inicial, tiempo) 
```


```python
xxA, yyA = resultado.T  # extraer posición y velocidad. 
```


```python
plt.plot(tiempo, xxA, c = 'r', label="Posicion") 
plt.plot(tiempo, yyA, c = 'k', label="Velocidad")
plt.legend(loc = 'best', prop={'size': 14})    
plt.xlabel('tiempo', fontsize = 14) 
plt.show()
```

### Amortiguamiento crítico

Si $\omega_0^2 = \Gamma^2$ se tiene movimiento aperiódico crítico (amortiguamiento crítico). 


```python
k = np.sqrt(.0625) # Constante del muelle
m = 1.0   # Masa
B = .5 # Constante de amortiguación
```


```python
omega0 = k/m
Gamma = B/(2*m)
```


```python
omega0**2, Gamma**2
```




    (0.0625, 0.0625)




```python
omega0**2 == Gamma**2
```




    True




```python
inicial = [0.6, 0.4] 
tiempo = np.arange(0, 20, .01)
resultado = odeint(armonico, inicial, tiempo) 
```


```python
xxC, yyC = resultado.T  # extraer posición y velocidad. 
```


```python
plt.plot(tiempo, xxC, c = 'r', label="Posicion") 
plt.plot(tiempo, yyC, c = 'k', label="Velocidad")
plt.legend(loc = 'best',prop={'size': 14})  
plt.xlabel('tiempo', fontsize = 14) 
plt.show()
```

En resumen, se tiene entonces: 


```python
fig, ((ax1, ax2, ax3), (ax4, ax5, ax6)) = plt.subplots(2, 3, sharex='col',
                                                       sharey='row',figsize =(10,6))
ax1.plot(tiempo, xx, c = 'k') 
ax1.set_title('Amortiguado', fontsize = 14)
ax1.set_ylabel('Posición', fontsize = 14)
ax2.plot(tiempo, xxA, c = 'b') 
ax2.set_title('Sobreamortiguado', fontsize = 14)
ax3.plot(tiempo, xxC, c = 'r') 
ax3.set_title('Crítico', fontsize = 16)
ax4.plot(tiempo, yy, c = 'k')
ax4.set_ylabel('Velocidad', fontsize = 14)
ax4.set_xlabel('tiempo', fontsize = 14)
ax5.plot(tiempo, yyA, c = 'b')
ax5.set_xlabel('tiempo', fontsize = 14)
ax6.plot(tiempo, yyC, c = 'r')
ax6.set_xlabel('tiempo', fontsize = 14)
plt.show()
```

> **Actividad**. ¿Cómo se ve el espacio fase para los diferentes casos así como para diferentes condiciones iniciales? 



<footer id="attribution" style="float:right; color:#808080; background:#fff;">
Created with Jupyter by Lázaro Alonso.
</footer>

300d5e5d525e22c0cd4b5af5cfca5294cf9d87d2
ipynb
Jupyter Notebook
SimulacionO2017/Modulo 1/Clase5_OsciladorAmortiguado.ipynb
jcmartinez67/SimulacionMatematica_O2017
085394f939c37717acbbe19a90cad661b279fbee
MIT
SimulacionO2017/Modulo 1/Clase5_OsciladorAmortiguado.ipynb
jcmartinez67/SimulacionMatematica_O2017
085394f939c37717acbbe19a90cad661b279fbee
MIT
SimulacionO2017/Modulo 1/Clase5_OsciladorAmortiguado.ipynb
jcmartinez67/SimulacionMatematica_O2017
085394f939c37717acbbe19a90cad661b279fbee
MIT
Qwen/Qwen-72B
1. YES
2. YES