Introduction

Each algorithm computes an approximation to a definite integral of the form,

where @math{w(x)} is a weight function (for general integrands @math{w(x)=1}). The user provides absolute and relative error bounds @math{(epsabs, epsrel)} which specify the following accuracy requirement,

where @math{RESULT} is the numerical approximation obtained by the algorithm. The algorithms attempt to estimate the absolute error @math{ABSERR = |RESULT - I|} in such a way that the following inequality holds,

The routines will fail to converge if the error bounds are too stringent, but always return the best approximation obtained up to that stage.

The algorithms in QUADPACK use a naming convention based on the following letters,

Q - quadrature routine

N - non-adaptive integrator
A - adaptive integrator

G - general integrand (user-defined)
W - weight function with integrand

S - singularities can be more readily integrated
P - points of special difficulty can be supplied
I - infinite range of integration
O - oscillatory weight function, cos or sin
F - Fourier integral
C - Cauchy principal value

The algorithms are built on pairs of quadrature rules, a higher order rule and a lower order rule. The higher order rule is used to compute the best approximation to an integral over a small range. The difference between the results of the higher order rule and the lower order rule gives an estimate of the error in the approximation.

The algorithms for general functions (without a weight function) are based on Gauss-Kronrod rules. A Gauss-Kronrod rule begins with a classical Gaussian quadrature rule of order @math{m}. This is extended with additional points between each of the abscissae to give a higher order Kronrod rule of order @math{2m+1}. The Kronrod rule is efficient because it reuses existing function evaluations from the Gaussian rule. The higher order Kronrod rule is used as the best approximation to the integral, and the difference between the two rules is used as an estimate of the error in the approximation.

For integrands with weight functions the algorithms use Clenshaw-Curtis quadrature rules. A Clenshaw-Curtis rule begins with an @math{n}-th order Chebyschev polynomial approximation to the integrand. This polynomial can be integrated exactly to give an approximation to the integral of the original function. The Chebyschev expansion can be extended to higher orders to improve the approximation. The presence of singularities (or other behavior) in the integrand can cause slow convergence in the Chebyschev approximation. The modified Clenshaw-Curtis rules used in QUADPACK separate out several common weight functions which cause slow convergence. These weight functions are integrated analytically against the Chebyschev polynomials to precompute modified Chebyschev moments. Combining the moments with the Chebyschev approximation to the function gives the desired integral. The use of analytic integration for the singular part of the function allows exact cancellations and substantially improves the overall convergence behavior of the integration.