Sorry for the basic question, but I was playing with Basix to approximate the integral of a toy 1D function with a normalizing constant of 1.5, but was getting 69.7 with Basix compared to 1.4998 with numpy.polynomial.legendre.leggauss()
. Below is the code - advice would be greatly appreciated.
import numpy as np
import basix
from basix import CellType
def test_1d_quadrature(library='basix'):
def f(xs):
assert np.all(xs <= 1.) and np.all(xs >= 0.)
ys = np.array([2 * x + 1. if x < 0.5 else -2 * x + 3. for x in xs])
return ys
if library == 'basix':
quad_points, quad_weights = basix.make_quadrature(CellType.interval, 100)
elif library == 'numpy':
quadrature_degree = 50
quad_points, quad_weights = np.polynomial.legendre.leggauss(quadrature_degree)
quad_points = (quad_points + 1) / 2 # Scale & shift [-1, 1] to [0, 1]
quad_weights = quad_weights / 2 # Divide the quadrature weights by 2 since volume halved.
integrate_fn = lambda f: np.sum(quad_weights * f(quad_points))
normalizing_constant = integrate_fn(f)
print('Actual normalizing constant: 1.5')
print('Quadrature normalizing constant: {}'.format(normalizing_constant))