Extrapolation-based Boundary Element Quadrature^{∗}

J. N. Lyness Mathematics and Computer Science Division Argonne National Laboratory 9700 South Cass Avenue Argonne, IL 60439

(In Rend. Semin. Mat. Univ. Pol, Torino: fasc. spec. 1991, pp. 189–203)

Preprint MCS-P208-0191

1 Introduction

The classical Euler Maclaurin summation formula expresses the discretization error made by the trapezoidal rule approximation to a finite integral as an asymptotic expansion in the mesh ratio. As such, it plays a significant role in the theory of numerical quadrature. In particular, it forms the basis for Richardson’s deferred approach to the limit, also known as Romberg integration and as quadrature by extrapolation. In this review, we draw attention to some more recent theory, based on the Euler Maclaurin summation formula (or expansion).

During the past 25 years, this theory has undergone significant development. Briefly, the expansion has been systematized to embrace all quadrature rules. It has been generalized to the N-dimensional cube and to the N-dimensional simplex. Moreover, a significant generalization has been made to cover functions having algebraic and logarithmic singularities at vertices. In this article we discuss some of these developments in the context of their possible use to analyze finite elements. In general, no proofs are provided, but references are given.

This article, and most of the available theory, is limited to the following class of problem (together with its higher-dimensional analogues). The integration regions are the square and the triangle. The integrand function is singular \underline{only at a vertex} of the integration region. Taking this as the origin, we treat only integrand functions of the type

f(x⃗) = ρ^{α} ln ρ Θ(θ)h(r)g(x.y).

Here x⃗ = (x, y) are Cartesian coordinates with the corresponding polar coordinates r, θ. The functions Θ(θ), h(r), and g(x, y) are regular, and rho is homogeneous about the origin. Examples of ρ(x⃗) include the following:

ρ = r
ρ = Ax + By,  A, B > 0
ρ = Ax^{3/2} + By^{3/2}.

The quadrature rules may be linear combinations of m_{i}-copies Q^{(m_{i})}f of any quadrature rule Qf = Q^{(1)}f = ∑ w_{i}f(x⃗_{i}).

∗This work was supported by the Applied Mathematical Sciences subprogram of the Office of Energy Research, U.S. Department of Energy, under Contract W-31-109-Eng-38.

2 The Classical One-Dimensional Background

In one dimension, we treat the interval [0, 1] and an integrand f(x) ∈ C^{p−1}[0, 1]. We denote the exact integral by

If = ∫_{0}^{1} f(x)dx,  (2.1)

and we consider any quadrature rule of the form

Qf = ∑_{j=1}^{ν} w_{j}f(x_{j}),  ∑_{j=1}^{ν} w_{j} = 1.  (2.2)

The m-copy version Q^{(m)} of this quadrature rule Q is obtained by subdividing the integration interval into m equal parts and applying the properly scaled version of Q to each. Specifically,

Q^{(m)}f = ∑_{k=0}^{m−1} ∑_{j=1}^{ν} \frac{w_{j}}{m} f (\frac{x_{j} + k}{m} ) .  (2.3)

Theorem 2.4. (The Euler Maclaurin summation formula). When f(x) ∈ C^{p}[0, 1],

Q^{(m)}f − If = ∑_{s=1}^{p−1} \frac{B_{s}}{m^{s}} + R_{p}^{(m)}(Q; f),  (2.4)

where B_{s} is independent of m and R_{p}^{(m)} = O(m^{−p}).

The important feature of this expansion is that the coefficients B_{s} = B_{s}(Q; f) are independent of m. In fact, they take the form

B_{s}(Q; f) = c_{s}(Q)If^{(s)},  (2.5)

where c_{s}(Q) is the result obtained by applying Q to the integrand function B_{s}(x)/s! and B_{s}(x) is the Bernoulli function. A simple integral representation for the remainder term is available.

We shall now demonstrate how the information that expansion (2.4) exists may be exploited to construct a quadrature technique or a new quadrature rule. In this demonstration, we choose Q to be the trapezoidal rule Qf = \frac{1}{2}(f(0) + f(1)) so that

Q^{(m)}f = \frac{1}{m} (\frac{1}{2}f(0) + ∑_{j=1}^{m−1} f(j/m) + \frac{1}{2}f(1)) .  (2.6)

Let us write down the expansion (2.4) with p = 6 for m = 1, 2 and 4. (We may omit B_{s} for s odd since Q is symmetric and so c_{s}(Q) = 0.) We find

Qf = If + B_{2} + B_{4} + R_{6}^{1}
Q^{(2)}f = If + B_{2}/2^{2} + B_{4}/2^{4} + B_{6}^{2}  (2.7)
Q^{(4)}f = If + B_{2}/4^{2} + B_{4}/4^{4} + R_{6}^{4}.

Extrapolation comprises carrying out the following ad hoc procedure. We remove the remainder terms and solve the resulting set of three equations in three unknowns. Note that, after we have removed the remainder terms, the unknowns are only approximations to If, B_{2}, and B_{4}. For classical reasons we denote by T_{0,2} this particular approximation to If. We find

T_{0,2} = \frac{1}{45}Q^{(1)}f − \frac{20}{45}Q^{(2)}f + \frac{64}{45}Q^{(4)}f,  (2.8)

which, using (2.6), we may reexpress in terms of function values

T_{0,2} = \frac{1}{90} \{7f(0) + 32f(1/4) + 12f(1/2) + 32f(3/4) + 7f(1)\} .  (2.9)

This is, in fact, the Newton Cotes equispaced rule of degree 5.

This straightforward example is given here to illustrate an important general point. Once the form (2.4) of the asymptotic expansion for Q^{(m)}f − If is known, one may construct an integration rule using only linear algebra. No further explicit information about moments or about interpolation polynomials is necessary. Any required information of that type must have already been contained implicitly in the expansion.

We note that (other than B_{s} = 0 for s odd) we do not need to know the form of B_{s} at all, merely that it is independent of m. Thus, starting with a rule Q, the only restriction being that Σw_{j} = 1, and using only Q^{(m)}f = If + ΣB_{s}/m^{s} + R_{p}(m), we can find a respectable quadrature rule of moderate polynomial degree.

We make this point strongly here because, in many elementary descriptions, the theory becomes enmeshed in algebraic detail and the reader can be left with the impression that serious restrictions are in place. The only serious restriction we know is that the integrand must be well behaved.

Note that this technique cannot produce all quadrature rules. One could, however, develop the theory of (constant weight finite interval) quadrature if one took advantage of detailed expressions for the term c_{s}(Q) which appears in (2.5) above. Such a development would be roughly equivalent to the standard development of numerical quadrature theory.

However, given a quadrature rule, one can sometimes find interesting properties or useful extensions simply using the Euler Maclaurin expansion. In one dimension, one discovers only already-known properties. The extension of these ideas to simplices and cubes with simple singularities at vertices can be helpful by providing new insight into standard intractable problems. This is the thrust of the rest of this paper.

3 Multidimensional Generalizations of Classical Results

The generalizations of the asymptotic expansion to the N-dimensional hypercube [0, 1)^{s} is straightforward. In this new context, Qf stands for any rule that integrates the constant function correctly, while Q^{(m)}f stands for the result of subdividing the hypercube into m^{N} subcubes each of side m^{−1} and applying the appropriately scaled version of Q to each. We find immediately that

Q^{(m)}f = If + ∑_{s=1}^{p−1} B_{s}/m^{s} + R_{p}^{(m)}(Q, f).  (3.1)

Here, B_{s} is an obvious generalization of (2.5), namely,

B_{s} = ∑_{s_{i}≥0^{Σs_{i}=s}} c_{s_{1},s_{2},...,s_{N}} (Q)If^{(s_{1},s_{2},...,s_{N})}.  (3.2)

However, the corresponding integral representation for R_{p} (which is O(m^{−p})) is a somewhat lengthy uncompromising expression.

Turning to the simplex, we encounter further complications. There are variant definitions for Q^{(m)}f. We shall proceed as follows. We treat the isosceles simplex

S : x_{i} ≥ 0; Σx_{i} < 1.  (3.3)

This forms a part of the hypercube [0, 1)^{N}. We define our rule Qf as if for this hypercube. Then Q^{(m)}f for the simplex is the result of applying the rule Q to a function θ that coincides with f within S and is zero outside S. (Certain conventions, not given here, apply on the boundary of S within H.) The reason for this somewhat complicated definition can be seen even in two dimensions. The subdivision of the triangle into four equal subtriangles by lines parallel to the edges of the triangle leaves the center triangle with a different orientation to the other three. If we simply defined a rule Q for the triangle S, we would leave ambiguity as to how to apply it in this center triangle. The approach given above provides a definite specification for this situation. The “rule” for the “other” triangle is the part of the rule for the square that lies in this “other” triangle.

An asymptotic expansion of precisely the same form as (3.1) appears to exist for the simplex. The coefficients B_{s} are, as before, independent of m but are given by an expression significantly more complicated than (3.2) above.

In the practice of higher-dimensional extrapolation, it is important to know whether any of these coefficients B_{s} vanish. In fact, certain cofactors c_{s_{1},s_{2},...,s_{n}}(Q) vanish when the rule Q enjoys certain symmetry and polynomial properties (SPP). This, in turn, causes terms B_{s} to vanish as specified in the following theorem.

Theorem 3.4 When Q is symmetric (about (\frac{1}{2}, \frac{1}{2}, …, \frac{1}{2})), then

B_{s} = 0,  ∀ s odd.  (3.4)

When Q is of polynomial degree d (for the hypercube), then

B_{s} = 0,  ∀ s ∈ [1, δ],  (3.5)

where δ = d for the hypercube and δ = d + 1 − N for the N-dimensional simplex.

Conditions (3.4) and (3.5) are termed SPP conditions. It is important to note that in this theorem, one inspects the behavior of Q over the \underline{hypercube} to determine its symmetry and degree. The result of this inspection determines which B_{s} may vanish when it is used over the \underline{simplex} and which B_{s} may vanish when it is used over the \underline{hypercube}. The possible properties of the truncated part of Q that may be applied to S may be of interest independently, but play no direct role in this theorem.

4 Singularities in One-Dimensional Quadrature

The error functional expansions given above are valid only when the integrand functions are sufficiently well behaved. In one dimension, for example, (2.4) is valid only when f(x) ∈ C^{p−1}[0, 1]. This is not simply a technical restriction. The expression for B_{s} contains a factor If^{(s)} that exists only when f^{(s)}(x) is integrable.

In 1961, Navot published a generalization of the Euler Maclaurin expansion which applies to the function f(x) = x^{α} with noninteger α (as well as with integer α). His original proof involved somewhat heavy algebraic manipulation. Successively simpler proofs appeared later. A slightly more general version of his result is the following theorem.

Theorem 4.1 Let f_{α}(x) = κx^{α}, where α ≠ −1; then

Q^{(m)}f_{α} − If_{α} = \frac{A_{α+1}}{m^{α+1}} + ∑_{s=1}^{p−1} \frac{B_{s}}{m^{s}} + O(m^{−p}),  (4.1)

where A_{α+1} and B_{s} are independent of m.

Simple analytic forms for these coefficients are known. A_{α+1} involves the generalized zeta function.

To specify these new coefficients B_{s}, we extend the notation If defined in (2.1) (or (4.2) below). When dealing with f_{α}, which is simply κx^{α}, we define If_{α}^{(s)} as a finite part integral, namely,

If_{α}^{(s)} = ∫_{0}^{1} f_{α}^{(s)}(x)dx,  α − s > −1,  (4.2)
If_{α}^{(s)} = − ∫_{1}^{∞} f_{α}^{(s)}(x)dx,  α − s < −1.  (4.3)

This is not defined for α − s = −1. For other values of α − s, one and only one of the integrals (4.2) or (4.3) converges and this one is used for the definition. In the next section this definition is extended to N dimensions with f_{α} any homogeneous function of order α.

The coefficient in (4.1) is

B_{s} = c_{s}(Q)If_{α}^{(s)}.  (4.4)

Thus, B_{s} is given by its standard formula (2.5) when that is meaningful but is modified when that is not meaningful. We note in passing that the theorem applies, as written, when α < −1, and so the conventional integral does not exist. The term If_{α} on the right of (4.1) is then simply a finite part integral as in (4.3) above. In fact, it is notationally convenient but initially bewildering to replace If_{α} in (4.1) by B_{0} and to include this in the sum on the right. Finally, we note that the only value of α for which Theorem 4.1 is invalid is α = −1. In that case the proper expansion includes other terms not given in (4.1).

Naturally, this expansion is not needed to integrate x^{α} whose integral is known to be 1/(α + 1). It is used to construct an expansion for the more sophisticated function x^{α}g(x), where g(x) ∈ C^{p}[0, 1]. Using a Taylor expansion (with remainder term), we set

f(x) = g(0)x^{α} + g′(0)x^{α+1} + \frac{g′′(0)}{2!}x^{α+2}  (4.5)

+ … + \frac{g^{(p−1)}(0)}{(p − 1)!}x^{α+p−1} + G_{p}(x)x^{α+p}.

One may apply Theorem 2.4 to the final term on the right. The other terms on the right are functions to which Theorem 4.1 applies (for the appropriate α).

Applying these theorems leads to the result

Q^{(m)}f − If ≃ ∑_{j=1} \frac{A_{α+1+j}}{m^{α+1+j}} + ∑_{s=1} \frac{B_{s}}{m^{s}} ,  α ≠ negative integer.  (4.6)

The structure of B_{s} here is not trivial. It is composed of a weighted sum of integrals over two possible intervals. However, it is independent of m, and each integral has a cofactor c_{s}(Q). Thus, these coefficients B_{s} also satisfy the SSP conditions (Theorem 3.4).

Note that, at this point, the restriction in Theorem 4.1 to the effect that α ≠ −1 has the effect of invalidating (4.6) for all negative integers α.

5 The Euler Maclaurin Expansion for a Vertex Singularity

The generalization to the N-dimensional hypercube of the result of the preceding section appeared in Lyness (1976a). Before stating this generalization, Theorem 5.10, we clarify some of the concepts on which it is based.

Definition 5.1 The N-dimensional function f(x⃗) is said to be homogeneous of degree α (about the origin) if

f(λx⃗) = λ^{α}f(x⃗),  ∀ λ > 0 and x⃗ ≠ 0.  (5.1)

We denote such a function by f_{α}(x⃗). We have anticipated this notation in Section 4. The one-dimensional function f_{α}(x) = κx^{α} is clearly homogeneous of degree α.

A monomial x^{p}y^{q} is homogeneous of degree p + q, and many properties relating to the polynomial degree of functions of monomials have direct analogues in the context of homogeneous degree. Thus, (f_{α})^{β} and f_{α}f_{β} are of degree αβ and α + β, respectively, and f_{α}(Mx⃗) when |M| ≠ 0 is also of degree α; ∂^{s}f_{α}/∂x^{s} is of degree α − s.

In more than one dimension, many more interesting functions are homogeneous. For example, r = (x^{2} +y^{2})^{1/2} is homogeneous of degree 1, and θ = arctan y/x is homogeneous of degree zero, as is any function ϕ(θ). Our basic result, (5.10) below, applies to a homogeneous function of specified degree.

We also require an N-dimensional generalization of the interval [a, b) when 0 ≤ a ≤ b.

Definition 5.2 L[a, b); {x⃗ | a ≤ max_{i} x_{i} < b; x_{i} ≥ 0}.

This is the region remaining when [0, a) is removed from [0, b). In two dimensions, it has the shape of the capital letter L.

We shall be integrating homogeneous functions over L-shaped regions. To this end we note first that

∫_{L[a,b)} f_{γ}(x⃗)d^{N}x = m^{−γ−N} ∫_{L[ma,mb)} f_{γ}(x⃗)d^{N}x.  (5.3)

This follows by elementary scaling. Next we express

L[1, ∞) = ⋃_{j=0}^{∞} L[m^{j}, m^{j+1}),  (5.4)

obtaining

∫_{L[1,∞)} f_{γ}(x⃗)d^{N}x = (1 + m^{γ+N} + m^{2(γ+N)} + …) ∫_{L[1,m)} f_{γ}(x⃗)d^{N}x,  (5.5)

which, so long as γ + N < 0, can be put in the form

∫_{L[1,m)} f_{γ}(x⃗)d^{N}x = (1 − m^{γ+N}) ∫_{L[1,∞)} f_{γ}(x⃗)d^{N}x.  (5.6)

When γ + N > 0, a similar derivation leads to

∫_{L[1,m)} f_{γ}(x⃗)d^{N}x = −(1 − m^{γ+N}) ∫_{[0,1)} f_{γ}(x⃗)d^{N}x.  (5.7)

These formulas are basic in establishing the main result (5.10) below. We note that the left-hand side of (5.7) appears at first sight to depend on m in an involved way. In fact, (5.7) shows that the dependence on m has a straightforward form; we write both (5.6) and (5.7) in the form

∫_{L[1,m)} f_{γ}(x⃗)d^{N}x = −(1 − m^{γ+N})If_{γ},  γ + N ≠ 0,  (5.8)

where we have used the following definition, which is an immediate generalization of the one given in (4.2) and (4.3) above.

Definition 5.9 When f_{γ} is of homogeneous degree γ ≠ −N,

If_{γ} = ∫_{[0,1)} f_{γ}(x⃗)d^{N}x  N + γ > 0,
= − ∫_{L[1,∞)} f_{γ}(x⃗)d^{N}x  N + γ < 0.  (5.9)

Theorem 5.10 Let f_{α}(x⃗) be an N-dimensional homogeneous function of degree α ≠ −N which has no singularity in L[1, 2). Let Qf be any integration rule for the hypercube [0, 1)^{N}. Then

Q^{(m)}f_{α} − If_{α} = \frac{A_{α+N}}{m^{α+N}} + \frac{C_{α+N} ln m}{m^{α+N}} + ∑_{s=1}^{p−1} \frac{B_{s}}{m^{s}} + O(m^{−p}),  (5.10)

the constants A_{j}, B_{j}, and C_{j} being independent of m.

In addition, the coefficients B_{j} and C_{j} have relatively simple forms:

B_{s} = ∑_{Σs_{i}=s} c_{s_{1},s_{2},...,s_{N}} (Q)If_{α}^{(s_{1},s_{2},...,s_{n})}  s ≠ α + N
B_{s} = 0  s = α + N
C_{α+N} = 0  α + N ≠ integer
C_{α+N} = ∑_{Σs_{i}=s} c_{s_{1},s_{2},...,s_{N}} (Q) ∫_{L[1,2)} f_{α}^{(s_{1},s_{2},...,s_{N})}d^{N}x/ ln 2 s = α + N.

The proof of this theorem is given in Lyness (1976a). An expression for A_{s} and an integral representation of the remainder term is given there. In fact, although (5.10) is an asymptotic expansion, no use of this is made in the derivation. It may be treated as an identity.

We close this section by making some minor and detailed observations about the formula.

1. The coefficients B_{s} and C_{s} involve the factors c_{s_{1},s_{2},...,s_{N}} (Q) in the same way as in the standard case (where f(x⃗) does not have a singularity). Thus, both B_{s} and C_{s} obey the SPP conditions (Theorem 3.4 above); that is, specified coefficients vanish when Q is symmetric and when Q is of specified polynomial degree. (In fact, C_{s} also satisfies Theorem 3.4 with δ replaced by d.)

2. The condition that f_{α}(x⃗) has no singularity in L[1, 2) is simply an alternative way of indicating that f_{α}(x⃗) has no singularity in any finite region R of the first octant 0 ≤ x_{i} which does not contain the origin.

3. For large α, this formula reduces to the standard N-dimensional Euler Maclaurin formula (3.1). While this is, of course, a convenient special case, the reader should bear in mind that (3.1) is used precisely in this context in the standard proof of (5.10).

4. We note that one can obtain the expansion for F (x⃗) = f_{α}(x⃗)g(x⃗), where g(x⃗) has loworder integrable derivatives, in precisely the same manner as in the one-dimensional case (see (4.6) above).

6 Further Generalizations of the Multidimensional Results

As mentioned above, when one sets into expansion (5.10) the forms of the constants A_{s}, B_{s}, and C_{s} and the integral representation for the remainder term, it becomes an identity. Because of this, it is not surprising to learn that the theorem is valid for all values of α and not only those for which If_{α} is a convergent integral. Intrinsically, it can be used to evaluate finite part integrals. For example, with N = 2 and α = −5/2, we would find an expansion with terms in m^{1/2}, m^{0}, m^{−1/2}, m^{−1}, …, and we might extrapolate to find the coefficient of m^{0} in this expansion. This is not the leading coefficient; this approach, akin to numerical differentiation, is not generally recommended by this author, but we mention it to clarify the underlying situation. When α is a negative integer and α ≤ −N, even this approach is extraordinarily difficult because the coefficient denoted by A_{0} is indistinguishable numerically from that denoted by B_{0}. The situation in one dimension has been clarified and is the subject of a forthcoming paper by Lyness and Monegato, but in higher dimensions it remains open.

Another consequence of the circumstance that this is an identity is that we may obtain a new identity by differentiating this one with respect to any incidental parameter. An obvious candidate is α. For example, when F (x, y) = (Ax + By)^{α}g(x, y), we find ∂F/∂α = F (x, y) ln(Ax + By); and it appears that, when α ≠ integer, we may obtain from the expansion for F an expansion for F (x, y) ln(Ax + By). When α ≠ integer, this derived expansion contains terms in m^{−(α+N+j)}, m^{−(α+n+j)} ln m, and m^{−s}.

At first sight this procedure might appear to be limited to noninteger α. This is because, as α passes through an integer value, one of the coefficients B_{s} in (5.10) changes abruptly as If_{α}^{(s)} changes its region of integration, and another coefficient C_{s} appears and disappears. Further investigation leads to the conclusion that Qf_{α} − If_{α} does not change abruptly and that these apparently abrupt changes result simply from the way the expansion is expressed. In Lyness (1976a), it is shown how to reformulate the expansion. It then appears that crude differentiation of form (5.10) gives, after all, a valid result. The reader may refer to that paper to see how Theorem 3.4 (about SSP conditions) can be applied in general in these cases. However, we recommend that potential users rederive any expansion they plan to use from (5.10) and carry out numerical experiments to confirm its validity.

Other generalizations of (5.10) are geometric in character. The formula as written applies equally to the simplex S : 0 ≤ x_{i} < 1; Σx_{i} < 1, though the coefficients B_{s} are different. This fact may be demonstrated geometrically. In two dimensions, for example, the square H may be considered to be the union of two triangles, S and S̄, with S defined as above. The Euler Maclaurin expansion for S is then the difference of the appropriate expansions for H and for S̄. Since f_{α} has no singularity in S̄, the expansion for S̄ contains only terms B_{s}/m^{s}. It follows that the expansion for S contains the same terms as the expansion for H. This brief synopsis of the two-dimensional derivation requires elaboration, and the N-dimensional result is significantly more difficult to derive; see Lyness (1978) and Lyness and Monegato (1980). It is again a result that is easy to state and is readily predictable, but is difficult to prove.

However, the next geometrical result is easy to prove and extraordinarily useful. The reader will be familiar with the circumstance that, by means of an affine transformation, one simplex may be transformed into another and that any quadrature rule of specified polynomial degree for one is transformed into one of the same degree for the other. The reader should also be aware that this remark refers only to situations where no weight function is involved. When a weight function is involved, the new rule applies with respect to a transformed and generally different weight function. One of the inconvenient aspects of the boundary element method seems to be that different Gaussian rules are required for differently shaped regions and only occasionally are such rules related.

Methods based on extrapolation do not suffer from this drawback. The underlying reason is that the affine transformation of f_{α}, a homogeneous function of degree α, is another homogeneous function of degree α. Furthermore, no weight function is involved in the formulation of the theory. To the author’s knowledge, this is the first time this possibility has been mentioned in the literature. Properly exploited, it allows the straightforward construction of powerful quadrature results for these regions and function singularities.

7 Concluding Remarks

This review article covers a small area of numerical quadrature. It is designed to treat only those aspects of quadrature by extrapolation that might be of interest in the theory of boundary integration methods. While one or two remarks are original, with few exceptions the mathematical theory described here has already been published in the open literature. The author hopes that this article will be helpful in informing the user as to what is available and where to look for further details.

Bibliography

C. T. H. Baker and G. S. Hodgson (1971). Asymptotic expansions for integration formulas in one or more dimensions, SIAM J. Numer. Anal., 8, 473–480.

F. L. Bauer, H. Rutishauser, and E. Stiefel (1963). New aspects in numerical quadrature, Proc. Symp. Appl. Math., 15, American Mathematical Society, Providence, R.I., 199–218.

R. Bulirsch (1964). Bemerkungen zur Romberg-Integration, Numer. Math., 6, 6–16.

R. Bulirsch and J. Stoer (1964). Fehlerabschätzungen und Extrapolation mit rationalen Funktionen bei Verfahren vom Richardson-Typus, Numer. Math., 6, 413–427.

E. de Doncker (1979). New Euler Maclaurin expansions and their application to quadrature over the s-dimensional simplex, Math. Comput., 33, 1003–1018.

E. de Doncker-Kapenga (1987). Asymptotic expansions and their applications in numerical integration, Numerical Integration, P. Keast and G. Fairweather (eds.), D. Reidel Publishing Co., 141–151.

L. Fox (1967). Romberg integration for a class of singular integrands, Comp. J., 10, 87–93.

L. Fox and L. Hayes (1970). On the definite integration of singular integrands, SIAM Rev., 12, 449–457.

F. R. de Hoog and R. Weiss (1973). Asymptotic expansions for product integration, Math. Comp., 27, 295–306.

D. C. Joyce (1971). Survey of extrapolation processes in numerical analysis, SIAM Rev., 13, 435–490.

S. Kumar (1983). Extrapolation to the limit, M.Sc. thesis, University of the South Pacific, Suva, Fiji.

D. B. Hunter (1973). The numerical evaluation of Cauchy principal values of integrals by Romberg integration, Num. Math., 21, 185–192.

J. N. Lyness (1971). The calculation of Fourier coefficients by the Mobius inversion of the Poisson summation formula – Part III. Functions having algebraic singularities, Math. Comp., 25, 483–494.

J. N. Lyness (1976a). An error functional expansion for n-dimensional quadrature with an integrand function singular at a point, Math. Comp. 30, 1–23.

J. N. Lyness (1976b). Applications of extrapolation techniques to multidimensional quadrature for some integrand functions with a singularity, J. Comp. Physics, 20, 346–364.

J. N. Lyness (1978). Quadrature over a simplex, SIAM J. Numer. Anal. 15, 122–133 and 870–877.

J. N. Lyness (1983). QUG2 – Integration over a triangle, Argonne National Laboratory technical memorandum MCS-TM-13.

J. N. Lyness (1985a). Extrapolation methods for multi-dimensional quadrature of singular and regular integrals, Argonne National Laboratory technical memorandum MCS-TM-45.

J. N. Lyness (1985b). The Euler Maclaurin expansion for the Cauchy principal value integral, Num. Math., 46, 611–622.

J. N. Lyness and L. M. Delves (1967). On numerical contour integration round a closed contour, Math. Comp., 21, 561–577.

J. N. Lyness and A. Genz (1980). On simplex trapezoidal rule families, SIAM J. Numer. Anal., 17, 126–147.

J. N. Lyness and B. J. J. McHugh (1970). On the remainder term in the N-dimensional Euler Maclaurin expansion, Numer. Math., 15, 333–344.

J. N. Lyness and G. Monegato (1980). Quadrature error functional expansion for the simplex when the integrand function has singularities at vertices, Math. Comp. 34, 213– 225.

J. N. Lyness and B. W. Ninham (1967). Numerical quadrature and asymptotic expansions, Math. Comp., 21, 162–178.

J. N. Lyness and K. K. Puri (1973). The Euler-Maclaurin expansion for the simplex, Math. Comp., 27, 273–293.

G. Monegato and J. N. Lyness (1979). On the numerical evaluation of a particular singular two-dimensional integral, Math. Comp., 33, 993–1002.

I. Navot (1961). An extension of the Euler-Maclaurin summation formula to functions with a branch singularity, J. Math. and Phys. 40, 271–276.

I. Navot (1962). A further extension of the Euler-Maclaurin summation formula, J. Math. and Phys., 41, 155–163.

W. Romberg (1955). Vereinfachte numerische Integration, Kgl. Norske Vid. Selsk. Forh., 38, 30–36.

A. Sidi (1979). Some properties of a generalization of the Richardson extrapolation process, J. Inst. Maths. Applics., 24, 327–346.

A. Sidi (1980). Extrapolation methods for oscillatory infinite integrals, J. Inst. Maths. Applics., 26, 1–20.

A. Sidi (1982). The numerical evaluation of very oscillatory infinite integrals by extrapolation, Math. Comp., 38, 517–529.

A. Sidi (1983). Euler-Maclaurin expansions for integrals over triangles and squares of functions having algebraic/logarithmic singularities along an edge, J. Approx. Th., 39, 39–53.

A. Sidi (1987). Extrapolation methods for divergent oscillatory infinite integrals that are defined in the sense of summability, J. Comp. Appl. Math., 17, 105–114.