# ¿Cómo funciona la suspensión de un auto?
 

<div>




</div>
 
> Una primer aproximación al modelo de la suspensión de un automovil es considerar el *oscilador armónico amortiguado*. 



Referencia: 
- https://es.wikipedia.org/wiki/Oscilador_arm%C3%B3nico#Oscilador_arm.C3.B3nico_amortiguado

Un **modelo** que describe el comportamiento del sistema mecánico anterior es

\begin{equation}
m\frac{d^2 x}{dt^2}=-c\frac{dx}{dt}-kx
\end{equation}

donde $c$ es la constante de amortiguamiento y $k$ es la constante de elasticidad. <font color=red> Revisar modelado </font>

Documentación de los paquetes que utilizaremos hoy.
- https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html
- https://docs.scipy.org/doc/scipy/reference/index.html

___
En `python` existe una función llamada <font color = blue>_odeint_</font> del paquete <font color = blue>_integrate_</font> de la libreria <font color = blue>_scipy_</font>, que permite integrar sistemas vectoriales de primer orden, del tipo

\begin{equation}
\frac{d\boldsymbol{y}}{dt} = \boldsymbol{f}(\boldsymbol{y},t); \qquad \text{ con }\quad \boldsymbol{y}\in\mathbb{R}^n,\quad \boldsymbol{f}:\mathbb{R}^n\times\mathbb{R}_{+}\to\mathbb{R}^n
\end{equation}

con condiciones iniciales $\boldsymbol{y}(0) = \boldsymbol{y}_{0}$. Notar que <font color=red> $\boldsymbol{y}$ representa un vector de $n$ componentes</font>.

Ahora, si nos fijamos bien, el modelo del *oscilador armónico amortiguado* que obtuvimos es una ecuación diferencial ordinaria (EDO) de segundo orden. No hay problema. La podemos convertir en un sistema de ecuaciones de primer orden de la siguiente manera:

1. Seleccionamos el vector $\boldsymbol{y}=\left[y_1\quad y_2\right]^T$, con $y_1=x$ y $y_2=\frac{dx}{dt}$.
2. Notamos que $\frac{dy_1}{dt}=\frac{dx}{dt}=y_2$ y $\frac{dy_2}{dt}=\frac{d^2x}{dt^2}=-\frac{c}{m}\frac{dx}{dt}-\frac{k}{m}x=-\frac{c}{m}y_2-\frac{k}{m}y_1$.
3. Entonces, el modelo de segundo orden lo podemos representar como el siguiente sistema vectorial de primer orden:

\begin{equation}
\frac{d\boldsymbol{y}}{dt}=\left[\begin{array}{c}\frac{dy_1}{dt} \\ \frac{dy_2}{dt}\end{array}\right]=\left[\begin{array}{c}y_2 \\ -\frac{k}{m}y_1-\frac{c}{m}y_2\end{array}\right]=\left[\begin{array}{cc}0 & 1 \\-\frac{k}{m} & -\frac{c}{m}\end{array}\right]\boldsymbol{y}.
\end{equation}


```python
# Primero importamos todas las librerias, paquetes y/o funciones que vamos a utlizar
from matplotlib import pyplot as plt
from scipy.integrate import odeint
import numpy as np
```


```python
# Ayuda función odeint
odeint?
```


```python
# Función f(y,t) que vamos a integrar
def amortiguado(y, t, k, m, c):
    y1 = y[0]
    y2 = y[1]
    return np.array([y2,
                     -k / m * y1 - c / m * y2])
```


```python
# Definimos los parámetros k, m y c
k = 3
m = 30
c = 5
# Condiciones iniciales
x0 = 0.5  # m
dx0 = 0   # m/s
y0 = [x0, dx0]
# Especificamos los puntos de tiempo donde queremos la solución
t = np.linspace(0, 100, 1000)
```


```python
# Solución numérica
y = odeint(func=amortiguado,
           y0=y0,
           t=t,
           args=(k, m, c))
```

¿Cómo entrega odeint las soluciones?


```python
# Averiguar la forma de solución
y
```




    array([[ 5.00000000e-01,  0.00000000e+00],
           [ 4.99750917e-01, -4.96265435e-03],
           [ 4.99009377e-01, -9.83828153e-03],
           ...,
           [ 4.12621841e-05,  3.30995434e-05],
           [ 4.45264114e-05,  3.21260247e-05],
           [ 4.76923498e-05,  3.11367079e-05]])




```python
# Mostrar la solución 
y.shape
```




    (1000, 2)



- $y$ es una matriz de n filas y 2 columnas.
- La primer columna de $y$ corresponde a $y_1$.
- La segunda columna de $y$ corresponde a $y_2$.

¿Cómo extraemos los resultados $y_1$ y $y_2$ independientemente?


```python
# Extraer y1 y y2
x = y[:, 0]
dx = y[:, 1]
```

### Para hacer participativamente...
- Graficar en una misma ventana $y_1$ vs. $t$ y $y_2$ vs. $t$... ¿qué pueden observar?


```python
# Gráfica
plt.figure(figsize=(6, 4))
plt.plot(t, x, label="Posición $y_1(t)=x(t)$")
plt.plot(t, dx, label="Velocidad $y_2(t)=x'(t)$")
plt.xlabel("Tiempo")
plt.legend()
```

- Graficar $y_2/\omega_0$ vs. $y_1$... ¿cómo se complementan estos gráficos? ¿conclusiones?


```python
# Gráfica
plt.figure(figsize=(6, 4))
plt.plot(x, dx)
plt.xlabel("Posición")
plt.ylabel("Velocidad")
```

## Dependiendo de los parámetros, 3 tipos de soluciones

Teníamos
\begin{equation}
m\frac{d^2 x}{dt^2} + c\frac{dx}{dt} + kx = 0
\end{equation}

si recordamos que $\omega_0 ^2 = \frac{k}{m}$ y definimos $\frac{c}{m}\equiv 2\Gamma$, tendremos

\begin{equation}
\frac{d^2 x}{dt^2} + 2\Gamma \frac{dx}{dt}+ \omega_0^2 x = 0
\end{equation}

<font color=blue>El comportamiento viene determinado por las raices de la ecuación característica. Ver en el tablero...</font>

### Subamortiguado

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


```python
k = 3
m = 30
c = 5
```


```python
w02 = k / m
G = c / (2 * m)
```


```python
w02, G**2
```




    (0.1, 0.006944444444444444)




```python
w02 > G**2
```




    True



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


```python
# Gráfica, de nuevo
plt.figure(figsize=(6, 4))
plt.plot(t, x, label="Posición $y_1(t)=x(t)$")
plt.plot(t, dx, label="Velocidad $y_2(t)=x'(t)$")
plt.xlabel("Tiempo")
plt.legend()
```


```python
plt.figure(figsize=(6, 4))
plt.plot(x, dx)
plt.xlabel("Posición")
plt.ylabel("Velocidad")
```

### Sobreamortiguado

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


```python
# Nuevas constantes 
k = .1 # Constante del muelle
m = 1.0   # Masa
c = 1 # Constante de amortiguación
```

Simular y graficar...


```python
w02 = k / m
G = c / (2 * m)
```


```python
w02, G
```




    (0.1, 0.5)




```python
w02 < G**2
```




    True




```python
# Simular
y = odeint(func=amortiguado,
           y0=y0,
           t=t,
           args=(k, m, c))
xs = y[:, 0]
dxs = y[:, 1]
```


```python
# Graficar
plt.figure(figsize=(6, 4))
plt.plot(t, xs, label="Posición $y_1(t)=x(t)$")
plt.plot(t, dxs, label="Velocidad $y_2(t)=x'(t)$")
plt.xlabel("Tiempo")
plt.legend()
```


```python
plt.figure(figsize=(6, 4))
plt.plot(xs, dxs)
plt.xlabel("Posición")
plt.ylabel("Velocidad")
```

### Amortiguamiento crítico

Si $\omega_0^2 = \Gamma^2$ se tiene movimiento *críticamente amortiguado*. 


```python
# Nuevas constantes
k = .0625 # Constante del muelle
m = 1.0   # Masa
c = .5 # Constante de amortiguación
```

Simular y graficar...


```python
w02 = k / m
G = c / (2 * m)
```


```python
w02, G**2
```




    (0.0625, 0.0625)




```python
w02 == G**2
```




    True




```python
# Simular
y = odeint(func=amortiguado,
           y0=y0,
           t=t,
           args=(k, m, c))
xc = y[:, 0]
dxc = y[:, 1]
```


```python
# Graficar
plt.figure(figsize=(6, 4))
plt.plot(t, xc, label="Posición $y_1(t)=x(t)$")
plt.plot(t, dxc, label="Velocidad $y_2(t)=x'(t)$")
plt.xlabel("Tiempo")
plt.legend()
```


```python
plt.figure(figsize=(6, 4))
plt.plot(xc, dxc)
plt.xlabel("Posición")
plt.ylabel("Velocidad")
```

En resumen, se tiene entonces: 


```python
tt = t
fig, ((ax1, ax2, ax3), (ax4, ax5, ax6)) = plt.subplots(2, 3, sharex='col',
                                                       sharey='row',figsize =(10,6))
ax1.plot(tt, x, c = 'k') 
ax1.set_title('Amortiguado', fontsize = 14)
ax1.set_ylabel('Posición', fontsize = 14)
ax2.plot(tt, xs, c = 'b') 
ax2.set_title('Sobreamortiguado', fontsize = 14)
ax3.plot(tt, xc, c = 'r') 
ax3.set_title('Crítico', fontsize = 16)
ax4.plot(tt, dx, c = 'k')
ax4.set_ylabel('Velocidad', fontsize = 14)
ax4.set_xlabel('tiempo', fontsize = 14)
ax5.plot(tt, dxs, c = 'b')
ax5.set_xlabel('tiempo', fontsize = 14)
ax6.plot(tt, dxc, c = 'r')
ax6.set_xlabel('tiempo', fontsize = 14)
plt.show()
```

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

> En un gráfico como el anterior, realizar gráficas del plano fase para los distintos movimientos y para cuatro conjuntos de condiciones iniciales distintas
- y0 = [1, 1]
- y0 = [1, -1]
- y0 = [-1, 1]
- y0 = [-1, -1]

Hacer lo anterior en un nuevo notebook de jupyter llamado Tarea7_ApellidoNombre.ipynb y subir en el espacio habilitado.



<footer id="attribution" style="float:right; color:#808080; background:#fff;">
Created with Jupyter by Lázaro Alonso. Modified by Esteban Jiménez Rodríguez.
</footer>

c507d37f77e9c734dcdb8e55c971048b6cca8c10
ipynb
Jupyter Notebook
Modulo3/Clase16_OsciladorAmortiguado.ipynb
Damaris02/simmatv2021
cfd3b49a50c543293b30cc797739ecdf7725d886
MIT
Modulo3/Clase16_OsciladorAmortiguado.ipynb
Damaris02/simmatv2021
cfd3b49a50c543293b30cc797739ecdf7725d886
MIT
Modulo3/Clase16_OsciladorAmortiguado.ipynb
Damaris02/simmatv2021
cfd3b49a50c543293b30cc797739ecdf7725d886
MIT
2021-05-25T21:51:37.000Z
2021-06-15T00:44:51.000Z
Qwen/Qwen-72B
1. YES
2. YES