The Subdivision Method

At each stage in the globally adaptive algorithm a selected subregion must be subdivided. A simple natural subdivision method is to cut at the midpoint of each edge, but this produces $2^n$ pieces. Another problem with this subdivision is that it is not as adaptive as the method we use because this subdivision is done without any analysis of the integrand, even though the error for the integral over a selected subregion is often due to irregularity of the integrand in only a small number of directions. Our subdivision strategy, which uses a division of the largest error subregion into at most four new pieces, and which takes account of differences in integrand behavior in different directions, allows the algorithm to proceed from one stage to the next in a controlled manner.

The subdivision procedure that we use, is a modified version of a procedures first described in [17] and further developed (for $n$ = 2 only) in [14]. Once a subregion has been selected for subdivision, the globally adaptive algorithm used by CUBPACK will recommend a subdivision into at most 2, 3 or 4 pieces, depending on the current progress of the integration. Our subdivision procedure then divides the subregion by cutting one, two or three edges of the selected subregion to produce a 2-division, 3-division or 4-division of the selected subregion, respectively. An $n$-simplex has $n(n+1)/2$ edge directions, and our algorithm chooses subdivision directions from these directions.

In order for the algorithm to be efficient, a method is needed for selecting good edges for subdivision, and therefore some measure of integrand irregularity is needed. A popular measure of integrand irregularity that has been successfully used with adaptive algorithms for hyper-rectangles is a fourth difference of the integrand. We follow this approach, using modifications of the methods described in [17] and [14]. Our simplex algorithm uses fourth differences centered at the centroid of the selected simplex. Let ${\bf v}_{k,0} ,\ {\bf v}_{k,1} ,\ ... ,\ {\bf v}_{k,n}$ be the vertices of the current largest error simplex subregion $T_k$, and let the edge directions be given by ${\bf d}_{i,j} = {\bf v}_{k,j}-{\bf v}_{k,i}$, for $0 \leq i < j \leq n$. Now define $f_{i,j} ( \alpha ) = f({\bf c}+ \alpha {\bf d}_{i,j}/(5(n+1)))$, with ${\bf c}= \sum_{i=0}^n {\bf v}_{k,i}/(n+1)$ (the centroid of $T_k$). Then a fourth difference operator for the $(i,j)^{th}$ direction is given by

$$D_{i,j}(f)\ =\ \vert\vert{\bf d}_{i,j}\vert\vert _1\ \vert ... ...j}(0)-4(f_{i,j}(2)+f_{i,j}(-2))+(f_{i,j}(4)+f_{i,j}(-4))\vert.$$

The scaling factor $\vert\vert{\bf d}_{i,j}\vert\vert _1$ is used to provide some bias for division of very long edges. All of the points where $f({\bf x})$ is computed for these differences lie within $T_k$. Our general algorithm is designed to allow for vector integrands ${\bf f}$. When ${\bf f}$ has more than one component, we define $D_{i,j} = \vert\vert D_{i,j}({\bf f})\vert\vert _1$.

The edges for subdivision are edges $(i,j)$ where $D_{i,j}$ is large. Let $D_{i_s,j_s} = \max_{i<j}(D_{i,j})$. If a 2-division has been recommended, then let ${\bf v}_c \ =\ ({\bf v}_{k,i_s}+{\bf v}_{k,j_s})/2$. Two new subregions are produced from $T_k$ that have the same vertices as $T_k$ except that the $i_s$ and $j_s$ vertices are respectively replaced by ${\bf v}_c$. For a 3- or 4-subdivision let $D_{i_t,j_t}=\max_{i<j, (i,j)\neq(i_s,j_s)}(D_{i,j})$. If a 4-division has been recommended, and $D_{i_t,j_t} > D_{i_s,j_s}/2$ (the two edges with largest $D$ values have similar $D$ values), then our algorithm first divides $T_k$ into two pieces using the algorithm for a 2-division and then halves each of the two new pieces by bisection of the $(i_t,j_t)$ edge for each piece. If either a 3-division has been recommended, or a 4-division is recommended but $D_{i_t,j_t} < D_{i_s,j_s}/2$, then our algorithm considers two possible 3-divisions. Let vertex index $l_s$ be defined by $D_{i_s,l_s}+D_{l_s,j_s}=\max_{l, l\neq i_s,l\neq j_s}(D_{i_s,l}+D_{l,j_s})$. The vertices indexed by $i_s$, $j_s$ and $l_s$ define a triangle. Now order these vertices, by exchanging $i_s$ and $j_s$ if necessary, so that $D_{i_s,j_s}\geq D_{j_s,l_s}\geq D_{i_s,l_s}$. If $D_{i_s,j_s}/8 \geq D_{j_s,l_s}$ (the largest $D$ value edge has $D$ significantly larger than the other $D$'s), then the edge $(i_s,j_s)$ is trisected, producing three new equal volume subregions of $T_k$. Otherwise, the edge $(i_s,j_s)$ is cut at the point ${\bf v}_c \ =\ (2{\bf v}_{k,i_s}+{\bf v}_{k,j_s})/3$, and two subregions of $T_k$ are produced. The subregion that has the original edge $(j_s,l_s)$ is then divided into two pieces by cutting that edge at the midpoint. The final result is three new equal volume subregions of $T_k$.

Procedure for Choice of Subdivision.

• Determine $D_{i_s,j_s} = \max_{i<j}(D_{i,j})$.
• If 2-division is recommended then
  bisect edge $(i_s,j_s)$ of $T_k$ to produce $T_k^1$ and $T_k^2$.
  else
  1. Determine $D_{i_t,j_t}=\max_{i<j, (i,j)\neq(i_s,j_s)}(D_{i,j})$.
  2. If 4-division is recommended and $D_{i_t,j_t} > D_{i_s,j_s}/2$ then
     i) bisect edge $(i_s,j_s)$ of $T_k$ to produce $T_k^1$ and $T_k^2$;
     ii) bisect edge $(i_t,j_t)$ of both $T_k^1$ and $T_k^2$ to produce four new subregions $T_k^{11}$, $T_k^{12}$, $T_k^{21}$ and $T_k^{22}$.
     else
     1. Determine $D_{i_s,l_s}+D_{l_s,j_s}=\max_{l, l\neq i_s,l\neq j_s}(D_{i_s,l}+D_{l,j_s})$.
     2. Reorder $i_s$, $j_s$, if necessary, so that $D_{i_s,j_s}\geq D_{j_s,l_s}\geq D_{i_s,l_s}$.
     3. If $D_{i_s,j_s}/8 \geq D_{j_s,l_s}$ then
        trisect edge $(i_s,j_s)$ of $T_k$ to produce new subregions $T_k^1$, $T_k^2$ and $T_k^3$.
        else
        i) cut edge $(i_s,j_s)$ of $T_k$ at $(2{\bf v}_{k,i_s}+{\bf v}_{k,j_s})/3$ to produce new subregions $T_k^1$ and $T_k^2$;
        ii) bisect edge $(j_s,l_s)$ of $T_k^2$ replacing $T_k^2$ with new subregions $T_k^{21}$ and $T_k^{22}$.
        end if.
     end if.
  end if.

The heuristic factors $1/2$ and $1/8$ used at steps 2 and 2c in our procedure were selected after extensive testing of the algorithm (see the next section). The cost of determining the subdivision of $T_k$ is $2n(n+1) + 1$ ${\bf f}$ values. The actual cost per new subregion is at most approximately $n(n+1)$ ${\bf f}$ values, and this is small compared to the cost of computing the local rules when the local rule has degree at least seven.