Intervalli diversi per la quadratura di Gauss-Legendre in numpy

Question

Come utilizzare il pacchetto NumPy numpy.polynomial.legendre.leggauss su intervalli diversi da [-1, 1]?

---

Il seguente esempio confronta scipy.integrate.quad al metodo di Gauss-Legendre sull'intervallo [-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) 

uscita:

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

Answer 1

Per change the interval, tradurre i valori x da [-1, 1] a [a, b] utilizzando, per esempio,

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

e poi scala la formula di quadratura di (b - a)/2:

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

Ecco una versione modificata del vostro esempio:

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) 

Esso stampa:

 
The QUADPACK solution: 1.0 with error: 1.11022302463e-14 
Gauss-Legendre solution: 1.0 
Difference between QUADPACK and Gauss-Legendre: 4.62963001269e-14 

Answer 2

quadpy (un piccolo progetto di miniera) come una sintassi più semplice per questo:

import numpy 
import quadpy 

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