Problem in Gaussian Quadrature with Hermite Polynomials

Question

I tried to make a code but it didn't give me what I expected, as I was expecting a value closer to the exact one.

import numpy as np
from scipy.special import roots_hermitenorm

def Gauss_hermite(func: callable, N: int) -> float:
    """
    Numerically computes the integral of a function using Gauss quadrature
    with Hermite polynomials.

    Parameters:
        func (callable): The function to be integrated.
        N (int): The number of intervals for integration.

    Returns:
        float: Approximate value of the integral of 'func' over the interval
            [a, b].
    """

    # Determine zeros and weights using scipy
    xx, ww = roots_hermitenorm(N)
    
    # Compute the integral
    integral = np.sum(ww * func(xx))
    
    return integral

# Example usage
result = Gauss_hermite(lambda x: np.exp(-x**2), 2046)
expected = np.sqrt(np.pi)

print(f"Result:   {result:.5f}")
print(f"Expected: {expected:.5f}")

which gives:

Result: 1.44720
Expected: 1.77245

Answer 1

Based on Gauss-Hermite Quadrature on Wikipedia, it looks like (conceptually) you want something like:

integral = np.sum(ww * func(xx)/np.exp(-xx**2/2))

The quadrature formula is for evaluating the integrand weighted by np.exp(-xx**2/2) (because the SciPy documentation says the polynomials are orthogonal with weight function np.exp(-x**2/2), not np.exp(-x**2)), so you need to undo that weighting.

This gives reasonable results for a lower-order polynomial (e.g. 64), but you'll run into numerical issues for orders like 2048. So really, rather than changing the weights, it's better to change your integrand by dividing it by np.exp(-x**2/2) analytically:

result = Gauss_hermite(lambda x: np.exp(-x**2/2), 2046)

If you have another integrand for which you can't remove the weighting so neatly, there are other tricks you can play to get around the numerical issues, or there are more appropriate quadrature rules to use, but that is a different question.

---

Change the weights:

import numpy as np
from scipy.special import roots_hermitenorm

def Gauss_hermite(func: callable, N: int) -> float:
    xx, ww = roots_hermitenorm(N)
    return np.sum(ww * func(xx)/np.exp(-xx**2/2))

res = Gauss_hermite(lambda x: np.exp(-x**2), 64)
print(res)  # 1.7724538509055154
np.testing.assert_allclose(res, np.sqrt(np.pi))  # passes

Change the integrand:

def Gauss_hermite(func: callable, N: int) -> float:
    xx, ww = roots_hermitenorm(N)
    return np.sum(ww * func(xx))

res = Gauss_hermite(lambda x: np.exp(-x**2/2), 2046)
print(res)  # 1.7724538509055427
np.testing.assert_equal(res, np.sqrt(np.pi))  # passes