Next: , Previous: QAWF adaptive integration for Fourier integrals, Up: Numerical Integration

17.11 CQUAD doubly-adaptive integration

cquad is a new doubly-adaptive general-purpose quadrature routine which can handle most types of singularities, non-numerical function values such as Inf or NaN, as well as some divergent integrals. It generally requires more function evaluations than the integration routines in quadpack, yet fails less often for difficult integrands.

The underlying algorithm uses a doubly-adaptive scheme in which Clenshaw-Curtis quadrature rules of increasing degree are used to compute the integral in each interval. The L_2-norm of the difference between the underlying interpolatory polynomials of two successive rules is used as an error estimate. The interval is subdivided if the difference between two successive rules is too large or a rule of maximum degree has been reached.

— Function: gsl_integration_cquad_workspace * gsl_integration_cquad_workspace_alloc (size_t n)

This function allocates a workspace sufficient to hold the data for n intervals. The number n is not the maximum number of intervals that will be evaluated. If the workspace is full, intervals with smaller error estimates will be discarded. A minimum of 3 intervals is required and or most functions, a workspace of size 100 is sufficient.

— Function: void gsl_integration_cquad_workspace_free (gsl_integration_cquad_workspace * w)

This function frees the memory associated with the workspace w.

— Function: int gsl_integration_cquad (const gsl_function * f, double a, double b, double epsabs, double epsrel, gsl_integration_cquad_workspace * workspace, double * result, double * abserr, size_t * nevals)

This function computes the integral of f over (a,b) within the desired absolute and relative error limits, epsabs and epsrel using the cquad algorithm. The function returns the final approximation, result, an estimate of the absolute error, abserr, and the number of function evaluations required, nevals.

The cquad algorithm divides the integration region into subintervals, and in each iteration, the subinterval with the largest estimated error is processed. The algorithm uses Clenshaw-Curits quadrature rules of degree 4, 8, 16 and 32 over 5, 9, 17 and 33 nodes respectively. Each interval is initialized with the lowest-degree rule. When an interval is processed, the next-higher degree rule is evaluated and an error estimate is computed based on the L_2-norm of the difference between the underlying interpolating polynomials of both rules. If the highest-degree rule has already been used, or the interpolatory polynomials differ significantly, the interval is bisected.

The subintervals and their results are stored in the memory provided by workspace. If the error estimate or the number of function evaluations is not needed, the pointers abserr and nevals can be set to NULL.