Diferentes intervalos para la cuadratura de Gauss-Legendre en números.

¿Cómo podemos usar el paquete NumPy numpy.polynomial.legendre.leggauss en intervalos distintos a [-1, 1] ?


El siguiente ejemplo compara scipy.integrate.quad con el método de Gauss-Legendre en el intervalo [-1, 1] .

 import numpy as np from scipy import integrate # Define function and interval a = -1. b = 1. f = lambda x: np.cos(x) # Gauss-Legendre (default interval is [-1, 1]) deg = 6 x, w = np.polynomial.legendre.leggauss(deg) gauss = sum(w * f(x)) # For comparison quad, quad_err = integrate.quad(f, a, b) print 'The QUADPACK solution: {0:.12} with error: {1:.12}'.format(quad, quad_err) print 'Gauss-Legendre solution: {0:.12}'.format(gauss) print 'Difference between QUADPACK and Gauss-Legendre: ', abs(gauss - quad) 

Salida:

 The QUADPACK solution: 1.68294196962 with error: 1.86844092378e-14 Gauss-Legendre solution: 1.68294196961 Difference between QUADPACK and Gauss-Legendre: 1.51301193796e-12 

Para cambiar el intervalo , traduzca los valores de x de [-1, 1] a [a, b] usando, por ejemplo,

 t = 0.5*(x + 1)*(b - a) + a 

y luego escalar la fórmula de cuadratura por (b – a) / 2:

 gauss = sum(w * f(t)) * 0.5*(b - a) 

Aquí hay una versión modificada de su ejemplo:

 import numpy as np from scipy import integrate # Define function and interval a = 0.0 b = np.pi/2 f = lambda x: np.cos(x) # Gauss-Legendre (default interval is [-1, 1]) deg = 6 x, w = np.polynomial.legendre.leggauss(deg) # Translate x values from the interval [-1, 1] to [a, b] t = 0.5*(x + 1)*(b - a) + a gauss = sum(w * f(t)) * 0.5*(b - a) # For comparison quad, quad_err = integrate.quad(f, a, b) print 'The QUADPACK solution: {0:.12} with error: {1:.12}'.format(quad, quad_err) print 'Gauss-Legendre solution: {0:.12}'.format(gauss) print 'Difference between QUADPACK and Gauss-Legendre: ', abs(gauss - quad) 

Se imprime:

 La solución QUADPACK: 1.0 con error: 1.11022302463e-14
 Solución de Gauss-Legendre: 1.0
 Diferencia entre QUADPACK y Gauss-Legendre: 4.62963001269e-14

quadpy (un pequeño proyecto mío) como una syntax más simple para esto:

 import numpy import quadpy out = quadpy.line_segment.integrate( numpy.cos, [1.1, 1.2], # the interval quadpy.line_segment.GaussLegendre(4) )