Gauss-Kronrod

Motivation

The \((2N+1)\)-point Gauss-Kronrod rule simultaneously produces a high-accuracy integral estimate and a reliable error bound in a single pass over \(2N+1\) function evaluations. It does so by augmenting an \(N\)-point Gauss-Legendre rule with \(N+1\) optimally placed Kronrod points: the original \(N\) Gauss nodes are reused, so no function evaluations are wasted.

This makes it the method of choice when both accuracy and a trustworthy error estimate are needed. It is the engine behind QUADPACK and SciPy’s scipy.integrate.quad.

Example

use integrate::gauss_kronrod::gauss_kronrod_rule;

let f = |x: f64| x.exp();

let a = 0.0_f64;
let b = 1.0_f64;
let n = 7_usize; // 15-point Gauss-Kronrod rule

match gauss_kronrod_rule(f, a, b, n) {
    Ok((integral, error)) => {
        println!("Integral:       {:.10}", integral);
        println!("Error estimate: {:.2e}", error);
    }
    Err(e) => println!("Error: {}", e),
}

Understanding Gauss-Kronrod rule

The Gauss-Kronrod rule is constructed in two stages.

Stage 1: Gauss-Legendre base rule

The \(N\)-point Gauss-Legendre rule approximates \(\int_{-1}^{1} f(x)\,dx\) by

\[ G_N(f) = \sum_{i=1}^{N} w_i^G f(x_i) \]

where \(x_1, \ldots, x_N\) are the zeros of the Legendre polynomial \(P_N\) and \(w_1^G, \ldots, w_N^G\) are positive weights chosen so that the formula is exact for all polynomials of degree at most \(2N - 1\).

Stage 2: Kronrod extension

The Stieltjes polynomial \(E_{N+1}\) of degree \(N+1\) is the unique monic polynomial satisfying

\[ \int_{-1}^{1} E_{N+1}(x) \cdot x^k \cdot P_N(x) \, dx = 0, \quad k = 0, 1, \ldots, N \]

Its \(N+1\) zeros \(t_1, \ldots, t_{N+1}\) interlace strictly with those of \(P_N\) and all lie in \((-1, 1)\). Adding them to the \(N\) existing Gauss-Legendre nodes gives \(2N+1\) distinct abscissas. Weights for all \(2N+1\) points can then be chosen so that the combined formula

\[ GK_{2N+1}(f) = \sum_{i=1}^{N} w_i^{GK} f(x_i) + \sum_{j=1}^{N+1} v_j^{GK} f(t_j) \]

is exact for all polynomials of degree at most \(3N + 1\) (when \(N\) is odd). Since the Gauss nodes are reused, the total cost is \(2N+1\) function evaluations.

Change of variables to \([a, b]\)

For an integral over a general interval \([a, b]\), the substitution \(x = \dfrac{(b-a)\,t + (b+a)}{2}\) gives

\[ \int_a^b f(x) \, dx = \frac{b-a}{2} \int_{-1}^{1} f\left(\frac{(b-a)\,t + (b+a)}{2}\right) dt \]

so the \((2N+1)\)-point Gauss-Kronrod approximation on \([a, b]\) is

\[ GK_{2N+1}(f, a, b) = \frac{b-a}{2} \cdot GK_{2N+1}(\tilde{f}), \quad \tilde{f}(t) = f\left(\frac{(b-a)\,t + (b+a)}{2}\right) \]

Error estimate

gauss_kronrod_rule returns both the integral and an error estimate:

• integral \(= GK_{2N+1}(f, a, b)\)
• error \(= \left\lvert GK_{2N+1}(f, a, b) - G_N(f, a, b) \right\rvert\)

The error is the absolute difference between the \((2N+1)\)-point Kronrod estimate and the \(N\)-point Gauss-Legendre estimate computed on the same set of function values. For smooth integrands this quantity is a reliable indicator of the absolute integration error.

Choosing \(N\)

The default \(N = 7\) (15-point rule) matches the default used by QUADPACK and SciPy on each subinterval and gives double-precision accuracy for most smooth integrands.

Limitations

The single-interval Gauss-Kronrod rule requires the integrand to be smooth on \([a, b]\). For integrands with singularities, discontinuities, or sharp peaks the error estimate may be unreliable and the accuracy may be poor. In such cases the interval should be split manually or an adaptive method should be used instead.

The rule is defined only on finite intervals \([a, b]\). For semi-infinite or doubly-infinite domains use the Gauss-Laguerre or Gauss-Hermite rules instead.

References

Node and weight computation ported from the original Fortran 77 by Robert Piessens and Maria Branders (C translation by John Burkardt, GNU LGPL):

\[ [1]: \text{Piessens R., Branders M., } \textit{A Note on the Optimal Addition of Abscissas to Quadrature Formulas of Gauss and Lobatto,} \\ \text{Mathematics of Computation, Vol. 28, No. 125, 1974, pp. 135\text{–}139.} \]