Integration of functions over geometries

No, that’s not the case. For polynomials Gaussian quadrature is exact, but for analytic functions Gaussian quadrature converges exponentially fast with the number of quadrature points n, so you can indeed get to nearly machine precision (nearly 16 digits) rapidly, and if you need to go to large n then FastGaussQuadrature.jl computes the quadrature rule in O(n) time.

A downside of FastGaussQuadrature.jl is that it just gives you an n-point rule, it doesn’t tell you what n should be for a given tolerance and a given function, and the rule is not a “nested” rule — it doesn’t give you an error estimate. So you are only your own in finding the correct n for your problem, e.g. by repeatedly doubling n until it converges. (Repeatedly increasing the order is also called p-adaptive quadrature.)

QuadGK is h-adaptive, in that it repeatedly subdivides the integration interval until a tolerance is achieved, and it uses an embedded rule to get an error estimate. This is convenient, and it is also fairly robust if your integrand is badly behaved in some localized region because it will only refine the quadrature around that region.

If your function is analytic, then using high-order quadrature is theoretically more efficient, but on the other hand you rarely have to go to very high order in practice, and the time for the computation of the rule is less important if you are doing integration over and over for many different functions and/or domains. So for smooth functions you can simply pass order=n for a larger n to QuadGK.jl to use a higher-order rule (which is computed once and cached). The computation of the rule is O(n^2), so it will be much slower than FastGaussQuadrature.jl for large n, but since in practice you rarely need n larger than a few hundred it doesn’t matter too much in my experience.

The other big advantage of FastGaussQuadrature.jl is that it has built-in rules for weighted quadrature with many different weight functions, which allow you to efficiently handle known endpoint singularities. See the discussion of weighted Gauss rules in the QuadGK manual, and integrands with singularities and discontinuities.

PS. Although QuadGK is not asymptotically exponentially convergent for a fixed order, unlike repeatedly increasing the order, it often gives the illusion of exponential convergence up to a given precision. See also this post: (Insanely-) Higher Order Derivatives - #17 by stevengj

5 Likes