# 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 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':
elif library == 'numpy':
quad_points = (quad_points + 1) / 2  # Scale & shift [-1, 1] to [0, 1]

normalizing_constant = integrate_fn(f)

print('Actual normalizing constant: 1.5')


Check out the shape of the quadrature arrays:

    print(f(quad_points).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