Quadrature over Curved Surfaces by Extrapolation^{∗}^{†}

J. N. Lyness

Mathematics and Computer Science Division Argonne National Laboratory 9700 South Cass Avenue Argonne, IL 60439

Abstract

In this paper we describe and justify a method for integrating over curved surfaces. This method does not require that the Jacobian be known explicitly. This is a natural extension of extrapolation (or Romberg integration) for planar squares or triangles.

1 Introduction

We treat the numerical integration over a curvilinear quadrilateral, σ, and over a curvilinear triangle, τ. These are embedded in a curved surface, which may be parameterized in the form

x = (x, y, z) = Φ(u) = Φ(u, v)  (1.1)
x = ϕ(u, v); y = ψ(u, v); z = χ(u, v).

The region σ (or τ) is that part of the curved surface that is mapped from the unit square S (or from the unit triangle T). That is, σ : Φ(u ∈ S) and τ : Φ(u ∈ T), where S and T are, respectively,

S : u ∈ [0, 1)^{2}  (1.2)
T : u + v < 1, u ≥ 0, v ≥ 0.

The method is designed for a situation in which a triangularization (see below) based on the mapping Φ is readily available; the functions ϕ, ψ, and χ and the Jacobian J (see below) need not be known explicitly. However, all these functions, known or unknown, are well behaved, and the surface in which σ and τ are embedded is smooth.

A situation in which this may occur is one in which the surface is defined implicitly:

H(x, y, z) = 0,

∗This work was supported by the Office of Scientific Computing, U.S. Department of Energy, under Contract W-31-109-Eng-38.

†This version (November 1993) has an appendix which replaces Section 6 of the previous version. (See the first paragraph of the Appendix.)

and ∇H is available. In some such cases, one may obtain accurate values of x iteratively, without having to know Φ explicitly. Georg (1991) discusses problems of this sort in the context of the boundary element method. He conjectures the existence of an asymptotic expansion, on which a convenient extrapolation scheme may be based. In this paper, we have been able to establish this expansion. This puts on a sound theoretical footing some of his innovative but partly heuristic work in this area.

We now describe a simple conventional method based on triangularization. Throughout this paper, m is a positive integer. The m-triangularization of the planar square [0, 1]^{2} comprises a partition into 2m^{2} distinct triangles using the lines

u = j/m;  v = j/m;  u + v = j/m  for all j.  (1.3)

The vertices of these triangles comprise a set of (m + 1)^{2} points (t_{j}, t_{k}), where j, k ∈ [0, m] and

t_{j} = t_{j}(m) = j/m.  (1.4)

We shall suppress the dependence on m in cases where no confusion is likely to arise.

The mapping function x = Φ(u) induces an m-triangularization of the curvilinear quadrilateral σ. The lines (1.3) are mapped into curves embedded in the curved surface, forming a set of 2m^{2} curvilinear triangles. The integral I_{σ}f may be approximated using as abscissas the (m + 1)^{2} distinct points that are vertices of these triangles. Let τ_{μ}, μ = 1, 2, ..., 2m^{2}, represent these curvilinear triangles; let f_{μ,1}, f_{μ,2}, f_{μ,3} represent the function values of f at the three vertices of τ_{μ}; and let \overline{Δ}_{μ} be the area of the planar triangle having these vertices. Then it is evident that

Q^{(m)}f = (1/3) ∑_{μ=1}^{2m^{2}} \overline{Δ}_{μ}(f_{μ,1} + f_{μ,2} + f_{μ,3})  (1.5)

is an approximation to I_{σ}f.

The principal result in this paper is that Q^{(m)}f has an h^{2}- error expansion. That is, setting h = 1/m, when f ∈ C^{(2p)}(σ) and Φ(u) is sufficiently well behaved, we have

Q^{(m)}f = I_{σ}f + B_{2}h^{2} + B_{4}h^{4} + .... + B_{2p}h^{2p} + O(h^{2p+1}),  (1.6)

where the coefficients B_{2q} are independent of h. In the course of obtaining this result, we uncover several other results of the same nature. Finally we obtain the result, given in Theorem 5.6 below, which is analogous to (1.6) when the region σ is replaced by the curvilinear triangle τ.

This theorem is in fact Conjecture 5.1 of Georg (1991), which is itself part of Conjecture 1 of Georg and Tausch (1993). This expansion may be used as a basis for h^{2}-extrapolation; see, for example, Bauer, Rutishauser and Stiefel (1963). One evaluates successively Q^{(m)}f for a sequence of values of m and, either by means of a Neville-Romberg T-table or by some other simple technique, one obtains a sophisticated result based on this sequence of possibly individually inadequate results. Some advantages and other features of this approach form part of the paper by Georg. The present paper is confined to one significant task: to establish the asymptotic expansion on which the underlying theory of his work may be based.

The author recently became aware that this same conjecture was established by Verlinden and Cools (1993) simultaneously and independently. Their paper appears in this issue of this journal.

As a preliminary, in Section 2, we describe a method based on extrapolation that one might use if the Jacobian function J in (2.2) below were readily available. The subsequent three sections are devoted to modifying this method so as to avoid the explicit evaluation of J. Section 3 is devoted to the surface, delineating the relation between \overline{Δ}_{μ} and J in some detail. Section 4 provides results for the curvilinear quadrilateral, and Section 5 widens these results to the curvilinear triangle. Section 6 provides details of some Euler-Maclaurin type expansions for the triangle.

2 Background Theory

In a recent article, Schwab and Wendland (1992) provide a broad survey and analysis of integration over curved surfaces required in the boundary element method. For the benefit of a wider audience, in this paper we start directly from the classical theory. The general theory of analytic integration over surfaces is treated in Courant (1956; see pp. 273 et seq.). We follow, as far as possible, the framework introduced there. We set

∫∫_{σ}f(x, y, z)dS = ∫∫_{S}f(ϕ(u, v), ψ(u, v), χ(u, v))J(u, v) dudv,  (2.1)

where

J^{2} = (ϕ_{u}ψ_{v} − ψ_{u}ϕ_{v})^{2} + (ψ_{u}χ_{v} − χ_{u}ψ_{v})^{2} + (χ_{u}ϕ_{v} − ϕ_{u}χ_{v})^{2}.  (2.2)

In this case the problem reduces to that of evaluating

Ig = ∫_{0}^{1}∫_{0}^{1}g(u, v)dudv,

where

g(u, v) = f(ϕ(u, v), ψ(u, v), χ(u, v))J(u, v).  (2.3)

There are many ways of integrating over a square. The classical approach (see, for example, Stroud (1971)) is by means of Gaussian formulas. A marginally less efficient but occasionally more convenient approach is by extrapolation (also known as Romberg integration), one version of which we now describe. Let

Rg = ∑_{i=1}^{ν} w_{i}g(u_{i}, v_{i})  (2.4)

be any quadrature rule for S that integrates a constant function correctly. We define the m^{2}-copy of this rule R as

R^{(m)}g = ∑_{i=1}^{ν} ∑_{j=0}^{m−1} ∑_{k=0}^{m−1} \frac{w_{i}}{m^{2}} g(\frac{(u_{i} + j)}{m}, \frac{(v_{i} + k)}{m}).  (2.5)

When g ∈ C^{p}(S) , the following minor generalization of the Euler-Maclaurin summation formula is valid:

R^{(m)}g − Ig = ∑_{s=1}^{p−1} \frac{B_{s}(R; g)}{m^{s}} + O(m^{−p}).  (2.6)

Note that the coefficients B_{s} are independent of m, the mesh ratio. Note also that, in general, this is a full expansion containing all inverse integer powers of m. However, when the rule R is symmetric about the point (1/2, 1/2), all terms B_{s} with odd s vanish. Further, when R is of polynomial degree d, we have also that B_{s} = 0 for all s ∈ [1, d]. Thus, in a context in which J(u, v) is readily available, an attractive way to carry out integration over a curvilinear quadrilateral is by transformation to a plane square and the use of extrapolation. A conventional rule to use in this context might be the product m-panel trapezoidal rule

R^{(m)}g = \frac{1}{m^{2}} ∑_{j=0}^{m}′′∑_{k=0}^{m}′′g(t_{j}, t_{k}).  (2.7)

Here, we can identify the positive integer m as the reciprocal of a step length h and define

t_{j} = t_{j}(m) = jh = j/m.  (2.8)

This is, of course, the m^{2}-copy of a four-point rule, which applies equal weight to each of the four vertices. This rule has an even error expansion and is used on occasion in the context of two-dimensional Romberg integration.

The rules with which we shall be concerned are more primitive than this. One of these is a rectangle rule defined by

R_{λ}g = g(0, 0).  (2.9)

Neither this nor its m-copy (given in (3.11) below), is symmetric, and the expansion (2.6) above is full. Such a rule would be rarely used in practice. However, it turns out to be convenient to develop the theory in terms of such primitive rules. We shall return to this rule and several variants in Section 3.

In the rest of this paper we deal with a context in which J(u, v) is not readily available. We describe how the simple approach given above can be modified to deal with this less transparent situation.

3 The m-Triangularization

In this section we are concerned with the relation between individual points of the mtriangularization of σ and the points that are the vertices of these triangles.

Definition 3.1 Let λ be a four-integer index (α, β, γ, δ). Let t_{i} = i/m and x_{j,k} = Φ(t_{j}, t_{k}). Then, depending on context, Δ_{j,k,λ}^{(m)} either denotes the triangle having vertices (t_{j}, t_{k}), (t_{j+α}, t_{k+β}),

(t_{j+γ}, t_{k+δ}), or denotes the area of this triangle; and \overline{Δ}_{j,k,λ}^{(m)} either denotes the planar triangle having vertices x_{j,k}, x_{j+α,k+β}, x_{j+γ,k+δ}, or denotes the area of that triangle.

Definition 3.2 In the above context,

det λ = αδ − βγ
−λ = (−α, −β, −γ, −δ).

When detλ = 0, the triangle Δ_{j,k,λ}^{(m)} has zero area. Moreover (α, β, γ, δ) and (γ, δ, α, β) refer to identical triangles. So, without any compromise, we may consider only indices λ for which detλ is a positive integer.

In this paper we shall use only six distinct indices λ, though much of the theory would allow a general assignment. These six refer to the distinct triangles of the triangularization that have x_{j,k} as a vertex. These indices are those that appear in Definition 3.6 below.

In this section we shall prove various results that are valid for all λ. We shall confine the details of the proofs to the case λ = (1, 0, 0, 1) so that \overline{Δ}_{j,k,λ}^{(m)} is the area of a planar triangle having vertices at x_{j,k}, x_{j+1,k}, x_{j,k+1} (the northeast elementary triangle).

Theorem 3.3 Let Φ(u) (defined in (1.1)) be C^{(p)}(S), and det λ = 1. Then

2\overline{Δ}_{j,k,λ}^{(m)} = h^{2}(δ_{0} + hδ_{1} + … + h^{p−1}δ_{p−1}) + O(h^{p+2}),  (3.1)

where

δ_{0} = J(t_{j}, t_{k})

and

δ_{s} = δ_{s,λ}(t_{j}, t_{k}),

where δ_{s,λ}(u, v) is a function of u and v having continuous partial derivatives of order p.

Proof. We treat only the case λ = (1, 0, 0, 1), and within this proof we abbreviate \overline{Δ}_{j,k,λ}^{(m)} to \overline{Δ}. The area of this triangle is one half of the absolute value of the triple vector product

| \begin{matrix}e_{1} & e_{2} & e_{3} \\ ϕ(t_{j+1}, t_{k}) − ϕ(t_{j}, t_{k}) & ψ(t_{j+1}, t_{k}) − ψ(t_{j}, t_{k}) & χ(t_{j+1}, t_{k}) − χ(t_{j}, t_{k}) \\ ϕ(t_{j}, t_{k+1}) − ϕ(t_{j}, t_{k}) & ψ(t_{j}, t_{k+1}) − ψ(t_{j}, t_{k}) & χ(t_{j}, t_{k+1}) − χ(t_{j}, t_{k}) \end{matrix} | .  (3.2)

We may expand this in the form

e_{1}D^{(1)} + e_{2}D^{(2)} + e_{3}D^{(3)},

giving

4\overline{Δ}^{2} = D^{(1)2} + D^{(2)2} + D^{(3)2}.

Since ϕ(u, v) and the other components of P are well behaved, we may use the expansion

ϕ(t_{j+1}, t_{k}) − ϕ(t_{j}, t_{k}) = hϕ_{u} + \frac{h^{2}}{2} ϕ_{uu} + …  (3.3)

and similar expansions for the other elements of D^{(3)}. This gives

D^{(3)} = h^{2} | \begin{matrix}ϕ_{u} + \frac{h}{2}ϕ_{uu} + … & ψ_{u} + \frac{h}{2}ψ_{uu} + … \\ ϕ_{v} + \frac{h}{2}ϕ_{vv} + … & ψ_{v} + \frac{h}{2}ϕ_{vv} + … \end{matrix} |

= h^{2}(D_{0}^{(3)} + hD_{1}^{(3)} + … + h^{p−1}D_{p−1}^{(3)}) + O(h^{p+2}),

where

D_{0}^{(3)} = | \begin{matrix}ϕ_{u}ψ_{u} \\ ϕ_{v}ψ_{v} \end{matrix} | = \frac{∂(ϕ, ψ)}{∂(u, v)}

and D_{s}^{(3)}, the coefficient of h^{s+2} in this expansion, is a sum of analogously constructed products of partial derivatives. It follows that 4Δ^{2} has an expansion in powers of h of the form

4\overline{Δ}^{2} = h^{4} (J^{2} + ∑_{s=1}^{p−1} h^{s}C_{s}) + O(h^{p+4}).

Finally, for sufficiently small h, we may take the square root of this to obtain

2\overline{Δ} = h^{2}J + ∑_{s=1}^{p−1} δ_{s}h^{s} + O(h^{p+2}).

The series manipulation is elementary; to establish the theorem, we need information about the continuity of the functions δ_{s,λ}(u, v).

Since the derivatives of ϕ, ψ, and χ of order p are continuous functions of u and v, it follows that so are the derivatives of D_{s}^{(3)} of order p − 1 − s; similarly, those of D_{s}^{(1)} and D_{s}^{(2)} have the same order of continuity. Then C_{s} has this order. Finally, so long as J(u, v) is bounded away from zero in the region S, we find that δ_{s} also has this order of continuity.

It is relatively straightforward to show that the effect of replacing λ by −λ in the above proof is the same as that of reversing the sign of h ; this gives

δ_{s,λ}(u, v) = (−1)^{s}δ_{s,−λ}(u, v).  (3.4)

It is no surprise that the first term in expansion (3.1) is h^{2}J. This is needed for the classical theory to survive. An expansion in powers of h is not unexpected. What is critical is that the coefficients δ_{s,λ}(u, v) are smooth functions of u and v. Although they occur in the form δ_{s,λ}(t_{j}, t_{k}) and t_{j}, t_{k} depend on m, the functions δ_{s,λ}(u, v) do not depend on m. Like J = δ_{0}(u, v), they depend on ϕ, ψ, and χ only. This circumstance is exploited later.

The condition det λ = 1 results in δ_{0,λ}(u, v) = J(u, v). If this were relaxed, the same theorem would hold with δ_{0,λ}(u, v) = J(u, v) det λ.

We now discuss the connection between the m-triangularization and the indices λ. The following definition is a first step in clarifying a somewhat involved situation.

Definition 3.4 Let B denote the square S or the triangle T defined in (1.2). Then δ_{j,k,λ}^{(m)}(B) =

1 or 0 depending on whether or not the triangle Δ_{j,k,λ}^{(m)} forms part of the m-triangularization (1.3) of B.

By inspection, one can verify that there are only six distinct values of index λ for which it is possible for any such triangle to be part of the triangularization. We refer to these six (which are precisely those that appear in Definition 3.6 below) as the triangularizationcompatible (T.C.) indices.

We now return to the construction of quadrature rules.

Definition 3.5 Let B denote S or T. We denote by R_{λ}^{(m)}(B) the rule that assigns a weight m^{−2} to all points (t_{j}, t_{k}) of the m-triangularization of B for which the triangle Δ_{j,k,λ}^{(m)} lies in B and is a member of the triangularization.

Definitions 3.4 and 3.5 will be applied in Section 5 in a context where B is the triangle T. In this section and in Section 4 we shall suppress the argument (S). We find that, for the six T.C. indices, the definition is satisfied as follows.

Definition 3.6 R_{λ}^{(m)} is the m^{2}-copy of R_{λ} given by

R_{(1,0,0,1)}g = g(0, 0)  (3.5)
R_{(0,1,−1,1)}g = g(1, 0)  (3.6)
R_{(−1,1,−1,0)}g = g(1, 0)  (3.7)
R_{(−1,0,0,−1)}g = g(1, 1)  (3.8)
R_{(0,−1,1,−1)}g = g(0, 1)  (3.9)
R_{(1,−1,1,0)}g = g(0, 1).  (3.10)

Clearly, for no index λ is R_{λ}^{(m)} the product trapezoidal rule (which is symmetric). It is the m^{2}-copy of a one-point rule that assigns full weight to one corner of the unit square. Specifically, in accordance with Definition 3.4,

R_{λ}^{(m)}g = (1/m^{2}) ∑_{j=0}^{m} ∑_{k=0}^{m} δ_{j,k,λ}^{(m)}g(t_{j}, t_{k}),  (3.11)

where δ_{j,k,λ}^{(m)} = 1 for only m^{2} of the (m + 1)^{2} specified abscissas in (3.11), and is zero otherwise. For example,

R_{(1,0,0,1)}^{(m)}g = \frac{1}{m^{2}} ∑_{j=0}^{m−1} ∑_{k=0}^{m−1} g(t_{j}, t_{k}).  (3.12)

However, since it is a copy rule, the full Euler-Maclaurin expansion (2.6) applies to this rule.

Corresponding to each R_{λ}^{(m)} we now define a quadrature rule for the curvilinear quadrilateral σ. This is obtained by first transforming the rule onto the curved surface and then replacing the Jacobian-dependent term by the area of an elementary triangle. The transformation gives

R_{λ}^{(m)}g = ∑_{j=0}^{m} ∑_{k=0}^{m} δ_{j,k,λ}^{(m)}h^{2}J(t_{j}, t_{k})f(x_{j,k}),  (3.13)

and the arbitrary replacement leads to the following definition.

Definition 3.7

Q_{λ}^{(m)}f = ∑_{j=0}^{m} ∑_{k=0}^{m} δ_{j,k,λ}^{(m)}2\overline{Δ}_{j,k,λ}^{(m)}f(x_{j,k}).  (3.14)

Although this appears to contain (m + 1)^{2} abscissas, in fact only m^{2} have nonzero weights. For example, in the case λ = (1, 0, 0, 1), abscissas with j = m and with k = m have zero weights.

Note carefully that, in spite of the notation, the rule in (3.14) is not an m^{2}-copy rule. It is obtained by modifying a term in a rule that is an m^{2}-copy rule in a different space. The expansion of Theorem 3.3 allows us to establish the following.

Theorem 3.8 Let Φ(u) (defined in (1.1)) be C^{(p)}(S), and let Q_{λ}^{(m)}f and R_{λ}^{(m)}g be given by (3.14) and (3.11), respectively. Then

Q_{λ}^{(m)}f = ∑_{s=0}^{p−1} h^{s}R_{λ}^{(m)}g_{s,λ} + O(h^{p}),  (3.15)

where

g_{s,λ}(u, v) = f(x(u, v))δ_{s,λ}(u, v)  (3.16)

and the functions δ_{s,λ} are defined in Theorem 3.3.

Proof. This is simply a matter of substituting for \overline{Δ}_{j,k,λ}^{(m)} in (3.14) the expression given in (3.1) and simplifying by using (3.11).

Note that g_{0,λ}(u, v) coincides with g(u, v). Note also that (3.15) is not an h-expansion, since R_{λ}^{(m)}g_{s} depends on h through m.

We have, as yet, made no commitment about the properties of f or g. To proceed, we need to be specific about the form of g_{s,λ}(u, v). In the rest of this paper we treat the case in which f(x) and consequently g_{s,λ}(u, v) are well behaved.

4 Rules for the Curvilinear Quadrilateral

Theorem 3.8 expresses Q_{λ}^{(m)}f, a rule over the curved surface σ, in an expansion, each term of which involves a product rectangle rule evaluation over the planar square S. As mentioned above, this is not an h-expansion.

We now restrict ourselves to an integrand g ∈ C^{(p)}(S). In this case we may make use of the Euler-Maclaurin expansion (2.6) applied to each rule R_{λ} in this expansion. This expansion is of the form

R_{λ}^{(m)}g = ∑_{j=0}^{p−1} h^{j}B_{j}(R_{λ}; g) + O(m^{−p}).  (4.1)

Here, we have set h = 1/m; bearing in mind that

B_{0}(R; g) = I(g) ≡ ∫_{0}^{1}∫_{0}^{1} g(u, v)dudv  (4.2)

and that B_{j}(R; g) is independent of m, we may substitute (4.1) into (3.15) to establish the following theorem.

Theorem 4.1 When Φ(u) (defined in (1.1)) and f(x(u, v)) are C^{(p)}(S), the quadrature rule Q_{λ}^{(m)}f of Definition 3.7 has an h-error expansion, namely,

Q_{λ}^{(m)}f − I_{σ}f = ∑_{w=1}^{p−1} \frac{B̃_{w}(λ; f)}{m^{w}} + O(h^{p}),  (4.3)

where

B̃_{w}(λ; f) = ∑_{j=0}^{w} B_{j}(R_{λ}; g_{w−j,λ}).  (4.4)

When λ is one of the six T.C. indices, the rule described by this theorem has the property that the points for function evaluation lie in the completion of σ.

Any one of these rules alone seems quite reasonable, but a little “lop-sided”. In the context of extrapolation, a rule having an h^{2}-error expansion would be preferable.

The rest of this section is devoted to deriving such a rule. It turns out (see Theorem 4.4 below) that the average of the six rules mentioned above has an h^{2}-error expansion.

Lemma 4.2 The coefficients in the Euler-Maclaurin expansion (4.1) satisfy

B_{j}(R_{λ}; g) = (−1)^{j}B_{j}(R_{−λ}; g).  (4.5)

Proof. The result may be established by direct evaluation of the coefficients in terms of Bernoulli functions. A more elegant approach is based on the circumstance that a symmetric rule has an m^{2}-error expansion (and an antisymmetric rule operator has an expansion involving only odd powers of m). It is readily verified that

Rg = \frac{1}{2}R_{λ}g + \frac{1}{2}R_{−λ}g  (4.6)

is symmetric. When j is odd, B_{j}(R; g) vanishes, leaving

B_{j}(R_{−λ}, g) = −B_{j}(R_{λ}, g) for all odd j.

A corresponding argument provides a corresponding result for even j, thus establishing (4.5).

Another simple proof using an integral representation of B_{j}(R_{λ}; g) is outlined in the appendix.

A feature of the foregoing theory, which has been specially built in, is that in many hexpansions involving λ, the effect of reversing the sign of h is the same as that of reversing the sign of λ. In particular, in view of (3.4), the functions defined in (3.16) satisfy

g_{s,λ}(u, v) = (−1)^{s}g_{s,−λ}(u, v).  (4.7)

We are now ready to consider the rule defined in (4.8) below. This rule applies to each fully interior point x_{i,j} a weight equal to the sum of the areas of two opposite triangles.

(For example, when λ = (1, 0, 0, 1), these are the one to the northeast and the one to the southwest.) In general, points on edges have weights involving only one of these triangles, the interior one. Two of the vertices have zero weight.

Theorem 4.3 Under the hypotheses of Theorem 4.1 the quadrature rule

Q^{(m)}f = \frac{1}{2}Q_{λ}^{(m)}f + \frac{1}{2}Q_{−λ}^{(m)}f  (4.8)

has an h^{2}-error expansion.

Proof. The proof relies entirely on Equations (4.5) and (4.7). These are used in (4.4) to show

B̃_{w}(λ; f) = ∑_{j=0}^{w} B_{j}(R_{λ}; g_{w−j,λ}) = ∑_{j=0}^{w}(−1)^{j}(−1)^{w−j}B_{j}(R_{−λ}; g_{w−j,−λ}) = (−1)^{w}B̃_{w}(−λ; f).

(4.9)

The error expansion for Q^{(m)}f in (4.8) is the average of two asymptotic series of the form (4.3), the second differing from the first only in that λ is replaced by −λ. Thus the coefficient of m^{−w} in this expansion is (1/2)(B̃_{w}(λ; f) + B̃_{w}(−λ; f)). In view of (4.9), this coefficient vanishes when w is odd. This establishes the theorem.

Our final result for the curvilinear quadrilateral is as follows.

Theorem 4.4 Under the hypotheses of Theorem 4.1 the quadrature rule

Q^{(m)}f = (1/6)(Q_{(1,0,0,1)}^{(m)}f + Q_{(0,1−,1,1)}^{(m)}f + Q_{(−1,1,−1,0)}^{(m)}f
+Q_{(−1,0,0,−1)}^{(m)}f + Q_{(0,−1,1,−1)}^{(m)}f + Q_{(1,−1,1,0)}^{(m)}f),  (4.10)

has an h^{2}-error expansion.

Proof. The proof is a simple corollary of Theorem 4.3. The quadrature rule here is obtained by taking the mean of three rules, each being of form (4.8); specifically, the three values of λ involved are those in the first three assignments in (4.10). Naturally, the mean of three h^{2}-error expansions is itself an h^{2}-error expansion.

This rule applies to each point x_{j,k} a contribution from each of the triangles inside σ of which it is a vertex. This contribution is one-third of the planar area of the triangle in question. Thus it could be reexpressed in terms of contributions from triangles. Doing this reveals the expression in (1.5). Thus this theorem is one of the principal results of this paper, foreshadowed in the introduction.

5 Rules for the Curvilinear Triangle

In this section we treat the numerical integration over the curvilinear triangle τ, which, we recall from Section 1, is embedded in the same curved surface (1.1) as is σ. Thus much of the theory in Sections 1–3 applies with only minor modification. The region τ is that part of the curved surface that is mapped from the unit triangle. That is, τ : Φ(u ∈ T), where T is

T : u + v < 1 u ≥ 0 v ≥ 0.  (5.1)

The same triangularization (1.3) is used, and so is the same notation introduced in Section 3 to describe individual triangles: m^{2} of the triangles in the m-triangularization of S form the m-triangularization of T.

We may construct a rule for T using the technique of Section 3. We use Definition 3.4 for δ_{j,k,λ}^{(m)}(T) and then define R_{λ}^{(m)}(T) in accordance with Definition 3.5. This gives

R_{λ}^{(m)}(T)g = (1/m^{2}) ∑_{j=0}^{m} ∑_{k=0}^{m} δ_{j,k,λ}^{(m)}(T)g(t_{j}, t_{k}).  (5.2)

Following the same steps as in Section 3, we are led to the correspondent of Definition 3.7.

Definition 5.1

Q_{λ}^{(m)}(τ)f = ∑_{j=0}^{m} ∑_{k=0}^{m} δ_{j,k,λ}^{(m)}(T)2\overline{Δ}_{j,k,λ}^{(m)}f(x_{j,k}).  (5.3)

One may verify that, when λ is one of the six T.C. indices, this employs either m(m + 1)/2 or m(m − 1)/2 function values. The result of Theorem 3.8 in this context is simply

Q_{λ}^{(m)}(τ)f = ∑_{s=0}^{p−1} h^{s}R_{λ}^{(m)}(T)g_{s,λ} + O(h^{p}),  (5.4)

where, as before,

g_{s,λ}(u, v) = f(x(u, v))δ_{s,λ}(u, v)  (5.5)

and the functions δ_{s,λ} are defined in Theorem 3.3.

In Section 4, we proceeded by invoking the two-dimensional Euler-Maclaurin expansion for the square. Less well known is the corresponding formula for the triangle. Besides the overall formula, we shall also require some details about the structure of the coefficients. The following theorem is a specialization of one given in Lyness and Puri (1973).

Theorem 5.2 When g ∈ C^{(p)}(T)

R_{λ}^{(m)}(T)g = ∑_{j=0}^{p−1} h^{j}B_{j}(R_{λ}; T; g) + O(m^{−p}),  (5.6)

where

B_{0}(R_{λ}; T; g) = I(T; g) ≡ ∫∫_{T} g(u, v)dudv.  (5.7)

An integral representation for B_{j}(R_{λ}; T; g) is given in the appendix. We now substitute (5.6) into (5.4) to obtain the following theorem.

Theorem 5.3 Under the hypotheses of Theorem 4.1 the quadrature rule Q_{λ}^{(m)}(τ)f of Definition 5.1 has an h-error expansion, namely,

Q_{λ}^{(m)}(τ)f − I_{τ} f = ∑_{w=1}^{p−1} \frac{B̃_{w}(λ; τ; f)}{m^{w}} + O(h^{p}),  (5.8)

where

B̃_{w}(λ; τ; f) = ∑_{j=0}^{w} B_{j}(R_{λ}; T; g_{w−j,λ}).  (5.9)

It remains to establish quadrature rules having h^{2}-error expansions. We follow precisely the same approach as in the preceding section for the square. We need the following result, which corresponds to Lemma 4.2 and relates coefficients in different Euler-Maclaurin expansions.

Lemma 5.4 The coefficients in the Euler-Maclaurin expansion (5.6) satisfy

B_{j}(R_{λ}; T; g) = (−1)^{j}B_{j}(R_{−λ}; T; g).  (5.10)

This is significantly more difficult to prove than the corresponding result for the square, and we defer this to the appendix. However, once this is established, the rest of the theory follows in a relatively straightforward manner. Corresponding to Theorem 4.3, we have the following theorem.

Theorem 5.5 Under the hypotheses of Theorem 4.1 the quadrature rule

Q^{(m)}(τ)f = \frac{1}{2}Q_{λ}^{(m)}(τ)f + \frac{1}{2}Q_{−λ}^{(m)}(τ)f  (5.11)

has an h^{2}-error expansion.

Proof. Exactly as in the proof of Theorem 4.3 we use (5.9), (5.10) and (4.7) to show that (1/2)(B̃_{w}(λ; τ; f) + B̃_{w}(−λ; τ; f)) vanishes when w is odd. Since this is the coefficient of m^{−w} in the error expansion of Q^{(m)}(τ)f, the theorem is established.

Our final result for the curvilinear triangle is as follows.

Theorem 5.6 Under the hypotheses of Theorem 4.1 the quadrature rule

Q^{(m)}(τ)f = (1/6)(Q_{(1,0,0,1)}^{(m)}(τ)f + Q_{(0,1,−1,1)}^{(m)}(τ)f + Q_{(−1,1,−1,0)}^{(m)}(τ)f
+Q_{(−1,0,0,−1)}^{(m)}(τ)f + Q_{(0,−1,1,−1)}^{(m)}(τ)f + Q_{(1,−1,1,0)}^{(m)}(τ)f)  (5.12)

has an h^{2}-error expansion.

This rule is simply the mean of three examples of the rule in Theorem 5.5. These three have the indices λ coinciding with the first three occurring in (5.12). The mean of three rules each of which has an h^{2}-error expansion also has an h^{2}-error expansion. As mentioned in Section 1, this theorem coincides with Conjecture 5.1 of Georg (1991) and with part of Conjecture 1 of Georg and Tausch (1993).

A Appendix: Notes on Euler-Maclaurin Expansions for Planar Regions

This appendix replaces Section 6 of a previous version of this paper. In that version, equations (6.5), (6.10), and (6.11) are incorrect. A correct example of (6.5) is given in (A.4). The other equations require modifications of the same nature.

The coefficients in the Euler-Maclaurin expansion (2.6) for the square, S, for any rule (2.5) (weights w_{i}; abscissas (u_{i}, v_{i})) have the following simple integral representation,

B_{j}(R; g) = ∑_{j_{i}≥0^{j_{1}+j_{2}=j}} c_{j_{1},j_{2}}(R) ∫_{S} g^{(j_{1},j_{2})}(u, v)dudv,  (A.1)

where

c_{j_{1},j_{2}}(R) = ∑_{i=1}^{ν} w_{i} \frac{B_{j_{1}}(u_{i})B_{j_{2}}(v_{i})}{j_{1}! j_{2}!}  (A.2)

is essentially the result of applying the rule R to the product of two Bernoulli polynomials.

For the rules R_{λ}, defined in Definition 3.6, the coefficients B_{j}(R_{λ}; g) in (4.1) are of the same form. Using the properties of the Bernoulli polynomials, we find

c_{j_{1},j_{2}}(R_{λ}) = \frac{B_{j_{1}}(ζ)B_{j_{2}}(η)}{j_{1}! j_{2}!},  (A.3)

where ζ, η are individually 0 or 1 and (ζ, η) occurs as argument of g in Equations (3.5) to (3.10).

Lemma 4.2 follows from (A.1) and (A.3) as a consequence of the symmetry properties of the Bernoulli polynomials, namely, B_{j}(0) = (−1)^{j}B_{j}(1).

The Euler-Maclaurin expansion for the triangle, T, is less known and is dealt with at length in Lyness and Puri (1973) and in Lyness (1978). In this case the coefficients have a more complicated structure than that in (A.1). However, in the special case of these trapezoidal-type rules R_{λ}, a simpler formula is available. For example, when λ = (1, 0, 0, 1),

B_{j}(R_{λ}; T; g) = ∑_{j_{i}≥0^{j_{1}+j_{2}=j}} c_{j_{1},j_{2}}(R_{λ}) ∫_{0}^{1} \frac{∂^{j_{1}}}{∂u^{j_{1}}} \{∫_{0}^{1−u} \frac{∂^{j_{2}}}{∂v^{j_{2}}} gdv\} du  (A.4)

with c_{j_{1},j_{2}}(R_{λ}) given by (A.3) above. When j_{1} = 0, the j_{1}, j_{2} term in the sum reduces to an integral over T, in direct analogy to (A.1) above. But otherwise the effect of differentiation with respect to the limit of an integral is to introduce additional terms.

In Section 5 the derivation of curvilinear rules Q(τ)f having h^{2}-error expansions rested on Lemma 5.4. We shall establish this lemma by showing two lemmas.

Lemma A.1 The operator (null rule) defined by

R_{λ−}^{(m)}(T)g = \frac{1}{2} \{R_{λ}^{(m)}(T)g − R_{−λ}^{(m)}(T)g\}  (A.5)

has an error expansion involving only odd powers of h = 1/m.

Lemma A.2 The rule defined by

R_{λ+}^{(m)}(T)g = \frac{1}{2} \{R_{λ}^{(m)}(T)g + R_{−λ}^{(m)}(T)g\}  (A.6)

has an h^{2}-error expansion.

Proof. We treat only the case λ = (1, 0, 0, 1). We recall that both R_{λ}^{(m)} and R_{−λ}^{(m)} apply a weight 1/m^{2} to every strictly interior point. In addition, R_{λ}^{(m)} applies this weight to all points on edges x = 0 and y = 0, except for points (1,0) and (0,1). On the other hand, R_{−λ}^{(m)} applies the same weight only to all points on the edge x + y = 1, except for the endpoints. Applying these assignments to the operators in (A.5) and (A.6), we find

R_{λ+}^{(m)}g = \frac{1}{m^{2}} ∑_{j>k>0^{j+k<m}} g (\frac{j}{m}, \frac{k}{m}) + \frac{1}{2m^{2}} ∑_{j=1}^{m−1} (g (\frac{j}{m}, 0) + g (0, \frac{j}{m}) + g (\frac{j}{m}, \frac{m − j}{m} ))

+\frac{1}{2m^{2}} g(0, 0),  (A.7)

while

R_{λ−}^{(m)}g = \frac{1}{2m^{2}} ∑_{j=1}^{m−1} (g (\frac{j}{m}, 0) + g (0, \frac{j}{m}) − g (\frac{j}{m}, \frac{m − j}{m} )) + \frac{1}{2m^{2}} g(0, 0)  (A.8)

which can be reexpressed in the form

2mR_{λ−}^{(m)}g = \frac{1}{m} ∑_{j=0}^{m}′′ (g (\frac{j}{m}, 0) + g (0, \frac{j}{m}) − g (\frac{j}{m}, \frac{m − j}{m} )) .  (A.9)

The right-hand side comprises three one-dimensional trapezoidal rule operators; since these are symmetric, each has an h^{2}-error expansion. Thus, the expansion for R_{λ−}^{(m)}g involves only odd inverse powers of m, establishing Lemma A.1.

In Lyness and Puri (1973), Section 7, rules for the triangle of the same general nature as (A.7) are discussed in detail. In particular, a rule denoted there by R_{s}^{(m)}g is defined that has an even expansion. This rule differs from R_{λ+}^{(m)}(T)g above only in that the weights attached to vertices are different. Specifically,

R_{s}^{(m)}g = R_{λ+}^{(m)}(T)g + \frac{1}{8m^{2}} (g(1, 0) + g(0, 1) − 2g(0, 0)).

These variant weights simply modify the coefficient of 1/m^{2} in the expansion; its even nature is, of course, not compromised This establishes Lemma A.2 above.

As a convenience to the reader, we note that Equations (3.5) through (3.10) may be abbreviated to

R_{(α,β,γ,δ)}g = g(ζ, η),

where, so long as det λ = 1, we have

ζ = 1 or 0 according as β > α or β < α η = 1 or 0 according as γ > δ or γ < δ.

Note that β + γ = 0; α ≠ β; γ ≠ δ; ζ − η = α + β − γ − δ.

References

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. Courant (1956). Differential and Integral Calculus, Volume 2, Interscience Publishers, Inc., New York. K. Georg (1991). Approximation of integrals for boundary element methods, SIAM J. Sci. Stat. Comput. 12, 443–453. K. Georg and J. Tausch (1993). Some error estimates for the numerical approximation of surface integrals, Math. Comp. (to appear). J. N. Lyness (1978). Quadrature over a simplex, SIAM J. Numer. Anal. 15, 122–133 and 870–877. J. N. Lyness and K. K. Puri (1973). The Euler-Maclaurin expansion for the simplex, Math. Comp. 27, 273–293. C. Schwab and W. L. Wendland (1992). On numerical cubatures of singular surface integrals in boundary element methods, Numer. Math. 62, 343–369. A. H. Stroud (1971). Approximate Calculation of Multiple Integrals, Prentice-Hall, Inc., Englewood Cliffs, New Jersey. P. Verlinden and R. Cools (1993). Proof of a conjectured asymptotic expansion for the approximation of surface integrals, Report TW189, KU Leuven 1993.