1D quadrature returning wrong result?

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))

Check out the shape of the quadrature arrays:

    print(f(quad_points).shape)
    print(quad_weights.shape)

which yields the following from Basix:

(51, 1)
(51,)

and the following from Numpy:

(50,)
(50,)

Since the arrays quad_weights and f(quad_points) do not have the same shape in the Basix version, they will be broadcast to a 51 \times 51 array when they are multiplied, which is not what you intend. If you flatten f(quad_points) before multiplying, you will obtain the correct result.

2 Likes

Don’t know how I missed that - thank you!

1 Like