Gaussian Quadrature

What Gaussian quadrature does

Unlike equally spaced rules, it selects the points to maximise polynomial exactness for a specified number of function evaluations.

For a weighted integral, the general $n$-point rule is

$$\int_a^b w(x)f(x)\,dx \approx \sum_{i=1}^{n} W_i f(x_i),$$

where $x_i$ are the quadrature nodes and $W_i$ are the quadrature weights. The word Gaussian refers to Carl Friedrich Gauss and the optimal-node construction, not to an assumption that the integrand or model inputs follow a normal distribution.

Why the nodes are different from equally spaced rules

Gaussian quadrature chooses nodes as the roots of an orthogonal polynomial associated with the interval and weight function. This construction uses both the node locations and their weights as degrees of freedom, allowing an $n$-point rule to integrate every polynomial of degree at most $2n-1$ exactly under its defining weight.

For Gauss–Legendre quadrature on $[-1,1]$, the weight function is $w(x)=1$ and the nodes are the roots of the Legendre polynomial $P_n(x)$. The approximation is

$$\int_{-1}^{1} f(x)\,dx \approx \sum_{i=1}^{n} W_i f(x_i).$$

The Gauss–Legendre weights can be written as

$$W_i=\frac{2}{\left(1-x_i^2\right)\left[P_n'(x_i)\right]^2}.$$

The exactness result applies to polynomials, but smooth non-polynomial functions are also often approximated accurately because they can be represented well by polynomials over a suitable interval. It does not mean that every smooth function is integrated exactly.

Common Gaussian quadrature families

The integration domain and any natural weight function determine the appropriate family of orthogonal polynomials. Choosing the family to match the mathematical form can avoid inefficient transformations and improve numerical behaviour.

The weight function is part of the quadrature definition and must not be counted twice. For example, Gauss–Hermite approximates an integral containing $e^{-x^2}$; a normal expectation requires the correct change of variables and normalising constant.

Mapping a finite interval to the standard domain

Gauss–Legendre nodes and weights are normally tabulated on $[-1,1]$. An affine transformation maps them to any finite interval $[a,b]$ without recomputing the standard nodes.

$$x_i^{[a,b]}=\frac{a+b}{2}+\frac{b-a}{2}x_i$$

and

$$\int_a^b f(x)\,dx =\frac{b-a}{2}\int_{-1}^{1}f\!\left(\frac{a+b}{2}+\frac{b-a}{2}z\right)dz \approx\frac{b-a}{2}\sum_{i=1}^{n}W_i f\!\left(x_i^{[a,b]}\right).$$

The factor $(b-a)/2$ is the Jacobian of the transformation. Omitting it is a common implementation error that changes the scale of the integral.

A three-point worked calculation

The three-point Gauss–Legendre rule provides a compact example with auditable values. Its nodes are $-\sqrt{3/5}$, $0$, and $\sqrt{3/5}$, with corresponding weights $5/9$, $8/9$, and $5/9$.

$$\int_{-1}^{1}f(x)\,dx \approx \frac{5}{9}f\!\left(-\sqrt{\frac35}\right) +\frac{8}{9}f(0) +\frac{5}{9}f\!\left(\sqrt{\frac35}\right).$$

For $f(x)=x^4$, the calculation is

$$\frac{5}{9}\left(\frac{9}{25}\right) +\frac{8}{9}(0) +\frac{5}{9}\left(\frac{9}{25}\right) =\frac{2}{5}.$$

This equals the analytical integral $\int_{-1}^{1}x^4\,dx=2/5$. The exact result is expected because a three-point Gaussian rule is exact through polynomial degree five.

A health-economic present-value example

Continuous-time health-economic models may integrate a flow of costs, health outcomes, hazards, or occupancy over time. Gaussian quadrature can evaluate such integrals efficiently when the integrand is smooth but repeated analytical integration is inconvenient.

Suppose a constant cost flow of £1,000 per year is discounted continuously at rate $r=0.03$ over five years. The present value is

$$PV=\int_0^5 1000e^{-0.03t}\,dt.$$

Using the three-point Gauss–Legendre nodes $z_i$ and weights $W_i$, transform each node as $t_i=2.5+2.5z_i$ and calculate

$$PV\approx 2.5\sum_{i=1}^{3}W_i\,1000e^{-0.03t_i} =4643.06745247.$$

The analytical value is

$$PV=1000\frac{1-e^{-0.15}}{0.03}=4643.06745250.$$

The absolute difference is approximately $2.6\times10^{-8}$ in currency units. This accuracy reflects the smooth integrand and should not be assumed for functions with discontinuities, sharp thresholds, or singular behaviour.

Error and convergence

For sufficiently smooth $f$, the Gauss–Legendre error depends on a derivative of order $2n$. The formal expression explains the high polynomial exactness, but it is rarely a convenient practical error estimate because the derivative may be unavailable or poorly bounded.

For some $\xi\in(-1,1)$, the $n$-point remainder has the form

$$R_n(f)=\frac{2^{2n+1}(n!)^4}{(2n+1)[(2n)!]^3}\,f^{(2n)}(\xi).$$

In practice, convergence is commonly assessed by repeating the calculation with a higher order and comparing results. Agreement between orders is evidence of numerical stability, not a proof of accuracy, so difficult integrands also require structural diagnostics or an independent method.

Composite and adaptive use

Increasing the order is not always the best response to a difficult integrand. Composite Gaussian quadrature divides the domain into subintervals and applies a fixed-order rule within each part, while adaptive schemes refine regions where estimated error is greatest.

A composite calculation can be written as

$$\int_a^b f(x)\,dx =\sum_{j=1}^{m}\int_{a_j}^{b_j}f(x)\,dx \approx\sum_{j=1}^{m}\frac{b_j-a_j}{2} \sum_{i=1}^{n}W_i f\!\left(\frac{a_j+b_j}{2}+\frac{b_j-a_j}{2}x_i\right).$$

Known discontinuities, kinks, treatment thresholds, or changes in hazard functions should usually become subinterval boundaries. Blindly raising $n$ across the whole interval can waste evaluations and may hide local numerical problems.

Expectations and multiple dimensions

Many model quantities are expectations rather than simple finite-interval integrals. Gaussian quadrature can approximate them by pairing a distribution-specific transformation or weight function with deterministic evaluations of the model.

For a continuous random variable $X$ with density $p(x)$,

$$\mathbb{E}[h(X)]=\int h(x)p(x)\,dx.$$

Independent dimensions can be combined using tensor-product rules. If each of $d$ dimensions uses $n$ nodes, however, the full grid requires $n^d$ evaluations, which is the curse of dimensionality.

$$\int\!\cdots\!\int f(x_1,\ldots,x_d)\,dx_1\cdots dx_d \approx\sum_{i_1=1}^{n}\cdots\sum_{i_d=1}^{n} W_{i_1}\cdots W_{i_d}f(x_{i_1},\ldots,x_{i_d}).$$

Sparse grids, dimension reduction, quasi-Monte Carlo methods, or ordinary Monte Carlo simulation may be preferable when dimensionality is high. Gaussian quadrature is deterministic and should not be described as a sampling distribution or a probabilistic sensitivity analysis by itself.

Choosing and implementing a rule

A reliable implementation starts with the integral's mathematical structure rather than a default order. The analyst should identify the domain, weight, smoothness, transformations, and required accuracy before selecting a quadrature family.

1. Define the integral, its units, and the decision output it contributes to.
2. Inspect the domain for infinite limits, singularities, discontinuities, and sharp changes.
3. Choose the Gaussian family whose domain and weight match the integral.
4. Transform the interval or random variable and retain every Jacobian and normalising factor.
5. Select an initial order and obtain validated nodes and weights.
6. Evaluate the integrand at each transformed node using consistent units and parameters.
7. Multiply each evaluation by its weight and sum the weighted values.
8. Increase the order or subdivide the interval and compare the resulting estimates.
9. Validate against an analytical case or an independent integration method.
10. Record the rule, order, transformation, tolerances, software, and convergence evidence.

In a spreadsheet, store one row per node with columns for the standard node, weight, transformed node, function value, and weighted contribution. The final integral should use an explicit sum of weighted contributions multiplied by the interval scale, which makes the Jacobian and every model evaluation auditable.

Validation checks for a model

Quadrature code should be tested independently of the health-economic model in which it is embedded. Simple functions with known integrals reveal incorrect nodes, weights, transformations, and normalisation before those errors affect a decision result.

• The weights for an $n$-point Gauss–Legendre rule on $[-1,1]$ should sum to $2$.
• The rule should reproduce monomials through degree $2n-1$ to numerical precision.
• Symmetric nodes and weights should cancel odd functions over a symmetric interval.
• Increasing the order or refining subintervals should stabilise the decision-relevant output.
• An independent adaptive integrator should give a compatible result within the required tolerance.
• Boundary cases and known discontinuities should be tested separately rather than averaged into a global check.

Validation must cover downstream results as well as the integral itself. A numerically small integration error may still matter when an incremental cost-effectiveness result is close to a decision threshold or when two strategies have nearly equal outcomes.

When Gaussian quadrature works poorly

Gaussian quadrature is most effective for smooth integrands that can be evaluated at arbitrary points. It is less suitable when only tabulated observations are available at fixed locations, because the method cannot freely choose its optimal nodes.

Performance may deteriorate with discontinuities, endpoint singularities, highly oscillatory functions, extremely narrow peaks, or expensive high-dimensional models. Appropriate remedies include splitting the interval, transforming the integrand, using a specialised weighted rule, applying adaptive integration, or selecting a different numerical method.

How to report the result

Transparent reporting allows another analyst to reproduce both the integral and its role in the model. Reporting only the final number conceals choices that can materially affect accuracy.

State the quadrature family, order, interval or distributional transformation, nodes-and-weights source, composite or adaptive strategy, stopping criterion, and validation method. Also report the model quantity being integrated, its units, any discounting convention, and whether numerical error was negligible relative to parameter uncertainty and the decision margin.