Taller de por Pybonacci y CAChemE se distribuye bajo una Licencia Creative Commons Atribución 4.0 Internacional.
SciPy es un conjunto de paquetes open-source para computación científica general (matématicas, ciencia e ingeniería). La mayoría hacen uso NumPy (biblioteca de Python para matrices multidimensionales), otros ofrecen interfaces a programas muy conocidos y utilizados escritos en Fortran, o C++. SciPy ofrece así un gran número de algoritmos científicos de alto nivel, algunos de los temas cubiertos son:
Antes de ponernos a programar, vale la pena comprobar si la funcionalidad que deseamos utilizar no está ya implementada en SciPy. Como programadores no profesionales, los científicos a menudo tendemos a reinventar la rueda, lo que conduce a errores, código no optimizado, difícil de compartir e imposible de mantener. Por el contrario, los programas de SciPy están optimizados y probados, y por lo tanto deben ser utilizados cuando sea posible.
Como siempre, para hacer el uso de estas funciones, importaremos lo paquetes que vayamos a necesitar:
import scipy as ...
Bueno, en realidad SciPy
nos ofrece muchas funciones y rara vez querremos usarlo todo en una misma sesión. Típicamente usaremos:
from scipy import ...
Debido a limitación de tiempo en este taller, vamos a centrarnos solo en el paquete de integración. Si quieres, puedes darle un vistazo rápido a Scientific Lectures: SciPy, high-level scientific computing donde puedes encontrar ejemplos más variados.
SciPy proporciona el paquete integrate
que proporciona algunas técnicas de integración tanto de funciones como de ecuaciones diferenciales. En primer lugar importémos el paquete y el resto de librerías que vayamos a necesitar.
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from scipy import integrate
Como se puede ver en la ayuda, si queremos realizar una integración numérica de una función de una variable, podemos utilizar el método de cuadratura de Gauss quad
(aunque también podemos usar trapz
, simps
...). La forma de acceder a ella tal y como hemos importado el paquete sería ejecutando integrate.quad
. Sin embargo, es más común importarla del siguiete modo:
from scipy.integrate import quad
¿Qué es lo primero que necesitamos hacer para integrar una función? Pues sí, la función... definamos una:
$$f(x) = x sin(x)$$def fun(x):
return x * np.sin(x)
Integraremos la función en el intervalo $[2, 9]$: $$ \int_{2}^{9} x sin(x) dx$$
x_inferior = 2
x_superior = 9
value, err = quad(fun, x_inferior, x_superior)
print('El resultado es: %f con un error de %e' %(value, err))
El resultado es: 6.870700 con un error de 2.864870e-13
Si lo deseamos, podemos hacer uso de lo que acabamos de aprender en la sección previa de este taller y representar la función e intervalo de integración que hemos tomado. Fíjate en como hemos hecho uso de $\LaTeX$ en el título del gráfico.
# Vector con valores de la variable independiente
x = np.linspace(0,10,100)
# Hacemos uso de la función para obtener los valores
# de la variable dependiente
y = fun(x)
# Representamos la figura
plt.plot(x,y, linewidth = 2)
# Añadimos un título
plt.title(u'$y = x sin(x)$', fontsize = 25)
# Creamos el intervalo para sombrear el área de integración
x_fill = np.linspace(2,9,100)
# Obtenemos los puntos sobre la curva
y_fill = fun(x_fill)
# Dibujamos el área sombreada
plt.fill_between(x_fill,y_fill, color = 'gray', alpha = 0.5)
# Representamos los resultados con una gradilla
plt.grid()
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}$.
SciPy ofrece dos maneras diferentes para resolver EDOs: Una API basada en la función odeint
y otra orientada a objetos basados en la clase ode
. Por lo general, odeint
es más sencilla de utilizar, por lo que la usaremos en este taller, pero la clase ode
ofrece un cierto nivel de control más preciso.
Para poder utilizar odeint
, debemos de importarlo desde el módulo scipy.integrate
y definir la función diferencial a integrar:
from scipy.integrate import odeint
def f(y, t):
return -y
Condiciones iniciales:
y0 = 1
Vector de tiempos donde realizamos la integración:
t = np.linspace(0, 3)
Integramos y representamos la solución:
sol = odeint(f, y0, t)
plt.plot(t, sol)
[<matplotlib.lines.Line2D at 0x7ff5e10>]
Tal y como podíamos esperar, la resolución de esta ecuación diferencial devuelve una exponencial negativa.
En el siguiente notebook, veremos un ejemplo aplicado de la resolución de un sistema de EDOs para modelar el comportamiento de una epidemia.