Se importa el módulo de integración de scipy
In [31]:
import scipy.integrate as integrate
import numpy as np
Tomemos por ejemplo la función $$f(x) = \frac{1}{x^2+a^2}$$ cuya integral entre $[0,\infty)$ es $\pi/(2a)$
Se define la f(x) como una función en python:
In [7]:
def functionA(x,a):
return 1.0/(x**2+a**2)
Existen varios métodos numéricos para obtener la integral. El más general es .quad, el cual acepta la función, el intervalo y argumentos para la función a integrar.
Evaluemos la integral por todos estos métodos y comparemos con el valor real
In [9]:
a = 1
integrate.quad(functionA, 0, np.inf, args=(a))
Out[9]:
La tupla retornada contiene el resultado de la integración y el error absoluto estimado. El valor exacto de la integral sería:
In [28]:
np.pi/(2*a)
Out[28]:
Algunas veces puede resultar más conveniente utilizar una función lambda para evitar tener que definir una función adicional. El mismo ejemplo anterior usando una función lambda sería así:
In [32]:
integrate.quad(lambda x: 1.0/(x**2+1), 0, np.inf)
Out[32]:
Hay casos en los que se tienen las muestras de la función (no necesariamente la función). En este caso son aplicables los métodos como la regla trapezoidal o la regla de simpson. Por ejemplo, consideremos la integral de $$ y = x^2 $$ Que sobre cualquier intervalo real es $$ \frac{1}{3}( x_f^3 - x_0^3) $$
Se crea un espacio lineal para representar 50 puntos de muestreo en el intervalo $[0,10]$:
In [14]:
x = np.linspace(0,10,50)
x.shape
Out[14]:
In [18]:
x[0], x[-1]
Out[18]:
Se obtienen los valores correspondientes de y
In [20]:
y = x**2
y[0], y[-1]
Out[20]:
Integración por el método trapezoidal
In [25]:
itra = integrate.trapz(y,x)
valorExacto = (10.**3-0**3)/3.
print(itra, valorExacto)
In [26]:
isim = integrate.simps(y,x)
print(isim, valorExacto)
Observamos que la regla de simpson arroja un error mucho menor
Módulo de integración de scipi https://docs.scipy.org/doc/scipy/reference/tutorial/integrate.html