Introduction

The problem considered in this paper is the numerical evaluation of integrals in the form

$$I[{\bf f}] \ =\ \int_T {\bf f}({\bf x}) dT,$$

where ${\bf x}$ is an $n$-vector, ${\bf f}$ is $l$-vector and $T$ is a collection of $m$ $n$-simplices. There has been only limited work on practical algorithms for this general problem. Most research has considered the case $n = 2$, where $T$ is a triangle. For this case there has been work on the development of integration rules (reviewed in the paper by Lyness and Cools [25]), and algorithms ([3,12,20,21,24] ). For the case $n = 3$, where $T$ is a tetrahedron, there has been more limited rule ([8,30]) and algorithm [4] development work. For $n > 3$, the main rule development papers are by Silvester [29], Grundmann and Möller [19], Keast [23] and Lyness and Genz [26]. Other rules are described and referenced in the books by Stroud [30] and Engels [13] and the papers by Cools and Rabinowitz [8] and Cools [9]. The only general algorithm widely available [28] is an automatic algorithm that uses Grundmann's and Möller's rules, although an experimental algorithm was developed by Kahaner and Wells [22]. A brief description of an earlier version of the algorithm described in this paper was given by Genz [17]. The purposes of this paper are to give a detailed description of that algorithm, along with several significant additions made since the original algorithm was described, and to describe the implementation and testing of the algorithm as a part of CUBPACK [10].

In the following sections the details of the algorithm for the general problem will be described. We employ a globally adaptive algorithm that has been used extensively so only a brief description of this general algorithm is given here.

The globally adaptive algorithm uses successive refinements or subdivisions of $T$, where each subdivision is used to provide a better approximation to $I[{\bf f}]$. These subdivisions are designed to dynamically concentrate the computational work in the subregions of $T$ where the integrand ${\bf f}({\bf x})$ is most irregular, and thus to adapt to the behavior of the integrand. The general structure of the globally adaptive algorithm consists of a sequence of stages. Each stage has the following five main steps:

1. Select a subregion with largest estimated error from the current set of subregions.
2. Divide the selected subregion.
3. Apply a local cubature rule to any new subregions.
4. Update the subregion set.
5. Update the global integral and error estimates, and check for termination.

The initial subregion set for the algorithm is the original collection of simplices $T$. The required input for such an algorithm is $T$, the integrand ${\bf f}({\bf x})$, a limit on the number of ${\bf f}$ values allowed, and a requested error tolerance $\epsilon$. The algorithm terminates when the estimated global error is less than $\epsilon$ or further subdivision would require too many ${\bf f}$ values.

Algorithms of this type have been extensively used for numerical cubature over hyper-rectangular regions (see for example [1]), and also for triangles (e.g. [3]) and tetrahedra [4]. Because procedures for the steps 1, 4 and 5 are the same for all regions, these have been implemented as part of the basic CUBPACK framework [6,10] and will not be discussed here. Our implementation allows the user to make use of specially designed high degree rules if $n = 2$ or $n = 3$. These rules and their error estimates are essentially the same as those used in the TOMS algorithms 706 [3] and 720 [4], and their inclusion in CUBPACK has already been described elsewhere [10]. The rest of this paper will focus on the choice of the local cubature rules, error estimation, subdivision strategy, testing and interface with the CUBPACK framework for the general simplex algorithm.