I'm trying to understand the local error estimate that Quadpack (and subsequently other libraries such as GSL, quadpack++, cubature, etc.) uses for it's general adaptive quadrature subroutine QAG. The integration procedure is a Gauss-Kronrod rule and Quadpack uses a heuristic error estimate that I struggle to make sense of.

The error estimate can be found on page 67 in the Quadpack book or on page 11 in this excellent (and free) review of Error Estimation in Adaptive Quadrature. Or in the Quadpack source code for that matter. It states:

Let $G_n[a, b]$ is the n-point Gauss quadrature rule of degree $2n−1$ and $K_{2n+1}[a, b]$ is the $2n+1$ point Gauss-Kronrod extension of degree $3n+1$ which is in turn used as the approximation to the integral. The local error estimate is:

$$\varepsilon_k = \tilde{I}_k \min \left\{ 1, \left(200\frac{\left|G_n[a_k, b_k]-K_{2n+1}[a_k, b_k]\right|}{ \tilde{I}_k} \right)^{3/2} \right\} $$


$$ \tilde{I}_k = \int_{a_k}^{b_k}\left| f(x) - \frac{K_{2n+1}[a_k, b_k]}{b_k-a_k} \right| dx$$

which is also evaluated numerically using the $K_{2n+1}[a,b]$ rule.

The power $3/2$ has been chosen experimentally in such a way that the error scales exponentially with a break-even point at approximately relative machine precision for single floating point arithmetic. This means that the error estimate is less pessimistic for small values of $\left|G_n[a_k, b_k]-K_{2n+1}[a_k, b_k]\right|$ and more reliable for large values.

Indeed this is verified by plotting (letting $\tilde{I}_k = 1$):


Furthermore Quadpack also evaluates $$ \hat{I}_k = \int_{a_k}^{b_k}| f(x) | dx$$ using the $K_{2n+1}[a,b]$ rule. Denote by $f_{min}$ and $\varepsilon_{mac}$ the minimum floating point number and the floating point machine precision respectively (these number depends on single or double precision).

Finally Quadpack checks if $$\hat{I}_k > \frac{f_{min}}{50 \varepsilon_{mac}}$$

holds true and in this case they set

$$\varepsilon_k = \max \left\{50 \varepsilon_{mac} \hat{I}_k, \varepsilon_k \right\}$$

which will be used as the final error estimation of the integral on $[a_k,b_k]$.

Question 1:

What is the reason for doing the last step with the absolute value integrand. I do not under stand the point of it at all.

Question 2:

Is there a particular reason why they choose the approximate single precision machine epsilon as the break even point, or is it just an arbitrary small number? Why not use the double precision machine epsilon?

Thanks in advance!

  • $\begingroup$ Why do you think, that only the single precision machine epsilon is used? The source code in the book has only a (single) real version. Implementations which support double and single precision on different machines (e.g. SLATEC) use calls to different functions to dynamically get epmach, i.e. epmach = d1mach(4) or epmach = r1mach(4) $\endgroup$ – gammatester Apr 13 '16 at 9:19
  • $\begingroup$ Yes, I agree with you but my question was regarding the power 1.5 which is experimentally chosen to be such that the "break even point" (see the graph) is approximately epmach = r1mach(4). The power is 1.5 both in the single precision and double precision implementations. $\endgroup$ – DoubleTrouble Apr 13 '16 at 13:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.