%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
Para integración numérica e integración de ecuaciones diferenciales, SciPy proporciona el paquete integrate
.
from scipy import integrate
Por ejemplo, supongamos que queremos integrar esta función:
Debemos definir una función en Python cuyo argumento sea la variable de integración:
def f(x):
return np.exp(-x ** 2)
Y ahora simplemente utilizamos la función quad
, que recibe como argumentos la función y los límites de integración.
integrate.quad(f, 0, 5)
(0.8862269254513955, 2.3183115159980698e-14)
Se devuelve el resultado de la integración y una estimación del error cometido.
Puede darse el caso en que nuestra función dependa de un cierto número de parámetros:
En este caso, debemos seguir respetando que el primer argumento es la variable de integración, pero a partir de ahí podemos añadir todos los argumentos que queramos:
def f(x, A, B):
return A * np.exp(-B * x ** 2)
A la hora de integrar esta función, tendremos que usar el parámetro args
para dar valores los argumentos extra de la función:
integrate.quad(f, 0, 5, args=(1.0, 1.0))
(0.8862269254513955, 2.3183115159980698e-14)
Dentro del paquete integrate
tenemos también funciones para resolver ecuaciones diferenciales ordinarias (EDOs), como es el caso de la función odeint
. Esta función puede resolver cualquier sistema del tipo:
Por ejemplo, supongamos que queremos resolver la ecuación diferencial:
Hemos de definir la función del sistema, las condiciones iniciales y el vector de tiempos donde realizaremos la integración.
def f(y, t):
return -y
y0 = 1.0
t = np.linspace(0, 3)
sol = integrate.odeint(f, y0, t)
plt.plot(t, sol)
[<matplotlib.lines.Line2D at 0x7ff8fee2dcd0>]
Y se obtiene la solución esperada: un decaimiento exponencial.
Para resolver ecuaciones de orden superior, habrá que realizar una reducción de orden:
def f(y, t):
return np.array([y[1], -y[0]])
t = np.linspace(0, 10)
y0 = np.array([1.0, 0.0])
sol = integrate.odeint(f, y0, t)
plt.plot(t, sol[:, 0], label='$y$')
plt.plot(t, sol[:, 1], '--k', label='$\dot{y}$')
plt.legend()
<matplotlib.legend.Legend at 0x7ff8fed267d0>
El paquete optimize
contiene unas cuantas funciones para optimización con y sin restricciones, minimización de funciones, múltiples algoritmos... Ahora nos vamos a concentrar en la búsqueda de ceros de funciones, utilizando la función root
.
from scipy import optimize
Por ejemplo, esta sería la función necesaria para resolver la ecuación de Kepler:
def F(E, e, M):
return E - e * np.sin(E) - M
La función root
recibe como argumento la función de la ecuación, la condición inicial y posibles argumentos extra de la función.
sol = optimize.root(F, 0.1, args=(0.016, 0.1))
sol.x
array([ 0.10162317])
El paquete optimize
trae también funciones para hacer ajustes. Para ver este tema explicado con detalle puedes leer el artículo Ajuste e interpolación unidimensionales básicos con SciPy y el notebook incluido en los materiales del curso.