Absolute and Relative Error Control in Composite Interpolatory Quadrature : the CIRQUE Algorithm

We introduce the CIRQUE algorithm, for approximating definite integrals of continuous, univariate, real-valued functions, using positive-coefficient composite interpolatory quadrature. CIRQUE estimates and controls absolute and/or relative error, without the need for a prior estimate of the magnitude of the integral. The limiting effects of roundoff error are catered for, and CIRQUE is able to provide estimates of error bounds as output. Moreover, if these bounds are deemed too large, it is a simple matter to rerun CIRQUE once to obtain an acceptable bound. We have demonstrated the algorithm using the Trapezium rule, Simpson’s rule and four-point Gauss-Legendre quadrature.


Introduction
Say we wish to find an approximation Q f to a definite integral I f of a continuous, univariate, real-valued function f (x) , such that where ε is a user-defined tolerance.In other words, we seek to control the relative error in the approximation.
Obviously, if we have a good estimate of the integral, I f ≈ I f , we can easily test whether or not the tolerance is satisfied by Q f , and if it is not, then we can, possibly, refine Q f .However, if we do not have I f , we cannot test the quality of Q f .One way around this problem is to use a numerical method of high order to compute I f , and a method of lower order to compute Q f ; such an approach requires the assumption that I f is suitably accurate without any means to verify such assumption and, furthermore, two different methods must be employed (Krommer and Ueberhuber, 1998).
So, the question is, can we design an algorithm that allows the relative error to be estimated and controlled without prior knowledge of I f , and without using more than one numerical method?Such an algorithm is the subject of this paper.Our algorithm, designated CIRQUE (an acronym formed from the initial letters of Relative Error Control in Interpolatory QUadrature), will also detect whether or not absolute error control is appropriate, instead of relative error control; we will provide a criterion for such a selection based on the magnitude of I f .Also, our algorithm will cater for any limitations arising from the presence of roundoff error.
In the next section, we briefly discuss concepts, terminology and notation relevant to the algorithm.These include interpolatory quadrature, error bounds in composite interpolatory quadrature, error control, roundoff error, and a distinction between relative and absolute error control based on computational efficiency.In Section 3 we describe the CIRQUE algorithm, and thereafter we present a few numerical examples.We also discuss a procedure for refining the algorithm's estimates of error bounds.In Section 6, we list the sequence of operations in CIRQUE, as a summary of the algorithm.

Relevant Concepts, Terminology and Notation
In this section, we provide appropriate background information.Most of the concepts used in this paper are wellestablished "classical" ideas and are drawn from numerous books on numerical analysis and computational integration.Such a list, which constitutes our general bibliography, includes Burden and Faires (2011), Davis andRabinowitz (1984), Engels (1980), Ghizetti and Ossiccini (1970), Isaacson and Keller (1994), Kincaid and Cheney (2002), Krommer and Ueberhuber (1998), Stroud (1974), and Stroud and Secrest (1966).

Interpolatory Quadrature
The integral of f (x) on [a, b] is denoted The wherein w i are the quadrature weights, and x i are the quadrature nodes.Clearly, Q f ; a, b is a linear combination of values of f sampled at n discrete nodes on [a, b], and is intended to be an approximation to I f ; a, b .The precise nature of the weights depends on the nature of the quadrature but, generally speaking, they are dependent on the length of [a, b].
In this paper, we restrict our considerations to a particular class of quadrature, known as interpolatory quadrature.We will not provide details of this class here (the reader is referred to the bibliography given earlier for extensive discussions of interpolatory quadrature), suffice it to say that the form of the absolute approximation error in interpolatory quadrature is ideally suited to our objectives, and we discuss this feature in the next section.Hence, from here on, Q f ; a, b refers to interpolatory quadrature.
If [a, b] is subdivided into N subintervals (denoted D j ), and Q f ; D j is applied on each subinterval, then the integral of f (x) on [a, b] is simply the sum of these individual quadratures.We denote such composite interpolatory quadrature by where the double sum notation indicates the sum of quadrature on each of N subintervals D j .
It is possible, and often convenient, to write Q f ; a, b in terms of a stepsize parameter h, as in where c i = w i /h.We choose to take this stepsize as the average separation of the nodes on [a, b] ; we will say more about h in the next section.We must make mention of the fact that in a feature that we will exploit later.

Error Bounds in Composite Interpolatory Quadrature
The absolute error in composite interpolatory quadrature is bounded as where A, θ and r are dependent on the number of nodes n in the underlying quadrature rule Q f ; D j , h (Isaacson and Keller, 1994).For example, the composite Simpson's rule has A = 1 180 , θ = 4 and r = 4 (r indicates the order of the method).If we impose a tolerance ε on the approximation, then the inequality allows the stepsize h to be determined, such that h is consistent with the maximum possible error and the imposed tolerance.Of course, if the maximum error is less than the imposed tolerance, then the actual error will also satisfy the tolerance.Once h has been determined, it is a simple matter to compute CQ f ; a, b, h .Indeed, h is given by We note here that if the underlying quadrature rule has nonuniform stepsize, such as a Gaussian rule, then h given in ( 9) is the average separation of the nodes and endpoints on [a, b] .For uniform stepsize rules, such as Newton-Cotes rules, h in ( 9) is the separation of any two adjacent nodes, because in such a rule the nodes are uniformly spaced.Note that If we choose to impose a tolerance ε on the relative error, we must demand so that ε I f ; a, b serves as an absolute tolerance, and h must be determined from The difficulty is immediately obvious: we do not know I f ; a, b .

Inclusion of Roundoff Error
Taking into account the presence of roundoff error, we write where w i and f (x i ) are exact, μ w i is the roundoff error in the computed value of w i , μ f i is the roundoff error in the computed value of f (x i ) , and the double sum notation has been explained previously.Obviously, it is now understood that CQ f ; a, b, h indicates an approximation computed with a finite-precision device. Thus, We define the second term on the RHS of (15) as the roundoff error RO in CQ f ; a, b, h , and so where μ is an upper bound on μ w i and μ f i (on our platform μ ∼ 10 −16 ), and in the last line we have ignored μ 2 .Now, if we use a quadrature rule in which all weights are positive (so that their sum is equal to the length of the interval on which they are defined), we then have (Burden and Faires, 2011;Isaacson and Keller, 1994) and if both max This last case is mentioned because it will be relevant later.
As regards absolute error control, roundoff error is incorporated as which gives If ε |RO| , h * is negative, imaginary or zero, none of which is admissible.We therefore require that ε > |RO| .
In other words, |RO| sets a lower limit on the tolerance that can be imposed.For relative error control, the condition becomes Once h * has been determined, we determine a new stepsize h using ( 9), consistent with an integer value of N.This new stepsize h is the actual stepsize used in computing the quadrature.
A very important point must be made here: in (20), h * clearly depends on |RO| so that, in (11), |RO| must be known in order to determine N.However, in ( 16), we see that |RO| is computed by summing over N subintervals; this is impossible if N is not known.Only if we use positive-coefficient quadrature can we eliminate the need to know N in computing |RO| .
We will emphasize this point again in Section 3.

Absolute and Relative Error Control
We give a criterion for choosing between absolute and relative error control.We have where B is an appropriate coefficient.Hence, for relative error control, and, for absolute error control, In the case I f ; a, b > 1, we use relative error control because of the larger stepsize.In the case I f ; a, b 1, we use absolute error control, also due to the larger stepsize.Using a larger stepsize implies greater efficiency, because fewer nodes are required on [a, b] .In other words, our criterion for implementing absolute or relative error control is based on considerations of computational efficiency.The choice is clearly dictated by the magnitude of I f ; a, b .

The CIRQUE Algorithm
We seek to estimate and control the error in a numerical approximation to We first make a change of variables In other words, we transform the integral to one defined on a unit interval.Any unit interval will do, provided f (mz + c) is bounded thereon.
So, we now have as the integral we must approximate.
Next, we define where m f (z) 1, then |g (z)| 1, so it is not necessary for M to be less than unity, since the purpose of M is merely to ensure that |g (z)| 1.Note that Stated otherwise, and so absolute error control is appropriate for approximating I g; a 2 , b 2 .
We demand where CQ g; a 2 , b 2 , h is any composite positive-coefficient interpolatory quadrature of g on [a 2 , b 2 ] , and ε g is a userdefined tolerance with ε g > 2μ (as stated earlier, 2μ is the maximal roundoff error associated with composite positivecoefficient interpolatory quadrature of a function with magnitude less than or equal to unity on a unit interval).We use ( 9) and ( 20) to find the largest stepsize h consistent with this inequality; it goes without saying that here we use max Now, we have Hence, The last inequality is due to the fact that the change of variable preserves both the integral and the quadrature (Krommer and Ueberhuber, 1998;Stroud, 1974).The LHS of the last inequality is the relative error in is an estimated upper bound on this relative error.Moreover, this estimate is good if ε g is small (because then the approximation CQ g; a 2 , b 2 , h is reliable).
Note the following: 1.If we find that MCQ g; a 2 , b 2 , h 1, then absolute error control is appropriate, and we can accept MCQ g; a 2 , b 2 , h as the approximation, consistent with the absolute tolerance Mε g .If Mε g is unacceptably large, then we simply redo the computation with a suitably smaller absolute tolerance ε g .

2.
Obviously, choosing ε g so that Mε g is of acceptable magnitude would avoid having to redo the computation, in the event that absolute error control is appropriate.Note that we can only know if absolute error control is appropriate once CQ g; a 2 , b 2 , h has been determined, and it may be that absolute error control is not appropriate -after all, CQ g; a 2 , b 2 , h is not known a priori.Nevertheless, we must be careful to ensure that ε g is not less than the maximal roundoff error 2μ.This might occur if M is very large, which probably means that relative error control would be relevant, anyway.
3. If we find that MCQ g; a 2 , b 2 , h > 1, then relative error control is appropriate.If the estimated bound ε g |CQ[g;a2,b2,h]| on the relative error is too large, we simply repeat the algorithm, with an appropriately reduced value of ε g (again, taking care that ε g 2μ).
4. Typically, we would choose ε g 2μ, but we should also ensure that ε g is small enough so that CQ g; a 2 , b 2 , h may be considered a reliable approximation.
We will say more about repeating CIRQUE with a modified ε g in the next section.

Numerical Examples
We demonstrate CIRQUE with two numerical examples.In both examples, we use μ = 2 −53 ≈ 10 −16 .In this section, we will use CQ g as shorthand for CQ g; a 2 , b 2 , h .
In Table 1, we show results for CIRQUE with Trapezium (A = 1 12 , θ = 2, r = 2), Simpson and four-point Gauss-Legendre quadrature (GL4, A = 0.00022, θ = 8, r = 8), for various tolerances.In this ] used for the composite quadrature; N f is the total number of evaluations of the integrand f (x) used in the composite quadrature; and N abs f is the number of evaluations of f (x) that would have been needed if absolute error control was implemented on [a 1 , b 1 ] , instead of relative error control.
For each ε g , |Δ R | UB is the same (to the indicated precision) for each method; this is simply due to the fact that, in each case, CQ g is sufficiently accurate.In all cases, the integral has large magnitude.Values of N and N f decrease with the order of the method, but increase with decreasing ε g .The values of N abs f show just how computationally expensive absolute error control can be, particularly for the low-order Trapezium rule.

The interpretation of |Δ
i.e. the exact value lies in a MCQ g |Δ R | UB -neighbourhood of MCQ g or, equivalently, a |Δ A | UB -neighbourhood of MCQ g , whatever the value of I f ; a 1 , b 1 might be.
Now, let us say we are unimpressed with the estimate |Δ R | UB = 3.2 × 10 −8 obtained for Simpson's rule with ε g = 10 −8 , and we would prefer a bound of 10 −9 .We simply repeat CIRQUE once using as the new tolerance.So, if the estimated upper bound on the relative error is not acceptable, it is a simple matter to correct it.The same process holds for the absolute error bound, as well.This is because both |Δ R | UB and |Δ A | UB are proportional to ε g .Generally, if we wish to modify |Δ R | UB and/or |Δ A | UB by a factor η, we must repeat CIRQUE with ε g modified according to subject, of course, to the condition that this modified value of ε g cannot be less than 2μ.This will yield estimated upper bounds of η |Δ R | UB and η |Δ A | UB .In the above example, η = 10 −9 3.2×10 −8 .2. We approximate 2π 0 sin xdx = 0 as an example of an integral for which a relative error cannot be computed.Transforming to the interval [0, 1], we find M = 2π.

Results are shown in
i.e. the exact value lies in a |Δ A | UB -neighbourhood of MCQ g , whatever I f ; a 1 , b 1 may be.The smallness of |Δ A | T in all cases is simply due to the high degree of antisymmetry present in the example, and the symmetrical node distribution in the three quadrature methods; if we did not know |Δ A | T , the only indication of the accuracy of the approximation is the |Δ A | UB -neighbourhood.

Bound Refinement
We see in where We

Summary of the CIRQUE Algorithm
It is worthwhile to summarize the CIRQUE algorithm, in the form of a list of the sequence of operations of the algorithm.
The order θ of the derivative f (θ) (x) must correlate with the underlying quadrature used by CIRQUE (e.g. for the Trapezium rule, θ = 2).CIRQUE then performs the following: 2. Determines M and normalizes the transformed integrand, so defining a new integrand g (z) .8. If the bounds in #5 or #6 (or #7, for that matter) are unacceptably large, CIRQUE is repeated with an appropriately modified tolerance.

Concluding Comments
We have developed an algorithm, designated CIRQUE, that computes a numerical approximation to a definite integral, using composite positive-coefficient interpolatory quadrature.CIRQUE is able to distinguish between the need for absolute and relative error control, and to implement such error control, without a priori knowledge of the magnitude of the integral.The criterion for choosing between absolute and relative error control is based on computational efficiency.Moreover, CIRQUE can provide an a posteriori estimate of the maximum error incurred in the quadrature process and, if such bound is unacceptably large, it is easy to rerun CIRQUE so as to achieve an acceptable bound.Roundoff error present in the quadrature process has been taken into account in the error control algorithm.The requirements of the integrand are that it is real-valued, univariate, and continuous on the interval of integration, and that the maximum magnitude of its relevant higher-order derivative is known or can be found.
Future research efforts should concern the development of CIRQUE to handle multivariate integrands, and the possible use of the algorithm in the context of adaptive quadrature.
where • • • indicates rounding up to the nearest integer, and ± refers to Q f ; a, b, h as being open (+) or closed (−).If it is left-open or right-open, then we use n in place of n ± 1.
has length (n ± 1) h (open (+) or closed (−) quadrature) or nh (left-or right-open quadrature), so that (n ± 1) hN or nhN is the length b − a of the interval of integration.
[a,b] | f (x)| and b − a are less than or equal to unity, we have |RO| 2μ.(18) so that M 1, always, and |g (z)| 1.If max [a 2 ,b 2 ] , h the nodes x i are not necessarily uniformly distributed on [a, b]; nor are the endpoints of [a, b] necessarily nodes.If the nodes are uniformly/not uniformly spaced, Q f ; a, b, h is said to have uniform/nonuniform stepsize.If a is a node, and b is not a node, then Q f ; a, b, h is said to be left-open; a similar definition holds for right-open quadrature.If neither endpoint is a node, Q f ; a, b, h is termed open.
table, |Δ R | UB is the estimated upper bound ε g |CQg| on the relative error; |Δ R | T is the actual upper bound on the relative error; |Δ A | UB is the estimated upper bound Mε g on the absolute error; |Δ A | T is the actual absolute error; N is the number of subintervals on [a 2 , b 2 Table 2, where symbols have the same meaning as in Table 1.The entries have the same character as those in Table 1, all exhibiting expected behaviour.We do not show |Δ R | UB or |Δ R | T because relative error control is meaningless in this example.The interpretation of |Δ A | UB is that Table 1 that |Δ A | T < |Δ A | UB , in some instances by two orders of magnitude.This implies that |Δ A | UB is not as tight an upper bound as we might like.The reason for this is that |Δ A | UB is, effectively, the upper bound consistent with the stepsize h * in (20).However, the stepsize actually used in CQ g; a 2 , b 2 , h is h in (9), in which b 2 −a 2 (n±1)h * has been rounded up to the nearest integer.Consequently, we have h h * .This gives see that |Δ| min and |Δ| max are quadrature approximation errors (in CQ g; a 2 , b 2 , h ) for the actual stepsize h, whereas |Δ A | UB is the quadrature error (in CQ g; a 2 , b 2 , h ) for h * plus a roundoff term.Now, it is easy to compute |Δ| min and |Δ| max , so we can present |Δ| min , |Δ| max as an interval indicating the minimum and maximum absolute error achievable with the value of the stepsize h used in CQ g; a 2 , b 2 , h .The actual absolute error lies within these bounds.Multiplying this interval by M gives error bounds on the approximation to MCQ g; a 2 , b 2 , h ≈ I f ; a 1 , b 1 , and dividing this interval by CQ g; a 2 , b 2 , h gives bounds on the relative error in the approximation.Since these bounds have been determined using the actual stepsize h, rather than h * , we expect them to be tighter bounds.For example, for f (x) = e x and ε g = 10 −4 , we findM |Δ| min , |Δ| max = [22, 441] M |Δ| min , |Δ| max = [2, 36](43)for Simpson and GL4 quadrature, respectively (we have rounded these numbers to nearest integer, for ease of presentation).In both cases,|Δ A | T ∈ M |Δ| min ,|Δ| max , and M |Δ| max < |Δ A | UB .For the relative error, we find |Δ| min , |Δ| max .Since 2μ is an upper bound on |RO| for the quadrature on [a 2 , b 2 ] , it could occur that |RO| = 0. Hence, we retain |Δ| min as the lower bound in (45).In summary, then, we can use |Δ| max to find tighter upper bounds on both relative and absolute errors in CIRQUE, and these bounds can be presented instead of |Δ R | UB and |Δ A | UB .Also, |Δ| min provides a lower bound on the accuracy, although this is not of primary interest.
R | UB .We have confirmed that this holds for all cases considered in the two examples (for f (x) = sin x, min g (θ) (x) = 0, so the very small values of |Δ A | T are accommodated).
composite positive-coefficient interpolatory quadrature CQ g; a 2 , b 2 , h to approximate the integral of g (z) on [a 2 , b 2 ] , with appropriate absolute error control, subject to tolerance ε g ., and the maximal absolute error as Mε g .6.If MCQ g; a 2 , b 2 , h 1, the maximal absolute error is estimated as Mε g .The relative error is not estimated in this case because such estimate could be unreliable, particularly if MCQ g; a 2 , b 2 , h is close to zero.7. The bounds in #5 and #6 could be replaced with the refined bounds |Δ| max |CQ[g;a2,b2,h]| and M |Δ| max .