# Solución de ecuaciones diferenciales

Dada la siguiente ecuación diferencial:

$$
\dot{x} = -x
$$

queremos obtener la respuesta del sistema que representa, es decir, los valores que va tomando $x$.

Si analizamos esta ecuación diferencial, podremos notar que la solución de este sistema es una función $\varphi(t)$, tal que cuando la derivemos obtengamos el negativo de esta misma función, es decir:

$$
\frac{d}{dt} \varphi(t) = -\varphi(t)
$$

y despues de un poco de pensar, podemos darnos cuenta de que la función que queremos es:

$$
\varphi(t) = e^{-t}
$$

Sin embargo muchas veces no tendremos funciones tan sencillas (ciertamente no es el caso en la robótica, donde usualmente tenemos ecuaciones diferenciales no lineales de orden $n$), por lo que en esta práctica veremos algunas estrategias para obtener soluciones a esta ecuación diferencial, tanto numéricas como simbolicas.

## Método de Euler

El [método de Euler](http://es.wikipedia.org/wiki/Método_de_Euler) para obtener el comportamiento de una ecuación diferencial, se basa en la intuición básica de la derivada; digamos que tenemos una ecuación diferencial general:

$$
\frac{dy}{dx} = y' = F(x, y)
$$

en donde $F(x, y)$ puede ser cualquier función que depende de $x$ y/o de $y$, entonces podemos dividir en pedazos el comportamiento de la gráfica de tal manera que solo calculemos un pequeño pedazo cada vez, aproximando el comportamiento de la ecuación diferencial, con el de una recta, cuya pendiente será la derivada:



«<a href="http://commons.wikimedia.org/wiki/File:M%C3%A9todo_de_Euler.jpg#/media/File:M%C3%A9todo_de_Euler.jpg">Método de Euler</a>» por <a href="//commons.wikimedia.org/w/index.php?title=User:Vero.delgado&amp;action=edit&amp;redlink=1" class="new" title="User:Vero.delgado (la página no existe)">Vero.delgado</a> - <span class="int-own-work" lang="es">Trabajo propio</span>. Disponible bajo la licencia <a href="http://creativecommons.org/licenses/by-sa/3.0" title="Creative Commons Attribution-Share Alike 3.0">CC BY-SA 3.0</a> vía <a href="//commons.wikimedia.org/wiki/">Wikimedia Commons</a>.

Esta recta que aproxima a la ecuación diferencial, podemos recordar que tiene una estructura:

$$
y = b + mx
$$

por lo que si sustituimos en $m$ la derivada y $b$ con el valor anterior de la ecuación diferencial, obtendremos algo como:

$$
\overbrace{y_{i+1}}^{\text{nuevo valor de }y} = \overbrace{y_i}^{\text{viejo valor de }y} + \overbrace{\frac{dy}{dx}}^{\text{pendiente}} \overbrace{\Delta x}^{\text{distancia en }x}
$$

pero conocemos el valor de $\frac{dy}{dx}$, es nuestra ecuación diferencial; por lo que podemos escribir esto como:

$$
y_{i+1} = y_i + F(x_i, y_i) \Delta x
$$

Resolvamos algunas iteraciones de nuestro sistema; empecemos haciendo 10 iteraciones a lo largo de 10 segundos, con condiciones iniciales $x(0) = 1$, eso quiere decir que:

$$
\begin{align}
\Delta t &= 1 \\
x(0) &= 1 \\
\dot{x}(0) &= 1
\end{align}
$$


```python
x0 = 1
Δt = 1 

# Para escribir simbolos griegos como Δ, tan solo tienes que escribir su nombre
# precedido de una diagonal (\Delta) y teclear tabulador una vez

F = lambda x : -x
x1 = x0 + F(x0)*Δt
x1
```


```python
x2 = x1 + F(x1)*Δt
x2
```

### Ejercicio

Crea codigo para una iteración mas con estos mismos parametros y despliega el resultado.


```python
x3 = # Escribe el codigo de tus calculos aqui
```


```python
from pruebas_2 import prueba_2_1
prueba_2_1(x0, x1, x2, x3, _)
```

Momento... que esta pasando? Resulta que este $\Delta t$ es demasiado grande, intentemos con 20 iteraciones:

$$
\begin{align}
\Delta t &= 0.5 \\
x(0) &= 1
\end{align}
$$


```python
x0 = 1
n = 20
Δt = 10/n

F = lambda x : -x
x1 = x0 + F(x0)*Δt
x1
```


```python
x2 = x1 + F(x1)*Δt
x2
```


```python
x3 = x2 + F(x2)*Δt
x3
```

Esto va a ser tardado, mejor digamosle a Python que es lo que tenemos que hacer, y que no nos moleste hasta que acabe, podemos usar un ciclo ```for``` y una lista para guardar todos los valores de la trayectoria:


```python
xs = [x0]
for t in range(20):
    xs.append(xs[-1] + F(xs[-1])*Δt)

xs
```

Ahora que tenemos estos valores, podemos graficar el comportamiento de este sistema, primero importamos la libreria ```matplotlib```:


```python
%matplotlib inline
from matplotlib.pyplot import plot
```

Mandamos a llamar la función ```plot```:


```python
plot(xs);
```

Sin embargo debido a que el periodo de integración que utilizamos es demasiado grande, la solución es bastante inexacta, podemos verlo al graficar contra la que sabemos es la solución de nuestro problema:


```python
from numpy import linspace, exp
```


```python
ts = linspace(0, 10, 20)
plot(xs)
plot(exp(-ts));
```

Si ahora utilizamos un numero de pedazos muy grande, podemos mejorar nuestra aproximación:


```python
xs = [x0]
n = 100
Δt = 10/n

for t in range(100):
    xs.append(xs[-1] + F(xs[-1])*Δt)
    
ts = linspace(0, 10, 100)
```


```python
plot(xs)
plot(exp(-ts));
```

## odeint

Este método funciona tan bien, que ya viene programado dentro de la libreria ```scipy```, por lo que solo tenemos que importar esta librería para utilizar este método.

Sin embargo debemos de tener cuidado al declarar la función $F(x, t)$. El primer argumento de la función se debe de referir al estado de la función, es decir $x$, y el segundo debe de ser la variable independiente, en nuestro caso el tiempo.


```python
from scipy.integrate import odeint
```


```python
F = lambda x, t : -x
```


```python
x0 = 1
ts = linspace(0, 10, 100)
xs = odeint(func=F, y0=x0, t=ts)
```


```python
plot(ts, xs);
```

### Ejercicio

Grafica el comportamiento de la siguiente ecuación diferencial.

$$
\dot{x}  = x^2 - 5 x + \frac{1}{2} \sin{x} - 2
$$

> Nota: Asegurate de impotar todas las librerias que puedas necesitar


```python
ts = # Escribe aqui el codigo que genera un arreglo de puntos equidistantes (linspace)
x0 = # Escribe el valor de la condicion inicial

# Importa las funciones de librerias que necesites aqui

G = lambda x, t: # Escribe aqui el codigo que describe los calculos que debe hacer la funcion

xs = # Escribe aqui el comando necesario para simular la ecuación diferencial

plot(ts, xs);
```


```python
from pruebas_2 import prueba_2_2
prueba_2_2(ts, xs)
```

## Sympy

Y por ultimo, hay veces en las que incluso podemos obtener una solución analítica de una ecuación diferencial, siempre y cuando cumpla ciertas condiciones de simplicidad.


```python
from sympy import var, Function, dsolve
from sympy.physics.mechanics import mlatex, mechanics_printing
mechanics_printing()
```


```python
var("t")
```


```python
x = Function("x")(t)
x, x.diff(t)
```


```python
solucion = dsolve(x.diff(t) + x, x)
solucion
```

### Ejercicio

Implementa el codigo necesario para obtener la solución analítica de la siguiente ecuación diferencial:

$$
\dot{x} = x^2 - 5x
$$


```python
# Declara la variable independiente de la ecuación diferencial
var("")

# Declara la variable dependiente de la ecuación diferencial
 = Function("")()

# Escribe la ecuación diferencial con el formato necesario (Ecuacion = 0)
# adentro de la función dsolve
sol = dsolve()
sol
```


```python
from pruebas_2 import prueba_2_3
prueba_2_3(sol)
```

## Solución a ecuaciones diferenciales de orden superior

Si ahora queremos obtener el comportamiento de una ecuacion diferencial de orden superior, como:

$$
\ddot{x} = -\dot{x} - x + 1
$$

Tenemos que convertirla en una ecuación diferencial de primer orden para poder resolverla numericamente, por lo que necesitaremos convertirla en una ecuación diferencial matricial, por lo que empezamos escribiendola junto con la identidad $\dot{x} = \dot{x}$ en un sistema de ecuaciones:

$$
\begin{align}
\dot{x} &= \dot{x} \\
\ddot{x} &= -\dot{x} - x + 1
\end{align}
$$

Si extraemos el operador derivada del lado izquierda, tenemos:

$$
$$
\begin{align}
\frac{d}{dt} x &= \dot{x} \\
\frac{d}{dt} \dot{x} &= -\dot{x} - x + 1
\end{align}
$$
$$

O bien, de manera matricial:

$$
\frac{d}{dt}
\begin{pmatrix}
x \\
\dot{x}
\end{pmatrix} =
\begin{pmatrix}
0 & 1 \\
-1 & -1
\end{pmatrix}
\begin{pmatrix}
x \\
\dot{x}
\end{pmatrix} +
\begin{pmatrix}
0 \\
1
\end{pmatrix}
$$

Esta ecuación ya _no_ es de segundo orden, es de hecho, de primer orden, sin embargo nuestra variable ha crecido a ser un vector de estados, por el momento le llamaremos $X$, asi pues, lo podemos escribir como:

$$
\frac{d}{dt} X = A X + B
$$

en donde:

$$
A = \begin{pmatrix}
0 & 1 \\
-1 & -1
\end{pmatrix} \quad \text{y} \quad B =
\begin{pmatrix}
0 \\
1
\end{pmatrix}
$$

y de manera similar, declarar una función para dar a ```odeint```.


```python
from numpy import matrix, array
```


```python
def F(X, t):
    A = matrix([[0, 1], [-1, -1]])
    B = matrix([[0], [1]])
    return array((A*matrix(X).T + B).T).tolist()[0]
```


```python
ts = linspace(0, 10, 100)
xs = odeint(func=F, y0=[0, 0], t=ts)
```


```python
plot(xs);
```

### Ejercicio

Implementa la solución de la siguiente ecuación diferencial, por medio de un modelo en representación de espacio de estados:

$$
\ddot{x} = -8\dot{x} - 15x + 1
$$

> Nota: Tomalo con calma y paso a paso
> * Empieza anotando la ecuación diferencial en tu cuaderno, junto a la misma identidad del ejemplo
> * Extrae la derivada del lado izquierdo, para que obtengas el _estado_ de tu sistema
> * Extrae las matrices A y B que corresponden a este sistema
> * Escribe el codigo necesario para representar estas matrices


```python
def G(X, t):
    A = # Escribe aqui el codigo para la matriz A
    B = # Escribe aqui el codigo para el vector B
    return array((A*matrix(X).T + B).T).tolist()[0]

ts = linspace(0, 10, 100) 
xs = odeint(func=G, y0=[0, 0], t=ts)
```


```python
plot(xs);
```


```python
from pruebas_2 import prueba_2_4
prueba_2_4(xs)
```

## Funciones de transferencia

Sin embargo, no es la manera mas facil de obtener la solución, tambien podemos aplicar una transformada de Laplace, y aplicar las funciones de la libreria de control para simular la función de transferencia de esta ecuación; al aplicar la transformada de Laplace, obtendremos:

$$
G(s) = \frac{1}{s^2 + s + 1}
$$


```python
from control import tf, step
```


```python
F = tf([0, 0, 1], [1, 1, 1])
```


```python
xs, ts = step(F)
```


```python
plot(ts, xs);
```

### Ejercicio

Modela matematicamente la ecuación diferencial del ejercicio anterior, usando una representación de función de transferencia.

> Nota: De nuevo, no desesperes, escribe tu ecuación diferencial y aplica la transformada de Laplaca tal como te enseñaron tus abuelos hace tantos años...


```python
G = tf([], []) # Escribe los coeficientes de la función de transferencia 

xs, ts = step(G)
plot(ts, xs);
```


```python
from pruebas_2 import prueba_2_5
prueba_2_5(ts, xs)
```

## Problemas

1. Modela matematicamente la suspensión de un automovil de masa $m = 1200 kg$, considera que los resortes de su suspensión tienen una constante $k = 15,000 \frac{N}{m}$ y un amortiguador con constante $c = 1,500 \frac{N s}{m}$.

2. Gráfica el comportamiento del sistema ante una fuerza $F = 1 N$.
3. Que tipo de comportamiento presenta este sistema? Estable, criticamente estable, inestable?

e790169bb82900e9e60477002a549f13249f7f1d
ipynb
Jupyter Notebook
Practicas/.ipynb_checkpoints/Practica 2 - Solucion de ecuaciones diferenciales-checkpoint.ipynb
robblack007/clase-dinamica-robot
f38cb358f2681e9c0dce979acbdcd81bf63bd59c
MIT
Practicas/.ipynb_checkpoints/Practica 2 - Solucion de ecuaciones diferenciales-checkpoint.ipynb
robblack007/clase-dinamica-robot
f38cb358f2681e9c0dce979acbdcd81bf63bd59c
MIT
2016-01-26T18:33:11.000Z
2016-05-30T23:58:07.000Z
Practicas/.ipynb_checkpoints/Practica 2 - Solucion de ecuaciones diferenciales-checkpoint.ipynb
robblack007/clase-dinamica-robot
f38cb358f2681e9c0dce979acbdcd81bf63bd59c
MIT
Qwen/Qwen-72B
1. YES
2. YES