16.1 Introduction

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

     I = \int_a^b f(x) w(x) dx

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

     |RESULT - I|  <= max(epsabs, epsrel |I|)

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

     |RESULT - I| <= ABSERR <= max(epsabs, epsrel |I|)

In short, the routines return the first approximation which has an absolute error smaller than epsabs or a relative error smaller than epsrel.

Note that this is an either-or constraint, not simultaneous. To compute to a specified absolute error, set epsrel to zero. To compute to a specified relative error, set epsabs to zero. 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.

• Integrands without weight functions
• Integrands with weight functions
• Integrands with singular weight functions