Implementation in CUBPACK and Testing

The algorithm described in the previous sections has been implemented in Fortran 95 as a part of CUBPACK. A standard call to use the algorithm is a call to the main CUBPACK integration routine CUBATR in the form:

USE PRECISION_MODEL
USE CUI
...
CALL CUBATR( $n$, $l$, ${\bf f}$, $m$, V, RGTYPE, VALUE, ERROR, INFORM,
FUNVLS, EPSABS, EPSREL, RESTART,
MINPTS, MAXPTS, KEY, JOB, TUNE )

CUBATR computes approximations in the vector VALUE to the vector integral

$$I[{\bf f}] = \int_T {\bf f}\ d x_1d x_2\cdots dx_n$$

where ${\bf f}$ is an $l$-vector of integrands and $T$ is collection of $m$ $n$-dimensional simplices with vertices given in a 3-dimensional array V. CUBATR attempts to compute VALUE with

$\vert$ I[$f_k]-$VALUE($k$) $\vert$ $\leq$ $\max$( EPSABS, EPSREL$\vert I[f_k]\vert$ ), for $k = 1, ..., l$,

where EPSABS and EPSREL are absolute and relative tolerances. A complete specification for all of the CUBATR parameters is given in [10]. In Appendix A we describe those that are specific to the simplex.

The simplex part of CUBPACK was tested using a combination of techniques. The primary testing used a modified version of the package developed by Genz [16]. This package was designed for testing cubature algorithms for hyper-rectangles. We use two applications of the package. In the first application, we used the same cubature problems generated by the package but divided the $n$-dimensional hyper-rectangular cubature region into $n!$ equal volume simplices in a standard way [26] and applied CUBATR to the resulting collection of pieces. Because of the increasingly larger numbers of simplices, we found that testing in this manner was feasible only for $n < 7$. In the second application of the package, we used the same cubature problems generated by the package but we transformed the standard unit simplex to the hyper-rectangle for each problem, using a constant Jacobian transformation described by Fang and Wang [15]: If the hyper-rectangle $H$ for a particular test problem has dimensions $[a_1,b_1]\times[a_2,b_2]\times \cdots \times [a_n,b_n]$, then the transformation from a point ${\bf x}\in T$ to a point ${\bf y}\in H$ is defined by $y_i = a_i + (b_i-a_i)((1-\sum_{j=i}^nx_j)/(1-\sum_{j=i+1}^nx_j))^i$, for $i=1,2, \dots, n$, with the Jacobian given by $n!\prod_{i=1}^n(b_i-a_i)$. The test families that we used are given in Table 3.

Table 3: Test Families for Hyper-Cubes

Test Family	Attribute
$f_1(x,y) = \cos(2\pi\beta_1+\sum_{i=1}^n\alpha_iy_i)$	Oscillatory
$f_2(x,y) = \prod_{i=1}^n(\alpha_i^{-2}+(y_i-\beta_i)^2)^{-1}$	Internal Peak
$f_3(x,y) = (1+\sum_{i=1}^n\alpha_iy_i)^{-(n+1)}$	Corner Peak
$f_4(x,y) = \exp(-\sum_{i=1}^n\alpha_i^{2}(y_i-\beta_i)^{2})$	Gaussian
$f_5(x,y) = \exp(-\sum_{i=1}^n\alpha_i\vert y_i-\beta_i\vert)$	$C_0$ Function

The integration region for each family is the unit $n$-cube $[0,1]^n$. To determine a particular test problem, the $\beta_i$ parameters are chosen uniformly random from $[0,1]$. Then, an additional set of parameters $\alpha'_i$ are also chosen uniformly random from $[0,1]$. These parameters are then all scaled by a constant $c$ with $\alpha_i = c\alpha'_i$, where $c$ is determined by the condition $n^{e_j}\sum_{i=1}^n\alpha_i = d_j$, with $e_j$ and $d_j$ fixed for each test family $j$. The tests that we carried out used $(e_1, e_2, e_3, e_4,e_5) = (1.5,2,2,1,2)$ and $(d_1, d_2, d_3, d_4, d_5) = (100,500,100,100,200)$. A series of tests were run using both applications of the package and results from the test package were used to determine what we believe are good values for the error estimation parameter $C_e$ as a function of the tuning parameter $C_t$. We also used the tests to tune the subdivision procedure.

Table 4: CUBATR Test Results for n = 7, KEY = 3 and TUNE = 0

${\bf f}$ Type	Estimated Digits	Actual Digits	Reliability	Incorrect Digits
Oscillatory	( 3.4, 3.7)	( 3.1, 3.5)	0.26	(0.2,0.3)
Product peak	( 2.6, 2.7)	( 2.7, 3.0)	0.84	(0.0,0.0)
Corner peak	( 3.3, 3.4)	( 3.2, 3.3)	0.36	(0.0,0.1)
Gaussian	( 2.2, 2.2)	( 2.8, 3.2)	0.98	(0.0,0.0)
C0 function	( 1.8, 2.0)	( 2.1, 2.5)	0.78	(0.0,0.0)

Tables 4 and 5 contain results from two sample test runs (after tuning) with the second application of the test package, with each run using 50 randomly chosen integrands from each of the five integrand types, for $n=7$ and KEY = 3 (degree 7). For each sample, CUBATR was called, with a maximum of 343000 ${\bf f}$ values allowed, and a requested relative accuracy set at $10^{-10}$ to force the maximum number of ${\bf f}$ values to be used. The entries in the tables of the form (a,b) are 97% confidence intervals for the respective sample medians. In these tables, the number of digits is defined to be $-\log_{10}$ of the relative error (estimated or actual). Incorrect digits are defined to the difference between estimated and actual digits, for those cases where the estimated number of correct digits is greater than the actual number of correct digits. The numbers in the reliability column give the fraction of the times when the number of estimated correct digits was less than the number of actual correct digits.

Table 5: CUBATR Test Results for n = 7, KEY = 3 and TUNE = 1

${\bf f}$ Type	Estimated Digits	Actual Digits	Reliability	Incorrect Digits
Oscillatory	( 2.6, 2.8)	( 3.2, 3.6)	1.00	(0.0,0.0)
Product peak	( 1.5, 1.7)	( 2.6, 2.9)	1.00	(0.0,0.0)
Corner peak	( 2.4, 2.5)	( 3.2, 3.4)	1.00	(0.0,0.0)
Gaussian	( 1.3, 1.4)	( 2.6, 2.9)	1.00	(0.0,0.0)
C0 function	( 0.7, 1.0)	( 2.0, 2.2)	1.00	(0.0,0.0)

The results in Tables 4 and 5 are typical of our tests for $1 < n \leq 8$. For high levels of reliability the user must be prepared to accept a more conservative error estimate. We also did extensive testing of the different rules with polynomials for all key values for $1 < n \leq 10$ to verify that we had a correct implementation. Test programs have been successfully run on SUN (f95), Linux (NAGWare f95) and Dec f90.

We did not complete any significant testing for singular integrals. The software is not expected to fail for singular integrals, but our experience has been that singular integrals can be more efficiently computed either by using appropriate transformations to remove the singularities, or by using extrapolation methods.