# Method selection for numeric quadrature

Several families of methods exist for numeric quadrature. If I have a specific class of integrands how do I select the ideal method?

What are the relevant questions to ask both about the integrand (e.g. is it smooth? does it have singularities?) and the computational problem (e.g. error tolerance, computational budget)?

How do answers to these questions rule out or promote the various families of methods? For simplicity lets consider just single or low-dimensional integrals.

For example the Wikipedia article on QUADPACK states that the fairly general QAGS routine "uses global adaptive quadrature based on 21-point Gauss–Kronrod quadrature within each subinterval, with acceleration by Peter Wynn's epsilon algorithm"

• Probably more specific information is needed to answer this properly. There is no one-size-fits-all criteria, gaussian quadrature often works well for very smooth problems whereas other quadratures may be used in the presence of mild singularities. But if you're periodic, then simple trapezoid might cut it. Mar 26 '13 at 20:43
• @Reid.Atcheson, I think you're answering the question right now. I'm not asking what is the best method, I'm asking what sorts of questions would you ask and what would those answers tell you? How does one approach these sorts of problems in general? Mar 26 '13 at 20:46

First of all, you need to ask yourself the question if you need an all-round quadrature routine that should take an integrand as a black box. If so, you cannot but go for adaptive quadrature where you hope that the adaptivity will catch "difficult" spots in the integrand. And that is one of the reasons Piessens et al. chose for a Gauss-Kronrod rule (this type of rule allows you to calculate an approximation of the integral and an estimate of the approximation error using the same function evaluations) of modest order applied in an adaptive scheme (with subdivision of the interval with the highest error) until the required tolerances are reached. The Wynn-epsilon algorithm allows to provide convergence acceleration and typically helps in the cases where there are end-point singularities.

But if you do know the "form" or "type" of your integrand, you can tailor your method to what you need so the computational cost is limited for the accuracy you need. So what you need to look at:

Integrand:

• Smoothness: can it be approximated (well) by a polynomial from an orthogonal polynomial family (if so, Gaussian quadrature will do well)
• Singularities: can the integral be split in integrals with only end-point-singularities (if so, the IMT-rule or double exponential quadrature will be good on each sub-interval)
• Computational cost for evaluation?
• Can the integrand be computed? Or is only limited point-wise data available?
• Highly oscillatory integrand: look for Levin-type methods.

When dealing with singularities, one typically prefers them to be at the end-point of the integrals (see IMT, double exponential). If this is not the case, one can resort to Clenshaw-Curtis integration where you capture the singularities in the weight-function. One typically defines forms of singularities like $|x-c|^{-\alpha}$ and establishes expressions for the weights of the quadrature as a function of $c$ and $\alpha$.

Integration interval: finite, semi-infinite or infinite. In case of semi-infinite or infinite intervals, can they be reduced to a finite interval by a variable transformation? If not, Laguerre or Hermite polynomials can be used in the Gaussian quadrature approach.

I don't have a reference for a real flow sheet for quadrature in general, but the QUADPACK book (not the Netlib manpages, but the real book) has a flow sheet to select the appropriate routine based on the integral you want to evaluate. The book also describes the choices in algorithms made by Piessens et al. for the different routines.

For low-dimensional integrals, one typically goes for nested one-dimensional quadrature. In the special case of two-dimensional integrals (cubature), there exist integration rules for different cases of integration domains. R. Cools has collected a large number of rules in his Encyclopedia of cubature formulas and is the main author of the Cubpack package. For high dimensional integrals, one typically resorts to Monte Carlo type methods. However, one needs typically a very large number of integrand evaluations to get reasonable accuracy. For low-dimensional integrals, approximation methods like quadrature/cubature/nested quadrature often out-perform these stochastic methods.

General interesting references:

1. Quadpack, Piessens, Robert; de Doncker-Kapenga, Elise; Überhuber, Christoph W.; Kahaner, David (1983). QUADPACK: A subroutine package for automatic integration. Springer-Verlag. ISBN 978-3-540-12553-2
2. Methods of Numerical Integration: Second Edition, Ph. Davis and Ph. Rabinowitz, 2007, Dover Books on Mathematics, ISBN 978-0486453392
• Nice response. Why would QUADPACK have chosen the 21 point Gauss-Kronrod method in particular? Why the magic number? Mar 26 '13 at 22:06
• @MRocklin: I guess it was a nice trade-off between accuracy and efficiency: you don't want to overkill your quadrature rule (costly) but you don't want it to be too weak neither (too much subdivisions in the adaptive part). To be complete: in the QAG routine, the user must specify the quadrature rule; in the QAGS (with extrapolation), the default is the 21 point rule but this can be overruled by using the extended calling routine QAGSE. Mar 27 '13 at 7:04
• @GertVdE Very nice response indeed. Can you elaborate on the use of Clenshaw-Curtis to capture mid-interval singularities, or provide references? I haven't heard it used in this way before, and couldn't find any details from a quick googling. Thank you! Apr 3 '13 at 19:41
• @OscarB: sorry for the long delay, was out without net access (ah the good life). See the Quadpack book §2.2.3.3 and further; Branders, Piessens, "An extension of Clenshaw-Curtis quadrature", 1975, J.Comp.Appl.Math., 1, 55-65; Piessens, Branders, "The evaluation and application of some modified moments", 1973, BIT, 13, 443-450; Piessens, Branders, "Computation of oscillating integrals", 1975, J.Comp.Appl.Math., 1, 153-164. If you do a literature search for "Robert Piessens" somewhere between 1972 and 1980, you will find a lot of interesting papers. Apr 22 '13 at 9:33