scipy.special.

roots_hermite#

scipy.special.roots_hermite(n, mu=False)[source]#

Gauss-Hermite (physicist’s) quadrature.

Compute the sample points and weights for Gauss-Hermite quadrature. The sample points are the roots of the nth degree Hermite polynomial, \(H_n(x)\). These sample points and weights correctly integrate polynomials of degree \(2n - 1\) or less over the interval \([-\infty, \infty]\) with weight function \(w(x) = e^{-x^2}\). See 22.2.14 in [AS] for details.

Parameters:

nint: quadrature order
mubool, optional: If True, return the sum of the weights, optional.

Returns:

xndarray: Sample points
wndarray: Weights
mufloat: Sum of the weights

See also

scipy.integrate.fixed_quad
numpy.polynomial.hermite.hermgauss
roots_hermitenorm

Notes

For small n up to 150 a modified version of the Golub-Welsch algorithm is used. Nodes are computed from the eigenvalue problem and improved by one step of a Newton iteration. The weights are computed from the well-known analytical formula.

For n larger than 150 an optimal asymptotic algorithm is applied which computes nodes and weights in a numerically stable manner. The algorithm has linear runtime making computation for very large n (several thousand or more) feasible.

References

[townsend.trogdon.olver-2014]

Townsend, A. and Trogdon, T. and Olver, S. (2014) Fast computation of Gauss quadrature nodes and weights on the whole real line. arXiv:1410.5286.

[townsend.trogdon.olver-2015]

Townsend, A. and Trogdon, T. and Olver, S. (2015) Fast computation of Gauss quadrature nodes and weights on the whole real line. IMA Journal of Numerical Analysis DOI:10.1093/imanum/drv002.

[AS]

Milton Abramowitz and Irene A. Stegun, eds. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. New York: Dover, 1972.

Examples

>>> from scipy.special import roots_hermite
>>> roots, weights = roots_hermite(5)
>>> roots
array([-2.02018287, -0.95857246,  0.        ,  0.95857246,  2.02018287])
>>> weights
array([0.01995324, 0.39361932, 0.94530872, 0.39361932, 0.01995324])

Verify that the values in roots are roots of the Hermite polynomial \(H_5(x)\).

>>> from scipy.special import eval_hermite
>>> eval_hermite(5, roots)
array([ 2.63775533e-13,  5.02429587e-15,  0.00000000e+00, -5.02429587e-15,
       -2.63775533e-13])

The values of \(H_5(x)\) evaluated at the roots are indeed rather close to zero. The increasing values for larger roots can be explained by the increasing derivative of the Hermite polynomial at the roots.

Verify that the sum of the weights equals the integral from \(-\infty\) to \(\infty\) of \(\exp(-x^2)\) which evaluates to \(\sqrt{\pi}\). There are two ways to obtain the sum of weights, both yielding the expected result within numerical precision.

>>> sum(weights)
np.float64(1.7724538509055157)
>>> roots, weights, sum_of_weights = roots_hermite(5, mu=True)
>>> sum_of_weights
np.float64(1.7724538509055159)

>>> from math import exp, pi, sqrt
>>> sqrt(pi)
1.7724538509055159

Roots and weights obtained from the Hermite polynomial \(H_n(x)\) are used in Gauss-Hermite quadrature where the integral from \(-\infty\) to \(\infty\) of \(f(x)\exp(-x^2)\) is evaluated. Roots and weights for order \(n\) are expected to yield the exact result for polynomials \(f(x)\) of a maximal order of \(2n-1\).

>>> f = lambda x: x**4
>>> weights @ f(roots)
np.float64(1.3293403881791352)

This result is indeed very close to the exact value of \(3\sqrt{\pi}/4\).

>>> 0.75*sqrt(pi)
1.329340388179137

In general, Gauss-Hermite quadrature will only yield an approximate value of the integral. Consider the integral from \(-\infty\) to \(\infty\) of \(\cos(x)\exp(-x^2)\) which evaluates to \(\sqrt{\pi}\exp(-1/4)\).

>>> import numpy as np
>>> weights @ np.cos(roots)
np.float64(1.3803900759356564)
>>> sqrt(pi) * exp(-0.25)
1.380388447043143

In order to improve the accuracy obtained from Gauss-Hermite quadrature, a large number of nodes can be chosen.

>>> roots, weights = roots_hermite(50)
>>> weights @ np.cos(roots)
np.float64(1.380388447043143)