Perhaps a little justification for the first point above is in order (pun not intended). Consider the standard problem of integrating f(x) = sin(x) over [-1, 1] using Gauss-Legendre quadrature. The standard L_2 norm error is plotted on the vertical axes, while the order is shown on the horizontal axes in the figure below.

Part of the issue here is down to my original sloppy implementation of the Golub-Welsch algorithm: using numpy's eigenvalue solver. If we replace this with scipy's we observe better results.

The amended lines of code in Parameter.py are given below:

def get_jacobi_eigenvectors(self, order=None):
"""
Computes the eigenvectors of the Jacobi matrix.
:param Parameter self:
An instance of the Parameter object.
:param int order:
Order of the recurrence coefficients.
"""
if order is None:
order = self.order + 1
JacobiMat = self.get_jacobi_matrix(order)
if order == 1:
V = [1.0]
else:
D, V = sc.linalg.eigh(JacobiMat)
idx = D.argsort()[::-1]
eigs = D[idx]
eigVecs = V[:, idx]
return eigVecs
def get_local_quadrature(self, order=None, ab=None):
if order is None:
order = self.order + 1
else:
order = order + 1
if ab is None:
# Get the recurrence coefficients & the jacobi matrix
JacobiMat = self.get_jacobi_matrix(order)
ab = self.get_recurrence_coefficients(order+1)
else:
ab = ab[0:order+1,:]
JacobiMat = self.get_jacobi_matrix(order, ab)
# If statement to handle the case where order = 1
if order == 1:
# Check to see whether upper and lower bound are defined:
if not self.lower or not self.upper:
p = np.asarray(self.distribution.mean).reshape((1,1))
else:
p = np.asarray((self.upper - self.lower)/(2.0) + self.lower).reshape((1,1))
w = [1.0]
else:
D, V = sc.linalg.eigh(JacobiMat)
idx = D.argsort()[::-1]
eigs = D[idx]
eigVecs = V[:, idx]
w = np.linspace(1,order+1,order) # create space for weights
p = np.ones((int(order),1))
for u in range(0, len(idx) ):
w[u] = float(ab[0,1]) * (eigVecs[0,idx[u]]**2) # replace weights with right value
p[u,0] = eigs[u]
return p, w

But this too is far from perfect! However, as pointed out in Hale and Townsend 2013, the only way to improve the accuracy with increasing orders is to address the instability that arises when using the standard Golub-Welsch algorithm.

The code to generate the first figure above can be found on EQ-live. The new feature pull request can be found here: