                TECHNION -- ISRAEL INSTITUTE OF TECHNOLOGY
                         DEPARTMENT OF MATHEMATICS

Multipoint Taylor formulæ

running title: Multipoint Taylor formulæ

                Shayne Waldron^{∗} (waldron@math.auckland.ac.nz)

Permanent address: Department of Mathematics, University of Auckland Private Bag 92019, Auckland, New Zealand

                         Technical Summary Report
                                March 1996

Summary

The main result of this paper is an abstract version of the Kowalewski – Ciarlet – Wagschal multipoint Taylor formula for representing the pointwise error in multivariate Lagrange interpolation. Several applications of this result are given in the paper. The most important of these is the construction of a multipoint Taylor error formula for a general finite element, together with the corresponding L_{p}–error bounds. Another application is the construction of a family of error formulæ for linear interpolation (indexed by real measures of unit mass) which includes some recently obtained formulæ. It is also shown how the problem of constructing an error formula for Lagrange interpolation from a D– invariant space of polynomials with the property that it involves only derivatives which annihilate the interpolating space can be reduced to the problem of finding such a formula for a ‘simpler’ one–point interpolation map.

AMS (MOS) Subject Classifications: primary 41A65, 41A80, 65D05, secondary 41A05, 41A10, 41A44, 41A55

Key Words: multipoint Taylor formula, Kowalewski’s remainder, Lagrange interpolation, linear interpolation on a triangle, finite elements, multivariate Hardy’s inequality, Taylor interpolation, one–point interpolation, least solution

∗ Partially supported by the Steenbock Professorship of Carl de Boor and the Israel Council for Higher Education.

1. Introduction

Overview

This paper is concerned with explicit representations of the linear map Q which is defined implicitly as the unique map which makes the following diagram commute

where U : X → Y, R : X → Z are linear maps with

ker U ⊂ ker R.

Maps Q of this type occur frequently in approximation theory and numerical analysis. For example, if U = D^{k+1} and R is a linear functional whose kernel (ker R) contains Π_{k} (the polynomials of degree ≤ k), then Q gives the Peano kernel representation

Rf = Q(D^{k+1}f).

There is particular interest in certain such maps Q whose norm provides the best constant in error bounds for various numerical schemes. The most notable of these are the finite element method and numerical differentiation and integration rules (see, e.g., Meinguet [M84] and the references therein).

The main result of this paper is that for certain maps of the form

R = A(1 − L)

part of the flexibility in Meinguet’s abstract method for representing Q (detailed below) can be used to satisfy conditions that ensure the representation for Q take a simple form. A precise statement of this result is given in Theorem 2.2.

The prototype of such representations is the multipoint Taylor formula for Lagrange interpolation (see Ciarlet and Wagschal [CW71]). In Section 3, this is discussed and then extended (in its unexpanded form) to all linear operators that reproduce Π_{k}.

In Section 4, this extension is discussed within the setting of finite elements. It is expanded in terms of the nodal variables to obtain the analogue of the original presentation of the multipoint Taylor formula.

In Section 5, L_{p}–error bounds for finite elements are obtained from these ‘expanded’ multipoint Taylor formulæ by using the multivariate Hardy’s inequality recently given by the author (and intended for precisely this type of situation).

In Section 6, to emphasise the wide applicability of the main ‘abstract result’ and to illustrate some of its more subtle points, it is used to obtain a family of formulæ for the error in linear interpolation which is indexed by real measures of unit mass. These include many of the known formulæ (some of which have only recently been obtained).

Finally, in Section 7, a general form of Taylor interpolation is introduced which allows the results of the previous sections to be extended. This ‘one–point interpolation from an arbitrary polynomial space’ has close connections with the least solution of de Boor and Ron [BR92] to the ‘polynomial interpolation problem’. It is shown how error formulæ for the least solution can be obtained from those for a simpler one–point interpolation problem. Representing the error in both of these maps involves some deep questions in the theory of multivariate Peano kernels.

Meinguet’s abstract method for representing Q

Given the situation above, i.e.,

QU = R  (1.1)

with R, U known linear maps, how can an explicit representation of Q be obtained? The approach most often used is to take

V : ran U → X

a right inverse for U, i.e., a map with

UV = 1,

and right multiply (1.1) by V to obtain

Q = RV.  (1.2)

This is stage (i) of what is referred to in Meinguet [M84] as the abstract method for representing Q. The prototypical example of (1.2) is the Peano kernel representation.

Stage (ii) of the abstract method is the observation that if the restriction of R to the range of V can be written as a finite sum

R|_{ran V} = ∑_{j} S_{j},

of what are referred to as standard linear mappings S_{j}, then (1.2) can be expanded as

Q = ∑_{j} S_{j}V.  (1.3)

Now assume that V is linear and Y, Z are normed linear spaces. In the final stage of the abstract method, stage (iii), the triangle inequality is applied to (1.3) to obtain

∥Q∥ ≤ ∑_{j} ∥S_{j}V∥.  (1.4)

If R, U are continuous and U is open, then Q is a continuous linear map. This was first observed by Sard [S48:Th.1] when X, Y, Z are Banach spaces, and is referred to as Sard’s factorisation theorem by Atteia [At92:Th.2.1,p98] who states it for topological vector spaces.

Meinguet does not suggest how best to choose V and {S_{j}}, other than that it should be done in such a way that bounding each ∥S_{j}V∥ is simpler than bounding ∥RV∥ directly. He gives many examples in [M75], …, [M84]. Earlier examples of the abstract method in the multivariate setting include Èzrohi [E57] and Sard [S63].

In each case known to the author the right inverse V is linear and comes via an associated linear projector P onto ker U as described below.

Proposition 1.5 (cf [M84]). Let U : X → Y be a linear map. There is a 1–1 (linear) correspondence between linear right inverses

V : ran U → X

for U and linear projectors

P : X → ker U

onto ker U, which is given by P = 1 − VU.

2. The main result

Next is the main result, that if

R = A(1 − L),  (2.1)

then it may be possible to choose V so that the representation of Q given by Meinguet’s abstract method takes a simple form. Later it will be shown how many practical estimates follow from this ‘abstract result’ for operators of the form (2.1), such as the error in the finite element method (see Examples 1 and 2 of Section 4). The reader might like to keep in mind the simple example

Rf = δ_{x}(1 − L)f = f(x) − Lf(x)

where L is the operator of Lagrange interpolation from Π_{k} and U = D^{k+1}, ker U = Π_{k}.

Theorem 2.2. Let Q be the linear map which is defined implicitly as the unique map that makes the following diagram commute where U : X → Y, A : X → Z, L : X → X are linear maps with

ker U ⊂ ker R.

If V : ran U → X is a right inverse of U, chosen so that

AVU = 0,  (2.3)

then

Q = −ALV.  (2.4)

Proof: Expanding (1.2) using R = A(1 − L) gives

Q = AV − ALV.  (2.5)

If AVU = 0 (which is equivalent to AV = 0), then (2.5) reduces to (2.4).

Above, the condition AVU = 0 was favoured over the equivalent one that AV = 0. This is because usually the right inverse V is a linear map which comes via an associated linear projector P := 1 − VU (see Proposition 1.5), and this choice makes it clear that (2.3) can be rewritten in the (often more convenient) form

A = AP.  (2.6)

At times, such as when U is the multivariate differentiation operator D^{k+1}, it is convenient to use Theorem 2.2 without any direct reference to U or V. The following corollary allows this.

Corollary 2.7. Let

R := A(1 − L) : X → Z

where A : X → Z, L : X → X are linear maps. If P : X → X is a linear projector with

ran P ⊂ ker R,  (2.8)

and

A = AP,  (2.9)

then

R = −AL(1 − P).  (2.10)

Proof: This corollary follows by taking U = 1 − P and choosing V to be the identity. However, it may be more instructive to give the direct proof that:

R = R(1 − P)  (since ran P ⊂ ker R)
= A(1 − L)(1 − P) = A(1 − P) − AL(1 − P)
= −AL(1 − P) (since A = AP).

Examples

As far as the author is aware, there are two examples of (what are effectively) representations of Q of the form (2.4). One is the multipoint Taylor formula for the error in Lagrange interpolation, and certain kinds of Hermite interpolation, from a space containing Π_{k}. This is discussed and then extended in the next section. The other is Gregory’s error formula for linear interpolation on a (standard) triangle (see next section).

3. Multipoint Taylor formulæ for linear operators that reproduce Π_{k}

In this section it is shown that the multipoint Taylor formula for the error in a Lagrange map (see below) is a natural example of equation (2.4) of Theorem 2.2, and it is extended (in its unexpanded form) to all linear operators that reproduce Π_{k}.

Lagrange maps

We say that Lagrange interpolation from a space of functions P to data at a finite set of points Θ ⊂ IR^{n} is correct if for each function f (defined at least on Θ) there is a unique p ∈ P with

p(v) = f(v), ∀v ∈ Θ.

The associated linear projector f ↦ p is called the Lagrange map (given by the space of interpolants P and points of interpolation Θ), and it is denoted by

L_{P,Θ} : f ↦ p.

The Lagrange form of a Lagrange map L_{P,Θ} is

L_{P,Θ}f = ∑_{v∈Θ} f(v) ℓ_{v},  (3.1)

where (3.1) uniquely defines

ℓ_{v} := ℓ_{v,P,Θ} ∈ P,

the Lagrange function for v ∈ Θ.

Lagrange maps with a space of interpolants that contains Π_{k} give rise to many commonly used finite elements, and are also frequently used to interpolate to scattered data. Particular examples of the latter, which have received much attention lately, are maps where the interpolants include radial basis functions, and the least solution to the polynomial interpolation problem (see Section 7).

The multipoint Taylor formula

The multipoint Taylor formula is a representation of the error

R := A(1 − L)

in a Lagrange map

L := L_{P,Θ}

for which Π_{k} ⊂ P, and where A is any of the derivatives

A := A_{g,x} := δ_{x}g(D) : f ↦ g(D)f(x), g ∈ Π_{k}.

Here, and whenever point evaluation at x is being thought of as a linear map on functions (and so should be applied from the left) it will be denoted by

δ_{x} : f ↦ f(x).

The differential operator induced by g ∈ Π_{k} is written g(D).

If P is a linear projector onto Π_{k} which interpolates A, i.e., for which

g(D)f(x) = g(D)Pf(x), ∀f,  (3.2)

then by Corollary 2.7

R = −AL(1 − P).  (3.3)

Here, the intent is to apply Theorem 2.2 with U = D^{k+1} (the k + 1 order multivariate differentiation operator) to obtain expansions of R(f) in terms of D^{k+1}f.

The linear functionals

{A_{g,x} : g ∈ Π_{k}}

are the interpolation conditions for the operator of Taylor interpolation from Π_{k} at x, which we denote T_{k,x}. This operator is a natural choice for P since it is the only linear projector (onto Π_{k}) that interpolates A = A_{g,x} for all g ∈ Π_{k}. For this choice the equation (3.3) becomes

g(D)(f − L_{P,Θ}f)(x) = −g(D)(L_{P,Θ}R_{k,x}f)(x), ∀f,  (3.4)

where

R_{k,x} := 1 − T_{k,x}

is the error (remainder) in Taylor interpolation from Π_{k} at x, which can be expressed as

R_{k,x}f(v) = \frac{1}{k!} ∫_{0}^{1}(1 − t)^{k}D_{v−x}^{k+1}f(x + t(v − x)) dt.  (3.5)

Here D_{y}f denotes the derivative of f in the direction y. For a given linear functional A, such as point evaluation at x, there are (for k > 0) many linear projectors satisfying (3.2). This is illustrated in Section 6 for linear interpolation, i.e., P = Π_{1}.

The careful reader will notice that the only property of L_{P,Θ} used so far is that it is a linear operator which reproduces Π_{k}, and so, in particular, (3.4) holds with L_{P,Θ} replaced by any such operator. This is a crucial observation, which leads to Theorem 3.15 and is discussed further in this section.

By expressing the Lagrange interpolant L_{P,Θ}(R_{k,x}f) in Lagrange form (3.1) the equation (3.4) can be expanded as

g(D)(f − L_{P,Θ}f)(x) = − ∑_{v∈Θ} R_{k,x}f(v) (g(D)ℓ_{v})(x).  (3.6)

This formula is referred to by Meinguet in [M84:p98] as the Kowalewski – Ciarlet – Wagschal multipoint Taylor formula out of deference to its origins which are as follows.

History of the (Kowalewski – Ciarlet – Wagschal) multipoint Taylor formula

The earliest occurrence of (3.6) is Kowalewski’s remainder (see [K32:p21-24], or Davis’s book [D75:p71]), which is the univariate case when the linear functional A is point evaluation at x.

In Ciarlet and Wagschal [CW71] the multivariate form of (3.6) was given for the special cases of interpolation by linear polynomials (Theorem 2) and by quadratic polynomials for derivatives A = δ_{x}g(D) up to order 1 (Theorem 4). The term multipoint Taylor formula was coined in that paper. The subsequent paper of Ciarlet and Raviart [CR72] proved the general formula (3.6). A discussion of the relationship of that proof with Theorem 2.2 and of the choice of the term ‘multipoint Taylor formula’ is given below.

The special case of formula (3.6) where A is point evaluation at x and Θ consists of the simplex points (dim Π_{k} equally spaced points on a simplex) was proved independently by Nicolaides in [N72:Th.5.1] and [N73:Th.2.1].

The Ciarlet–Raviart proof of the multipoint Taylor formula

To understand the proof of (3.6) given in [CR72] it is instructive to consider the following alternative simple proof found by the author (in trying to do so), which invited the generalisation to Theorem 2.2.

Since

f = T_{k,x}f + R_{k,x}f,

and L_{P,Θ} reproduces Π_{k}, i.e.,

L_{P,Θ}g = g, ∀g ∈ Π_{k},

it follows that

L_{P,Θ}f = L_{P,Θ}T_{k,x}f + L_{P,Θ}R_{k,x}f = T_{k,x}f + L_{P,Θ}R_{k,x}f.  (3.7)

In particular, since T_{k,x}f is the Taylor interpolant from Π_{k} to f at x, applying δ_{x}g(D) to each side of (3.7) gives

g(D)L_{P,Θ}f(x) = g(D)f(x) + (g(D)L_{P,Θ}R_{k,x}f)(x),  (3.8)

which (up to a rearrangement) is (3.4) and can be expanded to (3.6). Again, the careful reader will observe that (3.7) holds with L_{P,Θ} replaced by any linear operator that reproduces Π_{k}.

The proof of (3.6) given in [CR72] uses the m–th order differentiation operator, which (for simplicity) can be thought of as the map

D^{m} : f ↦ (D^{α}f : |α| = m)

to the sequence of partial derivatives of order m. There the Taylor formula

\begin{aligned}f(v) & = T_{k,x}f(v) + R_{k,x}f(v) \\ & = f(x) + D_{v−x}f(x) + ⋯ + \frac{1}{k!}D_{v−x}^{k}f(x) + R_{k,x}f(v) \end{aligned}  (3.9)

is used to expand the Lagrange form of L_{P,Θ}, and this expansion is differentiated to obtain

D^{m}L_{P,Θ}f(x) = ∑_{l=0}^{k} \frac{1}{l!} ∑_{v∈Θ} D_{v−x}^{l}f(x) (D^{m}ℓ_{v})(x) + ∑_{v∈Θ} R_{k,x}f(v) (D^{m}ℓ_{v})(x).  (3.10)

It is then proved that the double sum in (3.10) reduces to

∑_{l=0}^{k} \frac{1}{l!} ∑_{v∈Θ} D_{v−x}^{l}f(x) (D^{m}ℓ_{v})(x) = D^{m}f(x), 0 ≤ m ≤ k,  (3.11)

thereby proving (3.6), where D^{m}, 0 ≤ m ≤ k plays the role of g(D), g ∈ Π_{k}. With hindsight, since

∑_{l=0}^{k} \frac{1}{l!}D_{v−x}^{l}f(x) = T_{k,x}f(v),

it is seen that (3.11) is an expansion into Lagrange form of

D^{m}(L_{P,Θ}T_{k,x}f)(x) = D^{m}f(x), 0 ≤ m ≤ k,

an identity which holds because

L_{P,Θ}(T_{k,x}f) = T_{k,x}f,

and not because of the specific nature of the linear operator L_{P,Θ}. Instead, the proof of (3.11) given in [CR72] involves showing (by induction) what is described as ‘the (somehow unexpected) result’ that:

\frac{1}{l!} ∑_{v∈Θ} D_{v−x}^{l}f(x) (D^{m}ℓ_{v})(x) = \{ \begin{matrix}0, & 0 ≤ l < m; \\ D^{m}f(x), & l = m; \\ 0, & m < l ≤ k. \end{matrix}  (3.12)

Gregory’s error formula for linear interpolation on a triangle

For the multivariate case, in addition to (3.5), there are many expansions for the error in Taylor interpolation R_{k,x}f(v) in terms of D^{k+1}f that could be substituted into (3.6), each giving different sums of standard linear mappings (1.3) (and thus leading to possibly better norm bounds).

Gregory’s well known formulæ for the error in linear interpolation at

Θ := {(1, 0), (0, 1), (0, 0)}

is an example of this. In Gregory [G75] the formulæ (3.6) for remainders

f(x) − L_{Π_{1},Θ}f(x),  D^{(1,0)}(f − L_{Π_{1},Θ}f)(x),  D^{(0,1)}(f − L_{Π_{1},Θ}f)(x)  (3.13)

are expanded using a formula for R_{1,x}f(v) different from (3.5). To get an appreciation for the detail involved, the formula used, which is in the spirit of Sard [S63:Th.8,p163], is that for x, v ∈ T (the triangle with vertices Θ)

R_{1,x}f(v) = χ_{A_{1}}(v)[∫_{x_{1}}^{v_{1}}(v_{1} − s)D^{(2,0)}f(s, x_{2}) ds + ∫_{x_{2}}^{v_{2}} ∫_{x_{1}}^{v_{1}} D^{(1,1)}f(s, t) ds dt
+ ∫_{x_{2}}^{v_{2}}(v_{2} − t)D^{(0,2)}f(x_{1}, t) dt]
+ χ_{A_{2}}(v)[∫_{x_{1}}^{v_{1}} ∫_{x_{1}}^{t} D^{(2,0)}f(s, 1 − t) ds dt + ∫_{x_{1}}^{v_{1}} ∫_{x_{2}}^{1−s} D^{(1,1)}f(x_{1}, t) dt ds
+ ∫_{x_{1}}^{v_{1}} ∫_{1−s}^{v_{2}} D^{(1,1)}f(s, t) dt ds + ∫_{x_{2}}^{v_{2}}(v_{2} − t)D^{(0,2)}f(x_{1}, t) dt]
+ χ_{A_{3}}(v)[∫_{x_{1}}^{v_{1}}(v_{1} − s)D^{(2,0)}f(s, x_{2}) ds + ∫_{x_{2}}^{v_{2}} ∫_{1−t}^{v_{1}} D^{(1,1)}f(s, t) ds dt

+ ∫_{x_{2}}^{v_{2}} ∫_{x_{1}}^{1−t} D^{(1,1)}f(s, x_{2}) ds dt + ∫_{x_{2}}^{v_{2}} ∫_{x_{2}}^{s} D^{(0,2)}f(1 − s, t) dt ds],

(3.14)

where χ_{A_{i}} is the characteristic function of the set A_{i} which is defined by

A_{i} := \{ \begin{matrix}T \ (A_{2} ∪ A_{3}), & i = 1; \\ \{v ∈ T : v_{1} > 1 − x_{2}\}, & i = 2; \\ \{v ∈ T : v_{2} > 1 − x_{1}\}, & i = 3. \end{matrix}

Naturally, the formulæ for the remainders (3.13) so obtained are of a similar complexity to (3.14). Gregory was successful in obtaining useful L_{p}–error bounds from these formulæ.

The (unexpanded) multipoint Taylor formula for linear operators that reproduce Π_{k}

We end this section with the general result referred to previously.

Theorem 3.15. Suppose that L is a linear operator that reproduces Π_{k}. Then, for sufficiently smooth f

Lf = T_{k,x}f + LR_{k,x}f.  (3.16)

In particular, for g ∈ Π_{k}

g(D)(f − Lf)(x) = Q(D^{k+1}f),  (3.17)

where

Q(D^{k+1}f) := −g(D)(LR_{k,x}f)(x).  (3.18)

Proof: The proof is exactly as for (3.7) and (3.8). Since T_{k,x}f ∈ Π_{k} and L reproduces Π_{k},

Lf = LT_{k,x}f + LR_{k,x}f = T_{k,x}f + LR_{k,x}f.

Apply δ_{x}g(D) to each side of this, to obtain (after rearrangement) that

g(D)(f − Lf)(x) = −g(D)(LR_{k,x}f)(x).

From (3.5) it follows that −g(D)(LR_{k,x}f)(x) depends only on D^{k+1}f, as is claimed by the definition of (3.18).

Terminology

The term multipoint Taylor formula for (3.6) was adopted in [CW71], not because it involves the error R_{k,x}f(v) for many points v ∈ Θ, but rather because (3.6) can be rewritten as

g(D)f(x) = ∑_{v∈Θ} f(v) (g(D)ℓ_{v})(x) + error(D^{k+1}f),

which expresses the derivative g(D)f(x) in terms of the value of f at the multiple points v ∈ Θ plus an error term. Since a similar phenomenon will be observed for the linear projectors L considered in the next section, equation (3.17) will be referred to as the unexpanded multipoint Taylor formula for L. The adjective unexpanded is used since, as we have just seen, the error R_{k,x}f occurring in (3.18) can be expanded in many ways (as can L itself). These expansions are explored within the setting of finite elements in the next section.

4. Error formulæ for finite elements

In this section the expanded form of the multipoint Taylor formula (3.6) is extended to all finite elements, with the examples of the cubic Hermite n-simplex and Wilson’s brick being treated in detail.

Finite elements

The following (standard) definitions can be found in Brenner and Scott [BS94:Ch.3].

Definition 4.1. Let

(i) K ⊂ IR^{n} be a domain with piecewise smooth boundary (the element domain), (ii) P be a finite-dimensional space of functions on K (the shape functions), and (iii) N = {N_{i}} be a basis for P′ (the nodal variables). Then (K, P, N) is called a finite element.

It is assumed that the nodal variables {N_{i}} are defined on sufficiently smooth functions (e.g., those from some Sobolev space). For a finite element (K, P, N), the basis for P dual to {N_{i}} is called the nodal basis for P. It will be denoted by {ϕ_{i}}.

Definition 4.2. Given a finite element (K, P, N), its local interpolant (to f) is

I_{K}f := I_{(K,P,N)}f := ∑_{i} N_{i}(f) ϕ_{i}.  (4.3)

Observe that the map I_{K} defined by (4.3) is the linear projector onto P with interpolation conditions N, i.e., for which

N_{i}(I_{K}f) = N_{i}(f), ∀i, ∀f.

A typical example of a finite element is that (with its local interpolant) given by a Lagrange map L_{P,Θ}, i.e., a Lagrange finite element (see, e.g., Ciarlet and Lions [CL91:p95]). Here, the element domain is usually the convex hull of Θ, the shape functions are P, and the nodal variables are the point evaluations {δ_{v}}_{v∈Θ}. In this case the nodal basis consists of the Lagrange functions and (4.3) is exactly the Lagrange form (3.1) of L_{P,Θ}.

Expanding the (unexpanded) multipoint Taylor formula for a finite element in terms of the nodal basis

In Section 3, when L was a Lagrange map L_{P,Θ}, equation (3.17) was expanded in terms of the Lagrange functions (the nodal basis) to obtain the Kowalewski – Ciarlet – Wagschal multipoint Taylor formula (3.6). In exactly the same way, when L = I_{K} is any finite element it is possible to expand (3.17) in terms of the nodal basis {ϕ_{i}} as follows.

Theorem 4.4. Suppose that (K, P, N) is a finite element with Π_{k} ⊂ P ⊂ C^{k}. Then, for sufficiently smooth f and g ∈ Π_{k},

g(D)(f − I_{K}f) (x) = − ∑_{i} N_{i}(R_{k,x}f) (g(D)ϕ_{i})(x),  (4.5)

where N_{i}(R_{k,x}f) can be expressed in terms of D^{k+1}f using Lemma 4.11 (see below).

Proof: Since L := I_{K} is a linear projector that reproduces Π_{k}, Theorem 3.15 gives (3.17). Expanding I_{K}(R_{k,x}f) as it occurs in (3.17) using (4.3) gives (4.5). The expansion of N_{i}(R_{k,x}f) in terms of D^{k+1}f is discussed below.

For the reasons just mentioned, equation (4.5) (and expansions of it by Lemma 4.11) will be referred to as the expanded multipoint Taylor formula for the finite element (K, P, N) (or, equivalently, for its local interpolant I_{K}).

Using variations of the Ciarlet–Raviart proof of the multipoint Taylor formula (discussed earlier) the equivalent of (4.5) has been given for certain finite elements with nodal variables of the form

f ↦ D_{ξ_{1}} ⋯ D_{ξ_{d}}f(v),  0 ≤ d ≤ k,  (4.6)

which are often called Hermite finite elements. In Ciarlet and Raviart [CR72:Th.3] it is proved for P = Π_{k} and nodal variables with 0 ≤ d ≤ 1, with the result involving nodal variables with 2 ≤ d ≤ k being described as ‘long and cumbersome’ and hence omitted. In Gout [Go77:p414] the result is shown for Π_{k} ⊂ P (without restriction on 0 ≤ d ≤ k).

Expressing N_{i}(R_{k,x}f) in terms of D^{k+1}f

The majority of finite elements used in numerical modelling have nodal variables of the form

f ↦ q(D)f(v),  (4.7)

where q ∈ Π_{d}^{0} (the space of homogeneous polynomials of degree d), 0 ≤ d ≤ k. This is a (slightly) more general form than (4.6). Indeed, the very term nodal variable (degree of freedom is also used) comes from the fact that (4.6) and (4.7) involve the node v. A nodal variable of the form (4.7) will be said to have the usual form.

The exceptions given in Ciarlet and Lions monograph [CL91] are Wilson’s brick (page 214), and the Fraeijs de Veubeke triangle (page 309). In addition to nodal variables of the usual form (4.7), these involve some of the form

f ↦ ∫_{T} q(D)f(v) dv,  (4.8)

where T is the element domain K, or one of its faces, and dv is Lebesgue measure.

By considering the error in Taylor interpolation we can now give expressions for N_{i}(R_{k,x}f) in terms of D^{k+1}f for nodal variables N_{i} of either form. These will be presented using the following instance of the divided difference functional on IR^{n} (see (5.6) for the general definition). For 0 ≤ d ≤ k, let

∫_{[^{x,...,x}_{k+1−d},v]} f := \frac{1}{(k − d)!} ∫_{0}^{1}(1 − t)^{k−d}f(x + t(v − x)) dt,  (4.9)

which allows (3.5) to be rewritten as

R_{k−d,x}f(v) = ∫_{[^{x,...,x}_{k+1−d},v]} D_{v−x}^{k+1−d}f.  (4.10)

Lemma 4.11. Suppose that λ = N_{i} is a nodal variable of the usual form, i.e.,

λ(f) := q(D)f(v),

where q ∈ Π_{d}^{0}, 0 ≤ d ≤ k. Then, for sufficiently smooth f,

λ(R_{k,x}f) = ∫_{[^{x,...,x}_{k+1−d},v]} D_{v−x}^{k+1−d}q(D)f.  (4.12)

In particular, if f is C^{k+1} in a neighbourhood of the line segment from x to v, then there exists ξ belonging to the line segment from x to v for which

λ(R_{k,x}f) = \frac{1}{(k + 1 − d)!}(D_{v−x}^{k+1−d}q(D) f)(ξ).  (4.13)

If λ = N_{i} is a nodal variable which is of the form (4.8), i.e.,

λ(f) := ∫_{T} q(D)f (v) dv,

then, for sufficiently smooth f,

λ(R_{k,x}f) = ∫_{T} (∫_{[^{x,...,x}_{k+1−d},v]} D_{v−x}^{k+1−d}q(D)f)dv.  (4.14)

Proof: The proof relies upon the fact that

q(D)(T_{k,x}f) = T_{k−d,x}(q(D)f),  (7.5)

which is a special case of a more general result for the scale of mean value interpolations (see Waldron [W97_{1}:Prop.5.11]). Here is a direct proof. If f ∈ Π_{k}, then q(D)f ∈ Π_{k−d} and it is obvious that (7.5) holds. But both sides of (7.5) depend only on the derivatives up to order k of f at x, and so it holds for all (smooth) f.

From (7.5) it follows that

q(D)(R_{k,x}f) = R_{k−d,x}(q(D)f),

and hence, by (4.10),

λ(R_{k,x}f) = q(D)(R_{k,x}f)(v) = R_{k−d,x}(q(D)f) = ∫_{[^{x,...,x}_{k+1−d},v]} D_{v−x}^{k+1−d}q(D)f,

giving (4.12). Since

∫_{[^{x,...,x}_{k+1−d},v]} 1 = \frac{1}{(k + 1 − d)!},

the mean value theorem (for integrals) can be applied to (4.12) to obtain (4.13). Lastly, integrating (4.12) over v ∈ T gives (4.14).

Now we illustrate the use of Theorem 4.4 with three examples.

Example 1: The cubic Hermite n-simplex

Here we give the expanded multipoint Taylor error formula for the cubic Hermite nsimplex. This finite element was introduced (under the name approximation of type 3) in Ciarlet and Wagschal’s paper [CW71]. It is defined as follows.

Let K be an n-simplex in IR^{n} with vertices {a_{i}}_{i=1}^{n+1}, and let

a_{ijk} := (a_{i} + a_{j} + a_{k})/3, 1 ≤ i < j < k ≤ n + 1.

Let N be the set of linear functionals consisting of point evaluation at a_{i}, 1 ≤ i ≤ n + 1, and a_{ijk}, 1 ≤ i < j < k ≤ n + 1, together with

N_{ij} : f ↦ D_{a_{j}−a_{i}}f(a_{i}), 1 ≤ i, j ≤ n + 1, i ≠ j.

Then (K, Π_{3}, N) is a finite element called the cubic Hermite n-simplex.

a_{123} a_{1}  a_{2}  a_{3}

Fig. 4.1. The cubic Hermite triangle (dim Π_{3}(IR^{2}) = 10).

Its local interpolant is given by

Hf := I_{K}f = ∑_{i} f(a_{i}) ϕ_{i} + ∑_{i<j<k} f(a_{ijk}) ϕ_{ijk} + ∑_{i≠j} D_{a_{j}−a_{i}}f(a_{i}) ϕ_{ij},

where the nodal basis has the explicit representation

ϕ_{i} := −2λ_{i}^{3} + 3λ_{i}^{2} − 7λ_{i} ∑_{{j<k \atop j,k≠i}} λ_{j}λ_{k},

ϕ_{ijk} := 27λ_{i}λ_{j}λ_{k},

ϕ_{ij} := λ_{i}λ_{j}(2λ_{i} + λ_{j} − 1), with {λ_{i}} the barycentric coordinates corresponding to the points {a_{i}}.

Taking q(D) = D_{a_{j}−a_{i}} and v = a_{i} in Lemma 4.11 gives the expansion

N_{ij}(R_{3,x}f) = ∫_{[x,x,x,a_{i}]} D_{a_{i}−x}^{3}D_{a_{j}−a_{i}}f.

Thus, applying Theorem 4.4 to I_{K} = H (with k = 3) gives the error formula: for g ∈ Π_{3}

\begin{aligned}g(D)(f − Hf) (x) = & − ∑_{i}(∫_{[x,x,x,x,a_{i}]} D_{a_{i}−x}^{4}f)(g(D)ϕ_{i})(x) \\ & − ∑_{i<j<k}(∫_{[x,x,x,x,a_{ijk}]} D_{a_{ijk}−x}^{4}f)(g(D)ϕ_{ijk})(x) \\ & − ∑_{i≠j}(∫_{[x,x,x,a_{i}]} D_{a_{i}−x}^{3}D_{a_{j}−a_{i}}f)(g(D)ϕ_{ij})(x), \end{aligned}  (4.16)

which holds for all sufficiently smooth f.

In the next section it is shown how to obtain L_{p}-error bounds from this and other such formulæ.

Example 2: Wilson’s brick

The second example involves nodal variables which are not of the usual form. Let K be a rectangle in IR^{3} aligned in the coordinate directions e_{1}, e_{2}, e_{3}, with its length in these directions being 2h_{1}, 2h_{2}, 2h_{3} respectively, and its (eight) vertices being v ∈ Θ. Let

P := Π_{2} ⊕ span{( )^{1,1,1}},

where ( )^{1,1,1} : (x_{1}, x_{2}, x_{3}) ↦ x_{1}x_{2}x_{3}. Then, Wilson’s brick (see [WTDG73]) is the finite element (K, P, N), where the nodal variables N consist of the eight point evaluations {δ_{v} : v ∈ Θ} together with

N_{j} : f ↦ \frac{h_{j}^{2}}{h_{1}h_{2}h_{3}} ∫_{K} D_{j}^{2}f, 1 ≤ j ≤ 3.

Here D_{j}^{2}f denotes D_{e_{j}}^{2}f, the second partial derivative of f in the coordinate direction e_{j}.

2h_{1} 2h^{2} 2h_{3}  a_{1}  a_{2} a_{3} a_{4}  a_{5}  a_{6} a_{7} a_{8}

Fig. 4.2. Wilson’s brick, n = 3 (dim P = 11).

The local interpolant for Wilson’s brick is

Wf := I_{K}f := ∑_{v∈Θ} f(v) ϕ_{v} + ∑_{j}(\frac{h_{j}^{2}}{h_{1}h_{2}h_{3}} ∫_{K} D_{j}^{2}f) ϕ_{j}.

Formulæ for the nodal basis {ϕ_{v} : v ∈ Θ} ∪ {ϕ_{j} : 1 ≤ j ≤ 3} can be found in [CL91:p214].

Applying Theorem 4.4 to I_{K} = W (using (4.14) to express N_{j}(R_{2,x}f) in terms of D^{3}f) gives the following error formula: for g ∈ Π_{2}

\begin{aligned}g(D)(f − Wf) (x) & = − ∑_{v∈Θ}(∫_{[x,x,x,v]} D_{v−x}^{3}f) (g(D)ϕ_{v})(x) \\ & − ∑_{j} \frac{h_{j}^{2}}{h_{1}h_{2}h_{3}} ∫_{K}(∫_{[x,v]} D_{v−x}D_{j}^{2}f)dv (g(D)ϕ_{j})(x), \end{aligned}  (4.17)

for all sufficiently smooth f.

Example 3: Numerical differentiation and integration rules

Interpolatory (numerical) differentiation and integration rules based on point evaluations at v ∈ Θ and exact for Π_{k} are of the form

∑_{v∈Θ} w(v) f(v) = AL_{P,Θ}f,

where A is the derivative or integral to be approximated, L_{P,Θ} is a Lagrange map, and the weights w(v) are given by

w(v) = A(ℓ_{v}), v ∈ Θ.

The error in approximation by such a rule is

A(1 − L_{P,Θ})f,  (4.18)

which is of the form (2.1). For a differentiation rule, i.e.,

Af := g(D)f(x), g ∈ Π_{k},

this is the pointwise error in a Lagrange map, for which there is the Kowalewski – Ciarlet – Wagschal multipoint Taylor formula (3.6). For differentiation rules based on data that includes linear functionals other than point evaluations the expanded multipoint Taylor formula for a finite element (4.5) applies.

For an integration rule, i.e.,

Af := ∫_{Ω} f,

where Ω is some region of volume vol(Ω), it is also possible to apply Theorem 2.2. Let P be some linear projector onto Π_{k} satisfying

∫_{Ω} f = ∫_{Ω} Pf,  (4.19)

which is condition (2.9), and denote its error by E := 1 − P. Then, Corollary 2.7 gives the error formula

\begin{aligned}∫_{Ω} f − ∑_{v∈Θ} w(v) f(v) & = −AL_{P,Θ}Ef \\ & = −A(∑_{v∈Θ} E(v)ℓ_{v}) \\ & = − ∑_{v∈Θ} w(v)(Ef)(v). \end{aligned}  (4.20)

The hope is then to express the quantities (Ef)(v), v ∈ Θ occuring above in terms of D^{k+1}f. The practical difficulty with this scheme is finding (appropriate) projectors P which satisfy (4.19). At first glance it might seem the averaged Taylor interpolant

Pf := \frac{1}{vol(Ω)} ∫_{Ω} T_{k,x}f dx

would be a good choice, however this does not satisfy (4.19). The author has not pursued these questions further. There are a number of linear projectors P matching certain integrals, together with error representations Ef = Q(D^{k+1}f), the best known being the scale of mean value interpolations (see, e.g., [W97_{1}]), and so it would seem that at least in some situations such a theory might be feasible.

5. L_{p}–error bounds for finite elements

In this section, we use the multivariate Hardy’s inequality (recently introduced by the author) to obtain L_{p}–error bounds from the expanded multipoint Taylor formula for a finite element. To do this, the following facts about Sobolev spaces will be required.

Sobolev spaces

Let W_{p}^{k}(Ω) be the Sobolev space of functions defined on Ω (a bounded open set in IR^{n} with a Lipschitz boundary) with derivatives up to order k in L_{p}(Ω) equipped with the usual topology (see, e.g., Adams [Ad75]). It is convenient to include the condition that Ω have a Lipschitz boundary in the definition, so that Sobolev’s embedding theorem(s) can be applied. We will only need the following consequence of Sobolev’s embedding theorem(s) (which can be found in any text on Sobolev spaces). If k + 1 − d − n/p > 0, then

W_{p}^{k+1}(Ω) ⊂ C^{d}(Ω̄),  (5.1)

where Ω̄ is the closure of Ω.

To measure the size of its k–th order derivative, it is convenient to associate with each f ∈ W_{p}^{k}(Ω) the function |D^{k}f| ∈ L_{p}(Ω), given by the rule

|D^{k}f|(x) := sup_{{ξ_{1},...,ξ_{k}∈IR^{n} \atop ∥ξ_{i}∥≤1}} |D_{ξ_{1}} ⋯ D_{ξ_{k}}f(x)| = sup_{{u∈IR^{n} \atop ∥u∥=1}} |D_{u}^{k}f(x)|,  (5.2)

where the derivatives D_{ξ_{1}} ⋯ D_{ξ_{k}}f are computed from any (fixed) choice of representatives for the k–th order partial derivatives of f. The equality of the two suprema above follows from a classical result of Banach on the norm of a symmetric multilinear map. This definition of |D^{k}f| is consistent with its standard univariate interpretation. From (5.2), it is easy to see that |D^{k}f| is well–defined and satisfies

|D_{ξ_{1}} ⋯ D_{ξ_{k}}f| ≤ |D^{k}f| ∥ξ_{1}∥ ⋯ ∥ξ_{k}∥,  a.e.,  (5.3)

for ξ_{1}, …, ξ_{k} ∈ IR^{n}. The L_{p}(Ω)–norm of |D^{k}f| gives a seminorm on W_{p}^{k}(Ω)

f ↦ f _{k,p,Ω} := ∥ |D^{k}f| ∥_{L_{p}(Ω)}.  (5.4)

This coordinate–independent seminorm (5.4) is ideal for the analysis that follows because of (5.3). It is convenient to generalise (5.3) slightly. For q ∈ Π_{d}^{0}, let

∥q(D)∥ := max_{f∈Π_{d}^{0}\0} \frac{|q(D)f(0)|}{|D^{d}f|(0)},

which defines norm on the differential operators {q(D) : q ∈ Π_{d}^{0}}. It can easily be shown that: for f ∈ W_{p}^{k}(Ω)

|D_{ξ_{1}} ⋯ D_{ξ_{k−d}}q(D)f| ≤ ∥ξ_{1}∥ ⋯ ∥ξ_{k−d}∥ ∥q(D)∥ |D^{k}f|.  (5.5)

The multivariate Hardy’s inequality

The multivariate Hardy’s inequality involves the following linear functional called the divided difference functional on IR^{n} by Micchelli in [Mi80].

Definition. For Θ := [θ_{0}, …, θ_{k}] any sequence of k + 1 points in IR^{n}, let

f ↦ ∫_{Θ} f := ∫_{0}^{1} ∫_{0}^{s_{1}} ... ∫_{0}^{s_{k−1}} f(θ_{0} + s_{1}(θ_{1}−θ_{0}) + ⋯ + s_{k}(θ_{k}−θ_{k−1})) ds_{k} ⋯ ds_{2} ds_{1}, (5.6)

with the convention that ∫_{[ ]} f := 0.

This linear functional was introduced (with a change of variables) in (4.9), for the special case Θ = [x, …, x, v]. Let

(a)_{n} := (a)(a + 1)(a + 2) ⋯ (a + n − 1),

be the shifted factorial function, and #Θ denote the cardinality of Θ.

Multivariate Hardy’s inequality 5.7 ([W97_{2}:Th.3.2.1]). Let Θ be a finite sequence in IR^{n}, and let Ω be an open set in IR^{n} for which Ω̄ is starshaped with respect to Θ. If m − n/p > 0, then the rule

H_{m,Θ}f(x) := ∫_{[^{x,...,x}_{m},Θ]} f  (5.8)

induces a positive bounded linear map H_{m,Θ} : L_{p}(Ω) → L_{p}(Ω) with norm

∥H_{m,Θ}∥_{L_{p}(Ω)} ≤ \frac{1}{(m − 1)!(m − n/p)_{#Θ}} → ∞ as m − n/p → 0^{+}.  (5.9)

This upper bound for ∥H_{m,Θ}∥_{L_{p}(Ω)} is sharp when Θ involves only one point, i.e., when

Θ = [v, …, v],

and when p = ∞ (with the norm taken on only for the constant functions).

When n = 1, m = 1, Θ = [0] and Ω = (0, ∞) the inequality (5.9) is the well–known Hardy’s inequality: for p > 1

∥ x ↦ \frac{1}{x} ∫_{0}^{x} f ∥_{L_{p}(0,∞)} ≤ \frac{p}{p − 1}∥f∥_{L_{p}(0,∞)}, ∀f ∈ L_{p}(0, ∞),

and so (5.9) is referred to as the multivariate (form of) Hardy’s inequality.

L_{p}–error bounds from the multivariate Hardy’s inequality

Next, L_{p}–error bounds are obtained from the expanded multipoint Taylor formula for a finite element using the multivariate Hardy’s inequality. Let Ω be a bounded open subset of IR^{n} (with Lipschitz boundary), and

h_{v,Ω} := sup_{x∈Ω} ∥v − x∥ ≤ h := diam Ω (the diameter of Ω).

Lemma 5.10. Suppose that Ω̄ is starshaped with respect to v, and q ∈ Π_{d}^{0}, 0 ≤ d ≤ k. If k + 1 − d − n/p > 0 and f ∈ W_{p}^{k+1}(Ω), then

∥ x ↦ (q(D)R_{k,x}f)(v) ∥_{L_{p}(Ω)} ≤ \frac{(h_{v,Ω})^{k+1−d}∥q(D)∥}{(k − d)!(k + 1 − d − n/p)} f _{k+1,p,Ω}.  (5.11)

Proof: Since k +1−d−n/p > 0, it follows from Lemma 4.11 and the multivariate Hardy’s inequality 5.7 that

Sf : x ↦ (q(D)R_{k,x}f)(v)  (5.12)

defines a function in L_{p}(Ω). From (5.5) we obtain that

\begin{aligned}|D_{v−x}^{k+1−d}q(D)f| & ≤ ∥x − v∥^{k+1−d}∥q(D)∥ |D^{k+1}f| \\ & ≤ (h_{v,Ω})^{k+1−d}∥q(D)∥ |D^{k+1}f|, \end{aligned}

in L_{p}(Ω). Thus, by Lemma 4.11 and the multivariate Hardy’s inequality 5.7, we obtain

∥Sf∥_{L_{p}(Ω)} ≤ (h_{v,Ω})^{k+1−d}∥q(D)∥ ∥ x ↦ ∫_{[^{x,...,x}_{k+1−d},v]} |D^{k+1}f| ∥_{L_{p}(Ω)}

≤ (h_{v,Ω})^{k+1−d}∥q(D)∥\frac{1}{(k − d)!(k + 1 − d − n/p)}∥ |D^{k+1}f| ∥_{L_{p}(Ω)}

= \frac{(h_{v,Ω})^{k+1−d}∥q(D)∥}{(k − d)!(k + 1 − d − n/p)} f _{k+1,p,Ω},

which is (5.11).

We now extend the L_{p}–error bounds of Arcangeli and Gout [AG76] and Gout [Go77] to a general finite element and for a much wider class of seminorms, defined as follows. A family of seminorms | · |_{p,Ω} on W_{p}^{k}(Ω), 1 ≤ p ≤ ∞ is said to be of the usual form if either

|f|_{p,Ω} := ∥ (∥g_{j}(D)f∥_{L_{p}(Ω)} : j = 1, … m) ∥_{IR}^{m},  (5.13)

where the g_{j} ∈ Π_{k}(IR^{n}) are fixed and ∥ · ∥_{IR}^{m} is any norm on IR^{m}, or

| · |_{p,Ω} := · _{j,p,Ω},  0 ≤ j ≤ k.  (5.14)

Theorem 5.15. Suppose that (K, P, N) is a finite element with Π_{k} ⊂ P ⊂ C^{k}(Ω̄), and nodal variables {N_{i}} of the usual form

N_{i}(f) := q_{i}(D)f(v_{i}),

where q_{i} ∈ Π_{d_{i}}^{0}, 0 ≤ d_{i} ≤ k, and v_{i} ∈ IR^{n}. Let d := max_{i} d_{i}. If k + 1 − d − n/p > 0, then

|f − I_{K}f|_{p,Ω} ≤ ∑_{i} \frac{(h_{v_{i},Ω})^{k+1−d_{i}}∥q_{i}(D)∥ |ϕ_{i}|_{∞,Ω}}{(k − d_{i})!(k + 1 − d_{i} − n/p)} f _{k+1,p,Ω},  ∀f ∈ W_{p}^{k+1}(Ω), (5.16)

where | · |_{p,Ω} is any family of seminorms on W_{p}^{k}(Ω) of the usual form.

Proof: Let g_{j} ∈ Π_{k}. Since k + 1 − d − n/p > 0, Sobolev’s embedding theorem (5.1) implies that I_{K} is well defined on W_{p}^{k+1}(Ω). It follows from the multivariate Hardy’s inequality 5.7 and Theorem 4.4 that, for f ∈ W_{p}^{k+1}(Ω), the error satisfies

g_{j}(D)(f − I_{K}f) (x) = − ∑_{i}(q_{i}(D)R_{k,x}f)(v_{i}) (g_{j}(D)ϕ_{i})(x)

in L_{p}(Ω), giving

|g_{j}(D)(f − I_{K}f) (x)| ≤ ∑_{i} |(q_{i}(D)R_{k,x}f)(v_{i})| ∥g_{j}(D)ϕ_{i}∥_{L_{∞}(Ω)}.  (5.17)

Since ∥·∥_{L_{p}(Ω)} is a monotone norm, applying Lemma 5.10 to this (as a function of x) gives

∥g_{j}(D)(f − I_{K}f)∥_{L_{p}(Ω)} ≤ ∑_{i} \frac{(h_{v_{i},Ω})^{k+1−d_{i}}∥q_{i}(D)∥}{(k − d_{i})!(k + 1 − d_{i} − n/p)}∥g_{j}(D)ϕ_{i}∥_{L_{∞}(Ω)} f _{k+1,p,Ω}.

(5.18)

Taking the ∥ · ∥_{IR}^{m} norm of the inequality for m–vectors (indexed by j) which is given coordinatewise by (5.18) gives (5.16) for a family of seminorms of the form (5.13). Take the supremum of each side of (5.17) over the set of differential operators {D_{u}^{j} : ∥u∥ = 1} to obtain

|D^{j}(f − I_{K}f)|(x) ≤ ∑_{i} |(q_{i}(D)R_{k,x}f)(v_{i})| ϕ_{i} _{j,∞,Ω}, a.e. x.

Next, take the L_{p}(Ω)–norm of each side of this and apply Lemma 5.10 to get

f − I_{K}f _{j,p,Ω} ≤ ∑_{i} \frac{(h_{v_{i},Ω})^{k+1−d_{i}}∥q_{i}(D)∥ ϕ_{i} _{j,∞,Ω}}{(k − d_{i})!(k + 1 − d_{i} − n/p)} f _{k+1,p,Ω},

which is (5.16) for seminorms of the form (5.14).

Remark. If (K, P, N) involves nodal variables of the form (4.8), i.e.,

λ : f ↦ ∫_{T} q(D)f dv,

then the argument of Theorem 5.15 can be modified by using the integral form of Minkowski’s inequality to bring ∥ · ∥_{L_{p}(Ω)} inside the integral ∫_{T} (·) dv. This leads to estimates of the form

| x ↦ λ(R_{k,x}f)ϕ(x) |_{p,Ω} ≤ ∑_{i} \frac{h^{k+1−d}∥q(D)∥ vol(T)|ϕ|_{∞,Ω}}{(k − d)!(k + 1 − d − n/p)} f _{k+1,p,Ω}  (5.19)

in (5.16).

Now we illustrate Theorem 5.15 with some examples where for simplicity h_{v_{i},Ω} is replaced by the (possibly larger) constant h.

If I_{K} is a Lagrange map L_{P,Θ}, then (5.16) reduces to

|f − L_{P,Θ}f|_{p,Ω} ≤ \frac{1}{k!(k + 1 − n/p)} (∑_{v∈Θ} |ℓ_{v}|_{∞,Ω}) h^{k+1} f _{k+1,p,Ω},  (5.20)

which was proved by Arcangeli and Gout in [AG76:Th.1–1] when | · |_{p,Ω} is a seminorm of the form (5.14).

Following the work of [AG76], Gout [Go77] obtained Theorem 5.15 for finite elements with nodal variables of the form (4.6) and seminorms of the form (5.14). Both of these works are based on (the equivalent of) the multivariate Hardy’s inequality as used in the proof of Theorem 5.15.

For the cubic Hermite n–simplex (5.16) reduces to: for p > n/3 and f ∈ W_{p}^{4}(Ω)

|f − Hf|_{p,Ω} ≤ (∑_{i} \frac{|ϕ_{i}|_{∞,Ω}}{3!(4 − n/p)} + ∑_{i<j<k} \frac{|ϕ_{ijk}|_{∞,Ω}}{3!(4 − n/p)} + ∑_{i≠j} \frac{|ϕ_{ij}|_{∞,Ω}}{2!(3 − n/p)}) h^{4} f _{4,p,Ω},

where | · |_{p,Ω} is any of the usual seminorms on W_{p}^{3}(Ω). While for Wilsons’s brick, using (5.19), we obtain: for p > n and f ∈ W_{p}^{3}(Ω)

|f − Wf|_{p,Ω} ≤ (∑_{v∈Θ} \frac{|ϕ_{v}|_{∞,Ω}}{2!(3 − n/p)} + ∑_{j} \frac{|ϕ_{j}|_{∞,Ω}}{0!(1 − n/p)}) h^{3} f _{3,p,Ω},

where | · |_{p,Ω} is any of the usual seminorms on W_{p}^{2}(Ω).

6. Error formulæ for linear interpolation

In this section, Theorem 2.2 is used to construct a family of formulæ for pointwise error in linear interpolation. By linear interpolation we mean interpolation from the space of linear polynomials Π_{1} at Θ a set of n + 1 (affinely independent) points in IR^{n}, i.e., a Lagrange map

L_{Θ} := L_{Π_{1},Θ}.

If V is a linear right inverse for

U := D^{2},

and the corresponding linear projector P := 1 − VU onto Π_{1}(= ker D^{2}) satisfies

Pf(x) = f(x),  (6.1)

then by Theorem 2.2

f(x) − L_{Θ}f(x) = −(L_{Θ}V(D^{2}f))(x) = − ∑_{v∈Θ} V(D^{2}f)(v) ℓ_{v}(x).  (6.2)

Earlier, the choice P = T_{1,x} was made to obtain from (6.2) that

f(x) − L_{Θ}f(x) = − ∑_{v∈Θ} (∫_{[x,x,v]} D_{v−x}^{2}f) ℓ_{v}(x),  (6.3)

which is the expanded multipoint Taylor formula (4.5) for g = 1 (originally given in [CW71]).

In addition to P = T_{1,x}, there are many other linear projectors onto Π_{1} that satisfy (6.1). When a representation for V the corresponding right inverse of D^{2} is available, each of these leads to a different error formula (6.2). For example, one might take P to be Kergin interpolation at {x, ξ} with ξ ∈ IR^{n} (see, e.g., Waldron [W97_{1}]), for which

V(D^{2}f) = ∫_{[·,x,ξ]} D_{·−x}D_{·−ξ}f.  (6.4)

When ξ = x this is Taylor interpolation at x (compare (6.4) with (4.10) for k − d = 1). More generally, one can take a weighted average over ξ of such approximations to obtain the following family of formulæ of the type (6.2).

Theorem 6.5. Let μ be a real measure on IR^{n} of mass 1 (and finite variation). Then, the error in linear interpolation to Θ satisfies

f(x) − L_{Θ}f(x) = − ∑_{v∈Θ} ∫_{IR_{n}} (∫_{[v,x,ξ]} D_{v−x}D_{v−ξ}f) dμ(ξ) ℓ_{v}(x),  (6.6)

for sufficiently smooth f.

Proof: The convergence of the integrals of the operators used below (which can be interpreted pointwise) follows from the standard theory of integration.

Let K_{{x,ξ}} denote the operator of Kergin interpolation to {x, ξ} which has the representation

K_{{x,ξ}}f = f(x) + ∫_{[x,ξ]} D_{·−x}f.  (6.7)

Integrating the right hand side of (6.7) against the measure μ shows that

Pf := ∫ K_{{x,ξ}}f dμ(ξ)

defines a linear projector onto Π_{1} which satisfies (6.1). Similarly, the right inverses corresponding to K_{{x,ξ}} given by (6.4) can be integrated against μ to get

V(D^{2}f) = ∫_{IR_{n}} ∫_{[·,x,ξ]} D_{·−x}D_{·−ξ}f dμ(ξ),  (6.8)

the right inverse corresponding to P. Substitute (6.8) in (6.2) to obtain (6.6).

Certain choices of μ in Theorem 6.5 as a measure supported on a finite number of points give error formulæ that are known. Here are these examples which were obtained (some very recently) by a variety of methods. Example 1: Suppose that

μ = δ_{x},

i.e., point evaluation at x, then (6.6) gives the ‘classical’ multipoint Taylor formula (6.3). Example 2: Suppose that

μ = δ_{w},

i.e., point evaluation at w ∈ Θ, then (6.6) gives

f(x) − L_{Θ}f(x) = ∑_{{v∈Θ \atop v≠w}} (∫_{[x,w,v]} D_{x−v}D_{v−w}f) ℓ_{v}(x),  (6.9)

which is a special case of the Sauer–Xu formula [SX95:Cor.3.11] (also see de Boor [B96]). Example 3: Suppose that

μ = ∑_{w∈Θ} ℓ_{w}(x)δ_{w}

a weighted average of the point evaluations at w ∈ Θ, then (6.6) gives

f(x) − L_{Θ}f(x) = ∑_{v∈Θ} ∑_{w∈Θ} (∫_{[x,w,v]} D_{x−v}D_{v−w}f) ℓ_{v}(x)ℓ_{w}(x).  (6.10)

The ‘symmetrised’ form of (6.10) is

f(x) − L_{Θ}f(x) = ∑_{{\{v,w\}⊂Θ \atop v≠w}} (∫_{[x,w,v]} D_{w−v}D_{v−w}f) ℓ_{v}(x)ℓ_{w}(x),  (6.11)

which was given by Waldron in [W98]. Example 4: Suppose that

μ = \frac{1}{n + 1} ∑_{w∈Θ} δ_{w}

the average of the point evaluations at w ∈ Θ, then (6.6) gives

f(x) − L_{Θ}f(x) = \frac{1}{n + 1} ∑_{v∈Θ} ∑_{w∈Θ} (∫_{[x,v,w]} D_{v−x}D_{w−v}f) ℓ_{v}(x),  ()

which has the ‘symmetrised’ form

f(x) − L_{Θ}f(x) = \frac{1}{n + 1} ∑_{{\{v,w\}⊂Θ \atop v≠w}} ∫_{[x,w,v]} ℓ_{v}(x)D_{v−x}D_{w−v}f + ℓ_{w}(x)D_{w−x}D_{v−w}f. (6.13)

Both of these formulæ are new.

It is easy to see from these examples how many new formulæ for the error in linear interpolation can be obtained by choosing different μ in Theorem 6.5. This section was only intended to give an indication of how effectively Theorem 2.2 can be applied to obtain error formulæ. The generalisation of Theorem 6.5 to the error in interpolation from Π_{k}, where μ becomes a real measure on (IR^{n})^{k} of mass 1, is left to the reader. In the last section another promising application of Theorem 2.2 is outlined.

7. Error formulæ for Lagrange maps from those for one–point interpolation

In this section, we give a generalisation of Theorem 3.15 where Π_{k} is replaced by a polynomial space Q. It is shown that the problem of constructing an error formula for Lagrange interpolation from a D–invariant polynomial space P (such as a tensor product space, or more generally the least solution) that involves only derivatives which annihilate P can be reduced to the simpler problem of obtaining such a formula for ‘one–point’ (Taylor) interpolation from P at the origin (which has been solved in a number of cases). These results require the following facts about one–point interpolation.

One–point interpolation

We say that one–point interpolation from a (finite dimensional) space of polynomials Q to derivatives Q(D) at a point x ∈ IR^{n} is correct if for each function f (smooth enough at x) there is a unique q ∈ Q with

p(D)q(x) = p(D)f(x),  ∀p ∈ Q.  (7.1)

The associated linear projector f ↦ q is called the one–point interpolation map (from Q at x). It is denoted by

T_{Q,x} : f ↦ q,  (7.2)

and its error (remainder) by

R_{Q,x} : f ↦ f − q.  (7.3)

For Q = Π_{k} this one–point interpolation map is simply the Taylor map from Π_{k} at x, which was discussed in Section 3. There the (shorter) notations

T_{k,x} = T_{Π_{k},x}, R_{k,x} = R_{Π_{k},x}

were used.

One–point interpolation at the origin, which will be denoted by

T_{Q} := T_{Q,0} (with R_{Q} := R_{Q,0}),

is always correct (even if Q is not a homogeneous space of polynomials). To see this, and other properties of the set

OPI(Q) := {x ∈ IR^{n} : One–point interpolation from Q at x is correct},

we equip Π (the space of polynomials) with the inner product

⟨f, p⟩ := p(D)f(0),  (7.4)

which was used so effectively in de Boor and Ron [BR92] to study multivariate polynomial interpolation. With the polynomials thought of as an incomplete inner product space in this way we have the following.

Proposition 7.5. One–point interpolation from Q at x ∈ IR^{n} is correct for a.e. x. If Q is D–invariant (invariant under differentiation), then one–point interpolation is correct for all x ∈ IR^{n}, i.e.,

OPI(Q) = IR^{n},

and

T_{Q,x} = τ_{−x}T_{Q}τ_{x},  (7.6)

where τ_{x} is the translation operator

τ_{x} : f ↦ f(· + x).

Furthermore, one–point interpolation at the origin is always correct, and its restriction

T_{Q}|_{Π} : Π → Q

to the polynomials is the orthogonal projection onto Q.

Proof: The orthogonal projection of f ∈ Π onto Q is, by definition, the unique q ∈ Q for which

⟨f − q, Q⟩ = 0.  (7.7)

By (7.4) it is seen that (7.7) is precisely (7.1), i.e., the interpolation conditions for one–point interpolation from Q at the origin, which gives the last statement.

Let {q_{i}} be a basis for Q. One–point interpolation from Q at x is correct if and only if the Gramian matrix

G(x) := [q_{i}(D)q_{j}(x)]_{i,j}  (7.8)

is invertible. Thus, one–point interpolation from Q at x is not possible only when x is a zero of the polynomial

p : x ↦ det(G(x)).

But, we have shown that one–point interpolation at the origin is correct, so p(0) ≠ 0, and p is a nonzero polynomial, which therefore has zero set

Z(p) = IR^{n} \ OPI(Q)

of measure zero. In other words, one–point interpolation from Q is correct at a.e. x.

Suppose that Q is D–invariant. Then, Q is translation–invariant (see de Boor [B87]), and so

τ_{−x}T_{Q}τ_{x}  (7.9)

is a linear projector onto Q. Next, we show that (7.9) interpolates derivatives Q(D) at x, so that one–point interpolation from Q at x is correct and given by (7.6). This argument uses the facts that differentiations and translations commute, and

δ_{x}τ_{y} = δ_{x+y}.

Let p ∈ Q, then

δ_{x}p(D)(τ_{−x}T_{Q}τ_{x}f) = δ_{x}τ_{−x}p(D)T_{Q}τ_{x}f
= δ_{0}p(D)T_{Q}τ_{x}f
= δ_{0}p(D)τ_{x}f (T_{Q} interpolates δ_{0}p(D))
= δ_{0}τ_{x}p(D)f
= δ_{x}p(D)f,

so that τ_{−x}T_{Q}τ_{x} matches p(D)f(x).

The ‘almost everywhere’ in Proposition 7.5 is needed. Suppose that Q is the space spanned by the (univariate) polynomial

q_{1} : x ↦ 1 + ax, 0 ≠ a ∈ IR .

Then, the Gramian matrix (7.8) is

G(x) = [1 + a^{2} + ax],

and so one–point interpolation from Q to x ∈ IR is not correct when

x = −\frac{1 + a^{2}}{a} ,

and is correct for all other values of x.

More general multipoint Taylor formulæ for Lagrange maps and other linear projectors

We are now able to give a generalisation of Theorem 3.15 where Π_{k} is replaced by a (finite–dimensional) polynomial space Q.

Theorem 7.10. Suppose that L is a linear operator that reproduces a polynomial space Q. Then, for x ∈ OPI and sufficiently smooth f

Lf = T_{Q,x}f + LR_{Q,x}f.  (7.11)

In particular, for g ∈ Q and x ∈ OPI(Q)

g(D)(f − Lf)(x) = −g(G)(LR_{Q,x}f)(x).  (7.12)

If L is a Lagrange map L_{P,Θ}, then (7.12) can be expanded as

g(D)(f − L_{P,Θ}f)(x) = − ∑_{v∈Θ} R_{Q,x}f(v) (g(D)ℓ_{v})(x).  (7.13)

Proof: The proof is exactly as for Theorem 3.15 with Π_{k} replaced by Q and using the fact that

g(D)(T_{Q,x}f) = g(D)f(x), ∀x ∈ OPI(Q).

Proposition 7.5 gives the following useful corollary for Lagrange interpolation from a D–invariant space.

Corollary 7.14. Suppose that P is a (finite–dimensional) D–invariant space (necessarily a polynomial space) for which Lagrange interpolation at Θ ⊂ IR^{n} is correct, and L_{P,Θ} is the corresponding Lagrange map. Then, for g ∈ P and x ∈ IR^{n}

g(D)(f − L_{P,Θ}f)(x) = − ∑_{v∈Θ}(τ_{−x}R_{P} τ_{x}f)(v) (g(D)ℓ_{v})(x),  (7.15)

for sufficiently smooth f.

Proof: It is shown in de Boor [B87] that a finite–dimensional, D–invariant space is necessarily a polynomial space. By (7.6),

R_{P,x} = 1 − T_{P,x} = τ_{−x}(1 − T_{P})τ_{x} = τ_{−x}R_{P} τ_{x},

which is substituted into (7.13) to get (7.15), which holds for all x ∈ IR^{n} = OPI(P).

Thus, when P is D–invariant, to obtain an error formula for L_{P,Θ} that involves only derivatives that kill P, it is sufficient to do so for the ‘simpler’ problem of one–point interpolation at the origin. A particular D–invariant choice of great interest is the least solution of de Boor and Ron [BR92], which associates with any (finite) set of points Θ ⊂ IR^{n} a polynomial space, denoted by

P = Π_{Θ},

for which Lagrange interpolation at Θ is correct. In addition to being D–invariant the least space Π_{Θ} has many other desirable properties including constructability and being a minimal–degree solution, i.e.,

dim(Π_{Θ} ∩ Π_{k}) ≥ dim(P ∩ Π_{k}), ∀k,

for all polynomial spaces P for which Lagrange interpolation at Θ is correct.

References

[Ad75] Adams, R. A. (1975): Sobolev spaces. Academic Press, New York

[At92] Attéia, M. (1992): Hilbertian Kernels and Spline Functions. Elsevier Science Publishers, Amsterdam

[AG76] Arcangeli, R. and J. L. Gout (1976): Sur l’evaluation de l’erreur d’interpolation de Lagrange dans un ouvert de IR^{n}. Rev. Française Automat. Informat. Rech. Opér., Anal. Numer. 10(3), 5–27

[B87] Boor, C. de (1987): The polynomials in the linear span of integer translates of a compactly supported function. Constr. Approx. 3, 199-208

[B96] Boor, C. de (1996): On the Sauer-Xu formula for the error in multivariate polynomial interpolation. Math. Comp. 65, 1231–1234

[BR92] Boor, C. de and A. Ron (1992): The least solution for the polynomial interpolation problem. Math. Z. 210, 347–378

[BS94] Brenner, S. C. and L. R. Scott (1994): The mathematical theory of finite element methods. Springer-Verlag, New York

[CL91] Ciarlet, P.G. and J. L. Lions (1991): Handbook of Numerical Analysis, v. II, Finite Element Methods (Part 1). North-Holland, Amsterdam

[CW71] Ciarlet, P. G. and C. Wagschal (1971): Multipoint Taylor formulas and applications to the finite element method. Numer. Math. 17, 84–100

[CR72] Ciarlet, P. G. and P. A. Raviart (1972): General Lagrange and Hermite interpolation in IR^{N} with applications to finite element methods. Arch. Rational Mech. Anal. 46, 177–199

[D75] Davis, P. J. (1975): Interpolation and approximation. Dover, New York

[E57] Èzrohi, I. A. (1957): General forms of the remainder terms of linear formulæ in multidimensional approximate analysis II. Mat. Sb. (in Russian) 43(85), 9–28

[F74] Fraeijs de Veubeke, B. (1974): Variational principles and the patch test. Intern. J. Numer. Methods Eng. 8, 783–801

[Go77] Gout, J. L. (1977): Estimation de l’erreur d’interpolation d’Hermite dans IR^{n}. Numer. Math. 28, 407–429

[G75] Gregory, J. A. (1975): Error bounds for linear interpolation on triangles. In: J. Whiteman, ed, Mathematics of Finite Elements and Applications, 163–170. Academic Press, London

[K32] Kowalewski, G. (1932): Interpolation und genäherte Quadratur. B. G. Teubner, Berlin

[M75] Meinguet, J. (1975): Realistic estimates for generic constants in multivariate pointwise approximation. In: J. J. H. Miller, ed, Topics in Numerical Analysis II, 89–107. Academic Press, London

[M77] Meinguet, J. (1977): Structure et estimations de coefficients d’erreurs. Rev. Française Automat. Informat. Rech. Opér., Anal. Numer. 11, 355–368

[M78] Meinguet, J. (1978): A practical method for estimating approximation errors in Sobolev spaces. In: D. C. Handscomb, ed, Multivariate Approximation, 169–187. Academic Press, London [M79] Meinguet, J. (1979): A convolution approach to multivariate representation formulas.

In: W. Schempp and K. Zeller, eds, Multivariate Approximation Theory, 198–210. Birkhäuser, Basel

[M81] Meinguet, J. (1981): From Dirac distributions to multivariate representation formulas. In: Z. Ziegler, ed, Approximation Theory and Applications, 225–248. Academic Press, New York

[M82] Meinguet, J. (1982): Sharp “a priori” error bounds for polynomial approximation in Sobolev spaces. In: W. Schempp and K. Zeller, eds, Multivariate Approximation Theory II, 255–274. Birkhäuser, Basel

[M84] Meinguet, J. (1984): A practical method for obtaining a priori error bounds in pointwise and mean-square approximation problems. In: S. P. Singh, J. H. W. Burry, and B. Watson, eds, Approximation Theory and Spline Functions, 97–125. Reidel, Dordrecht

[Mi80] Micchelli, C. A. (1980): A constructive approach to Kergin interpolation in IR^{k}: multivariate B-splines and Lagrange interpolation. Rocky Mountain J. Math. 10, 485–497

[N72] Nicolaides, R. A. (1972): On a class of finite elements generated by Lagrange interpolation. SIAM J. Numer. Anal. 9(3), 435–445

[N73] Nicolaides, R. A. (1973): On a class of finite elements generated by Lagrange interpolation II. SIAM J. Numer. Anal. 10(1), 182–189

[S48] Sard, A. (1948): Integral representations of remainders. Duke Math. J. 15, 333–345 [S63] Sard, A. (1963): Linear approximation. Math. Survey 9, AMS,

[SX95] Sauer, T. and Yuan Xu (1995): On multivariate Lagrange interpolation. Math. Comp. 64, 1147–1170

[W97_{1}] Waldron, S. (1997): Integral error formulæ for the scale of mean value interpolations

which includes Kergin and Hakopian interpolation. Numer. Math. 77(1), 105–122 [W97_{2}] Waldron, S. (1997): A multivariate form of Hardy’s inequality and L_{p}-error bounds for

multivariate Lagrange interpolation schemes. SIAM J. Math. Anal. 28(1), 233–258 [W98] Waldron, S. (1998): The error in linear interpolation at vertices of a simplex. SIAM

J. Numer. Anal. xx, xxx–xxx

[WTDG73] Wilson, E. L., R. L. Taylor, W. P. Doherty, and J. Ghaboussi (1973): Incompatible displacement models. In: S. J. Fenves, N. Perrone, A. R. Robinson, W. C. Schnobrich, eds, Numerical and Computer Methods in Structural Mechanics, 43–57. Academic Press, New York