2015-10-31 13 views
5

Jak korzystać z pakietu NumPy numpy.polynomial.legendre.leggauss w odstępach innych niż [-1, 1]?Różne przedziały dla kwadratury Gaussa-Legendre w numpy


Poniższy przykład porównuje scipy.integrate.quad metody Gaussa-Legendre'a w przedziale [-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) 

wyjściowa:

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

Odpowiedz

5

do change the interval, tłumaczyć wartości X od [1, 1] do [a, b] za pomocą, na przykład,

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

a następnie skalę formuła kwadratury według (b - a)/2:

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

Oto zmodyfikowana wersja swoim przykładzie:

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) 

Drukuje:

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

quadpy (trochę projekt kopalni) jako prostszą składnię tego:

import numpy 
import quadpy 

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