Taller Python Científico - PyConES 2014

Parte III: SciPy

Juan Luis Cano, Kiko Correoso y Álex Sáez (@Pybonacci)

Fran Navarro, Zuria Bauer y Daniel Domene (@CAChemEorg)

Licencia Creative Commons
<span xmlns:dct="http://purl.org/dc/terms/" property="dct:title">Taller de </span> por <span xmlns:cc="http://creativecommons.org/ns#" property="cc:attributionName">Pybonacci y CAChemE</span> se distribuye bajo una Licencia Creative Commons Atribución 4.0 Internacional.

¿Qué es SciPy?

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.

Integración mediante métodos numéricos con SciPy/Python

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.

In [1]:
import numpy as np
import matplotlib.pyplot as plt

%matplotlib inline
In [2]:
from scipy import integrate

Cuadratura de Gauss

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:

In [3]:
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)$$
In [4]:
def fun(x):
    return x * np.sin(x)

Integraremos la función en el intervalo $[2, 9]$: $$ \int_{2}^{9} x sin(x) dx$$

In [5]:
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.

In [6]:
# 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()

De forma similar, SciPy dispone de funciones de integración numérica para integrales dobles y triples (dblquad y tplquad , respectivamente).

Ecuaciones diferenciales ordinarias

Para integrar EDOs vamos a usar la función odeint del paquete integrate, que permite integrar sistemas del tipo:

$$ \frac{d\mathbf{y}}{dt}=\mathbf{f}\left(\mathbf{y},t\right)$$

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:

In [7]:
from scipy.integrate import odeint
$$\frac{dy}{dt} + y = 0$$$$\frac{dy}{dt} = -y$$
In [8]:
def f(y, t):
    return -y
**¡Importante!**: Fíjate que la función del sistema recibe como primer argumento $\mathbf{y}$ (un array) y como segundo argumento el instante $t$ (un escalar). Esta convención va exactamente al revés que en MATLAB y nos equivocamos obtendremos errores o, lo que es peor, resultados incorrectos.

Condiciones iniciales:

In [9]:
y0 = 1

Vector de tiempos donde realizamos la integración:

In [10]:
t = np.linspace(0, 3)

Integramos y representamos la solución:

In [11]:
sol = odeint(f, y0, t)
plt.plot(t, sol)
Out[11]:
[<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.