Synopsis

#include <boost/math/quadrature/gauss_hermite.hpp>

namespace boost{ namespace math{ namespace quadrature{

template <class Real, unsigned Points, class Policy = boost::math::policies::policy<> >
struct gauss_hermite
{
   static const RandomAccessContainer& abscissa();
   static const RandomAccessContainer& weights();

   template <class F>
   static auto integrate(F f, Real* pL1 = nullptr)->decltype(std::declval<F>()(std::declval<Real>()));

   template <class F, class Norm>
   static auto integrate(F f, const decltype(f(Real(0)))& zero, Norm norm,
                         Real* pL1 = nullptr) -> decltype(f(Real(0)));
};

}}} // namespaces

Description

The gauss_hermite class template computes

using Points evaluations of f, where the abscissas /x_i/ are the zeros of the Hermite polynomial /H_Points. The rule is exact when /f is a polynomial of degree at most 2/Points/ - 1. Note that the weight function exp(-x²) is not part of f: pass only the factor that multiplies it.

Like gauss, this is a "one shot" non-adaptive routine. It obtains no estimate of the error, and is intended for smooth f that grow no faster than a polynomial, such as expectations of smooth functions of a normally distributed variable.

static const RandomAccessContainer& abscissa();
static const RandomAccessContainer& weights();

These provide direct access to the abscissas and weights. Only the nonnegative abscissas are stored, in increasing order; the rule is symmetric about zero.

template <class F>
static auto integrate(F f, Real* pL1 = nullptr)->decltype(std::declval<F>()(std::declval<Real>()));

Integrates f(x) exp(-x²) over the real line, and optionally sets *pL1 to the corresponding integral of |f|. If that is substantially larger than the absolute value of the result, then the sum was ill-conditioned. f may return a real or complex value.

template <class F, class Norm>
static auto integrate(F f, const decltype(f(Real(0)))& zero, Norm norm,
                      Real* pL1 = nullptr) -> decltype(f(Real(0)));

As for gauss, this overload handles vector- and matrix-valued integrands: zero is the additive identity of the returned type, and norm returns a nonnegative scalar convertible to Real.

Example

If X is normally distributed with mean μ and standard deviation σ, then substituting X = μ + √2 σ x gives

so, for example, E[cos X] = exp(-σ²/2) cos μ is computed by

using boost::math::quadrature::gauss_hermite;
using boost::math::constants::root_pi;
using boost::math::constants::root_two;
double mu = 0.5;
double sigma = 2;
auto g = [&](double x) { return std::cos(mu + root_two<double>() * sigma * x); };
double expectation = gauss_hermite<double, 30>::integrate(g) / root_pi<double>();

Choosing the number of points

gauss_hermite has pre-computed tables of abscissas and weights for 7, 10 and 15 points, at up to 110 decimal digits of precision. For other numbers of points, or higher precision types, the abscissas and weights are computed the first time they are needed and then cached. Each zero is isolated by bisection with a Sturm sequence, then refined with Halley's method, all evaluated with the three-term recurrence for orthonormal Hermite polynomials. In double precision the computed abscissas are within a few ulps of the correctly rounded values and the weights within about one hundred, for up to at least 200 points.

In double precision the recurrence overflows for more than 706 points (the outer weights underflow long before that). integrate then raises an evaluation_error through the Policy; if the policy does not throw, the result is a NaN.

As with gauss, defining BOOST_MATH_GAUSS_NO_COMPUTE_ON_DEMAND turns any combination of type and number of points without pre-computed values into a compile-time error. The program gauss_hermite_constants.cpp generates the pre-computed tables, and is easily modified for other numbers of points.

References