¿Te acuerdas de todos esos esquemas numéricos para integrar ecuaciones diferenciales ordinarias? Es bueno saber que existen y qué peculiaridades tiene cada uno, pero en este curso no queremos implementar esos esquemas: queremos resolver las ecuaciones. Los problemas de evolución están por todas partes en ingeniería y son de los más divertidos de programar.
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
Para integrar EDOs vamos a usar la función odeint
del paquete integrate
, que permite integrar sistemas del tipo:
con condiciones iniciales $\mathbf{y}(\mathbf{0}) = \mathbf{y_0}$.
from scipy.integrate import solve_ivp
Vamos a integrar primero una EDO elemental, cuya solución ya conocemos:
$$y' + y = 0$$$$f(y, t) = \frac{dy}{dt} = -y$$def f(t, y):
return np.array([-y])
Condiciones iniciales:
y0 = np.array([1])
tini = 0
tfin = 3
Integramos y representamos la solución:
sol = solve_ivp(f, (tini, tfin), y0)
plt.plot(sol.t, sol.y[0, :], 'o-')
[<matplotlib.lines.Line2D at 0x7fe3b2671828>]
Pero, ¿cómo se han seleccionado los puntos en los que se calcula la solución? El solver los ha calculado por nosotros. Si queremos tener control sobre estos puntos, podemos pasar de manera explícita el vector de tiempos:
time = np.linspace(tini, tfin, 30)
sol_2 = solve_ivp(f, (tini, tfin), y0, t_eval=time)
plt.plot(sol_2.t, sol_2.y[0, :], 'o-')
[<matplotlib.lines.Line2D at 0x7fe3b260d748>]
Probemos a pintar las dos soluciones anteriores, una encima de la otra:
plt.plot(sol_2.t, sol_2.y[0, :], 'o-')
plt.plot(sol.t, sol.y[0, :], 'o-')
[<matplotlib.lines.Line2D at 0x7fe3b25ea048>]
Podemos observar que a pesar de que en la primera se han usado muchos menos puntos, aquellos en los que se ha calculado la solución coinciden con el segundo resultado. Esto se debe a que, en realidad, el solver siempre da los pasos que considere necesarios para calcular la solución, pero sólo guarda los que nosotros le indicamos. Esto lo podemos ver del siguiente modo:
print(f"function evaluations in sol 1: {sol.nfev}")
print(f"function evaluations in sol 2: {sol_2.nfev}")
function evaluations in sol 1: 32 function evaluations in sol 2: 32
De hecho podemos usar la salida densa para obtener la solución en un punto cualquiera:
sol_3 = solve_ivp(f, (tini, tfin), y0, dense_output=True)
sol_3.sol(1.14567)
array([0.31820121])
t = np.linspace(tini, tfin, 45)
y = sol_3.sol(t)
plt.plot(t, y[0, :], 'o-')
[<matplotlib.lines.Line2D at 0x7fe3b2557c50>]
plt.plot(t, y[0, :], '^-')
plt.plot(sol_2.t, sol_2.y[0, :], 'x-')
plt.plot(sol.t, sol.y[0, :], 'o-')
[<matplotlib.lines.Line2D at 0x7fe3b24be710>]
Tendremos que acordarnos ahora de cómo reducir las ecuaciones de orden. De nuevo, vamos a probar con un ejemplo académico:
def f(t, y):
return np.array([y[1], -y[0]])
t0 = 0
t1 = 10
t = np.linspace(t0, t1)
y0 = np.array([1.0, 0.0])
sol = solve_ivp(f, (t0, t1), y0, t_eval=t)
plt.plot(t, sol.y[0, :], label='$y$')
plt.plot(t, sol.y[1, :], '--k', label='$\dot{y}$')
plt.legend()
<matplotlib.legend.Legend at 0x7fe3b24a5f98>
En nuestra edición anterior del curso de AeroPython puedes ver una aplicación muy interesante de lo que hemos visto hasta ahora al salto de Felix Baumgartner. ¡Aquí lo tienes!
$$\displaystyle m \frac{d^2 y}{d t^2} = -m g + D$$Las siguientes celdas contienen configuración del Notebook
Para visualizar y utlizar los enlaces a Twitter el notebook debe ejecutarse como seguro
File > Trusted Notebook
# Esta celda da el estilo al notebook
from IPython.core.display import HTML
css_file = '../styles/aeropython.css'
HTML(open(css_file, "r").read())