The Local Cubature Rules and Error Estimators

In this section we restrict our discussion to integrands with a single component $f({\bf x})$. For a vector integrand ${\bf f}({\bf x})$, the cubature rules discussed in this section can be applied to each of the components in ${\bf f}$. The local cubature rules that we use for the simplex algorithm take the form

$$B [f] \ =\ \sum_{i=1}^{N} w_i f( {\bf x}_i) \ \simeq \ \int_R f ( {\bf x}) d{\bf x},$$

for a chosen simplex subregion $R$. With these rules, the points ${\bf x}_i$ and weights $w_i$ are chosen to make the rule exact for all polynomials of degree $d$ or less, for some fixed $d$ and $n$. The rules that are used in practical calculations typically have degrees in the range 5-13, with $N$ a polynomial of degree $n$ in $d$. The points and weights are found for some standard $n$-simplex, then the linear structure of the rule allows an affine transformation, preserving the polynomial degree, to any other finite $n$-simplex. The use of rules that are exact for low-degree polynomials ensures rapid convergence if the subdivision is fine enough so that the integrand is reasonably smooth locally.

The local rules used by the simplex algorithm are from Grundmann and Möller [19]. Let $G_s$ denote a degree $2s+1$ Grundmann and Möller rule,

$$\begin{eqnarray*} G_s[f] &=& \sum_{i=0}^s2^{-2s}(-1)^i\frac{(2s+1+n-2i)^{2s+1}}{... ...}{2s+1+n-2i},...,\frac{2\beta_n+1}{2s+1+n-2i}} \right) \right). \end{eqnarray*}$$

The inner sum notation $\sum f(({\bf y}))$ denotes a sum over all unique permutations of the $(n+1)$-vector ${\bf y}$. Each $f$ value in the sum uses only the last $n$ components of ${\bf y}$. The rules as given are for an integral over the standard unit $n$-simplex defined by $\sum_{i=1}^n x_i \leq 1$ for all $x_i \geq 0$. A rule $G_s$ requires ${n+s+1}\choose{s}$ $f$ values, and includes embedded rules of degrees $2s-1$, $2s-3$, ..., 1.

The local error estimators used by our simplex cubature algorithm are ``Norwegian'' error estimators (see Berntsen and Espelid [2]). They require a sequence of null rules, which are combined to produce the final error estimate. Our simplex algorithm uses two sets of null rules. One set of null rules is constructed from the difference between a degree $2s+1$ Grundmann and Möller rule and a lower degree Grundmann and Möller rule. Define the degree $2i+1$ null rule $N_i$ by $N_i[f] = G_s[f] - G_i[f]$, for $0 \leq i < s$. The second set of null rules $\hat{N}_i[f]$ is defined by $\hat{N}_i[f] = G_s[f] - L_i[f]$, where $L_i$ is some degree $2i+1$ rule, different from $G_i[f]$, for $0 \leq i < s$. Our implementation allows $0 < s \leq 4$ and the $L$ rules implemented are: a) ($i = 3$) a degree 7 Mysovskikh [27] rule that requires $3(n+1)(n+2)/2$ additional points, b) ($i = 2$) a degree 5 Stroud rule [31] that requires $(n+1)(n+2)$ additional points, and c) and d) special ($i = 1$) degree 3 and ($i = 0$) degree 1 rules, respectively, that use subsets of the points of the degree 5 Stroud rule points.

The degree 7 Mysovskikh rule is described in a paper in Russian that is not easily accessible, so we give some details here of our implementation of this rule. The Mysovskikh rule structure is similar to that of the Grundmann and Möller degree 7 rule. Let $\{({\bf y}_i)\}$ be the sequence of $({\bf y})$'s for the $f(({\bf y}))$ values for the Grundmann and Möller rule, and $\{({\bf y}'_i)\}$ be the corresponding sequence for the Mysovskikh rule, and let the sequences $\{W_i\}$ and $\{W'_i\}$ be the corresponding sequences for the weights, the two rules are compared in Table 1.

Table 1: Grundmann and Möller and Mysovskikh Generators and Weights

$i$	${\bf y}_i$	$W_i$	${\bf y}'_i$	$W'_i$
1	$(\frac{1}{n+1},\frac{1}{n+1},\cdots,\frac{1}{n+1})$	$-\frac{\ (n+1)^7}{2^6(n+4)!3!}$	$= {\bf y}_1$	$W'_1$
2	$(\frac{3}{n+3},\frac{1}{n+3},\cdots,\frac{1}{n+3})$	$\ \frac{\ (n+3)^7}{2^6(n+5)!2!}$	$(1-n\alpha_1,\alpha_1,\cdots,\alpha_1)$	$W'_2$
3	$(\frac{5}{n+5},\frac{1}{n+5},\cdots,\frac{1}{n+5})$	$-\frac{\ (n+5)^7}{2^6(n+6)!}$	$(1-n\alpha_2,\alpha_2,\cdots,\alpha_2)$	$W'_3$
4	$(\frac{3}{n+5},\frac{3}{n+5},\frac{1}{n+5},\cdots,\frac{1}{n+5})$	$-\frac{\ (n+5)^7}{2^6(n+6)!}$	$= {\bf y}_4$	$= W_4$
5	$(\frac{7}{n+7},\frac{1}{n+7},\cdots,\frac{1}{n+7})$	$\ \frac{\ \ \ (n+7)^7}{2^6(n+7)!}$	$(1-n\alpha_3,\alpha_3,\cdots,\alpha_3)$	$W'_5$
6	$(\frac{5}{n+7},\frac{3}{n+7},\frac{1}{n+7},\cdots,\frac{1}{n+7})$	$\ \frac{\ \ \ (n+7)^7}{2^6(n+7)!}$	$(\frac{4}{n+7},\frac{4}{n+7},\frac{1}{n+7},\cdots,\frac{1}{n+7})$	$\ \frac{10(n+7)^7}{3^6(n+7)!}$
7	$(\frac{3}{n+7},\frac{3}{n+7},\frac{3}{n+7},\frac{1}{n+7},\cdots, \frac{1}{n+7})$	$\ \frac{\ \ \ (n+7)^7}{2^6(n+7)!}$	$= {\bf y}_7$	$= W_7$
8	-	-	$(\frac{11}{2(n+7)},\frac{5}{2(n+7)},\frac{1}{n+7},\cdots,\frac{1}{n+7})$	$\ \frac{2^6(n+7)^7}{3^8(n+7)!}$

Using $a_i = (n+1)\alpha_i-1$ for $i = 1, 2, 3$, Mysovskikh shows that the $a_i$'s are the zeros of the polynomial

$$\begin{eqnarray*} p(z) &=& -144( 142528 + n( 23073 - 115n ) )\\ & & -12( 66905... ...\\ & & -(n+7)(6386660+n(4411997+n(951821+n(61659-665n))))z^3 , \end{eqnarray*}$$

and

$$W'_2=\frac{U_7-(a_2+a_3)U_6+a_2a_3U_5}{a_1^5(a_1^2-(a_2+a_3)a_1+a_2a_3)},$$

$$W'_3=\frac{U_7-(a_1+a_3)U_6+a_1a_3U_5}{a_2^5(a_2^2-(a_1+a_3)a_2+a_1a_3)},$$

$$W'_5=\frac{U_7-(a_2+a_1)U_6+a_2a_1U_5}{a_3^5(a_3^2-(a_2+a_1)a_3+a_2a_1)},$$

where

$$U_5 = -\frac{6^3( 52212 - n( 6353 + n( 1934 - 27n ) ) )}{23328(n+6)!},$$

$$U_6 = \frac{6^4( 7884 - n( 1541 - 9n ) )}{23328(n+6)!},$$

$$U_7 = -\frac{6^5( 8292 - n( 1139 - 3n ) )}{23328(n+7)!}.$$

Mysovskikh's paper contains a table of $\alpha_i$'s and $W'_2$, $W'_3$ and $W'_5$ for $4 \leq n \leq 20$, accurate to eight decimal digits. We verified that the $a_i$ are real and distinct certainly for $n \leq 100$. In our implementation, the $a_i$'s are computed by the closed expression for the roots of a cubic derived by Cardan and Tartaglia in the 16th century. After finding accurate values for the $\alpha_i$'s and $W'_2$, $W'_3$ and $W'_5$, $W'_1$ is computed using

$$W'_1=\frac{1}{n!}-(n+1)(W'_2+W'_3+W'_5+ n(W'_4+W'_6 + 2W'_8 + (n-1)W'_7/3)/2).$$

The weights for the degree 5 Stroud rule are not given explicitly in Stroud's paper, but we derived some moderately simple expressions for these weights. The generators for this rule are $(\frac{1}{n+1},\frac{1}{n+1},\cdots,\frac{1}{n+1})$, and $(1-nr_i,r_i,\cdots,r_i)$, for $i=1,2$, and $((1-(n-1)u_i)/2,(1-(n-1))u_i/2,u_i,\cdots,u_i)$, for $i=1,2$, with respective weights $S_1$, $S_2$, $S_3$, $S_4$ and $S_5$, where $r_1=\frac{n+4-\sqrt{15}}{n^2+8n+1}$, $r_2=\frac{n+4+\sqrt{15}}{n^2+8n+1}$ $u_1=\frac{n+7+2\sqrt{15}}{n^2+14n-11}$, and $u_2=\frac{n+7-2\sqrt{15}}{n^2+14n-11}$. If we define $\lambda_i = 1-(n+1)r_i$ and $\delta_i = (1-(n+1)u_i)/2$ for $i=1,2$, then

$$S_2 = \frac{2(27-n)-\lambda_2(13-n)(n+5)}{\lambda_1^4(\lambd... ...\lambda_1(13-n)(n+5)}{\lambda_2^4(\lambda_2-\lambda_1)(n+5)!},$$

$$S_4 = \frac{2-\delta_2(n+5)}{\delta_1^4(\delta_1-\delta_2)(n... ...= \frac{2-\delta_1(n+5)}{\delta_2^4(\delta_2-\delta_1)(n+5)!},$$

$$S_1 = \frac{1}{n!}-(n+1)(S_2+S_3+ n(S_4+S_5)/2).$$

The special degree 3 rule uses the first three of the degree 5 Stroud rule generators, with weights

$$S'_2 = \frac{2-\lambda_2(n+3)}{\lambda_1^2(\lambda_1-\lambda... ...bda_2-\lambda_1)(n+3)!},\ S'_1 = \frac{1}{n!}-(n+1)(S_2+S_3).$$

Finally, the special degree 1 rule uses only the second of the degree 5 Stroud rule generators, with weight $S''_2 = \frac{1}{(n+1)!}$.

To test that the implementation of these rules was correct we checked for the exact (to machine precision) integration of monomials of the appropriate degrees. We were also able to accurately reproduce Mysovskikh's table.

The actual error estimators require the use of scaled orthogonal null rules, so we define the sequence $\{e_i\}_{i=1}^{2s}$ to be the set obtained by orthogonalizing and scaling (using Berntsen's and Espelid's procedure) the sequence $\hat{N}_{s-1}, N_{s-1},...,\hat{N}_0, N_0$. We also define related quantities $E_i$ and reduction factors $r_i$ required by the Norwegian algorithm, by $E_i = \sqrt{e_{2i}^2 + e_{2i-1}^2}$, for $i = 1, 2, ..., s$ and $r_i = E_i/E_{i+1}$, for $i = 1, 2, ..., s-1$. Finally, we define $r = \max(r_i)$, $\bar{E} = \max(E_i)$ and our final error estimate $E[f]$ by

$$E[f] = \left\{ \begin{array}{cl} C_e ( C_t\bar{E} + (1-C_t)E... ... r \geq 1 \\ rC_e E_s & \mbox{ if } r < 1 \end{array} \right.$$

where $C_e = 1 + 7 C_t$, and $0 \leq C_t \leq 1$. $C_t$ is a user selected tuning parameter. For $C_t = 1$, a very conservative error estimate is produced, and this is the most reliable case. For $C_t = 0$, a liberal error estimate is produced.

When the integrand is a vector, the error estimate $E[f]$ is applied to each component of ${\bf f}$, and the result is an $l$-vector of error estimates for the estimated integral of ${\bf f}$ over the selected subregion. In order to make decisions about possible further subdivisions, the globally adaptive algorithm requires a scalar measure of the error for that subregion. For our simplex algorithm, we use $\vert\vert E[{\bf f}]\vert\vert _\infty$ as an overall measure of the error for a selected subregion. This quantity is used for subsequent subdivision decisions by the region processor part of the globally adaptive algorithm. When the algorithm terminates, the error estimate returned for each component of ${\bf f}$ is the sum, taken over all of the subregions in the final subregion list, of the estimated subregion errors for that component.

In the Fortran 95 implementation of the algorithm the user is allowed to choose $s =$ 1, 2, 3 or 4 (degree $d =$ 3, 5, 7 or 9). The costs, in terms of ${\bf f}$ values, for computing the rules and error estimates are given in Table 2, for $1 < n \leq 10$.

Table 2: $\bf f$ Values Needed for Local Rules and Error Estimators
[2mm]

$d \backslash n$	2	3	4	5	6	7	8	9	10
3	7	9	11	13	15	17	19	21	23
5	16	23	31	40	50	61	73	86	100
7	32	49	86	126	176	237	310	396	496
9	65	114	201	315	470	675	940	1276	1695

The degree 3 and 5 rules are not recommended for most problems, because the error estimates are not as reliable as the error estimates for the degree 7 and 9 rules, which use more integrand values. Two special rules are also provided for $n = 2$ and $n = 3$. For $n = 2$, a 37 point degree 13 rule [3] is available, and for $n = 3$, a 43 point degree 8 rule [4] is available.