9.4: Cuadraturas Gaussianas

Anteriormente, hemos asumido en nuestro análisis de esquemas de integración numérica que los puntos de discretización están igualmente espaciados. Esta suposición no es estrictamente necesaria; por ejemplo, podemos modificar fácilmente las fórmulas de regla de punto medio y regla de trapecio para trabajar con puntos no equidistantes.

Sin embargo, si somos libres de elegir los puntos de discretización para realizar la integración numérica, y estos puntos no necesitan estar igualmente espaciados, entonces es posible explotar esta libertad para mejorar aún más la precisión de la integración numérica. La idea es optimizar la colocación de los puntos de discretización, para minimizar el error numérico resultante. Esta es la idea básica detrás del método de integración por cuadraturas gaussianas.

No discutiremos los detalles de este método de integración numérica. Para usarlo, puedes llamar a la función scipy.integrate.quad. Esta función utiliza una biblioteca numérica de bajo nivel llamada QUADPACK, que realiza la integración en cuadratura con cuadraturas adaptativas, es decir, determina automáticamente cuántos puntos de discretización deben usarse, y dónde deben ubicarse, para producir un resultado con el precisión numérica deseada.

Debido a que QUADPACK calcula los puntos de discretización por sí mismo, hay que pasar quad una función que represente el integrando, en lugar de una matriz de valores integrand como con trapz o simps. La forma estándar de llamar a la función es

t = quad(f, a, b)

que calcula la integral

\[t = \int_a^b f(x) dx.\]

El valor de retorno es una tupla de la forma (t, err), donde t es el valor de la integral y err es una estimación del error numérico. La función quad también acepta muchas otras entradas opcionales, que pueden usarse para especificar entradas adicionales para pasar a la función integrand, la tolerancia a errores, el número de subintervalos a usar, etc. Consulte la documentación para más detalles.

La elección de si realizar la integración numérica usando la regla de Simpson (simps) o cuadraturas gaussianas (quad) es situacional. Si ya conoces los valores de los integrands en un conjunto predeterminado de puntos de discretización (por ejemplo, a partir del resultado de un cálculo de diferencia finita), entonces usa simps. Si puedes definir una función que pueda calcular rápidamente el valor del integrando en cualquier punto, usa quad.

Aquí hay un ejemplo del uso de quad para calcular la integral\(\int_{0}^\infty \frac{dx}{x^2 + 1}\):

from scipy import *
from scipy.integrate import quad

def f(x):
    return 1.0 / (x*x+1)

integ, _ = quad(f, 0, 1000)
print(integ)

(Tenga en cuenta que quad devuelve dos valores; el primero es el valor calculado de la integral, y el otro es el error absoluto, que no nos interesa, por lo que lo tiramos en la variable “descartable” llamada _. Consulte la documentación para más detalles.) Al ejecutar el programa anterior se imprime el resultado\(1.569796\dots,\) que difiere del resultado exacto,\(\pi/2 = 1.570796\dots,\) por un error relativo de\(0.06\%\).

Aquí hay otro ejemplo, donde el integrando toma un parámetro adicional:\(\int_0^\infty x e^{-\lambda x}\, dx\):

from scipy import *
from scipy.integrate import quad

def f(x, lambd):
    return x * exp(-lambd * x)

integ, _ = quad(f, 0, 100, args=(0.5))
print(integ)

Al ejecutar el programa se imprime el resultado\(4.0\), el cual concuerda con el resultado exacto de\(1/\lambda^2\) para el valor elegido del parámetro\(\lambda =0.5\).