Quadrature Routines
JARLE BERNTSEN University of Bergen and
TERJEO. ESPELID University of Bergen
Anew i~lgorithm forestimating theerror in quadrature approximations is presented. Basedon the same integrand evaluations that we need for approximating the integral, one may, for many quadrature rules, compute a sequence ofnull rule approximations. These null rule approxima- tions are then used to produce an estimate of the local error. The algorithm allows us to take advantage of the degree of precision of the basic quadrature rule. In the experiments we show that the algorithm works satisfactorily for a selection of different quadrature rules on all test families of integrals.
Categories and Subject Descriptors: G.1.4 [Quadrature and Numerical Differentiation]:
Adaptive Quadrature
Genera~ Terms: Algorithm, Design, Reliability
Additional Key Words and Phrases: Adaptive quadrature, automatic quadrature, error estima- tion, r-ml lrules, quadrature rules, reliability, symmetric rules
1. INTRODUCTION
Automatic algorithms are now used widely for the numerical calculation of integrals. Since the first such algorithm was given by McKeeman [161 in 1963, many new and sophisticated algorithms, both adaptive and non-adap- tive, have been developed [4, 6, 9, 17, 19]. For a more complete reference see Davis and Rabinowitz [5, pp. 425-434].
In automatic quadrature algorithms the estimate of the true error in the approximation of the integral governs the decision on whether to return the current approximation and terminate or to continue. Both the efficiency and the reliability therefore depend heavily on the error estimating procedure. In many adaptive algorithms the local error estimate is simply taken as the absolute value of the difference between two quadrature approximations.
Testing has shown, Berntsen [1] and Berntsen et al. [3], that routines applying such a simple local error estimate may be very unreliable.
Authors’ Address: University of Bergen, Thormdhlensgate 55, N-5008 Bergen, Norway.
Permission to copy without fee all or part of this material is granted provided that the copies are not made or distributed for direct commercial advantage, the ACM copyright notice and the title of the publication and its date appear, and notice is given that copying is by permission of the Association for Computing Machinery. To copy otherwise, or to republish, requires a fee and/or specific permission.
@ 1991 0098-3500/91/0600-0233 $01.50
ACM TransactIons on Mathematical Software, Vol. 17, No 2, June 1991, Pages 233-252
234 . J. Berntsen and T O. Espelid
More sophisticated local error estimating algorithms have been suggested by several authors: Laurie [11, 121, de Boor [61, and Piessens et al. [17].
Probably the most successful algorithm so far, see [1] and Espelid and S@evik [8], is QAG in QUADPACK, [171. This algorithm, using a Gauss- Kronrod rule as a basic rule, has a heuristic local error estimating algorithm developed especially for this type of rule. Data from experiments, showing the performance of the usual error estimating procedure, has been used to construct this local error estimating algorithm.
We will, in this paper, focus on the local error estimation procedure to be used in an adaptive algorithm for a one-dimensional integral over a finite interval H = [a, b]. A new local error estimation algorithm, which can be used on a selection of different basic integration rules in principally the same adaptive algorithm, will be developed. In the development of this new procedure both reliability y and economy are taken into account. The authors experience with adaptive multidimensional integration has been of great help in this work, and we have applied many of the ideas presented in this paper in multidimensional numerical integration [2, 4, 7].
We define the integral to be computed by I[f] =
/b f(x)dx.
a
An adaptive quadrature algorithm has input a, b, f and an error tolerance e.
The outpuf will normally be an estimate of the integral, Q, and an error estimate E. If possible, the algorithm computes Q with E s ~. This does not ensure that I 1 – Q I < ~, but it does suggest that it is likely. Basic elements in an adaptive quadrature algorithm include:
–An interval collection organized in some data structure, e.g., a heap. At a certain stage we may have M intervals, say HA, k = 1,...,M.
—A quadrature rule, Q, to produce a local estimate to the integral over each interval in the collection.
—A local error estimate procedure to estimate the error in each of the local approximations produced by Q. This error estimate, Ek, is usually con- structed using one or more null rules.
–A strategy for choosing the next interval to be processed. In a globally adaptive algorithm this is the interval with the greatest error estimate.
—A strategy for how to divide such an interval. Usually an interval is bisected.
Suppose that the local estimates are Q~ and E~ for the local integrals and errors respectively. Then the global estimates can be computed
(1)
(2)
ACM TransactIons on Mathematical Software, Vol. 17, No 2, June 1991
We give here the basic structure in such an algorithm:
A Globally Adaptive Quadrature Algorithm Initialize:
Control:
Process intervals:
Update:
Initialize the ~nter~al colle~tion and put M = 1; Produce QI and El;APut Q = QI and E = El;
wlhile E > cdo begin
Pick an interval from the collection; say interval Hk; Bisect this interval;
compute: Q(l), ~~1), Q~2J and fi~z);
q = Q + (QE) + Qf’) - Q,);
E = E+ (E~l) + E@ – E,);
Let these two new inte~vals replace interval H~ in the coLlection and put M = M + 1;
end
In the next section we give some general definitions and theorems concerning quadrature rules and null rules. In Section 3 we give a useful connection between divided differences and null rules. We develop the new error esti- mating algorithm in detail for symmetric quadrature rules and symmetric null rules only and give some special results concerning such rules in Section 4. The new error estimating algorithm is then described in Section 5 and in Section 6 we present the error estimating algorithm used in QAG in QUAD- PACK in our framework. In Section 7 we give the test results and in Section 8 some summarizing comments.
2. QUADRATURE RULES AND NULL RULES
Given n + 1 distinct points %,, i = O, . . . . n, and a quadrature rule based on these points
n
(3) xlandw,, i= O, l,..., n, are the rule’s nodes and weights respectively. By a simple translation this rule may be used on any of the local intervals produced by the adaptive algorithm.
Definition 1. A quadrature rule Q[ f] has degree d if it integrates exactly all polynomials of degree < d and fails to integrate exactly ~(x) = .xd+ 1.
We will only consider quadrature rules of at least degree O. The following theorem contains some well known results
THEOREM 1. An interpolator quadrature rule based on n + 1 distinct points has degree at least n. A quadrature rule based on n -t- 1 distinct points of degree d > n is unique and therefore has to be interpolator. A quadrature rule based on n + 1points has degree at most 2 n + 1(Gauss).
The term null rule was first used in 1965 by J. N. Lyness [13]. The following definition of a null rule is useful in this context.
Definition 2. A rule
n
N[f] = ,po U, f(x,).
ACM TransactIons on Mathematical Software, Vol
(4)
17, No. 2, June 1991
236 . J Berntsen and T. O. Espelid
is a null rule iff it has at least one nonzero weight, and in addition
‘&
L,=o.C=()
A null rule is furthermore said to have degree d if it integrates to zero all polynomials of degree s d and fails to do so with f(x) = Xd+l.
Suppose that Q and N have degrees d~ >0 and d~ respectively. Then
Q,[f] = Q[f] +
AJJ[f] = fo(LU, + W)f(XJ
is a quadrature rule of degree d > min (d~, d~). Thus we see that any null rule can be written as the difference between two different quadrature rules of degrees ? O, e.g.,
hN[f] = Q.[f] - Q[f]>
and the scaling of a null rule is similar to replacing the rule of reference Q~.
THEOREM 2. A null rule based on n + 1 points has degree of precision at
mostn – 1.
This is obvious from the fact that if a null rule has degree of precision at least n, then two different quadrature rules of precision at least n have to exist contradicting the uniqueness theorem.
Suppose that the n + 1 nodes are fixed and that the unique interpolator rule Q of degree d is chosen as the quadrature rule in the adaptive algo- rithm. In order to estimate the error in the approximation produced by Q we may use a sequence of rules of reference Ql, Qz, . . . of decreasing degrees d>dl>dz >.. .such that
NJ[f] = Q[f] - QJ[f]
is a null rule of degree d~. The Q~’s may be constructed by repeatedly removing points from the initial set of nodes checking that the interpolatory rule based on the remaining nodes actually has lower degree than the previous rule (this will be the case if the last removed node was used (weight # O) by the previous rule in the sequence). To produce rules of reference by removing nodes has been suggested by a number of authors.
Alternatively we may construct the null rules NJ directly by setting up the moment equations. Thus, in view of Theorem 2, we have
[’
X.X.n—11 xlxln—11 X2X2n—l..,1 ...““”. . . Xn—lx,n11[)
Unu~U1 = o. (5)The coefficient matrix in (5) has dimension n x (n + 1) and rank n. There- fore the null space of this matrix has dimension 1, giving a solution to (5)
ACM TransactIons on Mathematical Software, Vol 17, No 2, June 1991
which is unique except for a scaling factor. The system of equations
[
1 1 1 ... 1
N
u~X. xl X2 ““” Xn U1
. . . = o. (6)
X;–m x;–’” x;–m . . . x;”-”’ u,
has a coefficient matrix of dimension (n – m + 1) x (n + 1). It has rank n – m + 1 and the dimension of its null space is therefore m. Denote this null space S~, for m = 1, . . . . n. obviously SI c Sz c . . . c s.. Each of these null spaces is a linear vector space. A vector u, which is a solution of (6), is then in S~. Suppose that v is a vector in this null space too, then we may define an inner product between null rules as follows
(l-VU,
NV) =
(U, v) = f .?4,.,, E=() and we have a norm of a null rule and a ruleWe may furthermore define two null rules as orthogonal when (Aru, IVv) = O.
The idea of using inner products in connection with null rules was introduced by Espelid [7]. Note that S ~ ~ ~ (1 S; has dimension one and that the vector in this space will be unique, except for a scaling factor, and orthogonal to all vectors in S~.
Definition 3. Define a null rule N and an interpolator rule, based on the same set of points, as equally strong if the null rule is scaled such that II Arll, == IIQll,. A sequence of null rules NJ, j = 1,2,..., m, based on the same set of points, is said to be equally strong if each null rule in the sequence and the interpolator rule Q are equally strong.
THEOREM 3. Given a set of n -I- 1 distinct points. Then there exists a unique
interpo!atory rule Q and a unique (except for a – 1 scaling) sequence of orthogonal, equally strong null rules Nl, Nz, . . ., N. of decreasing degrees n–l, n–2 , . ...0. The subsequence Nl, N2, . . . . N~ forms an orthogonal basis for S~.
The last statement is obvious from the fact that the whole sequence is an orthogonal basis for S. and given N e S~ then
.
Since each null rule in the sequence Nm + ~, ., ., N. is orthogonal to any null rule in S~thenh, =0, i= m+l, . . ..n.
ACM Transactions on Mathematxal Software, Vol. 17, No. 2, June 1991
238 . J Berntsen and T. O. Espelid
In order to actually construct this equally strong sequence of null rules we may find Iil first solving (5). Suppose that we have found IVl, IV2, . . ., IVn. 1 then N may be found solving (6) in such a way that we know that N f sm. ~ and next use Gram–Schmidt’s orthogonalization process and scaling to find lvm .
Suppose now that the points { x,} are in an interval H~ of length h. If Q is the interpolatory rule based on these points then Q[ll = h and IIQ IIz = O(h).
Any equally strong null rule will therefore have II N l\~ = 0(h). In the next section we will find a relation between null rules and divided differences. We will use this relation to find to which order in h we can expect I N[ f] \ to be when used on the interval Hk. This observation will be used as a guide in the construction of the error estimating procedure based on such a sequence of equally strong null rules.
3. DIVIDED DIFFERENCES AND NULL RULES
In this section we will show the connection between null rules and divided differences. Let ~[ XO, xl, . . . . x ~] be a divided difference for the function f based on the set of distinct points { XO, xl, ., Xn} in the interval Hk of length h. The following formula is well known, for example see Isaacson and Keller [101,
(9)
If ~ is sufficiently smooth, then
f[xo>%j. ... xm] =f@@ ($n)lm!j (lo)
for a value of ~~ EHk. Now, (9) shows that a divided difference is a linear combination of function values and (10) shows that this linear combination of function values gives the value zero for all polynomials up to degree m – 1 and the value one for f(x) = x‘. Therefore the divided difference given in (9) is a null rule of degree m – 1. Let Q be the interpolatory rule based on x,, i=O, l7 ...> rz for n > m. Then the following null rule and Q is equally strong
n-m+, [f] = Cmzhm+’f[xo, xI,. . . . Xm]
N (11)
with the constant _c~, independent of h, chosen such that II ~.-=+l II~ = IIQll ~.
The null rules {N. _,+l} for i = m, m + 1,....n span Sn.l. Let N be any null rule of degree m – 1 then there exists a set of real numbers {h, ] such that
(12)
By applying (10) to each rule in this basis we find that each rule, N._ ~+ ~, in the sequence of equally strong and orthogonal null rules discussed in the
ACM Transactions on Mathematical Software, Vol 17, No 2, June 1991
previous section may be written
~.m+,[f] = f
h,h’+’cLfqgL)/i!N
L=m
(13)
where both c, and & are independent of h and all ( ,’s are in Hk.
THEOIREM 4. Given a set of n + 1 points in an interval Hh of length h. Let Qbe th~ interpolator rule based on these points and let N._ ~~, be the equally strong null rule of degree m – 1 in the unique (except for a – I scaling) sequence of orthogonal and equally strong null rules. If h is sufficiently small and f is sufficiently smooth then IN~_~+l[ f ] I = O(hrn+l).
Finally we remark that instead of solving (6), we may use the divided difference formulation to pick a null rule which is in S~ and not in Sm. ~ directly and then only the orthogonalization arid the scaling remain to be done. Similarly, NI is simply proportional to n XO, xl, . . . . x.].
4. SYMMETRIC RULES
The algorithm we are going to develop will be based on symmetric quadra- ture rules and symmetric null rules. The symmetry implies that the results of applying a quadrature rule and a null rule are independent of the direction of integration.
The points in an 2n + 1 point symmetric quadrature rule, defined on [–l,ll, maybe numbered xi for i= –n, – n+ 1,. ... n with
xo=o<x1<x2<. ..<2cn <1,
where .x_z= —xl, i = 1,2, ..., n. Symmetric rules have weights w_, = w, andu_, =u,, i=l,2, ..., n, for the quadrature rule and null rule respec- tively. Thus, on the interval [ – 1, 1], we have
l[f] = /1
f(x)dx= ~ w, f(xz) = Q[f],–1 ~=—~
N[f] = fi utf(x,).
L=—n
Furthermore we will consider, as quadrature rules, interpolator rules only.
This means that the quadrature rule’s degree of precision, d, will be at least 2n + 1.From Section 2 we know that it is possible to construct a sequence of 2n orthogonal and equally strong null rules of decreasing degrees 2n – 1,2n –2 , . ...0. Any symmetric null rule, N, has to have odd degree since
N[ X 2 n - 1] = O for m = 1,2,. . . due to symmetry. In order to construct a
sequence of n symmetric, orthogonal and equally strong null rules NI, NZ,. ... N. of degrees 2n – 1,2 n – 3, . . . . 1 respectively, we may solve
ACM Transactions on Mathematical Software, Vol 17, No 2, June 1991
240 . J Berntsen and T. O. Espelid
the following system of equations
/ 1/2 1 1 . . . 1
~ )[1 U.
o x: x; ““” x: u,
= Om,
. . . (14)
o xlZ(n–m) X22(n–m) . . . x2(n–m)
n Un
form =1,2, ..., n. This ( n – m + 1) x ( n + 1) coefficient matrix has rank n – m + 1and the dimension of its null space is therefore m. Let us denote this null space S~, for m = 1,2, . . . . n. Obviously SI c S2 c “ “ “ c S.. A vector u, which is a solution of (14), is then in S~ and to each such vector we have a symmetric null rule of degree at least 2( n – m) + 1. The inner product defined between null rules applies in the symmetric case too, thus we define two symmetric null rules as orthogonal if ( IVU, lVv) = O, where both u= u, and u., = v, fori=l,2, ..., n. Defining the inner product between tw~ vectors u and v in t$~ to be
(U, v) = u,+ 2 ~ UL”L L=l
gives us ( NU, Nv) = (u, v) for any to vectors u and v in Sm and their associated symmetric null rules. We may now give a theorem similar to Theorem 3 in the symmetric case.
THEOREM 5. Given a set of 2 n + 1 distinct points, symmetric in [– 1, 1].
Then there exists a unique (except for a – 1 scaling) sequence of symmetric, orthogonal and equally strong null rules Nl, Nz, . . . . N. of decreasing degrees 2n – l,2n – 3, . . . .1.The subsequence Nl, N2, . . . . N~ forms an orthogonal basis for all symmetric null rules of degree at least 2( n – m) + 1.
Instead of solving (14) we may, as in Section 2, use the divided difference f[x-n+n-l, X-n+m,. ~, Xn-m+, 1 to construct a symmetric null rule of degree 2(n – m) + 1.Given NI, N2, . . ., N~_l then use the Gram–Schmidt’s or- thogonalization process and scaling to find N~. Note that Theorem 4 implies that I N~[ f] I = 0(h2(n-m)+3).
5. A LOCAL ERROR ESTIMATING ALGORITHM
In this section we will develop a new local error estimating algorithm. This algorithm consists of several new elements, each one with its own motiva- tion. Let us start this section with summarizing our knowledge about error estimation in adaptive quadrature:
Asymptotic / nonasymptotic. The classical way to estimate the local error is based on the knowledge that asymptotically, that is for f sufficiently smooth and h sufficiently small, we have I 1 – Q I < I N[ f] 1.As we know, using I N[ f] I as a local error estimate in an adaptive algorithm may give an unreliable routine and it is therefore important to test whether we, in a particular interval, have a situation where the asymptotic theory may be
ACM TransactIons on Mathematical Software, Vol 17, No 2, June 1991
applied or not. We will design such a test using a sequence of local error
~+ ~ of decreasing asymptotic order in h. Thus we estimates El, Ez, . . . . E
may test if the computed values correspond to the expected order of this sequence. In the asymptotic case we will like to use an optimistic error estimate based both on the knowledge of the null rule’s and the basic rule’s degree of precision. In the nonasymptotic case we will base our error estimate on statistical considerations. The choice of value of K is, in addition, a question of reliability versus economy.
—Phase-factor. As Lyness and Kaganove [14] have pointed out, phase factor effects may ruin \ iV[ f] I as an estimate of the true error. This is true both in the asymptotic and the nonasymptotic case. In order to reduce such phase factor effects we construct each element in the sequence {Ej} by two linearly independent (in fact they are orthogonal) null rules.
–Rounding-noise. Local error estimation focus on estimating the local trun- cation error as reliable and economical as possible. Rounding errors may of course influence both the estimate of the local integral and the evaluation of all null rules. We will, as do Piessens et al. [17], define a certain noise level. All local error estimates below this noise level are redefined to be on the noise level.
Suppose we have a symmetric set of 2n + 1 distinct points. Given the sequence of symmetric, orthogonal and equally strong null rules. For a given subinterval and function f we may compute the estimates
eJ=NJ[ f], j=l,2, . . ..n.
Now pick el and ez, then it is easy to prove that
iVl and Nz span Sz (or rather the vectors associated with these null rules span $Z ). Therefore El is the greatest error an equally strong null rule from f$z may give. The idea of using El as an error estimate instead of I el I was introduced by Espelid [7]. El is less sensitive to accidental small values in
\ el I (phase factors, [14]). In these one-dimensional problems we obviously loose some degree of precision while doing this. Hopefully this price is worth paying in order to reduce phase factor effects. Now let us define
1/2
E~ = (e~~.1 + e~~) , j= ’1,2,..., [n/2].
Thus we replace the sequence { I e~I } by a new sequence {Ej}. This last sequence is based on null rules of degree of precision 2n – 3,2 n – 7, . ...2 n
– 4[n/2] + 1. However, now the sequence of null rules to be used are not fixed, but the null spaces in which we find the null rules ~re fixed: e.g., El is based on a null rule in Sz, Ez is based on a null rule in SA n ~~ , and so on.
We do not have to use our particular choice lVI and Nz to compute El since any orthogonal pair of null rules in Sz will do. This is true for all estimates in the sequence EJ. Each element E~ is constructed with phase factors in
ACM Transactions on Mathematical Software, Vol 17, No. 2, June 1991.
242 . J. Berntsenand T. O. Espelid
mind and is based on independent null rules and may therefore be considered as independent estimates of the same local error. Define the local error
Eo=[Q[f]
-I[f]l.
Let us use the K + 1 first elements in the sequence EJ, where ~ + 1 s [n/21.
Suppose furthermore that the local interval, say H~ of length h, is suffi- ciently small and the function f sufficiently smooth in this interval. Then we have, based on the previous section, that
E7 = 0(h2fn-2~1+3), j=l,2, . . .. K+l. (15)
In addition we have, based on the same assumptions as (15), a similar expression for the true error
E. =
o(h~+z),
(16)where the value of d, depending on the basic rule in question, will be in the range 2n + 1 s d < 4n + 1.The lower bound is achieved for, e.g., the Clen- shaw-Curtis rules and the upper bound is achieved for the Gauss rules. (15) and (16) imply that, for h sufficiently small and f sufficiently smooth in that interval, we have
EO<<EI << E2 << .“” << E~~l. (17) Let us define the reduction factors
r] = EJIEJ+I, j=l,2, . . ..K.
and
r= max r~.
J= 1,2, ,K
In view of (17) we see that a necessary test on whether h is small enough and
f
sufficiently smooth, is to check that r < rC,,t,C.l for a heuristic value of r ,,LtLC.l < 1. If this test is passed, we may apply an optimistic error estimate based on (15) and (16)E=cr”E1. (18)
Note that we have r = 0( h4) and observing that r. = EO /E1 is of order O(hd-(z.-s) ) we may choose a value of a in the range 1 s u s (d – (2 n – 3))/4 depending on the degree of optimism we want to put into this algo- rithm. We will use the term strong test on asymptotic behavior on the test r < rC,Lt,C,l. The test r ~r,tl,~l < r < 1 we will denote the weak test on asymp- totic behavior and when this test is passed we will suggest to damp the optimism, using an error estimate as follows
In order to give a continuous error estimating function we choose
The last test to consider is 1< r which we denote the test on nonasymptotic behavior. In this case we cannot use (15) and (16) as a guide in the construc- tion of the error estimate. We prefer in this case, following Espelid [71, to use a statistical argument, although we are well aware that the statistical assumptions, on which this argument is based, do not hold. Define ~ = ]/h,
ACM TransactIons on Mathematical Software, Vol 17, No 2,June 1991
where 1 is the value of the integral over a particular interval of length h. We may consider ~ the expected value of the function f over this interval. Let us define A f, = f( x,) – ~. Apply the basic rule Q on the function f(x) – ~ to get (now x, means the transformed abscissas in the new interval and the weights WZare the correctly scaled weights with x w, = h).
Q[f-f] = ~ w,Af(x,) = Q[f] -I[f].
~=—~
Assu~ming that the values A f( x,) are independent and each one is uniformly distributed in an interval [ -A f, A f], then the expected value of Q[ f] - 1[ t]
is O and the standard deviation proportional to IIQ II~. For any null rule IV the same argument may be applied, e.g.,
2J[f-~] = ~ u, Af(x,)=lJ[~] –N[t] =AT[f].
t=—~
Here the standard deviation will be proportional to IIN II~. By scaling the null rule such that IIQ II~ = II N II~, we therefore have the same standard deviation in the two expressions. Therefore the elements in the set E, will have the same expected size through this scaling. In the nonasymptotic region we consider any ordering of the elements of this set as equally probable. Therefore we may, by simply bad luck, have r s 1 even though we are in the nonasymptotic region. However, the probability that this will occur is then 1/( K + 1)!. In addition since we distinguish between weak and strong asymptotic behavior, we assume that we can usuaIly maintain the reliability in spite of such failures. We suggest the following error estimate in the nonasymptotic case
E=lo max EJ . (20)
J= 1,2,. ,K+-1
The constant 10 is chosen as a reasonable guarding constant to take care of the rare, but possible, situation that EO ;> max E~. Having chosen this constant equal to 10, and in order to have as smooth an error estimating function as possible, we choose
c’ = c rC:;t~Cal= 10, giving
c = lorl–”
CTttlcal
We summarize this section by giving a local error estimating algorithm.
A Local Error Estimating Algorithm Compute:
lVonasymptotic:
Weak asymptotic:
Strong asymptotic:
e =N~[f], j=l,2, . . ..2K+ 2;
~J = (e;,-1 + e~J)’’2~~,~, ~,;l..., K + 1;
r] = E1jEl~l, J = ~, * Y r=max J–1,22_ ~,Kr;
if r> lthen E= ~Omax~.l,z, ,K+l EJ
elseifArC,,~cC.l< r then E = 10 rEl else E = 10r&j$.aLraE1
endif
ACM Transactions on Mathematical Software, Vol 17, No. 2, June 1991
244 . J. Berntsen and T. O Espelid
Note: we may get r >1, even though we are in the asymptotic region, simply because the precision of the actual computer may influence the computations.
If the correct values of ~, and E, are very small, then they both may consist mainly of noise from the computations. Therefore r > 1 is possible even though we are in the region of asymptotic behavior.
treat this noise situation is discussed in Section 7.
6. ERROR ESTIMATION IN QAG IN QUADPACK
In order to illustrate the difference between the error
presented in this paper and the one used QAG, we present the last one in the framework used in this paper. The general purpose algorithm QAG is based on a 2n + 1point Gauss–Kronrod quadrature rule, Q, and an n point Gauss rule, Q’, as a rule of reference. In the error estimation they use, in addition,
Iiow to-discover and
estimating algorithm
Two independent error estimates are now computed
el=Q[f] -Q’[f]=lV, [f], (21)
2=Q[l&~l] = ~ w,l~(.x,)-;l. (22)
L=—n
As indicated (21) is a null rule of degree 2 n – 1 and it turns out that II~1 IIZ = IIQ II~ to three digits when n = 10. z given in (22) integrates constants to O and we have that E = 0( h2), if f has a bounded first deriva- tive in this local interval of length h. This operator therefore resembles a null rule of degree O. Define
r= lel l/Z.
The error estimating algorithm in QAG in QUADPACK may now be de- scribed as follows
The Local Error Estimating Algorithm in QAG
Compute: ed = Nl[ f];
f= Q[fl/h;
Z= Q[lf–~l];
r= lel l/i?;
Nona.symptotic: if r >A1/200 then E = ? Asymptotic: else E = ci-1’2 I ell
endif
Note that we have, for h sufficiently small and f sufficiently smooth, r = 0(h2” - l). Therefore we can expect, for h small enough, r < 1/200 (this value is chosen based on tuning the error estimate on a selection of test problems) and thus asymptotically that E = 0( h3’+ ~). Now the true error 11 = I Q[f ] – 1[ f ] I = 0( h3’ + 3). Therefore we may state that this estimate is slightly cautious by a factor h5/2. The constant is chosen to be c = 2003’2 in QAG. In principle the two algorithms have certain similarities. Both try to decide whether there is reason to believe that we are in the asymptotic region
ACM Transactions on Mathematical Software, Vol 17, No 2, June 1991
or not. Both use an optimistic and pessimistic error estimation depending on that decision. Phase factors may influence the QAG error estimation more than the error estimation suggested in Section 5, since they trust one null rule only in the optimistic case. ? is not likely to be effected by this phase factor due to the use of the absolute values.
7. NUMERICAL EXPERIMENTS.
Automatic routines for computing one dimensional integrals all have a FORTRAN calling sequence of the type
CALL INTGRL(A,B,F,EPS,MAXVLS,VALUE,ACCUR, . ..)
for an integral 1[ l’] over the interval from A to Il. The user supplies a requested relative error EPS and a maximum allowed number of integrand evaluations MAXVLS, and the routine returns an approximation to the integml, VALUE, with estimated absolute error ACCUR.
We want to test:
—the subroutine DQAG (this is the double precision version of QAG) from QUADPACK with routines DQK21 and DQK61 for estimating the integral and the error over each subregion. DQK21 applies a 21 point Gauss-Kron - rod rule. DQK61 applies a 61 point Gauss–Kronrod rule.
— six lmodifications of DQAG. DQK21 is replaced by routines that apply 21 point rules of type Gauss, Lobatto, Gauss -Kronrod and Clenshaw-Curtis respectively. DQK61 is replaced by routines that apply a 61 point Gauss rule, and a 61 point Gauss-Kronrod rule respectively.
The error estimating procedure is described in Section 5, and we have chosen
= 1/4,
‘Cr ZtLCd
a= (d- (2n - 3))/4, K=3,
where d = 41 for the 21 point Gauss rule, d = 39 for the Lobatto rule, d = 31 for the Gauss-Kronrod rule, d = 21 for the Clenshaw-Curtis rule, d = 121 for the 61 point Gauss rule and d = 91 for the 61 point Gauss- Kronrod rule.
In order to test whether we have reached the noise level or not, (see Section 5) we have also included a noise test similar to the one used in QAG [171, following our local error estimating algorithm:
Compute: noise = 50* Cn*RESABS;
The noise test: if El < noise and Eg < noise then E = noise
qndif
E = max( E, noise);
where En denotes the relative machine accuracy and RESABS is an approxi- mation to 1[ I ~ I]. The test on round off in the global approximations Q in DQAG is removed, and in addition to the replacement of DQK21 and DQK61, this is the only modification to DQAG that we have done in the last six routines.
ACM Transactions on Mathematical Software, Vol 17, No. 2, June 1991.
246 . J. Berntsen and T. O. Espelid
The testing technique that we have used is due to Lyness and Kaganove
[151. This technique is based on a selection of test families of integrands.
Each test family has a special attribute, a difficulty parameter and one or
more random parameters. The test families used in our experiments are
given below, and they are picked from Berntsen [11, Lyness and Kaganove
[151 and Sorevik [181.
Test families Attributes
1 ];(]x–hl)~’dx Singularity
2 /;fz(~) dx Discontinuous
{
o if x<A
where ~z( x) =
exp(az x) otherwise 3
/{w&a31
X– Al)d~4 /l 10”4/((x – A)z + 102”4) dx 5 /; ~f=llO”’/((x – ~,)z + 102”5) dx 6 /:2B(x – h)COS(B(X – h)2) dx
where B = 10”’/max(h2, (1 – A)z) In our experiments we have chosen the 1,...,6, to be (numbered from family 1 to 6):
CO function One peak Four peaks
Nonlinear oscillatory
difficulty parameters a,, i =
a = (-0.5,0.5,2.0, -4.0, - 2.0,3.0)
—
The random parameters, h (or X,, i = 1, . . . ,4, for test family 5), are picked randomly from the region of integration.
MAXVLS is in all experiments set equal to 50000. EPS = 10-1,...,10-7.
(For test family 1 we stop at 10-5.) For these values of EPS and for each test family we have asked all routines to compute the integrals for 1000 samples of random parameters, and in all cases the routines report that VALUE satisfies the error request. The experiments are run on an Alliant FX/8 (1 processor) in double precision and have taken approximately 22 hours of CPU time.
In Figures 1-6 we report on the performance of the different routines for each test family and for each error request. For each routine the four numbers from left to right are
(1) the average number of integrand evaluations;
(2) the average number of correct digits in VALUE;
(3) the number of cases (out of 1000) where the actual error is greater than the error request;
(4) the average number of wrong digits in the cases described above.
The number of wrong digits is computed by
Wrong digits = I – loglo( EPS) + log10(6.C,) 1,
where ~.Ct is the actual relative error in VALUE. For the last six subroutines we let the name of the quadrature rule define the name of the routine. The numbers in the parentheses indicate the number of evaluation points used by
ACM TransactIons on Mathematical Software, Vol 17, No 2, June 1991
m DQAG(21) Gauss-Kronrod(21)
~ 233 1.9 6 0.1 228 1.7 18 0.1
~()-2
519 2.9 19 0.5 502 2.7 49 0.2
~o-3 794 3.9 39 0.6 757 3.6 52 0.3
10-4
1388 5.5 9 1.5 1367 5.6 0 0.010-5 1981 6.6 0 0.0 1946 6.6 0 0.0
EPS Gauss(21) Lobatto(21) Clenshaw-Curtis (21)
~ 224 1.7 20 0.1 229 1.8 22 0.1 225 1.7 32 0.1
~()-2
495 2.7 53 0.2 503 2.7 61 0.3 497 2.7 89 0.3
~(y-3
762 3.6 -51 0.3 772 3.6 61 0.3 807 3.6 79 0.3
10-4
1125 5.0 0 0.0 1865 5.9 23 0.9 1929 6.1 20 0.910-5 1704 6.0 0 0.0 2473 7.0 0 0.0 2552 7.2 0 0.0
Fig. 1. Test family 1 (singularity)
DQAG(21). . . Gauss-Kronrod~
155 2.7 1 1.0
290 3.7 1 2.0
424 4.7 3 1.1
561 5.7 12 1.0
695 7.2 18 1.4
823 10.9 22 2.1 885 14.7 28 2.5
Gauss(21)
180 2.8 1 1.0
322 3.9 1 2.0
456 4.8 7 0.7
590 5.8 20 0.9
721 7.2 28 1.5
837 10.7 36 2.1 876 14.5 43 2.7
184 2.9 1 ‘ 1:0
323 3.9 1 2.0
458 4.8 3 1.1
595 5.9 12 1.0
729 7.5 18 1.4
844 11.1 22 2.1
885 14.7 28 2.5
Lobatto(21) Clenshaw-Curtis( 21)
189 2.9 0 0.0 182 2.8 0 0.0
326 3.9 0 0.0 325 3.8 0 0.0
466 4.8 0 0.0 462 4.8 0 0.0
606 5.9 0 0.0 603 5.8 0 0.0
741 6.9 0 0.0 740 6.8 0
0.0
870 8.0 0 0.0 876 8.0 0 0.0
983 9.2 0 0.0 1016 9.7 0 0.0
Fig. 2. Test family 2 (discontinuous).
the quadrature rule. For the routines applying 61 point quadrature rules, only the results for test family 6 are given.
8. COMMENTS AND CONCLUSIONS
We will first comment on the performance of the routines on the different test families of integrands.
— The singular problems. The integrals are difficult and all routines fail occasionally (from O to 9 per cent of the cases depending on the routine and
ACMTransactions on Mathematical Software, Vol. 17, No 2, June 1991
248 . J Berntsen and 1, 0. Espelid
EPS DQAG(21) Gauss- Kronrod(21)
10-1 68 3.9 0 0.0 39 3.6 0 0.0
~()-2 127 4.8 0 0.0 80 4.2 0 0.0
IO-3 192 5.8 2 0.1 146 5.1 4 0.2
10-4 262 6.8 8 0.3 217 6.1 2 0.1
10-5
326 7.7 12 0.5 283 7.1 1 0.310-6 396 8.7 13 0.6 350 8.1 2 l.o
~()-7
466 9.7 19 0.8 424 9.1 8 1.0EPS Gauss(21) Lobatto(21) Clenshaw-Curtis( 21)
1o-1 39 3.6 0 0.0 39 3.5 0 0.0 38 3.5 0 0.0
IO-2
79 4.2 0 0.0 78 4.2 0 0.0 77 4.1 0 0.010-3 143 5.1 1 0.2 148 5.2 6 0.2 146 5.1 7 0,2
10-4 214 6.1 0 0.0 218 6.2 2 0.1 216 6.1 5 0.1
10-5 280 7.2 1 0.3 285 7.1 1 0.4 284 7.1 1 0.4
10-6 349 8.1 5 0.9 356 8.1 0 0.0 356 8.1 2 0.3
10-7 419 9.1 13 1.1 430 9.2 3 0.1 426 9.1 3 0.3
F1g 4 Test family 4 (one peak) Fig 3. Test family 3 (Co function)
DQAG(21) Gauss-Kronrod(21)
482 7.3 7 1.0 464 6.1 0 0.0
521 8.8 2 2.0 493 7.6 0 0.0
552 10.2 0 0.0 528 9.0 0 0,0
577 11.4 0 0.0 557 10.5 0 0.0
602 12.3 0 0.0 574 11.2
0 0.0
634 13.0 0 0.0 606 12.2 0 0.0
682 13.4 0 0.0 642 12.8 0 0.0
Gauss(21) Lobatto(21)
461 6.2 2 0.8 466 6.0 3 1.0
490 7.8 0 0.0 496 7.6 0 0.0
523 9.2 0 0.0 530 9.1 0 0.0
550 10.7 0 0.0 554 10.5 0 0.0
566 11.5 0 0.0 572 11.4 0 0.0
585 12.1 0 0.0 595 12.2 0 0.0
613 12.9 0 0.0 623 12.9 0 0.0
497 6.9 0 0.0
533 7.9 0 0.0
566 8.9 0 0.0
598 9,9 0 0.0
641 10.9 0 0.0
701 11.7 0 0.0
ACM Transactions on Mathematical Software, Vol. 17, No 2, June 1991
-EFS DQAG(21)
ml 407 6.2 2 0.4
10-2 480 7.8 3 0.9
10-3 522 9.7 0 0.0
10-4 550 10.8 1 0,0 10-5 627 11.8 1 1.0 10-6 695 12.9 0 0.0
sl%%%F-
10-1
327 4.4 0 0.010-2 405 6.6 0 0,0
10-3 491 8.6 0 0.0
10-4 524 10.3 0 0.0 10-5 579 11.4 0 0.0 10-6 641 12.1 0 0.0 10-7 679 12.8 0 0.0
Gauss-Kronrod(21)
326 4.3 0 0.0
409 6,5 0 0,0
494 8.4 0 0.0
531 10.3 0 0.0 591 11.3 0 0.0 668 12.1 0 0.0 718 13.2 0 0.0
Lobatto(21)
326 4.2 0 0.0
412 6.4 0 0.0
495 8.3 0 0.0
.529 10.0 0 0.0 587 11.0 0 0.0 652 11.8 0 0.0 689 12.5 0 0.0 ~
Clenshaw-Curtis( 21)
327 4.1 1 0.1
413 6.2 2 0.3
494 7.6 0 0.0
545 9.0 0 0.0
622 9.9 0 0.0
720 10.9 0 0.0
773 11.8 0 0.0
Fig. 5. Test family 5 (four peaks).
the error request). The routines using Gauss and Gauss -Kronrod rules with our error estimate fail more frequently than DQAG. On the other hand the average numbers of wrong digits are greater for the latter routine, and the routine applying the Gauss rule is somewhat more effi- cient. The routines applying closed rules (Lobatto and Clenshaw-Curtis) are more expensive and also less reliable than the routines applying open rules (Gauss and Gauss –Kronrod).
—The discontinuous problems. The routines applying closed rules are supe- rior to the routines applying open rules with respect to reliability. The disposition of evaluation points seems to govern the reliability. DQAG and the routine applying the Gauss- Kronrod rule in combination with our error estimate thus have exactly the same numbers of failures.
—The Co problems. The routines applying our error estimate are all more reliable than DQAG, and at the same time more efficient. The routine applying the Lobatto rule seems to be somewhat more reliable than all other routines.
—The peak problems. On these problems all routines are very reliable. When we request only a few digits of accuracy, the routines tend to fail occasion- ally. However, when we combine the Gauss –Kronrod rule with our error estimate no failures are reported. The routines using our error estimate are also more efficient than DQAG (except for the routine applying Clen- shaw– Curtis in a few cases). We also note how the efficiency depends on the polynomial degree of the quadrature rule, the routine applying the Gauss rule being clearly most efficient.
ACMTransactions on Mathematical Software, Vol. 17, No. 2, June 1991.
250 . J Berntsen and T. O. Espehd
DQAG(21)
3519 11.1 26 1.0 3860 11.6 4 2.0 4208 11.8 0 0.0 4574 11.9 0 0.0 4985 11.9 0 0.0 5403 11.9 0 0.0 5852 12.0 0 0.0
Gauss(21) 3821 11.7 6 1.0 4169 11.9 0 0.0 4428 11.9 0 0.0 4635 11.9 0 0.0 4849 11.9 0 0.0 5069 11.9 0 0,0 5293 11.9 0 0.0
DQAG(61)
2439 11.5 16 1.0
2548 11.6 2 2.0
2641 11.7 1 3.0
2727 11.7 1 4.0
2815 11.7 0 0.0 2919 11.7 0 0.0 3022 11,8 0 0.0
Gauss-Kronrod(21 ) 3846 11.6 1 1.0 4206 11.8 0 0.0 4513 11.9 0 0.0 4813 11,9 0 0.0 5121 11.9 0 0.0 5445 11.9 0 0.0 5796 11.9 0 0.0
Lobatto(21) 3860 11.8 0 0.0 4206 11.8 0 0.0 4471 11.8 0 0.0 4696 11.9 0 0.0 4923 11.9 0 0.0 5157 11.9 0 0.0 5399 11.9 0 0.0 Gauss-Kronrod(61 ) 2393 11.6 5 1.0 2479 11.6 0 0.0 2557 11.7 0 0.0
2627 11.7 0 0.0
2692. 11.7 0 0.0 2774 11.7 0 0.0 3042 11.7 0 0.0
Clenshaw-Curtis( 21)
3874 7.1 0 0.0
4289 8.5 0 0.0
4755 9.5 0 0.0
5270 10.3 0 0.0 5831 11.2 0 0.0 5465 11.9 0 0.0 7165 12,1 0 0.0
Gauss(61) 2379 11.6 3 1.0 2432 11.6 0 0.0 2456 11.6 0 0.0 2470 11.6 0 0.0 2480 11.6 0 0.0 2500 11.6 0 0.0 2676 11.6 0 0.0
Fig 6 Test family 6 (nonlinear oscillatory)
— The nonlinear oscillatory problems. Also for these problems all routines are very reliable. As for the peak problems there are a few failures when we request only a few digits of accuracy. However, the numbers of failures are clearly less for the routines applying our error estimate. Again the efficiency of the routines increase with the degree of the quadrature rules, and by using routines that apply 61 point rules the numbers of integrand evaluations needed are almost halved.
From these more detailed considerations we conclude that our local error estimating algorithm seems to work satisfactorily with all the used quadra- ture rules. We should note that we have used the same heuristics for all rules. In addition, we have used as high a value of the parameter a as possible in each case (in agreement with the basic rule’s degree of precision).
In final implementations other choices can be made in order to achieve balance between efficiency and reliability.
ACMTransactions on Mathematical Software, Vol 17, No 2, June 1991.
We have also seen that this new local error estimating algorithm gives us the opportunity to benefit from a high polynomial degree of the quadrature rule. Therefore, we see no reason for not using Gauss rules instead of Gauss- Kronrod rules in automatic adaptive quadrature routines.
On discontinuous problems the closed rules give clearly more reliable software. The routine applying the Lobatto rule is in addition very efficient on all the other problems except for the singular problems. It seems therefore natural to include a routine based on the Lobatto rules in an adaptive software package.
The use of this new error estimating algorithm involves the computation of 2 K + 2 inner products in order to compute the necessary null rule approxi- mations. With the usual error estimate only one such inner product is needed However, the same integrand evaluations are used, and for most integral computations it is believed that the integrand evaluations dominate the total computation time. The inner product evaluations may also, on many computers, be vectorized.
We believe that the main achievement of this paper is that we, by combin- ing theoretical considerations and practical experience, have produced an error estimating algorithm which can be used with different underlying rules and seems to work satisfactorily with all. With the new algorithm one may benefit from using quadrature rules of high polynomial degree, and the tests indicate that the reliability of the routines applying the new error estimate compare favorably with the reliability of DQAG from QUADPACK.
ACKNOWLEDGMENTS
The authors want to thank James N. Lyness for helpful discussions and especially for pointing out that the theory concerning null rules should be given for a general set of n + 1 points and not only for the symmetric case.
Furthermore the authors want to thank the referees for many valuable suggestions and corrections.
REFERENCES
1. BERNTSEN) J. A test of some well known quadrature routines. Reports in Informatics 20, Dept. of Informatics, Univ. of Bergen, 1986.
2. BERNTSEN, J. Practical error estimation in adaptive multidimensional quadrature routines.
J. Comput. Appl. Math. 25 (1989), 327-340.
3. BmmmsEN, J., ESPELID, T. O., AND GENZ, A. A test of ADMINT. Reports in Informatics 31, Dept. of Informatics, Univ. of Bergen, 1988.
4. BERNTSEN, J., ESPELID, T. O., AND GENZ, A. An adaptive algorithm for the approximate calculation of multiple integrals. To appear in ACM Trans. Math. SoftuJ.
5. DAVIS, P. J., AND RABINOWITZ, P. Methods of Numer~cal Integration. Academic Press, 1984.
6. DE BOOR, C. On writing an automatic integration algorithm. In Mathematical Software., J.
R. Rice, Ed., Academic Press, 1971.
7. ESPELID, T. O. Integration rules, null rules and error estimation. Reports in Informatics 33, Dept. of Informatics, Univ. of Bergen, 1988.
8. ESPELID, T. O., AND S6REVIK, T. A discussion of a new error estimate for adaptive quadrature. BIT. 29 (1989), 283-294.
ACM Transactions on Mathematical Software, Vol. 17, No 2, June 1991.
252 . J Berntsen and T. O. Espelld
9 GENZ, A C,, AND M~LIKj A. A. An adaptive algorlthm for numerical integration over an N-dimensional rectangular region. J. Comput Appl Math. 6 (1980), 295-302.
10 ISAACSON, E., AND KELLER, H B. Analyszs of Numerzcal Methods North Oxford Academic, 1966
11. LAURIE, D. P Sharper error estimate in adaptive quadrature. BIT 23 (1983), 258-261 12. LAURIE, D. P. Practical error estimation in numerical integration J. Cornput. Appl.
Math. 12 & 13 (1985), 425-431.
13. LYNESS, J. N. Symmetric mtegratlon rules for hypercubes III Construction of integration rules using null rules Math. Comput. 19 (1965), 625-637
14.LYNESS, J. N., AND KA~ANOVE, J. J. Comments on the nature of automatic quadrature routines. ACM Trans. Math Softw 2, 1 (1976), 65-81.
15, LYNESS, J, N., AND KAC,ANOVE, J. J. A technique for comparing automatic quadrature routines Comput J 20 (1977), 170-177.
16.MCKEEMAN, W. M Certification of algorithm 145. Adaptive numerical integration by S,mpson’s rule. Commun. ACM 6 (1963), 167-168
17 PIESSENS, R., DE DONCIIER-KAPENGA, E., UBERHUBER, C W , AND KAHANER, D K QUAD- PACK, a subroutine package for automatic integration. Series m Cornputatzonal Math 1 Springer-Verlag, 1983.
18 SOREVIK, T. Reliable and efficient algorithms for adaptive quadrature Tech. Rep. Thesis for the degree Doctor Scientiarum, Dept. of Informatics, Univ. of Bergen, 1988.
19 VAN DOOREN, P., AND DE RIDDER, L. An adaptive algorithm for numerical integration over an N-dimensional cube J Comput. Appl Math. 2 (1976), 207-217
Received May 1989; revised January 1990; accepted January 1990
ACM TransactIons on Mathematical Software, Vol 17, No 2, June 1991