We are interested in producing estimates of the error (
)
when R[f] is a degree 2m+1 approximation.
It is normally assumed (see Genz and Malik [17]) that | E[f] | is bounded by an expression
for some null rule N of degree less than 2m+1.
Numerical experiments have shown (Berntsen [3, 5] and [15]) that a more reliable
routine has to use a more sophisticated error estimation procedure.
The error estimates produced by the algorithm are based on the approximations
delivered by the null rules (),
with a common scaling given by
when 
is the region
of integration.
Lyness and Kaganove [20] have shown
how ``phase factors'' may cause severe problems for the
error estimation. In order to avoid the problems with ``phase factors''
we therefore use pairs of linearly independent null rules 
and
and select for each subregion the error estimate
1
will then be the greatest error estimate produced
by a null rule in the
space spanned by and with norm satisfying ()
(Espelid [14]).
Furthermore the null rules 
, for i=1,2,3, are of degree 2m-1, 2m-3
and 2m-5 respectively.
In [14] it is shown that the number of values of 
that
have to be considered when computing the maximum in (7), is equal to the
number of generators of the integration rule applied, and therefore the work
involved when computing is almost negligible.
In order to check the asymptotic assumption behind ()
we also check the inequalities([5, 6])
for suitable constants 
.
If the tests () are passed, we choose
otherwise we choose
where 
and 
are heuristic constants that are selected
to achieve a reasonable balance between reliability and efficiency.
In globally adaptive quadrature algorithms based on bisection of subregions, we build up binary trees of subregions as the computation proceeds. The final estimates of the integral and the error are based on the corresponding estimates from the subregions represented by the leaves of such a binary tree. Going from one level in the tree to the next, all information from the ``old'' function evaluations are often neglected. At least this is the case for the QUADPACK (Piessens et al. [23]) routine QAG and multidimensional ADAPT like routines ([17] and [25]). Any error estimate based on a finite number of integrand evaluations may be unreliable, and even though the actual error is very great, the estimated error may be negligible. The result of this may therefore be that a problem spot discovered through a large error estimate on one level of the tree, may be ``lost'' on the next level.
In order to try to avoid this problem we use two-level
error estimates. Let 
be an approximation to the integral over a
given region H, and 
for j = 1, 2, be two approximations over
the two children subregions we get by bisecting H.
We then get a two level error estimate by the expression
This error estimate is combined with the local error estimates 
and 
to produce the final estimated errors, 
and 
, over each of the two new
regions (see Sørevik [24] and [5]),
We summarize the error estimation procedure as follows:
For each of the children subregions we first compute
and then we impose the test
if 
and
then
else
endif,
where and are local error estimates
over the children subregions.
A two level error estimate 
is computed by
The final estimates
are then computed.
The table below shows the values of the heuristic constants used by the
different sets of rules in DCUHRE, the FORTRAN implementation of this algorithm.
All heuristic constants are selected to achieve a reasonable level of
reliability
without increasing the cost too much.
Our technical report [9] provides numerical evidence on how DCUHRE
( note that the original name of this routine was ADMINT)
performs with these choices of heuristic constants.