Quadrature Error Expansions – Part II: The Full Corner Singularity

J. N. Lyness^{1} Mathematics and Computer Science Division Argonne National Laboratory, Argonne, IL 60439, U.S.A. lyness@mcs.anl.gov

and

E. de Doncker Department of Computer Science Western Michigan University Kalamazoo, MI 49008 elise@cs.wmich.edu

Abstract (Summary):

We continue the work of Part I, treating in detail the theory of numerical quadrature over a square [0, 1)^{2} using an m^{2} copy, Q^{(m)}, of a one-point quadrature rule. As before, we determine the nature of an asymptotic expansion for the quadrature error functional Q^{(m)}F − IF in inverse powers of m and related functions, valid for specified classes of the integrand function F. The extreme case treated here is one in which the integrand function has a full-corner algebraic singularity. This has the form x^{λ}y^{μ}r_{ρ}(x, y). Here λ, μ, and ρ need not be integer, and r_{ρ} is (x^{2} + y^{2})^{\frac{1}{2}} or some other similar homogeneous function. The error expansion forms the theoretic basis for the use of extrapolation for this kind of integrand.

Mathematics Subject Classifications (1991). 65D32, 65B15

1 Introduction

The successful practice of extrapolation in quadrature requires that it be based on the expansion of the error functional Q^{(m)}f −If that corresponds to the nature of the integrand function. Romberg integration is based on the Euler-Maclaurin expansion and so is efficient for well-behaved functions. The construction of error functional expansions for other integrand functions is no trivial task. Recently, however, significant progress has been reported. The theory has been extended to integrand functions having end-point algebraic and logarithmic singularities (Navot [5] [6]), to multidimensional hypercubes and simplices, and to some of the more sophisticated multidimensional algebraic and logarithmic singularities (Lyness [2], Sidi [7]).

1This 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.

The present paper is a contribution to this latter theory. In Part I (Lyness and de Doncker-Kapenga [3]) we assembled some known results relating to two-dimensional quadrature. These were mainly minor extensions of one-dimensional results of Navot and other results of Sidi. There we established the notational framework that we use here extensively to handle a more intransigent singularity, which we term the full corner singularity. To keep this part self-contained, we have redefined here all terms used. The reader should refer to Part I for motivation and discussion of these terms.

We denote by

Q_{αβ}^{(m)}G = \frac{1}{m^{2}} ∑_{j=0}^{m−1} ∑_{k=0}^{m−1} G(\frac{α + j}{m} , \frac{β + k}{m}) (1)

the m^{2} copy of the one-point quadrature rule Q_{αβ}G = G(α, β), which approximates the integral IG of G(x, y) over the unit square

H_{2} : 0 ≤ x < 1, 0 ≤ y < 1. (2)

We assume that (α, β) ∈ H_{2}, but this assumption is not essential. For some classes of function, expansions for the error functional Q_{αβ}^{(m)} −IG are known. The most familiar error functional expansion is the Euler-Maclaurin expansion. A two-dimensional version is given as Theorem 1 below.

Theorem 1. When G(x, y) and its partial derivatives of order p are integrable over H_{2}, then

Q_{αβ}G = IG + ∑_{s=1}^{p−1} \frac{B_{s}(H_{2}, Q_{αβ}, G)}{m^{s}} + R_{p}(m; H_{2}, Q_{αβ}, G), (3)

where R_{p} = O(m^{−p}).

As is well known, the coefficients can be expressed in terms of double integrals. We use a concise notation for these, which was introduced in Part I and is useful subsequently. We set

Definition 1 (I.3.19)

t_{j,k}^{[θ^{1},θ^{2}]}(m, F)_{a c}^{b d} = ∫_{a}^{b} dx ∫_{c}^{d} dy h_{j}^{[θ^{1}]}(α, mx) h_{k}^{[θ^{2}]}(β, my) F^{(j,k)}(x, y). (4)

Here θ^{1} and θ^{2} take the values 1 or 0 only. The quadrature rule parameters enter through the kernel functions

\begin{aligned}h_{s}^{[1]}(β, t) & = & c_{s}(β) = B_{s}(β)/(s!), & s ≥ 0, \\ h_{s}^{[0]}(β, t) & = & (B_{s}(β) − \overline{B}_{s}(β − t))/(s!), & s ≥ 1, \end{aligned}

where B_{s}(x) and \overline{B}_{s}(x) are the sth Bernoulli polynomial and its periodic extension. In terms of these, we have

B_{s}(H_{2}, Q_{αβ}, G) = ∑_{k=0}^{s} t_{k,s−k}^{[1,1]}(1, G)_{0 0}^{1 1} (5)

and

R_{p}(m; H_{2}, Q_{αβ}, G) = \frac{1}{m^{p}} ∑_{k=0}^{p} t_{k,p−k}^{[..]}(1, G)_{0 0}^{1 1}. (6)

Here the superscript [..], which occurs in most remainder terms in this paper, is simply an abbreviation for an index pair [θ_{k,p}^{1}, θ_{p−k,p}^{2}] each member of which is 1 or 0. Possible choices are specified in Definition 3.1 of Part I; one choice is

[θ_{k,i+k}^{1}, θ_{i,i+k}^{2}] = \{ \begin{matrix}[0, 1] & k − i ≥ 2, \\ [0, 0] & k − i = 1, 0, \\ [1, 0] & k − i ≤ −1. \end{matrix}

The actual choice plays no essential role in the present argument.

Theorem 2 below is a two-dimensional generalization of a theorem first proved by Navot [5] [6] . It deals with an integrand function with a singularity on an edge of the square. Theorem 2 was proved by Sidi [7]. The conditions stated below are minor variants of those used by him.

Theorem 2. Let

Γ_{λ}(x, y) = x^{λ}G(x, y), λ > −1,

where G^{(r,s)}(x, y) is continuous over the square when 0 ≤ r ≤ p and 0 ≤ s ≤ ⌊p − λ⌋, p ≥ 1. Then

Q_{αβ}^{(r,s)}Γ_{λ} = ∑_{s=0}^{p−1} \frac{B_{s}(H_{2}, Q_{αβ}, Γ_{λ})}{m^{s}} + ∑_{t=0}^{⌊p−λ⌋−1} \frac{E_{λ+1+t}(H_{2}, Q_{αβ}, Γ_{λ})}{m^{λ+1+t}}
+R_{p}(m; H_{2}, Q_{αβ}, Γ_{λ}) + R̃_{p}^{(1)}(m; H_{2}, Q_{αβ}, Γ_{λ}),  (7)

with R_{p}^{(1)} =O(m^{−p}) and R̃_{p}^{(1)} =O(m^{−p−{λ}}).

For later convenience we have divided the remainder term into two parts of different orders, given explicitly in (11) and (12) below.

In both of the formulas (3) and (7), the coefficients B_{s} and E_{λ+1+t} are independent of m but depend on α and β. In (7), B_{s} is not generally given by the same formula (5) as before, both because it is not valid and also because, if applied, it does not converge. In Part I we showed how to circumvent this difficulty by a technique akin to subtracting out the singularity. The subtraction functions used are ϕ_{w}^{1}(x, y) and ϕ_{w}^{2}(x, y). The first is the function given in (I.3.29) and designed so that, when 0 < a < b < ∞, the function ϕ_{⌊k−λ⌋}^{1}^{(k,q)}(x, y) is integrable over [1, ∞) × [a, b], and

F^{(k,q)}(x, y) − ϕ_{⌊k−λ⌋}^{1}^{(k,q)}(x, y) is integrable over [0, 1] × [a, b]. The second is analogous, the roles of x and y being reversed.

Definition 2

\begin{matrix}ϕ_{w}^{1}(x, y) & = & ∑_{ℓ=0}^{w−1} F_{λ+ℓ}^{1}(x, y) & = & ∑_{ℓ=0}^{w−1} x^{λ+ℓ} G^{(ℓ,0)}(0, y)/(ℓ!), \\ ϕ_{w}^{2}(x, y) & = & ∑_{ℓ=0}^{w−1} F_{μ+ℓ}^{2}(x, y) & = & ∑_{ℓ=0}^{w−1} y^{μ+ℓ} G^{(0,ℓ)}(x, 0)/(ℓ!). \end{matrix} (8)

Note that F_{λ+ℓ}^{1}(x, y) is of the form assumed by Theorem I.3.2 (I.3.13), as long as G^{(ℓ,0)}(0, y) has a sufficient number of integrable derivatives.

Simple integral representations for all the coefficients in Theorem 2 may be stated in terms of the integrals in Definition 1. These are

B_{s}(H_{2}, Q_{αβ}, Γ_{λ}) = ∑_{k=0}^{s} {t_{k,s−k}^{[1,1]}(1, ϕ_{⌊k−λ⌋}^{1})_{∞ 0}^{1 1} + t_{k,s−k}^{[1,1]}(1, F − ϕ_{⌊k−λ⌋}^{1})_{0 0}^{1 1}} , (9)

E_{λ+1+t}^{(1)}(H_{2}, Q_{αβ}, Γ_{λ}) = ∑_{ℓ=0}^{t} U_{t−ℓ}^{[1]}(1, F_{λ+ℓ}^{1})_{0}^{1}, (10)

R_{p}(H_{2}, Q_{αβ}, Γ_{λ}) = \frac{1}{m^{p}} ∑_{k=0}^{p} \{t_{k,p−k}^{[. .]}(1, ϕ_{⌊k−λ⌋}^{1})_{∞ 0}^{1 1} + t_{k,s−k}^{[. .]}(1, Γ_{λ} − ϕ_{⌊k−λ⌋}^{1})_{0 0}^{1 1}\} , (11)

R̃_{p}^{(1)}(H_{2}, Q_{αβ}, Γ_{λ}) = \frac{1}{m^{p+λ}} ∑_{ℓ=0}^{⌊p−λ⌋−1} θ_{⌊λ+ℓ+1⌋,p}^{1} U_{p−⌊λ+ℓ+1⌋}(1, F_{λ+ℓ}^{1})_{0}^{1}. (12)

Here we have set

U_{k}^{[θ^{2}]}(m, F_{λ′})_{c}^{d} = t_{\overline{s},k}^{[0,θ^{2}]}(m, F_{λ′})_{0 c}^{1 d} − ∑_{s=\overline{s}}^{\overline{p}−1} t_{s,k}^{[1,θ^{2}]}(m, F_{λ′})_{a_{s−λ′−1}c}^{1  d} − t_{\overline{p},k}^{[0,θ^{2}]}(m, F_{λ′})_{∞ c}^{1 d},  (13)

where \overline{s} and \overline{p} can take any values satisfying 0 ≤ \overline{s} < λ′ + 1 < \overline{p}, a_{ν} = 0, 2, ∞ when ν < 0, ν = 0 or ν > 0 respectively, and F_{λ′} is of the form x^{λ′}G(y), with G(y) and its first p derivatives integrable. Note that in context, c and d are finite and the explicit choice of ν = s − λ′ − 1 in the integration limit a_{ν} ensures that the integral exists. When λ′ is integer, the terms on the right of (13) are zero when the order of the x-derivative exceeds λ′. In this case, U_{k}^{[θ^{2}]}(m, F_{λ′})_{c}^{d} is also zero. In the next section, we shall use expression (13) with the function F_{λ′+ℓ}^{1}(x, y) as integrand. There, when c and d are omitted, it is assumed that d = 1 and c = a_{k−(μ+ρ−ℓ)−1}; that is, c is chosen to make the integral

converge in the y direction. A similar function, V_{k}^{[θ^{1}]}(m, F_{μ})_{a}^{b}, is defined in (I.3.21a). This may be obtained from (13) by interchanging the roles of x and y.

We note that the validity of Theorem 2 above can also be established under the less strict condition of continuity of the partial derivatives G^{(r,s)}(x, y) whenever 0 ≤ r + s ≤ ⌊p − λ⌋.

The principal purpose of Part I was not to prove these theorems, both of which are well known, but to establish a unified notation and, using this, to establish explicit integral representations both for the coefficients B_{s} and E_{λ+1+t} and for the remainder terms R_{p}. In this part, we use these integral forms to obtain corresponding expansions for functions having more sophisticated singularities.

The reader interested in verifying our subsequent results in detail should consult Part I to obtain the details of some of the conventions we employ. However, to follow the thrust of the developments, the reader need only bear in mind that (4) is merely an abbreviated notation for a two-dimensional integral, the important aspects of which are the region [a, b] × [c, d] and the factor F^{(j,k)}(x, y) in the integrand function.

Our fundamental result, obtained in Sections 2 and 3, is for the integrand function

F (x, y) = x^{λ}y^{μ}r_{ρ}(x, y),

where r_{ρ} is homogeneous of order ρ (see Definition 1 below). We note that when μ and ρ/2 are integer, the integrand function is of the correct form for Sidi’s theorem to apply; while when λ and μ are integer, the expansion is also known (Lyness [2] [1]). Our derivation of the expansion for general λ, μ, and ρ is by no means straightforward. It requires nearly all the results in Part I and constitutes an extensive generalization of the approach in [2].

2 The Full-Corner Algebraic Singularity – Elements of the Expansion

In this section and in Section 3, we treat the integrand function

F (x, y) = x^{λ}y^{μ}r_{ρ}(x, y), λ > −1, μ > −1, λ + μ + ρ > −2, (14)

where r_{ρ}(x, y) is homogeneous of degree ρ and has no singularity in H_{2} other than possibly at the origin.

Definition 3. A function f(x, y) is homogeneous of degree ρ about the origin if

f(mx, my) = m^{ρ}f(x, y) for all m > 0, (x, y) ≠ (0, 0).

We shall treat homogeneous functions r_{ρ}(x, y) whose derivatives of all orders exist in H_{2} = [0, 1)^{2} except possibly at the origin.

Simple examples of r_{ρ}(x, y) include

r_{ρ}(x, y) = r^{ρ}, (γx + ηy)^{ρ}, with γ, η > 0, and r^{ρ}h(θ),

where r^{2} = x^{2} + y^{2}, θ = arctan(y/x), and h(θ) is analytic in θ for 0 ≤ θ ≤ π/2.

The error expansion for Q_{αβ}F − IF given at the end of this section, in (25–27), is not in Euler-Maclaurin form (in inverse powers of m and with coefficients independent of m). The additional work in Section 3 is needed to obtain the standard Euler-Maclaurin form of the expansion.

For the subsequent development, we subdivide the plane ℜ^{2} into squares of unit side.

Definition 4

H_{2}(k_{1}, k_{2}) : k_{1} ≤ x < k_{1} + 1; k_{2} ≤ y < k_{2} + 1.

H_{2}(0, 0) coincides with H_{2} defined in (2) above. We adopt the following approach. Let us rescale the problem of integration over H_{2} using an m^{2} copy of Q to one of integration over [0, m)×[0, m). In this context the quadrature rule used over each of the m^{2} squares [k_{1}, k_{1} + 1) × [k_{2}, k_{2} + 1) is a translated version of the rule Q over H_{2}. We deal with each square individually and assemble the results. It is convenient to define four different rectangles.

Definition 5

\begin{aligned}R(1, 1) & = & [0, 1) × [0, 1) & = & H_{2}(0, 0) = H_{2}, \\ R(1, m) & = & [0, 1) × [1, m) & = & ∪_{k_{2}=1}^{m−1}H_{2}(0, k_{2}), \\ R(m, 1) & = & [1, m) × [0, 1) & = & ∪_{k_{1}=1}^{m−1}H_{2}(k_{1}, 0), \\ R(m, m) & = & [1, m) × [1, m) & = & ∪_{k_{1}=1}^{m−1} ∪_{k_{2}=1}^{m−1} H_{2}(k_{1}, k_{2}). \end{aligned}

These comprise 1, m − 1, m − 1, and (m − 1)^{2} unit squares, respectively.

We now present the calculation of an expansion for Q_{αβ}^{(m)}F. Since F (x, y) given by (14) is homogeneous of degree λ + μ + ρ, we may set

F (\frac{α + k_{1}}{m} , \frac{β + k_{2}}{m}) = m^{−(λ+μ+ρ)}F (α + k_{1}, β + k_{2}).

Then, using (1), we may scale the problem as follows:

Q_{αβ}^{(m)}F = \frac{1}{m^{2}} ∑_{k_{1}=0}^{m−1} ∑_{k_{2}=0}^{m−1} F (\frac{α + k_{1}}{m} , \frac{β + k_{2}}{m})

= m^{−(λ+μ+ρ+2)} ∑_{k_{1}=0}^{m−1} ∑_{k_{2}=0}^{m−1} F (α + k_{1}, β + k_{2})  (15)

= m^{−(λ+μ+ρ+2)} ∑_{k_{1}=0}^{m−1} ∑_{k_{2}=0}^{m−1} Q_{α,β}(H_{2}(k_{1}, k_{2}))F.

This expresses Q_{αβ}^{(m)}(H_{2}, F) as a sum of one-point rules applied to each of m^{2} unit squares.

The integrand function has no singularity in R(m, m). In all the (m − 1)^{2} squares in R(m, m) we employ the standard two-dimensional Euler-Maclaurin expansion (3) above. In each of the m − 1 squares of R(1, m), the integrand has an x^{λ} singularity on the left-hand edge. Here we employ expansion (7) above, which explicitly treats this singularity. The m − 1 squares of R(m, 1) are treated by using an analogous expansion. The anomalous contribution from R(1, 1) is already in the correct form to contribute to the term in m^{−(λ+μ+ρ+2)} of the expansion.

We now treat the terms in (16) corresponding to the squares in R(1, m). As mentioned above, we apply (7) to each of the squares to obtain a result of the form

∑_{k_{2}=1}^{m−1} Q_{αβ}(H_{2}(0, k_{2}))F = ∑_{s=0}^{p−1} B̃_{s}(1, m) + ∑_{t=0}^{⌊p−λ⌋−1} Ẽ_{λ+1+t}^{(1)}(1, m) + R̃_{p}(1, m) + R̃_{p}^{(1)}(1, m). (16)

To proceed, we need the explicit forms for the coefficients in (7) given in (9–12). Inserting these into (7) and collecting like terms, we find

B̃_{s}(1, m) = ∑_{k=0}^{s} {t_{k,s−k}^{[1,1]}(1, ϕ_{⌊k−λ⌋}^{1})_{∞ 1}^{1 m} + t_{k,s−k}^{[1,1]}(1, F − ϕ_{⌊−λ⌋}^{1})_{0 1}^{1 m}} , (17)

Ẽ_{λ+1+t}^{(1)}(1, m) = ∑_{ℓ=0}^{t} U_{t−ℓ}^{[1]}(1, F_{λ+ℓ}^{1})_{1}^{m}, (18)

R̃_{p}(1, m) = ∑_{k=0}^{p} {t_{k,p−k}^{[. .]}(1, ϕ_{⌊k−λ⌋}^{1})_{∞ 1}^{1 m} + t_{k,s−k}^{[. .]}(1, F − ϕ_{⌊k−λ⌋}^{1})_{0 1}^{1 m}} , (19)

R̃_{p}^{(1)}(1, m) = ∑_{ℓ=0}^{⌊p−λ⌋−1} θ_{⌊λ+ℓ+1⌋,p}^{1} U_{p−⌊λ+ℓ+1⌋}(1, F_{λ+ℓ}^{1})_{1}^{m}. (20)

We note that these are obtained from (9–12) by simply setting m = 1 in (11) and (12) and replacing the integration limits (0, 1) in the y variable by (1, m). The justification for this procedure lies partly in that the kernel functions h_{s}^{[0]}(β, y) are periodic with period 1. Note also that the m dependence does not occur in the integrand functions, but only in the integration limits, as we are adding together m − 1 contributions each of which is m independent.

Using the same procedure for the set of terms in (16) corresponding to squares lying in R(m, 1) gives a result of the same character as (16), namely,

∑_{k_{1}=1}^{m−1} Q_{αβ}(H_{2}(k_{1}, 0))F = ∑_{s=0}^{p−1} B̃_{s}(m, 1) + ∑_{t=0}^{⌊p−μ⌋−1} Ẽ_{μ+1+t}^{(2)}(m, 1) + R̃_{p}(m, 1) + R̃_{p}^{(2)}(m, 1). (21)

The formulas for Ẽ_{μ+1+t}^{2} and R̃_{p}^{(2)} correspond to (18) and (20); μ replaces λ, θ^{2} replaces θ^{1}, and a V coefficient replaces a U coefficient. The formulas for B̃_{s}(m, 1) and R̃_{p}(m, 1) will be needed again in (26–27) below.

An analogous procedure applied to the set of squares lying in R(m, m) gives a simpler result; using the standard two-dimensional Euler-Maclaurin expansion (3), we find

∑_{k_{1}=1}^{m−1} ∑_{k_{2}=1}^{m−1} Q(H_{2}(k_{1}, k_{2}))F = ∑_{s=0}^{p−1} B̃_{s}(m, m) + R̃_{p}(m, m), (22)

with

B̃_{s}(m, m) = ∑_{k=0}^{s} t_{k,s−k}^{[1,1]}(1, F)_{1 1}^{m m} (23)

and

R̃_{p}(m, m) = ∑_{k=0}^{p} t_{k,p−k}^{[. .]}(1, F)_{1 1}^{m m}. (24)

Thus, substituting (16), (21), and (22) into the right-hand side of (16), we obtain

Q_{αβ}^{(m)}(H_{2})F = \frac{1}{m^{λ+μ+ρ+2}} \{Q_{αβ}(H_{2})F + ∑_{s=0}^{p−1} B̃_{s} + ∑_{t=0}^{⌊p−λ⌋−1} Ẽ_{λ+1+t}^{(1)}(1, m)

+ ∑_{t=0}^{⌊p−μ⌋−1} Ẽ_{μ+1+t}^{(2)}(m, 1) + R̃_{p} + R̃_{p}^{(1)}(1, m) + R̃_{p}^{(2)}(m, 1)},  (25)

where

B̃_{s} = B̃_{s}(1, m) + B̃_{s}(m, m) + B̃_{s}(m, 1) (26)

and

R̃_{p} = R̃_{p}(1, m) + R̃_{p}(m, m) + R̃_{p}(m, 1). (27)

Explicit integral representations now have been given for all of the above coefficients that carry a tilde and have the argument (m, m) or (1, m). These appear in (17–20), (23), and (24). Integral representations for coefficients with the argument (m, 1) have not been given explicitly but can be obtained from (17–20) by interchanging the roles of x and y.

3 Error Expansion in Euler-Maclaurin Form

The error expansion (25) involves many terms, defined by integral representations having various integrands over several different integration regions. These terms depend on m through the region of integration. The present section is devoted to the somewhat lengthy task of re-expressing the coefficients as sums of inverse powers of m and collecting like terms. The result is Theorem 3 below.

We treat first the contribution of the two terms (18) and (20) to the sum in (25). This contribution is given in (31) below, the relevant m independent coefficients being defined in (32–35). All these involve the U coefficient defined in (13). To this end we establish the following lemma.

Lemma 1

\begin{matrix}U_{t−ℓ}^{[1]}(1, F_{λ+ℓ}^{1})_{1}^{m} & = & (m^{μ+ρ+1−t} − 1)U_{t−ℓ}^{[1]}(1, F_{λ+ℓ}^{1}), μ + ρ + 1 − t ≠ 0, \\ & = & − log_{2} m U_{t−l}^{[1]}(1, F_{λ+ℓ}^{1})_{2}^{1}, μ + ρ + 1 − t = 0; \end{matrix} (28)

\begin{matrix}U_{p−⌊λ⌋−ℓ−1}^{[0]}(1, F_{λ+ℓ}^{1})_{1}^{m} & = & −U_{p−⌊λ⌋−ℓ−1}^{[0]}(1, F_{λ+ℓ}^{1})_{∞}^{1} \\ & & +m^{−p+[λ]+μ+ρ+2} U_{p−⌊λ⌋−ℓ−1}^{[0]}(m, F_{λ+ℓ}^{1})_{∞}^{1}, \\ & & p > λ + μ + ρ + 2. \end{matrix} (29)

Proof. By its definition (13), the U coefficient is a sum of terms each of the form t_{s,t−ℓ}^{[θ^{1},θ^{2}]}(1, F_{λ+ℓ}^{1})_{a_{ν}1}^{1 m}, with various values of s and θ^{1}. By definition (4) when θ^{2} = 1, the integral in the y coordinate in each of these integrals is of the form

∫_{1}^{m} F_{λ+ℓ}^{1}^{(s,t−ℓ)}(x, y)dy,

where F_{λ+ℓ}^{1}(x, y) is homogeneous of degree μ+ρ−ℓ in y. The key part of this proof is to note that every element of the subtraction function is homogeneous of degree λ + μ + ρ and that element F_{λ+ℓ}^{1}(x, y) is of the form x^{λ+ℓ}y^{μ+ρ−ℓ}. To establish (28), we need establish the same result for the integral above. This task is relatively simple. We may apply the first part of Lemma I.2.3 with λ replaced by μ + ρ − ℓ and s replaced by t − ℓ. This produces the result stated above for each t coefficient separately, and applying (13) again yields (28) as written. A similar treatment, using (I.2.15), yields (29).

□

Note that in the situation where μ + ρ − ℓ integer ≥ 0 and μ + ρ < t, both sides of (28) reduce to zero. In particular, when μ + ρ − t + 1 = 1, both sides of the second part of (28) are zero unless μ + ρ − ℓ = −1 and t − ℓ = 0.

Substituting (28) and (29) into (18) and (20) gives

Ẽ_{λ+1+t}^{(1)}(1, m) = ∑_{ℓ=0}^{t}(m^{μ+ρ+1−t} − 1)U_{t−ℓ}^{[1]}(1, F_{λ+ℓ}^{1}), μ + ρ + 1 − t ≠ 0,
= − log_{2} m ∑_{ℓ=0}^{t} U_{t−ℓ}^{[1]}(1, F_{λ+ℓ}^{1})_{2}^{1} = − log_{2} m U_{0}^{[1]}(1, F_{λ+ℓ}^{1})_{2}^{1},
μ + ρ + 1 − t = 0,

R̃_{λ+1+t}^{(1)}(1, m) =

− ∑_{ℓ=0}^{[p−λ]−1} θ_{[λ+ℓ+1],p}^{1} (U_{p−⌊λ⌋−ℓ−1}^{[0]}(1, F_{λ+ℓ}^{1})_{∞}^{1} − m^{−p+[λ]+μ+ρ+2} U_{p−⌊λ⌋−ℓ−1}^{[0]}(m, F_{λ+ℓ}^{1})_{∞}^{1}).

If we introduce

t^{∗} = μ + ρ + 1, (30)

the contribution of the above terms to the sum in (25) is

\frac{1}{m^{λ+1+t^{∗}}} \{ ∑_{t=0}^{[p−λ]−1} Ẽ_{λ+1+t}^{(1)}(1, m) + R̃_{p}^{(1)}(1, m)\} = ∑_{t=0}^{[p−λ]−1} \frac{E_{λ+1+t}^{(1)}(H_{2}, Q_{αβ}, F)}{m^{λ+1+t}}

+\frac{A_{λ+1+t^{∗}}^{(1)}(H_{2}, Q_{αβ}, F) + C_{λ+1+t^{∗}}^{(1)}(H_{2}, Q_{αβ}, F) log_{2} m}{m^{λ+1+t^{∗}}} + R_{p}^{(1)}(H_{2}, Q_{αβ}, F), (31)

where

E_{λ+1+t}^{(1)}(H_{2}, Q_{αβ}, F) = (1 − δ_{t−t^{∗}}) ∑_{ℓ=0}^{t} U_{t−ℓ}^{[1]}(1, F_{λ+ℓ}^{1}), (32)

R_{p}^{(1)}(m; H_{2}, Q_{αβ}, F) = m^{−p−⌊λ⌋} ∑_{ℓ=0}^{⌊p−λ⌋−1} θ_{⌊λ+ℓ+1⌋,p}^{1} U_{p−⌊λ⌋−ℓ−1}^{[0]}(m, F_{λ+ℓ}^{1})_{∞}^{1}, (33)

A_{λ+1+t^{∗}}^{(1)}(H_{2}, Q_{αβ}, F) = − ∑_{t=0}^{⌊p−λ⌋−1} E_{λ+1+t}^{(1)}(H_{2}, Q_{αβ}, F) − R_{p}^{(1)}(1; H_{2}, Q_{αβ}, F), (34)

C_{λ+1+t^{∗}}^{(1)}(H_{2}, Q_{αβ}, F) = −δ_{{t^{∗}}} ∑_{ℓ=0}^{t^{∗}} U_{t^{∗}−ℓ}^{[1]}(1, F_{λ+ℓ}^{1})_{2}^{1}
= −δ_{{t^{∗}}}U_{0}^{[1]}(1, F_{λ+ℓ}^{1})_{2}^{1}.  (35)

Here, δ_{0} = 1 and δ_{α} = 0 when α ≠ 0.

We note that the expansion (31) has a different form depending on whether t^{∗} = μ + ρ + 1 is an integer. Generally it is not an integer. The log_{2} m term is then absent from the expansion; and the expansion contains terms in m^{−(λ+1+t)}, t = 0, 1, ..., and a single term in m^{−(λ+1+t^{∗})}. However, when t^{∗} is an integer, one of the set of terms in m^{−(λ+1+t)} (the one with t = t^{∗}) apparently coincides in form with the single term in m^{−(λ+1+t^{∗})}. However, the factors (1 − δ_{t−t^{∗}}) and δ_{{t^{∗}}} in (32) and (35) resolve the situation by removing one of these terms and including instead a term in log_{2} m/m^{(λ+1+t^{∗})}. This is a familiar phenomenon.

At this point we have re-expressed the terms Ẽ_{λ+1+t}^{(1)}(1, m) and R̃_{p}^{(1)}(1, m) in (25) in the required form in (31). The terms Ẽ_{μ+1+t}^{(2)}(m, 1) and R̃_{p}^{(2)}(m, 1) can be expressed in an analogous manner. The term Q_{αβ}(H_{2})F is already independent of m. It remains then to deal with the B̃_{s} and R̃_{p}. The derivation for R̃_{p} resembles closely that for B̃_{s} and is not given completely. The derivation for the terms B̃_{s} is significantly more complicated than the one just given for Ẽ and R̃_{p}^{(1)}. The result is given in (48–51) below.

According to (27), R̃_{p} is the sum of terms given by (24), (19), and a formula analogous to (19). These involve in total many different integrals, each of which is expressed in terms of (4). For manipulative convenience we introduce more symbols.

Definition 6

\begin{matrix}T_{k,ℓ}^{[θ^{1},θ^{2}]}(μ, F, Φ, Ψ)_{m} = t_{k,ℓ}^{[θ^{1},θ^{2}]}(μ, Φ)_{∞ 1}^{1 m} + t_{k,ℓ}^{[θ^{1},θ^{2}]}(μ, F − Φ)_{0 1}^{1 m} \\ +t_{k,ℓ}^{[θ^{1},θ^{2}]}(μ, F)_{1 1}^{m m} + t_{k,ℓ}^{[θ^{1},θ^{2}]}(μ, Ψ)_{1 ∞}^{m 1} + t_{k,ℓ}^{[θ^{1},θ^{2}]}(μ, F − Ψ)_{1 0}^{m 1}. \end{matrix} (36)

Here F, Φ, and Ψ are functions of x and y. Using this notation, we find (from (17), (19), (23), (24), (26), and (27)) that

B̃_{s} = ∑_{k=0}^{s} T_{k,s−k}^{[1,1]}(1, F, ϕ_{⌊k−λ⌋}^{1}, ϕ_{⌊s−k−μ⌋}^{2})_{m} (37)

and

R̃_{p} = ∑_{k=0}^{p} T_{k,p−k}^{[. .]}(1, F, ϕ_{⌊k−λ⌋}^{1})_{m}, (38)

where F, ϕ_{w}^{1}, and ϕ_{w}^{2} are given by (14) and (8).

To simplify (36), we first generalize it.

Definition 8

I_{[a,b)}(F, Φ, Ψ) = ∫_{a}^{b} dx ∫_{0}^{a} dy (F − Ψ) + ∫_{0}^{a} dx ∫_{a}^{b} dy(F − Φ)
+ ∫_{a}^{b} dx ∫_{a}^{b} dy (F − Φ − Ψ) + ∫_{a}^{b} dx ∫_{b}^{∞} dy (−Ψ)
+ ∫_{b}^{∞} dx ∫_{a}^{b} dy (−Φ).

It is straightforward to verify that

I_{[a,b)}(F, Φ, Ψ) + I_{[b,c)}(F, Φ, Ψ) = I_{[a,c)}(F, Φ, Ψ) (39)

and

T_{k,ℓ}^{[1,1]}(1, F, Φ, Ψ)_{m} = I_{[1,m)}(F^{(k,ℓ)}, Φ^{(k,ℓ)}, Ψ^{(k,ℓ)}).

The reason for introducing I_{[a,b)} is in part that the dependence of I_{[1,m)}(F^{(k,ℓ)}, Ψ^{(k,ℓ)}, Ψ^{(k,ℓ)}) on m, which occurs in the integration limits, can be isolated easily, as shown in the following lemma.

Lemma 2. If F_{δ}, Φ_{δ}, and Ψ_{δ} are homogeneous about the origin of degree δ, and if I_{[a,b)}(F_{δ}, Φ_{δ}, Ψ_{δ}) exists for finite a and b, then the integrals occurring on the right below also exist and

\begin{matrix}I_{[1,m)}(F_{δ}, Φ_{δ}, Ψ_{δ}) & = & (m^{δ+2} − 1)I_{[0,1)}(F_{δ}, Φ_{δ}, Ψ_{δ}), & δ > −2, \\ & = & log_{2} m I_{[1,2)}(F_{δ}, Φ_{δ}, Ψ_{δ}), & δ = −2, \\ & = & −(m^{δ+2} − 1)I_{[1,∞)}(F_{δ}, Φ_{δ}, Ψ_{δ}), & δ < −2. \end{matrix} (40)

Proof. The key to the proof is the identity

I_{[a,b)}(F_{δ}, Φ_{δ}, Ψ_{δ}) = m^{−δ−2}I_{[ma,mb)}(F_{δ}, Φ_{δ}, Ψ_{δ}). (41)

This follows as a simple exercise using the homogeneity property of the integrand functions. In view of (41) we may define

I_{[1,∞)}(F_{δ}, Φ_{δ}, Ψ_{δ}) = lim_{n→∞} ∑_{k=0}^{n} I_{[m^{k},m^{k+1})}(F_{δ}, Φ_{δ}, Ψ_{δ}). (42)

So long as the limits exist, we use (41) with m replaced by m^{k} and with [a, b) = [1, m) to re-express (42) in the form

I_{[1,∞)}(F_{δ}, Φ_{δ}, Ψ_{δ}) = lim_{n→∞} ∑_{k=0}^{n} m^{k(δ+2)}I_{[1,m)}(F_{δ}, Φ_{δ}, Ψ_{δ}). (43)

When δ < −2, the limit on the right exists, and the third relation given in the lemma follows immediately. The first relation may be established in a similar manner. When δ = −2, one may use (41) to show that I_{[a,b)}(F_{−2}, Φ_{−2}, Ψ_{−2}) is a function of b/a, say χ(b/a), and use (39) to show that this function has the property χ(x_{1}) + χ(x_{2}) = χ(x_{1}x_{2}). The only nontrivial function with this property is a multiple of log x, and this leads to the second relation in (40).

□

While Lemma 2 will be used to isolate the m dependence in B̃_{s} given by (37), the following lemma will allow us to deal with R̃_{p} given by (38).

Lemma 3. Let h(x⃗) = h(x, y) be bounded for all x, y and periodic with period 1. Let δ < −2.

Then

I_{[1,m)} (h(x⃗)F_{δ} , h(x⃗)Φ_{δ} , h(x⃗)Ψ_{δ}) = I_{[1,∞)} (h(x⃗)F_{δ} , h(x⃗)Φ_{δ} , h_{(}x⃗)Ψ_{δ})  (44)
−m^{δ+2} I_{[1,∞)} (h(mx⃗)F_{δ} , h(mx⃗)Φ_{δ} , h(mx⃗)Ψ_{δ}).

By direct application of Lemma 2 and Lemma 3 with the functions used in (17) and (19), and with the corresponding subtraction functions used along y = 0, we obtain

T_{k,s−k}^{[1,1]}(1, F, ϕ_{⌊k−λ⌋}^{1}, ϕ_{⌊s−k−μ⌋}^{2})_{m}

\begin{matrix}= (m^{λ+μ+ρ−s+2} − 1)T_{k,s−k}^{[1,1]}(1, F, ϕ_{⌊k−λ⌋}^{1}, ϕ_{⌊s−k−μ⌋}^{2})_{0}, & λ + μ + ρ − s > −2, \\ = − log_{2} m T_{k,s−k}^{[1,1]}(1, F, ϕ_{⌊k−λ⌋}^{1}, ϕ_{⌊s−k−μ⌋}^{2})_{2}, & λ + μ + ρ − s = −2, \\ = −(m^{λ+μ+ρ−s+2} − 1)T_{k,s−k}^{[1,1]}(1, F, ϕ_{⌊k−λ⌋}^{1}, ϕ_{⌊s−k−μ⌋}^{2})_{∞}, & λ + μ + ρ − s < −2; \end{matrix}

T_{k,p−k}^{[. .]}(1, F, ϕ_{⌊k−λ⌋}^{1}, ϕ_{⌊p−k−μ⌋}^{2})_{m} = T_{k,p−k}^{[. .]}(1, F, ϕ_{⌊k−λ⌋}^{1}, ϕ_{[p−k−μ]}^{2})_{∞}  (45)
−m^{λ+μ+ρ−p+2}T_{k,p−k}^{[. .]}(1, F, ϕ_{⌊k−λ⌋}^{1}, ϕ_{[p−k−μ]}^{2})_{∞}.

This puts us in a position to re-specify the contribution of B̃_{s} and R̃_{p} to the sum in (25). Setting

s^{∗} = λ + μ + ρ + 2, (46)

we obtain

\frac{1}{m^{λ+μ+ρ+2}} \{∑_{s=0}^{p−1} B̃_{s} + R̃_{p}\} = ∑_{s=0}^{p−1} \frac{B_{s}(H_{2}, Q_{αβ}, F)}{m^{s}} + \frac{A_{s^{∗}}^{(1,2)}(H_{2}, Q_{αβ}, F)}{m^{s∗}}  (47)

+\frac{C_{s^{∗}}^{(1,2)}(H_{2}, Q_{αβ}, F) log_{2} m}{m^{s∗}} + R_{p}^{(1,2)}(m; H_{2}, Q_{αβ}, F),

where

B_{s}(H_{2}, Q_{αβ}, F) = (1 − δ_{s−s∗}) ∑_{k=0}^{s} T_{k,s−k}^{[1,1]}(1, F, ϕ_{⌊k−λ⌋}^{1}, ϕ_{⌊s−k−μ⌋}^{2})_{a_{s−s^{∗}}} , (48)

R_{p}^{(1,2)}(m; H_{2}, Q_{αβ}, F) = −\frac{1}{m^{p}} ∑_{k=0}^{p} T_{k,p−k}^{[. .]}(1, F, ϕ_{⌊k−λ⌋}^{1}, ϕ_{[p−k−μ]}^{2})_{∞}, (49)

A_{s^{∗}}^{(1,2)}(H_{2}, Q_{αβ}, F) = − ∑_{s=0}^{p−1} B_{s}(H_{2}, Q_{αβ}, F) − R_{p}^{(1,2)}(1; H_{2}, Q_{αβ}, F), (50)

and

C_{s^{∗}}^{(1,2)}(H_{2}, Q_{αβ}, F) = −δ_{{s^{∗}}} ∑_{k=0}^{s^{∗}} T_{k,s^{∗}−k}^{[1,1]}(1, F, ϕ_{⌊k−λ⌋}^{1}, ϕ_{[s^{∗}−k−μ]}^{2})_{2}. (51)

Replacing the corresponding part of expansion (25) with (48) and taking (31) into account, we are now ready to reformulate (25) as an expansion in inverse powers of m. This leads us to the principal result of this paper.

Theorem 3. Let F (x, y) = x^{λ}y^{μ}r_{ρ}(x, y), λ, μ > −1, λ + μ + ρ > −2, and p > λ + μ + ρ + 2. Then

Q_{αβ}^{(m)}F =

∑_{s=0}^{p−1} \frac{B_{s}(H_{2}, Q_{αβ}, F)}{m^{s}} + \frac{A_{λ+μ+ρ+2}(H_{2}, Q_{αβ}, F) + C_{λ+μ+ρ+2}(H_{2}, Q_{αβ}, F) log_{2} m}{m^{λ+μ+ρ+2}}

+ ∑_{t=0}^{[p−λ]−1} \frac{E_{λ+1+t}^{(1)}(H_{2}, Q_{αβ}, F)}{m^{λ+1+t}} + ∑_{t=0}^{[p−μ]−1} \frac{E_{μ+1+t}^{(2)}(H_{2}, Q_{αβ}, F)}{m^{μ+1+t}} + R_{p}(m; H_{2}, Q_{αβ}, F),  (52)

where the coefficients B, A, C, E^{(1)}, and E^{(2)} are independent of m and

R_{p} = O(m^{−p}).
□

Two of the coefficients, B_{s} and E^{(1)}, have already been given explicitly (by (48) and (32), respectively). The rest of the coefficients and the remainder comprise

\begin{matrix}A_{λ+μ+ρ+2}(H_{2}, Q_{αβ}, F) & = & Q_{αβ}(H_{2}, Q_{αβ})F + A^{(1)} + A^{(1,2)} + A^{(2)}, \\ C_{λ+μ+ρ+2}(H_{2}, Q_{αβ}, F) & = & C^{(1)} + C^{(1,2)} + C^{(2)}, \\ R_{p}(m; H_{2}, Q_{αβ}, F) & = & R_{p}^{(1)} + R_{p}^{(1,2)} + R_{p}^{(2)}, \end{matrix} (53)

and E^{(2)}, respectively. The coefficients with superscript (1,2) are given explicitly by (49–51), the coefficients with superscript (1) by (32–35), and the coefficients with superscript (2) by formulas analogous to those with superscript (1).

Note furthermore that Theorem 3 is stated for a value of p > λ + μ + ρ + 2. However, for any p̄, 0 ≤ p̄ ≤ p, one could truncate the summations listed in (52) to retain terms whose order exceeds O(1/mp̄) only, and absorb all terms of orders ≤ O(1/mp̄) in a remainder term R̄p̄ = O(1/mp̄). Notationally, the case p̄ = 0 reduces to R̄_{0}(m; H_{2}, Q_{αβ}, F) = Q_{αβ}^{(m)}(H_{2})F.

We close this section with a remark about the conditions under which some of these terms in (52) vanish. The term of the form (C_{r} log m)/m^{r} occurs only when λ, μ, and ρ are such that one of s^{∗} or t^{∗} defined in (46) or (30) or w^{∗} = λ + ρ + 1 is integer. Detailed reference to these formulas reveals that C_{r} = 0 unless the same subscript r occurs in one of the three other summations in (52). (See the remarks just after equation (35).)

We remind the reader that the corresponding expansion for the m^{2} copy of a standard quadrature rule

Qf = ∑_{j=1}^{ν} ω_{j}f(α_{j}, β_{j}) (54)

is simply the weighted sum of m^{2} copies of ν distinct one-point quadrature rules Q_{α_{j}β_{j}}. Consequently we have treated here only the one-point rule. The reader should bear in mind that while assembling the expansion for the rule (54) is straightforward, there is a chance of some terms disappearing. For example, when Qf is symmetric, the terms

B_{s}(H_{2}, Q, f) = ∑_{j=1}^{ν} w_{j}B_{s}(H_{2}, Q_{α_{j}β_{j}}, f) (55)

occurring in (3) and (7) are zero for all odd s.

In addition, if the rule Q in (54) above is of polynomial degree d, B_{s} given by (55) vanishes for s ∈ [1, d]. In fact, the contribution to C_{s} denoted by C^{(1,2)} also satisfies these conditions. However, the other contributions to C_{s} do not. In the special case in which λ and μ are both integers, C_{s} satisfies these conditions. This is the case treated in Lyness [2].

4 Applications

These asymptotic expansions have a natural application in cubature by extrapolation, particularly in the context of the finite element. We have carried out extensive numerical experiments, both to verify the expansions numerically and to assess their utility in the cubature context.

We arrange this extrapolation so as to calculate the analogues of the terms that would occur in the classical Romberg table. We choose an increasing sequence of mesh ratios m_{0}, m_{1},... and assume an expansion containing the non-zero terms in (52). As in the T -table, the element T_{k,p} is an approximation to If based on p + 1 approximations Q^{(m_{j})}f with j ∈ [k, p + k]. Specifically, it is the solution for B_{0} of the set of p + 1 linear equations obtained by discarding all but the first p + 1 terms of (52). In standard Romberg integration, one may use the Neville algorithm to calculate elements of the table recursively. Here we have simply used a linear equation solver.

Naturally, the cost of cubature by extrapolation depends on the number of terms in the expansion that have to be eliminated. Any unnecessary term included is likely to increase the cost significantly. But any necessary term omitted slows down the convergence to a rate commensurate with that term.

Because of this cost pattern, we have to be concerned that, in any particular case, all displayed terms in (52) are needed. When we consider the special cases f(x, y) = x^{λ}, and y^{ν} and r^{ρ} separately, we see that any umbrella expression such as (52) has to include either these terms or analogous sequences of terms of a similar nature and number. Nevertheless, we have carried out many numerical examples to satisfy ourselves that at least the early terms are usually present. In doing so we found extensive numerical evidence that, under some readily recognizable circumstances, some do not occur. We express this as a conjecture.

Conjecture. When λ + μ + ρ = −1, the coefficient C_{1} occurring in (52) is zero, unless all of λ, μ, and ρ are integers.

We have no proof. Since r(x, y) is a general homogeneous function, there may be a class of these functions that we have overlooked for which the conjecture is not valid. However, since each additional term in the expansion adds significantly to the expense, we feel obliged to mention this possible economy.

The following example is included to give the reader some feeling for the large difference in cost that may be experienced in the same problem using variant mesh ratio sequences and expansions. It is anecdotal in nature and is typical of our experience.

Numerical Example. f(x, y) = (xyr^{4})^{−0.2}. Using the mid-square rule (α = β = 1/2), we see from (52) that the principal term in the error expansion is of form (log m)/m^{0.8}, and the other terms are m^{−j} with j = 0.8, 1.8, 2, 2.8, 3.8, 4, ... .

(i) The proper investigation of this expansion requires a highly accurate, reliable numerical approximation for the integral If. Making the substitution x = X^{5} and y = Y^{5} produces an easier integral over the same square. The new integrand function F (X, Y) = 25X^{3}Y^{3}/(X^{10} + Y^{10})^{0.4} is homogeneous of degree 2 and has a much simpler error expansion. This can be evaluated by extrapolation or by an adaptive quadrature routine.

(ii) Returning to the original example, our first numerical task was to verify the expansion numerically. Using the correct expansion, and the geometric sequence 1, 2, 4, 8, 16, … of mesh ratios, we can obtain nine-figure accuracy after nine iterations. The last approximation uses 256^{2} = 65, 536 function values, the total being about 87, 381. We can do better than this using the conventional F-sequence 1, 2, 3, 4, 6, 8, 12, … of mesh ratios which includes 1, 2, 3 and the double of any member already present. Here, we obtain ten figures after twelve iterations, the final one using 64^{2} = 4, 096 function values, the total being about 8, 530. The harmonic sequence 1, 2, 3, 4, 5, 6 … became unstable before we reached nine figures; while we cannot recommend them, for the record we note that we can find near-harmonic sequences that obtain nine-figure accuracy in this example using between 1200 and 2000 function values. In all cases, we found a T-table of the expected form, and no evidence that there could be a missing term in the expansion or that there was an unnecessary term. Nevertheless, several numerical experiments of the type described in (iii) and (iv) below were carried out.

(iii) We introduced an extra term B_{1}/m into the expansion. Using the same geometric or F-sequences as before gave results less accurate by about two decimal places. When we extended the calculation to obtain the same accuracy as before, we needed one extra extrapolation using the geometric sequence and two using the F-sequence. In both cases this additional extrapolation involved quadrupling the cost.

(iv) We omitted the term B_{2}/m^{2} from the expansion. This procedure essentially destroyed the calculation. Terminating at the same point as before, we found five or six correct figures in place of nine or ten. Our estimate of the cost of obtaining nine or ten ran into the millions.

Our conclusion is that this sort of integration shares many of the features of the classical twodimensional analogue of Romberg integration. So long as the proper expansion is used, the accuracy pattern in the T-table is much the same.

5 Concluding Remarks

In this paper, we have managed to avoid the use of asymptotic theorems. Naturally we have stated the order of remainder terms, which of course is vital to the purpose for which these expansions are obtained. However, in all cases, we have provided completely specified integral representations for the remainder terms. Thus, the asymptotic expansions have the status of identities. This can be of significant help when the theory comes to be extended, as it often happens that elements of some remainder term contribute to earlier terms in a more developed expansion.

The expansion for an integrand function having a singularity of the form

f(x, y) = x^{λ}y^{μ} (56)

is much simpler to derive. One obtains with little difficulty almost the same expansion as that obtained by setting ρ = 0 in (52). The difference is that in the corresponding expansion for (56) the logarithmic term is omitted, whether or not λ + μ is an integer. This term appears in (52) with a non-zero coefficient when λ + μ + ρ is an integer. This difference arises because, in two dimensions, there are homogeneous functions of degree zero that are not constant. These are functions of θ, the second polar coordinate. For this wider class of function, the extra logarithmic term is required. While even then this coefficient is generally zero, it is non-zero when λ + μ is an integer; this situation includes all the simple cases in which λ and μ are both integers.

The result of Section 3 may be readily extended and generalized in several ways. For example, an obvious extension is to the integrand function

f(x, y) = x^{λ}y^{μ}r_{ρ}(x, y)g(x, y)

where g(x, y) is regular, the other terms being subject to the restrictions of Theorem 3. This is

Q_{αβ}^{(m)}f ∼ ∑_{s=0} \frac{B_{s}}{m^{s}} + ∑_{j=0} \frac{A_{λ+μ+ρ+j+2} + C_{λ+μ+ρ+j+2} log_{2} m}{m^{λ+μ+ρ+j+2}}  (57)

+ ∑_{t=0} \frac{E_{λ+1+t}^{(1)}}{m^{λ+1+t}} + ∑_{t=0} \frac{E_{μ+1+t}^{(2)}}{m^{μ+1+t}} .

A derivation of integral representations for the expansion coefficients and the remainder term is not given in this paper. It is trivial to write these down in the case in which g(x, y) is a polynomial.

Extensions to integration over a triangle, rather than a square, follow the lines of a similar generalization described in Lyness and Monegato [4] and lead to an expansion of precisely the same form as (57) above. Other extensions to integrand functions having logarithmic singularities are clearly possible, but require detailed justification along the lines of Section 5 of Lyness [2]. Our understanding is that expansions of this type are needed to handle elements required in some of the more recent applications of the boundary element method.

Most of this work was carried out during a three-month visit by one author to Argonne National Laboratory in 1979. Our original proof was longer than the present one, and less general. Since then, we have searched at length for a less extended and detailed proof of these results, with only limited success. We have decided to discontinue this search temporarily; we have presented this proof in its present form in the expectation that others, armed with a clear description of and confidence in the results, may be more successful than we were in condensing it.

References

[1] Lyness, J. N. (1976): Applications of extrapolation techniques to multidimensional quadrature of some integrand functions with a singularity. Journal of Computational Physics 20, 346–364

[2] Lyness, J. N. (1976): An error functional expansion for N-dimensional quadrature with an integrand function singular at a point. Mathematics of Computation 30, 1–23

[3] Lyness, J. N., de Doncker-Kapenga, E. (1987): On quadrature error expansions. Part I. Journal of Computational and Applied Mathematics 17, 131–149

[4] Lyness, J. N., Monegato, G. (1980): Quadrature error expansions for the simplex when the integrand has singularities at vertices. Mathematics of Computation 34, 213–225

[5] Navot, I. (1961): An extension of the Euler-Maclaurin summation formula to functions with a branch singularity. Journal of Mathematics and Physics 40, 271–276

[6] Navot, I. (1962): A further extension of the Euler-Maclaurin summation formula. Journal of Mathematics and Physics 41, 155–163

[7] Sidi, A. (1983): Euler-Maclaurin expansions for integrals over triangles and squares of functions having algebraic/logarithmic singularities along an edge. Journal of Approximation Theory 39, 39–53 Appendix

Remarks and Errata On Quadrature Error Expansions. Part I. JCAM 17 131-149 (1987)

• P134. Sentence after Equation (2.13). Replace p ≤ λ by p ≤ λ + 1 (or λ integer).

• P135. Third equation of Lemma 2.3. Note that both sides reduce to 0 for s = λ + 1.

• P136. First sentence after the proof of Theorem 2.4. Replace p > λ + 1 by p ≥ λ + 1.

• P137. Equation (2.27). Replace subscript λ + 1 by λ + l.

• P138. Line before (2.38). Replace final (0) by (1).

• P139. Equation (2.41). Replace exponent μ + t + 1 by μ + t − 1. Equation (2.46). Replace -2 in the upper bound on the second summation by -1. Replace p > 1 by p ≥ 1.

• P142. Definition 3.1. Replace p > 1 by p ≥ 1. Sentence after Equation (3.9). Replace (3.4) by (3.5).

• P143. Theorem 3.2. Add the condition p ≥ 1.

• P144. Line after Equation (3.21a). Replace ρ by p̄.

• P145. Theorem 3.4 and Theorem 3.5. Add the condition p ≥ 1. Equation (3.31). Replace p as a subscript of Q in the first term on the right by β.

• P145. Before Theorem 3.5. Add the following sentence: Note that the following theorem was proved by Sidi (Journal of Approximation Theory 39, 1983, pp. 39-53), under slightly different conditions.

• P147. First sentence. Replace “leads” by “lead”.

• P147. Add before Section 4: Theorem 3.5 can be proved valid under the less strict condition of continuity on the derivatives G^{(r,s)}(x, y) for 0 ≤ r + s ≤ [p − λ].

• P148. The second displayed equation in Section 5 should contain on the right-hand side the additional term

∑_{i=1}^{2} ∑_{j=3}^{4} ∑_{t≥0} \frac{A_{λ_{i}+λ_{j}+t+2}}{m^{λ_{i}+λ_{j}+t+2}} .

• P136 (first sentence of Theorem 2.5), P139 (top sentence and first sentence of Theorem 2.8), P145 (Theorem 3.5). Replace the requirement “integrable” on the indicated derivatives by “continuous”.