                             Technical Report
                              17 August 2001

Mean value interpolation for points in general position

Shayne Waldron

Department of Mathematics, University of Auckland, Private Bag 92019, Auckland, New Zealand e–mail: waldron@math.auckland.ac.nz (http:www.math.auckland.ac.nz/˜ waldron)

ABSTRACT

The reduced interpolation conditions for mean value interpolation at a multiset of points in IR^{d} are given explicitly, i.e., the interpolation conditions are expressed using (tangential) derivatives of the lowest possible orders averaged over the smallest possible (polytope) regions. For points in general position the result takes a particularly simple form. Here we give the ‘Lagrange’ form of the extension of the mean value interpolant to less smooth functions, and give an application of it to the finite element method. For points which are not in general position the statement and proof of the result are more technical. These are closely related to an interesting binomial identity which counts the number of ‘allowable’ submultisets of a finite multiset of points in IR^{d}. For example, if the multiset consists of r points in general position each repeated m times, then the identity gives

\binom{mr − 1 + d}{d} = ∑_{j=1}^{d} \binom{r}{j} ∑_{k_{1}=1}^{m} ⋯ ∑_{k_{j}=1}^{m} \binom{k_{1} + ⋯ + k_{j} − 1 + d − j}{d − j},

which reduces to the Chu–Vandermonde convolution identity when m = 1. From our construction we are able to define multivariate analogues of Hermite–Fejér interpolation. We also outline the numerical calculation of the reduced functionals and the corresponding mean value interpolants.

Key Words: mean value interpolation, Kergin interpolation, Hakopian interpolation, multivariate interpolation, finite element method, binomial identity, Chu–Vandermonde identity, simplex spline, Kergin–Fejér interpolation

AMS (MOS) Subject Classifications: primary 05A19, 41A05, 65D07, secondary 41A10, 41A63,

1. Introduction

Kergin interpolation K_{Θ} is a natural multivariate extension of Hermite interpolation, much studied since its introduction in [K80]. It is the first of the mean value interpolation operators

H_{Θ}^{(m)} : C^{n−m}(IR^{d}) → Π_{n−m}(IR^{d}), 0 ≤ m ≤ d − 1

which depend on the n + 1 > m (possibly coincident) ‘interpolation’ points Θ ⊂ IR^{d}. The last of this family H_{Θ}^{(d−1)} is known as the Hakopian map (see [H81]). They are most easily defined using the linear functional

f ↦ ∫_{[θ_{0},...,θ_{k}]} f := \frac{1}{k! vol_{k}(S)} ∫_{S} f ∘ A,  (1.1)

where S is any k–simplex in IR^{s} with (k–dimensional) volume vol_{k}(S), and A : IR^{s} → IR^{d} is any affine map taking the k + 1 vertices of S onto the points θ_{0}, …, θ_{k}. The change of variables formula shows that (1.1) does not depend on the particular choice of S and A, e.g., one can take

∫_{[θ_{0},...,θ_{k}]} 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}.

Then the linear projector H_{Θ}^{(m)} can be defined via its Newton form (see Goodman [G83])

H_{Θ}^{(m)}f(x) := m! ∑_{j=m+1}^{n+1} ∑_{{Ψ⊂Θ_{j−1} \atop #Ψ=m}} ∫_{Θ_{j}} D_{x−Θ_{j−1}\Ψ}f.  (1.2)

Here Ψ is an m–subsequence of Θ_{j−1}, which is the first j − 1 points of Θ = [θ_{0}, …, θ_{n}], and D_{x−Θ_{j−1}\Ψ}f is the (j − 1 − m)–th order derivative

D_{x−Θ_{j−1}\Ψ}f := ( ∏_{{θ∈Θ_{j−1} \atop θ∉Ψ}} D_{x−θ})f, D_{y}f := y_{1}D_{1}f + ⋯ + y_{d}D_{d}f.

Desirable properties of mean value interpolation (cf [BHS93:Ch.12,13]) include:

\underline{Interpolation}: The following mean values are matched

∫_{Ψ} H_{Θ}^{(m)}f = ∫_{Ψ} f, Ψ ⊂ Θ, #Ψ = m + 1.

These are just the function values f(θ_{i}) in the case of Kergin interpolation.

\underline{Symmetry in Θ}: H_{Θ}^{(m)} does not depend on the ordering of the points in Θ.

\underline{PDE Faithfulness}: For any homogeneous constant coefficient partial differential operators p(D) and q(D) of degree k ≤ n − m,

p(D)f = q(D)g ⟹ p(D)(H_{Θ}^{(m)}f) = q(D)(H_{Θ}^{(m)}g).

\underline{Continuity and Coalescence}: For each f ∈ C^{n−m}(IR^{d}), the map

IR^{d×(n+1)} → Π_{n−m}(IR^{d}) : Θ ↦ H_{Θ}^{(m)}f

is continuous. Moreover, if Θ = [a, …, a], then substitution into (1.2) shows that H_{Θ}^{(m)}f is the degree n − m Taylor polynomial of f centred at a.

Because of the coalescence property, the mean value interpolant is defined in general only for C^{n−m} functions. However, for points in general position, we will show that H_{Θ}^{(m)} can be extended to C^{d−1−m} functions. In Theorem 3.1 we give an explicit description of the interpolation conditions of this continuous extension. These ‘reduced’ functionals are decomposed into their dependence on derivatives of lowest possible orders (averaged over simplices with vertices from Θ which have the lowest possible dimension). Since these ‘nodal variables’ are of low smoothness, they invite use in finite elements. This application is considered in Section 4. In Section 5, we describe the extension of H_{Θ}^{(m)} to functions with the least smoothness possible for an arbitrary multiset Θ. This more technical result is closely related to an interesting binomial identity which counts the number of ‘allowable’ submultisets of Θ. For example, if the multiset consists of r points in general position each repeated m times, then the identity gives

\binom{mr − 1 + d}{d} = ∑_{j=1}^{d} \binom{r}{j} ∑_{k_{1}=1}^{m} ⋯ ∑_{k_{j}=1}^{m} \binom{k_{1} + ⋯ + k_{j} − 1 + d − j}{d − j},  (1.3)

which reduces to the Chu–Vandermonde convolution identity when m = 1. From this construction we are able to define multivariate analogues of Hermite–Fejér interpolation. In Section 7, we outline the numerical calculation of the reduced functionals and the corresponding mean value interpolants.

The emphasis of this paper is on the explicit construction of the so called reduced functionals, and applications. Already, in [K80:Remark 4.8] it was known that Kergin interpolation H_{Θ}^{(0)} had a continuous extension to C^{d−1} when the points Θ were in general position in IR^{d}. Not until [BC97] were the reduced functionals described in the bivariate case (d = 2) and a formula for the continuous extension given. A complete discussion of the relationship of the results presented here with the existing literature is given at the end of Section 3.

2. Preliminary results

Here we give a version of the extension [M80:Th.6], which is in the spirit of the following manifestation of the fundamental theorem of calculus

∫_{[Θ,v,w]} D_{v−w}f = ∫_{[Θ,v]} f − ∫_{[Θ,w]} f.  (2.1)

This is used to ‘reduce’ the interpolation conditions of the mean value interpolant.

As usual, let Π_{k}(V) and Π_{k}^{0} := Π_{k}^{0}(V) denote the polynomials defined on a subspace V ⊂ IR^{d} that are of degree ≤ k and homogeneous of degree k, respectively. Let aff(Ψ) denote the affine hull of Ψ ⊂ IR^{d}, and aff(Ψ)^{⊥} the linear space of vectors in IR^{d} orthogonal to it. Then the derivatives of order k which are (in directions) orthogonal to aff(Ψ) can be described by

{q(D) : q ∈ Π_{k}^{0}(aff(Ψ)^{⊥})} = span{D_{n_{1}}D_{n_{2}} ⋯ D_{n_{k}} : n_{1}, …, n_{k} ⊥ aff Ψ},

where q(D) is the differential operator induced by q. If Ψ has affine dimension s, then this space has dimension

dim(Π_{k}^{0}(aff(Ψ)^{⊥})) = dim(Π_{k}^{0}(IR^{d−s})) = ( \begin{matrix}k + d − s − 1 \\ d − s − 1 \end{matrix} ) .  (2.2)

Lemma 2.3 (cf [M80:Th.6]). Suppose Θ is a sequence of two or more points in IR^{d}. Then the linear functional

λ : f ↦ ∫_{Θ} h(D)f, h ∈ Π_{k}^{0}(IR^{d}), k ≥ 1

belongs to A + B where

A := span{f ↦ ∫_{Ψ} p(D)f : Ψ ⊂ Θ, #Ψ = #Θ − 1, p ∈ Π_{k−1}^{0}(IR^{d})},
B := span{f ↦ ∫_{Θ} q(D)f : q ∈ Π_{k}^{0}(aff(Θ)^{⊥})}.

Proof: It is sufficient to assume that h(D) = D_{ξ_{1}}D_{ξ_{2}} ⋯ D_{ξ_{k}} (a spanning set). Express ξ_{1} as

ξ_{1} = y + n_{1}, y ∥ aff(Θ), n_{1} ∈ aff(Θ)^{⊥}.

Since y is parallel to the affine hull of Θ := [θ_{1}, …, θ_{r}], it can be expressed in the form

y = ∑_{i=1}^{r} μ_{i}θ_{i}, μ_{i} ∈ IR, ∑_{i=1}^{r} μ_{i} = 0.

Thus, by [M80:Th.6], we obtain that

∫_{Θ} h(D)f = − ∑_{i=1}^{r} μ_{i} ∫_{Θ\θ_{i}} D_{ξ_{2}} ⋯ D_{ξ_{k}}f + ∫_{Θ} D_{n_{1}}D_{ξ_{2}} ⋯ D_{ξ_{k}}f.  (2.4)

Hence λ can be expressed as the sum of an element of A and the linear functional

f ↦ ∫_{Θ} D_{n_{1}}D_{ξ_{2}} ⋯ D_{ξ_{k}}f.

Repeating this process (on the above functional) an additional k − 1 times (for ξ_{2}, …, ξ_{k}) shows that λ is the sum of an element of A and a linear functional of the form

f ↦ ∫_{Θ} D_{n_{1}}D_{n_{2}} ⋯ D_{n_{k}}f, n_{1}, …, n_{k} ⊥ aff Θ,

which is an element of B.

We will also need the following combinatorial result.

Lemma 2.5 (Chu–Vandermonde convolution). For 0 ≤ m ≤ d − 1, n ≥ m,

dim(Π_{n−m}(IR^{d})) = ( \begin{matrix}n − m + d \\ d \end{matrix} ) = ∑_{j=m+1}^{d} ( \begin{matrix}d − m − 1 \\ d − j \end{matrix} ) ( \begin{matrix}n + 1 \\ j \end{matrix} ) .  (2.6)

Proof: The right hand side counts the number of ways of choosing d objects from n − m + d after having split them into two groups of sizes d − m − 1 and n + 1, and so equals the left hand side.

3. The extension of H_{Θ}^{(m)} to C^{d−1−m} for points in general position

A sequence of points is said to be in general position in IR^{d} if every subsequence has the maximum affine dimension that its cardinality allows, e.g., points are in general position in IR^{2} if no three of them are collinear.

We now show that if Θ is in general position in IR^{d}, then the mean value interpolation operators can be extended from C^{n−m} to C^{d−1−m} functions. The argument given below can be viewed as the multivariate analogue of taking the Newton form for Hermite interpolation at points in general position in IR (no points repeated) and using the divided difference identity (cf (2.1) and Lemma 2.3) to obtain the Lagrange form.

A finite rank linear projector P : X → X is uniquely determined by its range V, and the range of its dual Λ, called its interpolation conditions. That is P f is the unique element of V matching Λ, i.e.,

λ(P f) = λ(f), ∀λ ∈ Λ,

where dim(V) = dim(Λ). Denote by conv(Θ) the convex hull of Θ.

Theorem 3.1 (Reduced functionals). Suppose the n + 1 points of Θ are in general position in IR^{d}. Then the mean value interpolation map has a continuous extension

H_{Θ}^{(m)} : C^{d−1−m}(conv(Θ)) → Π_{n−m}(IR^{d}),

which is uniquely determined by the interpolation conditions

Λ := ⨁_{j=m+1}^{d} ⨁_{{Ψ⊂Θ \atop #Ψ=j}} Λ_{Ψ}^{(m)},  (3.2)

where

Λ_{Ψ}^{(m)} := span{f ↦ ∫_{Ψ} q(D)f : q ∈ Π_{j−m−1}^{0}(aff(Ψ)^{⊥})}.  (3.3)

Proof: First we show the interpolation conditions of H_{Θ}^{(m)} are contained in Λ. By using (for simplicity) the basis {D^{α} : α ∈ ZZ_{+}^{d}, |α| ≤ n − m} for the space of constant coefficient differential operators of order ≤ n − m, the Newton form (1.2) can be expanded

H_{Θ}^{(m)}f(x) = ∑_{j=m+1}^{n+1} ∑_{|α|=j−1−m} p_{α}(x)λ_{α}(f),  λ_{α} : f ↦ ∫_{Θ_{j}} D^{α}f,  (3.4)

where p_{α} ∈ Π_{n−m}(IR^{d}). Since H_{Θ}^{(m)} is a linear projector onto Π_{n−m}(IR^{d}), and the number of terms in (3.4) is

∑_{k=0}^{n−m} #{D^{α} : |α| = k} = dim(Π_{n−m}(IR^{d})),

the polynomials (p_{α}) are a basis for Π_{n−m}(IR^{d}) with dual basis (λ_{α}), and (λ_{α}) is a basis for the interpolation conditions. Since the j points of Θ_{j} are in general position,

dim(Π_{k}^{0}(aff(Ψ)^{⊥})) = 0, Ψ ⊂ Θ_{j}, #Ψ ≥ d + 1,

and so j − 1 − m applications of Lemma 2.3 to λ_{α} : f ↦ ∫_{Θ_{j}} D^{α}f gives

λ_{α} ∈ span{f ↦ ∫_{Ψ} q(D)f : Ψ ⊂ Θ_{j}, m + 1 ≤ #Ψ ≤ d, q ∈ Π_{#Ψ−m−1}^{0}(aff(Ψ)^{⊥})} ⊂ Λ.

Since λ_{α} ∈ Λ is continuous on C^{d−1−m}(conv(Θ)), the desired extension of H_{Θ}^{(m)} is given by

H_{Θ}^{(m)}f := ∑_{j=m+1}^{n+1} ∑_{|α|=j−1−m} p_{α}λ_{α}(f), ∀f ∈ C^{d−1−m}(conv(Θ)),

and it is unique, since (the restriction of) C^{n−m}(IR^{d}) is dense in C^{d−1−m}(conv(Θ)). It remains only to show that the sum in (3.2) is direct. By (2.2) and Lemma 2.5,

dim(Λ) ≤ ∑_{j=m+1}^{d} ∑_{{Ψ⊂Θ \atop #Ψ=j}} dim(Λ_{Ψ}^{(m)}) ≤ ∑_{j=m+1}^{d} ( \begin{matrix}d − m − 1 \\ d − j \end{matrix} ) ( \begin{matrix}n + 1 \\ j \end{matrix} ) = dim(Π_{n−m}(IR^{d})).

Since dim(Λ) = dim(Π_{n−m}(IR^{d})), there must be equality throughout above, and so by a dimension count the sum is direct.

The Lagrange form

For each Ψ, let B_{Ψ} = B_{Ψ}^{(m)} be some basis for Π_{#Ψ−m−1}^{0}(aff(Ψ)^{⊥}). Then the extension can be represented in the Lagrange form

H_{Θ}^{(m)}f(x) = ∑_{j=m+1}^{d} ∑_{{Ψ⊂Θ \atop #Ψ=j}} ∑_{q∈B_{Ψ}} p_{Ψ,q}(x) ∫_{Ψ} q(D)f,  ∀f ∈ C^{d−1−m}(conv(Θ)), (3.5)

where (p_{Ψ,q}) is the basis of Π_{n−m}(IR^{d}) which is dual to the basis of interpolation conditions

{f ↦ ∫_{Ψ} q(D)f : Ψ ⊂ Θ, m + 1 ≤ #Ψ ≤ d, q ∈ B_{Ψ}}.

We refer to (3.5) as the Lagrange form, as opposed to the Newton form (1.2), since it involves derivatives of f of the lowest possible orders. In [BW01] it is shown that the ‘Lagrange polynomials’ p_{Ψ,q} are multi–ridge functions of d + m + 1 − #Ψ variables.

f ↦ ∫_{[a,b]} D_{n}f  conditions^{\binom{n+1}{3}} conditions 2\binom{n+1}{2}  f ↦ ∫_{[a,b,c]} D_{n}^{2}f f ↦ f(a) n + 1 conditions

Fig. 3.1. The \binom{n+3}{3} reduced interpolation conditions for the Lagrange form of Kergin interpolation at n + 1 points in general position in IR^{3}.

Examples

In IR, points in general position are simply distinct points, and so H_{Θ}^{(0)} is Lagrange interpolation at Θ, with (1.2) giving the classical Newton form, via the Hermite–Genocchi formula for divided differences

[θ_{0}, …, θ_{k}]f = ∫_{[θ_{0},...,θ_{k}]} D^{k}f.  (3.6)

In (5.6) all the Ψ consist of a single point, with Π_{#Ψ−0−1}^{0}(aff(Ψ)^{⊥}) = Π_{0}(IR^{d}), for which it is natural to choose the basis B_{Ψ} = {1}. Then (3.5) is the classical Lagrange form.

In IR^{2}, the extension of Kergin interpolation K_{Θ} := H_{Θ}^{(0)} is defined on C^{1} functions, and the subsets Ψ in (3.2) consist of either one or two points. For one point subsets we have function evaluation, while for subsets Ψ = {θ_{s}, θ_{t}}, the average over the line segment from θ_{s} to θ_{t} of the derivative in the direction (θ_{t} − θ_{s})^{⊥} orthogonal to it is matched. By choosing the bases B_{Ψ} = {1} and B_{Ψ}(D) = {D_{(θ_{t}−θ_{s})^{⊥}}} in each case respectively, the Lagrange form (3.5) for Θ = (θ_{0}, …, θ_{n}) becomes

K_{Θ}f = ∑_{r=0}^{n} p_{r}f(θ_{r}) + ∑_{0≤s≤t≤n} p_{st} ∫_{[θ_{s},θ_{t}]} D_{(θ_{t}−θ_{s})^{⊥}}f,  ∀f ∈ C^{1}(conv Θ).  (3.7)

This extension, together with some beautiful formulæ (involving the identification of IR^{2} with C) for the Lagrange polynomials p_{r}, p_{st} ∈ Π_{n}(IR^{2}), was recently obtained by Bos and Calvi [BC97:Th.2.3].

In IR^{3}, the extension K_{Θ} is defined for C^{2} functions. In addition to reduced interpolation conditions like those in (3.7), the data ∫_{[a,b,c]} D_{n}^{2}f is matched, where n is a normal to the plane through a, b, c ∈ Θ (see Fig. 3.1). For four points these 20 reduced functionals were derived in Lorentz [L92:Th.13.2.4] (also see the applications below).

For the Hakopian map H_{Θ}^{(d−1)}, only subsets Ψ with d elements are taken in (3.2), and so the original presentation (which required the points be in general position) is obtained.

In Cavaretta, Goodman, Micchelli and Sharma [CGMS83] the extension of Theorem 3.1 is given in the cases m = d − 2 (Theorem 2) and m = d − 3 (Theorem 3), together with some explicit formulæ for the corresponding Lagrange polynomials (Theorems 4 and 5).

4. An application to finite elements

The extended mean value interpolation maps of Theorem 3.1 give rise to a family of finite elements. First we recall some definitions (see, e.g., Brenner and Scott [BS94]).

Definition 4.1. Let

(i) K ⊂ IR^{d} 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. 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.

If Θ ⊂ IR^{d} is a sequence of n + 1 points in general position, n ≥ d, then K := conv(Θ) is an element domain, and (3.5) is the local interpolant of a finite element with shape functions P = Π_{n−m}(IR^{d}), and nodal variables given by the interpolation conditions (3.2). We now state this result and give some examples.

Theorem 4.4 (Finite elements). Suppose Θ ⊂ IR^{d} is a sequence of n + 1 points in general position, n ≥ d > m ≥ 0. Let K := conv(Θ), P := Π_{n−m}(IR^{d}), and N be the nodal variables given by

f ↦ ∫_{Ψ} q(D)f,  Ψ ⊂ Θ, m + 1 ≤ #Ψ ≤ d, q ∈ B_{Ψ},  (4.5)

where B_{Ψ} is some basis for Π_{#Ψ−m−1}^{0}(aff(Ψ)^{⊥}). Then (K, P, N) is a finite element with its local interpolant given by (3.5).

For the general theory of the error in such a local interpolant see [C78:Ch.3]. This effectively says that the reproduction of the polynomials Π_{n−m}(IR^{d}) gives approximation order h^{m−n+1} where h is the diameter of the element domain K. Estimates of the constants involved here can be found using an error formula for the local interpolant, involving only the (n − m + 1)–th derivative of f, which was given in [W97_{1}:Th.5.12], and applying to it a multivariate version of Hardy’s inequality to get L_{p}–estimates for p > d (see [W97_{2}:Th.4.3.1]).

Examples

The simplest bivariate example is Θ = (θ_{0}, θ_{1}, θ_{2}) the vertices of a triangle K, and m = 0 (Kergin). Here P = Π_{2}(IR^{2}), and the data matched by the nodal variables (4.5) is

f(θ_{0}), f(θ_{1}), f(θ_{2}), ∫_{θ_{0}}^{θ_{1}} D_{(θ_{1}−θ_{0})^{⊥}}f, ∫_{θ_{1}}^{θ_{2}} D_{(θ_{2}−θ_{1})^{⊥}}f, ∫_{θ_{0}}^{θ_{2}} D_{(θ_{2}−θ_{0})^{⊥}}f,  (4.6)

with ∫_{a}^{b} denoting the line integral over the segment from a to b. This finite element is given in Lorentz [L92:Th.13.2.2,Cor.13.2.3]. For Θ = (θ_{0}, θ_{1}, θ_{2}, θ_{3}) the vertices of a convex quadrilateral K, P = Π_{3}(IR^{2}), and the 10 nodal variables match the values

f(θ_{r}), 0 ≤ r ≤ 3, ∫_{θ_{s}}^{θ_{t}} D_{(θ_{t}−θ_{s})^{⊥}}f, 0 ≤ s < t ≤ 3,

and so forth. See Bos and Calvi [BC97:Th.2.3] for a description of the local interpolant.

In IR^{3}, the first example is Θ = (θ_{0}, θ_{1}, θ_{2}, θ_{3}) the vertices of a tetrahedron K. Here P = Π_{3}(IR^{3}), and the data matched by the nodal variables is the 4 function values f(θ_{r}), 0 ≤ r ≤ 3, the 12 values of the form

∫_{θ_{s}}^{θ_{t}} D_{(θ_{t}−θ_{s})^{⊥}}f, 0 ≤ s < t ≤ 3,

where D_{(θ_{t}−θ_{s})^{⊥}}f now runs over some basis for the 2–dimensional space of derivatives in directions orthogonal to the line through θ_{t} and θ_{s}, and the 4 values

∫_{conv(θ_{r},θ_{s},θ_{t})} D_{n}^{2}f, 0 ≤ r < s < t ≤ 3,

where the integration above is over the triangle with vertices {θ_{r}, θ_{s}, θ_{t}}, and n is some vector orthogonal to this triangle. This finite element is given in [L92:Th.13.2.4,Cor.13.2.5].

It might be hoped to obtain more practicable finite elements by replacing the averages occuring in the nodal variables by averages over sets of smaller affine dimension (ideally point evaluations). For example, Morley’s finite element (see Ciarlet [C78:p.374]) can be obtained by replacing (4.6) by

f(θ_{0}), f(θ_{1}), f(θ_{2}), D_{(θ_{1}−θ_{0})^{⊥}}f(\frac{θ_{0} + θ_{1}}{2} ), D_{(θ_{2}−θ_{1})^{⊥}}f(\frac{θ_{1} + θ_{2}}{2} ), D_{(θ_{2}−θ_{0})^{⊥}}f(\frac{θ_{0} + θ_{2}}{2} ).

The author is currently unaware of any general results in this direction.

5. The reduced functionals for H_{Θ}^{(m)} where Θ is any multiset

Here we find the reduced functionals for H_{Θ}^{(m)} where Θ is any sequence (multiset) of n + 1 points in IR^{d}. By considering the basis (λ_{α}) of (3.4), it follows the interpolation conditions of H_{Θ}^{(m)} are contained in

∑_{{Φ⊂Θ \atop #Φ≥m+1}} L^{(m)}(Φ), L^{(m)}(Φ) := span\{f ↦ ∫_{Φ} q(D)f : q ∈ Π_{i−1−m}^{0}(IR^{d}), i = #Φ\}.

As before, our plan is to ‘reduce’ the spaces L_{Θ}^{(m)}(Φ) as much as possible, so that we obtain a direct sum decomposition of reduced functionals. Since now the points need not be in general position this is more delicate, requiring the introduction of certain spaces of simplex splines to describe the averaging of the tangential derivatives. We are also unable to appeal to the Chu-Vandermonde convolution to show the sum we obtain is direct, though once proved we show that it implies a most interesting generalisation of the Chu-Vandermonde convolution identity. The result for Hakopian interpolation (m = d − 1) has already been obtained in [H84], although there the directness of the sum is not proved.

The (affine) dimension of a submultiset Ψ ⊂ Θ is defined to be the dimension of the affine hull of the underlying set, i.e., dim(Ψ) := dim(aff(Ψ)). For example, if Ψ consists of two distinct points (giving a line), each repeated a number of times, then its dimension is 1. The multiplicity of a point in multiset is its repetition number. Let S ∩^{∗} Θ denote the intersection of an affine space S ⊂ IR^{d} with a multiset Θ ⊂ IR^{d} counting multiplicities. For example, Θ ∩^{∗} IR^{d} is the multiset Θ, whilst Θ ∩ IR^{d} is the underlying set.

The simplex spline corresponding to a finite multiset (or sequence) Θ := [θ_{0}, …, θ_{k}] of points in IR^{d} (sometimes called knots) is defined to be the representer of the distribution

f ↦ ∫_{Θ} f, ∫_{[]} f := 0,

defined by (1.1). It is a nonnegative piecewise polynomial function M(·|Θ) defined on aff(Θ), supported on conv(Θ), with the normalisation

∫_{aff(Θ)} M(·|Θ)f = ∫_{conv(Θ)} M(·|Θ)f = k! ∫_{Θ} f,  (5.1)

which ensures it integrates to 1. Properties of these splines have been studied in great detail (see, e.g., [BHS93]). We are interested only in a particular spline space S_{i}(Ψ), which occurs in our construction of the reduced functionals. This space of functions on aff(Ψ) is defined as

S_{i}(Ψ) := span{M(·|Φ) : Φ ∈ A_{i}(Ψ)}, i = 1, …, #Ψ,

where

A_{i}(Ψ) := {Φ ⊂ Ψ : #Φ = i, dim(Φ) = dim(Ψ)}.

Note S_{i}(Ψ) = 0 if i + 1 < dim(Ψ). When A_{i}(Ψ) = {Φ ⊂ Ψ : #Φ = i}, as is the case for Ψ in general position in aff(Ψ), Dahmen and Micchelli [DM83:Th.1] have shown that

dim(S_{i}(Ψ)) = \binom{#Ψ − i + dim(Ψ)}{dim(Ψ)}.

Fig. 5.1. Simplex splines for points equally spaced on a circle in IR^{2}: 3 knots (discontinous piecewise constant function), 4 knots (continuous piecewise linear function), 5 knots (C^{1} piecewise quadratic function), and 6 knots (C^{2} piecewise cubic function).

Fig. 5.2. Simplex splines for 5 points in IR^{2} which are not in general position (piecewise quadratics). In the first 3 knots are collinear, and the spline is only C^{0} across this line. In the second 4 knots are collinear, and the spline is discontinuous across this line.

Lemma 5.2. Let Ψ consist of j affinely independent points with multiplicities m_{1}, …, m_{j}. Then for i = j, …, #Ψ

dim(S_{i}(Ψ)) = #A_{i}(Ψ) = c_{i}(m_{1}, …, m_{j}),  (5.3)

where

c_{i}(m_{1}, …, m_{j}) := coefficient of t^{i} in ∏_{k=1}^{j}(t + t^{2} + ⋯ + t^{m_{k}}).  (5.4)

Proof: The submultisets in A_{i}(Ψ) consist of those Φ of size i which contain each of the j points at least once. There is not a general formula for c_{i}(m_{1}, …, m_{j}), the number of these, but it can be calculated from the generating function (5.4), or by using the principle of inclusion and exclusion.

Let ξ_{1}, …, ξ_{j} be the barycentric coordinates of a point in aff(Ψ) relative to the j points, and k_{1}, …, k_{j} be the multiplicities of these points as they occur in Φ ∈ A_{i}(Ψ). Then the restriction of M(·|Φ) to the simplex T := conv(Ψ) is a scalar multiple of

ξ_{1}^{k_{1}}ξ_{2}^{k_{2}} ⋯ ξ_{j}^{k_{j}} (k_{1} + ⋯ + k_{j} = i).

It is well known that these polynomials of degree i which appear in the definition of the Bernstein operator (of degree i for the simplex T) are linearly independent. Hence the M(·|Φ) are linearly independent, and so dim(S_{i}(Ψ)) = #A_{i}(Ψ).

Theorem 5.5 (Reduced functionals). Let Θ be any sequence (multiset) of n+1 points in IR^{d}. Then Λ_{Θ}^{(m)} the interpolation conditions for H_{Θ}^{(m)} can be written as the direct sum

Λ_{Θ}^{(m)} = ⨁_{j=1}^{d} ⨁_{{Ψ⊂Θ,#Ψ≥m+1 \atop dim(Ψ)=j−1} aff(Ψ)∩^{∗}Θ=Ψ} ⨁_{i=m+1}^{#Ψ} Λ_{Ψ,i}^{(m)},  (5.6)

where

Λ_{Ψ,i}^{(m)} := span{f ↦ ∫_{Φ} q(D)f : q ∈ Π_{i−1−m}^{0}(aff(Ψ)^{⊥}), Φ ∈ A_{i}(Ψ)},  (5.7)

and

dim(Λ_{Ψ,i}^{(m)}) = \binom{i − 1 − m + d − j}{d − j} dim(S_{i}(Ψ)),  j − 1 := dim(Ψ).  (5.8)

Proof: Suppose λ ∈ L^{(m)}(Φ), Φ ⊂ Θ, #Φ = i ≥ m + 1, then, by Lemma 2.3,

λ ∈ Λ_{Ψ,i}^{(m)} + ∑_{{Ξ⊂Φ \atop #Ξ=i−1}} L^{(m)}(Ξ),  Ψ := aff(Φ) ∩^{∗} Θ.  (5.9)

Note #Ψ ≥ i and Λ_{Ψ,i}^{(m)} = 0 if dim(Ψ) = d. Since Λ_{Θ}^{(m)} ⊂ ∑_{{Φ⊂Θ \atop #Φ≥m+1}} L^{(m)}(Φ), repeated application of (5.9) gives

Λ_{Θ}^{(m)} ⊂ ∑_{j=1}^{d} ∑_{{Ψ⊂Θ,#Ψ≥m+1 \atop dim(Ψ)=j−1} aff(Ψ)∩^{∗}Θ=Ψ} ∑_{i=m+1}^{#Ψ} Λ_{Ψ,i}^{(m)}.  (5.10)

The reverse inclusion follows from (5.7), and so it suffices to show the sum in (5.10) is direct, i.e.,

λ ∈ Λ_{Ψ_{1},i_{1}}^{(m)} ∩ Λ_{Ψ_{2},i_{2}}^{(m)}, (Ψ_{1}, i_{1}) ≠ (Ψ_{2}, i_{2}) ⟹ λ = 0.

If aff(Ψ_{1}) ∩ aff(Ψ_{2}) = ϕ, then by considering the support of the distribution λ, we have

supp(λ) ⊂ aff(Ψ_{1}) ∩ aff(Ψ_{2}) = ϕ ⟹ λ = 0.

Hence, it remains to consider the case aff(Ψ_{1}) ∩ aff(Ψ_{2}) ≠ ϕ. For this, our proof relies on the identification of Λ_{Ψ,i}^{(m)} with a space of functions F_{Ψ,i}^{(m)}. Since IR^{d} = aff(Ψ)^{⊥} ⊕ aff(Ψ), the tensor product

F_{Ψ,i}^{(m)} := Π_{i−1−m}^{0}(aff(Ψ)^{⊥}) ⊗ S_{i}(Ψ)
:= span{f = p ⊗ s : p ∈ Π_{i−1−m}^{0}(aff(Ψ)^{⊥}), s ∈ S_{i}(Ψ)},

(5.11)

with

f(x + y) := (p ⊗ s)(x, y) = p(x)s(y), ∀x ∈ aff(Ψ)^{⊥}, ∀y ∈ aff(Ψ),

is a space of functions on IR^{d}. Note, by (5.1), that

Λ_{Ψ,i}^{(m)} = span{f ↦ ∫_{aff(Ψ)}(p(D)f)s : p ∈ Π_{i−1−m}^{0}(aff(Ψ)^{⊥}), s ∈ S_{i}(Ψ)}.  (5.12)

Let P_{Ψ,i}^{(m)} : Λ_{Ψ,i}^{(m)} → F_{Ψ,i}^{(m)} be the linear map λ ↦ p ⊗ s, where λ(f) = ∫_{aff(Ψ)}(p(D)f)s. This map is clearly onto, and it is 1–1 since

p ⊗ s = 0, s ≠ 0 ⟹ p(x)s(y) = 0 ⟹ p(x) = 0, ∀x ∈ aff(Ψ)^{⊥} ⟹ p = 0.

Hence the linear spaces Λ_{Ψ,i}^{(m)} and F_{Ψ,i}^{(m)} are isomorphic, and in particular,

dim(Λ_{Ψ,i}^{(m)}) = dim(F_{Ψ,i}^{(m)}) = dim(Π_{i−1}^{0}(aff(Ψ)^{⊥})) dim(S_{i}(Ψ))

= \binom{i − 1 − m + d − j}{d − j} dim(S_{i}(Ψ)),

which is (5.8). In view these isomorphisms, it is sufficient to show that the sum

F_{Ψ_{1},i_{1}}^{(m)} + F_{Ψ_{2},i_{2}}^{(m)},  (Ψ_{1}, i_{1}) ≠ (Ψ_{2}, i_{2}), aff(Ψ_{1}) ∩ aff(Ψ_{2}) ≠ ϕ,  (5.13)

is direct, i.e., that

f ∈ F_{Ψ_{1},i_{1}}^{(m)} ∩ F_{Ψ_{2},i_{2}}^{(m)}, ⟹ f = 0.

This we do in two cases. Note, by (5.11), that

supp(f) = supp(f|_{conv(Ψ)}) + aff(Ψ)^{⊥},  ∀f ∈ F_{Ψ,i}^{(m)},  (5.14)

and if L is a line in IR^{d}

L is orthogonal to aff(Ψ) ⟹ f|_{L} ∈ Π_{i−1−m}^{0}(L).  (5.15)

Case 1: aff(Ψ_{1})^{⊥} = aff(Ψ_{2})^{⊥}, i.e., aff(Ψ_{1}) and aff(Ψ_{2}) are parallel. The condition aff(Ψ_{1}) ∩ aff(Ψ_{2}) ≠ ϕ implies aff(Ψ_{1}) = aff(Ψ_{2}), i_{1} ≠ i_{2}. By (5.15), the restriction of f to each line orthogonal to aff(Ψ_{1}) is a homogeneous polynomial of degree i_{1} − 1 − m and i_{2} − 1 − m, i.e., the zero polynomial (i_{1} ≠ i_{2}). Since aff(Ψ_{1}) ≠ IR^{d} (dim(Ψ_{1}) = j − 1 < d), this implies f = 0.

Case 2: aff(Ψ_{1})^{⊥} ≠ aff(Ψ_{2})^{⊥}. If B is a ball in IR^{d} containing Ψ_{1} and Ψ_{2}, then

supp(f|_{conv(Ψ_{i})}) ⊂ supp(f) ∩ B =: K, i = 1, 2.

Since aff(Ψ_{i})^{⊥} is a linear subspace, using (5.14), this implies

supp(f) = supp(f|_{conv(Ψ_{i})}) + aff(Ψ_{i})^{⊥} ⊂ K + aff(Ψ_{i})^{⊥} ⊂ supp(f) + aff(Ψ_{i})^{⊥} = supp(f),

giving

K + aff(Ψ_{1})^{⊥} = K + aff(Ψ_{2})^{⊥}.  (5.16)

If f ≠ 0, then K ≠ ϕ is compact, and (5.16) implies

aff(Ψ_{1})^{⊥} = aff(Ψ_{2})^{⊥},

which contradicts our assumption.

This completes the proof.

Example 1 (Points in general position). For Θ in general position, aff(Ψ) ∩^{∗} Θ = Ψ, dim(Ψ) = #Ψ + 1, A_{i}(Ψ) = 0 if i < #Ψ and A_{i}(Ψ) = {Ψ} if i = #Ψ, so (5.6) becomes

Λ_{Θ}^{(m)} = ⨁_{j=m+1}^{d} ⨁_{{Ψ⊂Θ \atop #Ψ=j}} Λ_{Ψ,j}^{(m)},  (5.17)

where Λ_{Ψ,j}^{(m)} = Λ_{Ψ}^{(m)} of (3.3), and we recover Theorem 3.1.

Example 2 (All points coincident). If all points are the same, i.e., Θ = [a, …, a], only Ψ = Θ satisfies aff(Ψ) ∩^{∗} Θ = Ψ. Since aff(Ψ)^{⊥} = IR^{d}, and ∫_{Φ} f, Φ ⊂ Θ, #Φ ≥ 1 is a multiple of f(a), the direct sum (5.6) becomes

⨁_{i=m+1}^{n+1} Λ_{Θ,i}^{(m)}, Λ_{Θ,i}^{(m)} = span{f ↦ (q(D)f)(a) : q ∈ Π_{i−1−m}^{0}(IR^{d})},

i.e., the interpolation conditions for Taylor interpolation of degree n − m at a.

Example 3 (Hermite interpolation). Let Θ consist of the points x_{1}, …, x_{r} in IR repeated with multiplicities m_{1}, …, m_{r}. Then the Ψ in (5.6) consist of a point x_{k} repeated m_{k} times, giving

Λ_{Θ}^{(0)} = ⨁_{k=1}^{r} ⨁_{i=1}^{m_{k}} Λ_{[^{x_{k},...,x_{k}}_{m_{k}}],i} = ⨁_{k=1}^{r} ⨁_{i=1}^{m_{k}} span{f ↦ D^{i−1}f(x_{k})},  (5.18)

which is the usual presentation of the interpolation conditions for Hermite interpolation. By Hermite–Genocchi (3.6), the interpolation conditions (λ_{α}) in (3.4) are equivalent to

Λ_{Θ}^{(0)} = span{f ↦ [θ_{1}, …, θ_{r}]f : [θ_{1}, …, θ_{r}] ⊂ Θ, r = 1, …, n},

where [θ_{1}, …, θ_{r}]f is the divided difference of f at the sequence of points θ_{1}, …, θ_{r}.

Example 4 (Points on a circle in IR^{d}). Let Θ consist of n+1 points in general position in some plane in IR^{d}, d ≥ 2, e.g., points on some circle. Then the Ψ in (5.6) are those of size 1, 2 and n + 1. For the Ψ of sizes j = 1, 2, the corresponding (reduced) interpolation conditions are

f(θ_{r}), \binom{n + 1}{1} conditions, \} j = 1

∫_{[θ_{r},θ_{t}]} D_{(θ_{t}−θ_{r})^{⊥}}f, \binom{n + 1}{2}(d − 1) conditions, \} j = 2

where θ_{r}, θ_{t} ∈ Θ, θ_{r} ≠ θ_{t}. The only other Ψ occuring is Ψ = Θ, for j = 3. When Θ consists of n + 1 points in general position in an s–dimensional subspace of IR^{d}, it is possible to give an explicit subset of A_{i}(Θ) whose corresponding splines are a basis of S_{i}(Θ). This s–generating sequence of [DM83:Th.4.2] is defined as follows

B_{i}(Θ) := {[θ_{1}, …, θ_{i−s}] ∪ Φ : Φ ⊂ [θ_{i−s+1}, …, θ_{n+1}], #Φ = s},

where Θ = [θ_{1}, …, θ_{i−s}, θ_{i−s+1}, …, θ_{n+1}]. In particular, for our points on a circle (s = 2) the remaining interpolation conditions (for Ψ = Θ) can be expressed as: for i = 3, …, n+1

∫_{[θ_{1},...,θ_{i−2},θ_{r},θ_{t}]} D_{⊥}^{i−1}f, \binom{i + d − 4}{d − 3}\binom{n + 3 − i}{2} conditions\} j = 3,

where i − 2 < r < t ≤ n + 1.

Example 5 (Hakopian interpolation). For Hakopian interpolation (m = d − 1) with aff(Θ) = d the result was obtained in [H84]. In the notation used there our sum ⊕_{j=1}^{d} is described in two parts. For ⊕_{j=2}^{d} the set of Ψ (dim Ψ = j−1) in the second sum corresponds to Hakopian’s set L_{p}, 1 ≤ p = j − 1 ≤ d − 1 of all p–dimensional hyperplanes spanned by points from Θ, and for a fixed L_{p}^{q} ∈ L_{p} the B_{α,L_{p}^{q}∩{x^{0},...,x^{r}}}^{p} are bases of simplex splines for our S_{i}(Ψ), i = d, …, #ψ. For ⊕_{j=1} only the points a with multiplicity r := #Ψ ≥ d, Ψ := {a} ∩^{∗} Θ contribute interpolation conditions, namely

Λ_{Ψ,i}^{(d−1)} = span{f ↦ (q(D)f)(a) : q ∈ Π_{i−d}(IR^{d})}, i = d, …, r

which are clearly the same as those given by Hakopian. Hakopian states than an inductive argument shows the number of interpolation conditions equals dim(Π_{n−d+1}(IR^{d})) = \binom{n+1}{d}, and so the sum is direct, but it is not included.

Example 6 (Extensions of H_{Θ}^{(m)}). The linear functionals in (5.6) are called reduced interpolation conditions because the derivatives in them are tangential to the polytopes they are averaged over, and so the order of these derivatives cannot be decreased any further by application of Lemma 2.3. As in Theorem 3.1 the mean value interpolation operator H_{Θ}^{(m)} can be extended continuously to any functions these reduced functionals are defined on. In particular, there is an extension to C^{r}(IR^{d}) → Π_{n−m}(IR^{d}), where

r := max_{{Ψ⊂Θ,#Ψ≥m+1 \atop dim(Ψ)<d} aff(Ψ)∩^{∗}Θ=Ψ} (#Ψ − 1 − m).

6. Chu–Vandermonde convolution identities

Since the sum (5.6) is direct we obtain the following dimension count.

Corollary 6.1 (Chu–Vandermonde). For any multiset Θ of n + 1 points in IR^{d} and 0 ≤ m ≤ d − 1,

\binom{n − m + d}{d} = ∑_{j=1}^{d} ∑_{{Ψ⊂Θ,#Ψ≥m+1 \atop dim(Ψ)=j−1} aff(Ψ)∩^{∗}Θ=Ψ} ∑_{i=m+1}^{#Ψ} \binom{i − 1 − m + d − j}{d − j} dim(S_{i}(Ψ)),  (6.2)

where

1 ≤ dim(S_{i}(Ψ)) ≤ \binom{#Ψ − i + j − 1}{j − 1}.  (6.3)

The upper bound in (6.3) is attained when A_{i}(Ψ) = {Φ ⊂ Ψ : #Φ = i}. This happens if the underlying point set of Ψ is in general position in the (j − 1)–dimensional affine space that it spans, and the number of repeated points Ψ \ (aff(Ψ) ∩ Ψ) is ≤ i − j.

Proof: Since (3.2) is direct, dim(Λ_{Θ}^{(m)}) = dim(Π_{n−m}(IR^{d})) = \binom{n−m+d}{d} equals the sum of the dimensions of the Λ_{Ψ,i}^{(m)} which are given by (5.8). The estimate (6.3) follows from the proof of Lemma 4.1 in Dahmen and Micchelli [DM83], where it is shown

dim(S_{i}(Θ)) = dim span{x ↦ ∏_{θ∈Φ}(1 + θ^{∗}x) : Φ ⊂ Θ, #Φ = n − i}

≤ dim(Π_{n−i}(IR^{s})) = \binom{n − i + s}{s},

together with the asserted equality. Furthermore, if the stated condition holds, then every Φ ⊂ Ψ of size i contains j affinely independent points, and so belongs to A_{i}(Ψ), i.e., A_{i}(Ψ) = {Φ ⊂ Ψ : #Φ = i}.

When the points Θ are in general position we obtain from (5.17) the Chu-Vandermonde convolution (2.6). It is possible to interpret (6.2) as a way of counting the ‘allowable’ submultisets of a finite multset of points Θ in IR^{d}. We will illustrate this for m = 0.

Definition 6.4. A submultiset Φ ⊂ Θ is said to be allowable if

(a) Θ ∩ conv(Φ) ⊂ Φ, i.e., all points of Θ in conv(Φ) are in Φ.

(b) points of Θ in the relative interior of conv(Φ) are repeated in Φ with full multiplicity. For short we will write Φ ⊂^{∗} Θ.

Theorem 6.5 (Allowable submultisets). Suppose Θ is a multiset of n + 1 points in IR^{d}, and either

(a) The points Θ ∩ IR^{d} are in general position in an s–dimensional affine subspace of IR^{d}, where s = 0, s = 1 or s = d.

(b) IR^{d} = IR or IR^{d} = IR^{2}. then

\binom{n + d}{d} = ∑_{i=1}^{n+1} ∑_{{Φ⊂^{∗}Θ \atop #Φ=i}} \binom{i + d − dim(Φ) − 2}{d − dim(Φ) − 1},  (6.6)

where the sum in (6.6) is over all allowable submultisets (of size i).

Proof: Every allowable Φ ⊂ Θ which contributes to the sum (6.6) belongs to a unique Ψ with aff(Ψ) ∩^{∗} Θ = Ψ and dim(Ψ) = j − 1 < d, namely Ψ := aff(Φ) ∩^{∗} Θ. Hence

∑_{i=1}^{n+1} ∑_{{Φ⊂^{∗}Θ \atop #Φ=i}} \binom{i + d − dim(Φ) − 2}{d − dim(Φ) − 1}

= ∑_{j=1}^{d} ∑_{{Ψ⊂Θ \atop dim(Ψ)=j−1} aff(Ψ)∩^{∗}Θ=Ψ} ∑_{i=j}^{#Ψ} \binom{i − 1 + d − j}{d − j}#\{Φ ∈ A_{i}(Ψ) : Φ is allowable\}.

Thus, (6.2) with m = 0 implies (6.6), provided

#{Φ ∈ A_{i}(Ψ) : Φ is allowable} = dim(S_{i}(Ψ)),  i = j, …, #Ψ.  (6.7)

We now show (6.7) holds in the cases (a) and (b).

Case (a): Suppose (a) holds. Then the possible values of dim(Ψ) are 0, 1, …, s. Those Ψ with aff(Ψ)∩^{∗} Θ = Ψ, and dim(Ψ)+1 = j = 1, 2, …, s, consist of j affinely independent points repeated with full multiplicity. Fix i = j, …, #Ψ. Then each Φ ∈ A_{i}(Ψ) is allowable, and so Lemma 5.2 gives

{Φ ∈ A_{i}(Ψ) : Φ is allowable} = #A_{i}(Ψ) = dim(S_{i}(Ψ)),

which is (6.7). If dim(Ψ) = s, then there is exactly one Ψ with aff(Ψ) ∩^{∗} Θ = Ψ, namely Ψ = Θ, and so (6.7) is only required to hold when s < d. For s = 1, by taking a B–spline basis (see, e.g., [BHS93]) for Θ ordered θ_{1} ≤ θ_{2} ≤ ⋯ ≤ θ_{n+1} (as they appear on the line) it follows

#{Φ ∈ A_{i}(Θ) : Φ is allowable} = #{(θ_{j}, …, θ_{j+i−1}) : θ_{j} < θ_{j+i−1}, j = 1, …, n + 2 − i}
= dim(S_{i}(Θ)), i = 2, …, n + 1,

as required. For s = 0, Θ is a single point repeated, and (6.7) holds because the one submultiset of each size i is allowable.

Case (b): The case of IR follows from (a), take d = 1, and s = 0, 1. For IR^{2}, it is necessary to show that (6.7) holds for j − 1 = 0, 1, which again follows from (a).

Here are some examples of the binomial identities given by Theorem 6.5.

Example 1 (Points in general position). If the points are in general position, then every submultiset of size i up to d is allowable, with dim(Φ) = i − 1. Those of size d + 1 or larger have dim(Φ) = d, and so contribute zero to the sum. Hence

\binom{n + d}{d} = ∑_{i=1}^{d} \binom{n + 1}{i}\binom{i − 1 + (d − i + 1) − 1}{(d − i + 1) − 1} = ∑_{i=1}^{d} \binom{n + 1}{i}\binom{d − 1}{d − i},  (6.8)

which is the Chu–Vandermonde identity.

Example 2 (All points coincident). If all the points coincide, then there is exactly one allowable submultiset of size i, giving the binomial identity

\binom{n + d}{d} = ∑_{i=1}^{n+1} \binom{i + d − 2}{d − 1}.  (6.9)

Example 3 (All points on a line in IR^{d}, with no repetitions). If there are n + 1 points in IR^{d} all distinct and lying on a line, then there are exactly n + 2 − i allowable subsets of size i (consecutive points), which leads to

\binom{n + d}{d} = (n + 1) + ∑_{i=2}^{n+1}(n + 2 − i)\binom{i + d − 3}{d − 2},  (6.10)

and

\binom{n + d}{d} = ∑_{i=1}^{n+1}(n + 2 − i)\binom{i + d − 3}{d − 2},  d ≥ 2.  (6.11)

Example 4 (All points on a line in IR^{d}, with repetitions). Suppose that the points in Θ lie on a line in IR^{d}, say in the order x_{1}, …, x_{r}, occurring with multiplicities m_{1}, …, m_{r}, where m_{1} + ⋯ + m_{r} = n + 1. Then the allowable submultisets are either one point repeated a number of times, or start with a point x_{i} (repeated 1, …, m_{i} times) contain all points x_{i+1}, x_{i+2}, …, x_{j−1} with full multiplicity and the point x_{j} (repeated 1, …, m_{j} times). Using (6.9) to count the submultisets of one point repeated, we obtain the identity

\begin{aligned}\binom{n + d}{d} = & ∑_{i=1}^{r} \binom{m_{i} − 1 + d}{d} \\ & + ∑_{1≤i<j≤r} ∑_{k_{1}=1}^{m_{i}} ∑_{k_{2}=1}^{m_{j}} \binom{k_{1} + m_{i+1} + m_{i+2} + ⋯ + m_{j−1} + k_{2} − 1 + d − 2}{d − 2}. \end{aligned}

(6.12)

When all the points occur with multiplicity 1 (Example 3), this reduces to

\binom{n + d}{d} = (n + 1) + ∑_{1≤i<j≤n+1} \binom{j − i + d − 2}{d − 2}.  (6.13)

For points in IR^{2}, we get

\binom{n + 2}{2} = \binom{m_{1} + ⋯ + m_{r} + 1}{2} = ∑_{i=1}^{r} \binom{m_{i} + 1}{2} + ∑_{1≤i<j≤r} m_{i}m_{j},  (6.14)

and for two points occuring with multiplicities m and n + 1 − m respectively, we get

\binom{n + d}{d} = \binom{m − 1 + d}{d} + \binom{n − m + d}{d} + ∑_{i=1}^{m} ∑_{j=1}^{n+1−m} \binom{i + j − 1 + d − 2}{d − 2}. (6.15)

Example 5 (Points in general position in IR^{d}, each repeated m times). Suppose that Θ consists of r points in general position in IR^{d}, each repeated m times. Then the allowable submultisets which contribute to the sum consist of up to d points, each repeated with some multiplicity (≤ m). We count these as follows. The i_{1} = 1, …, d distinct points can be chosen in \binom{r}{i_{1}} ways. Of these, i_{2} of them occur with multiplicity at least two, this selection can be done in \binom{i_{1}}{i_{2}} ways, and i_{3} of these i_{2} points occur with multiplicity three or more, and so forth, giving

\begin{matrix}( & {mr − 1 + d \atop d} & ) \\ = ∑_{i_{1}=1}^{d} \binom{r}{i_{1}} ∑_{i_{2}=0}^{i_{1}} \binom{i_{1}}{i_{2}} ∑_{i_{3}=0}^{i_{2}} \binom{i_{2}}{i_{3}} ⋯ ∑_{i_{m}=0}^{i_{m−1}} \binom{i_{m−1}}{i_{m}}\binom{d + i_{2} + i_{3} + … + i_{m} − 1}{d − i_{1}}. \end{matrix} (6.16)

Example 6 (Points in general position in IR^{d}, with repetitions). Suppose that Θ consists of r points in general position in IR^{d}, repeated with multiplicities m_{1}, …, m_{r}, where m_{1} +⋯+m_{r} = n+1. Then the allowable submultisets which contribute to the sum consist of up to d points, each repeated with some multiplicity. This gives the identity

\binom{n + d}{d} = ∑_{j=1}^{d} ∑_{1≤i_{1}<i_{2}<···<i_{j}≤r} ∑_{k_{1}=1}^{m_{i_{1}}} ⋯ ∑_{k_{j}=1}^{m_{i_{j}}} \binom{k_{1} + ⋯ + k_{j} − 1 + d − j}{d − j}.  (6.17)

In the case of the previous example (m_{1}, …, m_{r} = m), this gives

\binom{mr − 1 + d}{d} = ∑_{j=1}^{d} \binom{r}{j} ∑_{k_{1}=1}^{m} ⋯ ∑_{k_{j}=1}^{m} \binom{k_{1} + ⋯ + k_{j} − 1 + d − j}{d − j},  (6.18)

which is (1.3). Compare this with (6.16), which was obtained by counting the allowable submultisets in a different way.

Example 7 (Points in general position in an affine subspace of IR^{d}). Suppose that Θ consists of n + 1 > s (distinct) points which are in general position in an s–dimensional affine subspace of IR^{d}. Then the submultisets Ψ which contribute to the summation are all those of size j up to s, and Ψ = Θ. This gives

\begin{matrix}\binom{n + d}{d} = ∑_{j=1}^{s} \binom{n + 1}{j}\binom{d − 1}{d − j} + ∑_{i=s+1}^{n+1} \binom{i + d − s − 2}{d − s − 1} & \binom{n + 1 − i + s}{s}, \\ n + 1 > s, 1 ≤ s ≤ d. \end{matrix} (6.19)

In the special cases Example 1 and 3 (s = d and s = 1), this reduces to (6.8) and (6.10).

Example 8 (Different allowable submultisets). For d ≥ 3, the first example where (6.6) does not hold is 4 distinct points in IR^{3} which lie in a plane, of which no three are collinear. In this case \binom{n+d}{d} = 20. All submultisets of sizes 1, 2 are allowable, contributing 4 × 1 + 6 × 2 = 16 to the sum. If none of the points is in the convex hull of the other three (take the vertices of a square for example), then there are 4 allowable submultisets of size 3 and 1 allowable submultiset of size 4, giving total sum of 16+4×1+1×1 ≠ 20. It is possible to modify the definition of ‘allowable’ to include such configurations. For example, from the proof of Theorem 6.5, it follows one can take any definition of ‘allowable’ for which (6.7) holds.

7. Kergin–Fejèr interpolation in IR^{d}

The Hermite–Fejèr interpolation operator at distinct points x_{1}, x_{2}, …, x_{r} in IR is a modified form of Hermite interpolation at these points repeated with multiplicity 2, where the function values at x_{1}, …, x_{r} are matched, but the remaining reduced interpolation conditions from (5.18) are taken to be zero (i.e., the derivative of the interpolant is zero at the points x_{1}, …, x_{n}). Let K_{Θ} := H_{Θ}^{(0)} be the Kergin interpolation operator.

For any Θ, the linear projector K_{Θ} : C^{n}(IR^{d}) → Π_{n}(IR^{d}) can be written as

K_{Θ} = ⨁_{j=1}^{d} ⨁_{{Ψ⊂Θ \atop dim(Ψ)=j−1} aff(Ψ)∩^{∗}Θ=Ψ} ⨁_{i=1}^{#Ψ} K_{Θ,Ψ,i},  K_{Θ,Ψ,i} : C^{n}(IR^{d}) → Π_{Θ,Ψ,i},  (7.1)

where K_{Θ,Ψ,i} is the linear projector defined by the interpolation conditions Λ_{Ψ,i} := Λ_{Ψ,i}^{(0)} of (5.7), and the space of interpolants

Π_{Θ,Ψ,i} := {p ∈ Π_{n}(IR^{d}) : p ∈ ker(Λ_{Φ,k}), ∀(Φ, k) ≠ (Ψ, i)},

This ‘direct sum’ satisfies

K_{Θ,Ψ_{1},i_{1}}K_{Θ,Ψ_{2},i_{2}} = 0, (Ψ_{1}, i_{1}) ≠ (Ψ_{2}, i_{2}).

The polynomial spaces Π_{Θ,Ψ,i} are dual to the Λ_{Ψ,i}, and the dual of (5.6) is

Π_{n}(IR^{d}) = ⨁_{j=1}^{d} ⨁_{{Ψ⊂Θ \atop dim(Ψ)=j−1} aff(Ψ)∩^{∗}Θ=Ψ} ⨁_{i=1}^{#Ψ} Π_{Θ,Ψ,i}.  (7.2)

In the above terms, the Hermite–Fejèr interpolation operator is

F_{Θ} := ⨁_{k=1}^{r} K_{Θ,[x_{k},x_{k}],1},  Θ := {x_{1}, x_{1}, x_{2}, x_{2}, …, x_{r}, x_{r}}.  (7.3)

By analogy, we define Kergin–Fejèr interpolation in IR^{d} by (7.3), when x_{1}, …, x_{r} are distinct points in IR^{d}. This operator has a continuous extension to functions defined at the points x_{1}, …, x_{r}, and maps into an r–dimensional subspace of Π_{2r−1}(IR^{d}). In order to compute this operator it is necessary to find the reduced interpolation conditions for Kergin interpolation at Θ, which we now do for points in general position in IR^{2} and IR^{3}.

Example 1 (Kergin–Fejèr interpolation in IR^{2}). The Ψ in the direct sum (7.1) have affine dimension 0 or 1. Those of dimension 0 are of the form Ψ = [x_{k}, x_{k}], and give the matching of function values and first order derivatives at x_{k}, i.e., f(x_{k}) and ∇f(x_{k}). Those of dimension 1 (j = 2) are of the form Ψ = [x_{s}, x_{s}, x_{t}, x_{t}], with

A_{2}(Ψ) = {[x_{s}, x_{t}]}, A_{3}(Ψ) = {[x_{s}, x_{s}, x_{t}], [x_{s}, x_{t}, x_{t}]}, A_{4}(Ψ) = {[x_{s}, x_{s}, x_{t}, x_{t}]}.

With (x_{t} − x_{s})^{⊥} retaining its previous definition, we can now write down the reduced interpolation conditions (as applied to f)

f(x_{k}), ∇f(x_{k}), 3n conditions,
∫_{[x_{s},x_{t}]} D_{(x_{t}−x_{s})^{⊥}}f, s ≠ t, \frac{1}{2}n(n − 1) conditions,
∫_{[x_{s},x_{s},x_{t}]} D_{(x_{t}−x_{s})^{⊥}}^{2}f, s ≠ t, n(n − 1) conditions,
∫_{[x_{s},x_{s},x_{t},x_{t}]} D_{(x_{t}−x_{s})^{⊥}}^{3}f, s ≠ t, \frac{1}{2}n(n − 1) conditions,

giving a total of dim(Π_{2r−1}(IR^{2})) = \binom{2r+1}{2} linearly independent conditions.

Example 2 (Kergin–Fejèr interpolation in IR^{3}). The Ψ in (7.1) are of dimensions 0, 1, 2, and are of the forms

Ψ = [x_{s}, x_{s}], Ψ = [x_{s}, x_{s}, x_{t}, x_{t}], Ψ = [x_{s}, x_{s}, x_{t}, x_{t}, x_{u}, x_{u}],

respectively (where s, t, u are distinct). Using D_{⊥} to indicate a directional derivative in a direction orthogonal to the affine hull of the points that occur in the integral, we can express the reduced interpolation conditions (applied to f) as

f(x_{s}), ∇f(x_{s}),  4r conditions,  } j = 1

\begin{matrix}∫_{[x_{s},x_{t}]} D_{⊥}f,  2\binom{r}{2} conditions, \\ ∫_{[x_{s},x_{s},x_{t}]} D_{⊥}^{2}f,  6\binom{r}{2} conditions, \\ ∫_{[x_{s},x_{s},x_{t},x_{t}]} D_{⊥}^{3}f,  4\binom{r}{2} conditions, \end{matrix}  \} j = 2

\begin{matrix}∫_{[x_{s},x_{t},x_{u}]} D_{⊥}^{2}f,  \binom{r}{3} conditions, \\ ∫_{[x_{s},x_{s},x_{t},x_{u}]} D_{⊥}^{3}f,  3\binom{r}{3} conditions, \\ ∫_{[x_{s},x_{s},x_{t},x_{t},x_{u}]} D_{⊥}^{4}f,  3\binom{r}{3} conditions, \\ ∫_{[x_{s},x_{s},x_{t},x_{t},x_{u},x_{u}]} D_{⊥}^{5}f,  \binom{r}{3} conditions, \end{matrix} \} j = 3

The author has not investigated to what extent the convergence properties (for large numbers of nodes) of Kergin–Fejèr interpolation in IR^{2} and IR^{3} are similar those of the well understood case of Hermite–Fejèr interpolation in IR. Any results in this direction would be of interest.

8. Computational considerations

In this final section, we outline how the subspaces Λ_{Ψ,i}, their duals Π_{Θ,Ψ,i}, and the related polynomial interpolants K_{Θ,Ψ,i}f of (7.1) can be computed numerically. Similar considerations apply to the other mean value interpolation operators.

Identification of the Ψ

The first step is to identify the submultisets Ψ in the direct sum decomposition (3.2), if they have not already been given.

In exact arithmetic this can be done simply by, at the j–th stage, taking the Ψ of dimension j − 1 and adding new points to them (as many as possible) to get the Ψ of dimensional j. These calculations should be arranged efficiently. First determine the Ψ of dimension 0, i.e., find the multiplicities of each of the distinct points. At the j–th stage, take each Ψ of dimension j − 1 and add to it a new point θ (up to full multiplicity) which gives a j–dimensional affine subspace different from one already found at this stage (i.e., the j–dimensional multiset Ψ ∪ {θ, …, θ} is not contained some Ψ̃ of dimension j already obtained), then add to Ψ∪{θ, …, θ} as many of the remaining points as is possible without changing its dimension, to obtain a new Ψ̃ of dimension j. Stop as soon as all the Ψ have been found, i.e., when j = d or Ψ = Θ.

If computations are only done to machine precision, then some tolerance, say tol, should be used when deciding whether or not a point lies in the affine hull of some set of points. The issues involved here are perhaps best illustrated by an example. Consider two points Θ := {θ_{1}, θ_{2}} in the plane which are close, but not equal. Though the map (θ_{1}, θ_{2}) ↦ K_{[θ_{1},θ_{2}]}f is continuous for f ∈ C^{1}(IR^{2}), the computation of K_{[θ_{1},θ_{2}]}f from the reduced interpolation conditions suffers a loss of significance as θ_{1}, θ_{2} → θ. Here the Λ_{Ψ,i} are 1-dimensional

Λ_{[θ_{i}],1} = span{f ↦ f(θ_{i})}, i = 1, 2, Λ_{[θ_{1},θ_{2}],2} = span{f ↦ ∫_{[θ_{1},θ_{2}]} D_{(θ_{2}−θ_{1})^{⊥}}f}, with corresponding dual polynomial spaces

Π_{Θ,[θ_{1}],1} = span{(θ_{2} − θ_{1})^{T} (· − θ_{2})}, Π_{Θ,[θ_{2}],1} = span{(θ_{2} − θ_{1})^{T} (· − θ_{1})},

Π_{Θ,[θ_{1},θ_{2}],2} = span\{(θ_{2} − θ_{1})^{⊥T} (· − \frac{θ_{2} + θ_{1}}{2} )\}.

This gives the Kergin interpolant

\begin{aligned}K_{[θ_{1},θ_{2}]}f(x) = f(θ_{1}) & \frac{(θ_{1} − θ_{2})^{T} (x − θ_{2})}{∥θ_{2} − θ_{1}∥^{2}} + f(θ_{2})\frac{(θ_{2} − θ_{1})^{T} (x − θ_{1})}{∥θ_{2} − θ_{1}∥^{2}} \\ & + ∫_{[θ_{1},θ_{2}]} D_{(θ_{2}−θ_{1})^{⊥}}f \frac{(θ_{2} − θ_{1})^{⊥T} (2x − (θ_{1} + θ_{2}))}{2∥(θ_{2} − θ_{1})^{⊥}∥^{2}}. \end{aligned}  (8.1)

There is a loss of significance in the computed values of each of the terms in (8.1) when θ_{1}, θ_{2} are close to θ. Here the calculated value of K_{[θ_{1},θ_{2}]}f should be close to K_{[θ,θ]}f (the tangent plane to f at θ), which can be computed from f(θ), ∇f(θ) without loss of significance. This behaviour accurately reflects the general situation.

In general, any points θ ∈ Θ at a distance less than tol from the affine hull of a Ψ being computed should be included in Ψ (without increasing its affine dimension), and the Kergin interpolant calculated from the corresponding reduced interpolation conditions. To do this, translate the points in Ψ (and θ) so that the affine subspace they span becomes a linear subspace (subtract one of the points, or better their average, from each of them). Form the matrix A = [Ψ, θ] with columns given by the points in Ψ and θ, and find a QR– factorisation A = QR. The absolute value of the last diagonal entry of R is the distance of θ from the affine subspace aff(Ψ), and if this is < tol, then θ should be included in Ψ. If the translation of subtracting to first point in Ψ is adopted, then the Q and R can simply be updated one column at a time as each point is added to Ψ. Alternatively, one could use a standard package such as MATLAB’s rank(A,tol) to compare the numerical rank of A = [Ψ] and A = [Ψ, θ].

Finding the interpolant K_{Θ}f, Λ_{Ψ,i} and Π_{Θ,Ψ,i}

Now suppose the Ψ in the direct sum have been obtained, either exactly or numerically (as just discussed). To find the Kergin interpolant K_{Θ}f, or the interpolants K_{Θ,Ψ,i}f (where some of the interpolation conditions are zero) it is necessary to solve a linear system for the coefficients {a_{k}} of the interpolant p = ∑_{k} a_{k}p_{k} relative to some choice of basis {p_{k}} for Π_{n−1}(IR^{d}), i.e.,

∑_{k} a_{k}λ(p_{k}) = λ(p) = λ(f),  (8.2)

for λ from some spanning set for the reduced interpolation conditions. One could take {p_{k}} to be the basis of shifted monomials with centre the average of the points in Θ. Choosing a spanning set of λ ∈ Λ_{Ψ,i}, as defined by (5.7), can be done as follows. Let N := {u_{1}, u_{2}, …, u_{d−j+1}} be an orthonormal basis for aff(Ψ)^{⊥}, and choose the p(D) in (5.12) to be of the form p(D) = D_{y_{1}}D_{y_{2}} ⋯ D_{y_{i−1}}, where y_{1}, …, y_{i−1} ∈ N. For each such choice, taking s from a basis for S_{i}(Ψ) leads to a basis of λ for Λ_{Ψ,i}. In the absence of a known basis for S_{i}(Ψ), it is natural to take s from the spanning set {M(·|Φ) : Φ ∈ A_{i}(Θ)}, which in general leads to a system (8.2) of more equations than unknowns. If solved numerically, this system should be solved using the least squares method. Even if a basis for S_{i}(Ψ) is available, as is outlined in some special cases in Example 16, the author suggests finding the least squares solution of the (possibly) larger system of equations (which is symmetric in the points Ψ, and likely to give a more accurate answer).

The dual polynomial spaces Π_{Θ,Ψ,i} can be found as follows. Order the Λ_{Ψ,i} as Λ_{1}, Λ_{2}, …, and suppose that a spanning set λ_{1}^{(k)}, … λ_{m_{k}}^{(k)} is given for each Λ_{k}, as above. Let Λ : Π_{n−1}(IR^{d}) → IR^{m}, m := ∑_{k} m_{k}, be the data map given by

Λ(f) := [ Λ_{1}(f) ; Λ_{2}(f) ; ⋮ ] , Λ_{k}(f) := [ λ_{1}^{(k)}(f) ; ⋮ ; λ_{m_{k}}^{(k)}(f) ] ,

and P : IR^{ℓ} → Π_{n−1}(IR^{d}), ℓ := \binom{n−1+d}{d}, be the basis map (a_{k}) ↦ ∑_{k} a_{k}p_{k}. Apply elementary column operations to the m × ℓ matrix A := ΛP to obtain a block diagonal form

AM = (ΛP)M = [ \begin{matrix}C_{1} & 0 & ⋯ & 0 \\ 0 & C_{2} & ⋯ & 0 \\ ⋮ & ⋮ & & ⋮ \\ 0 & 0 & ⋯ & C_{n} \end{matrix} ] , M ∈ IR^{ℓ×ℓ} invertible,

where the blocks C_{k} are of size m_{k} × ℓ_{k}, ℓ_{k} := dim(Λ_{k}), and have linearly independent columns. This can be done by using Gauss elimination (with pivoting) to find a column echelon form for A, or using Gauss elimination by segments applied to the ‘segmented matrix’

A^{T} = [(Λ_{1}P)^{T} , (Λ_{2}P)^{T} , …] ∈ IR^{ℓ×m}

(as detailed in de Boor [B94]). The sequence of ℓ_{k} polynomials

P (Me_{i}), k = ℓ_{1} + ⋯ + ℓ_{k−1} + 1, …, ℓ_{1} + ⋯ + ℓ_{k}

is a basis for Λ_{k} (above Me_{i} is the i–th column of the matrix M).

A typical use of the dual polynomial spaces is in the computation of the Kergin–Fejèr interpolant F_{Θ}(f) from the point evaluations of f at {x_{1}, …, x_{r}}, which requires a basis for the r–dimensional space

⨁_{k=1}^{r} Π_{Θ,[x_{k},x_{k}],1} ⊂ Π_{2r−1}(IR^{d}), Θ := {x_{1}, x_{1}, …, x_{r}, x_{r}}.

Code for the computations outlined in this section can be obtained from the author.

Acknowledgement

Thanks to Thomas Grandine (Seattle) for providing the data for the plots of simplex splines in Figures 5.1 and 5.2 from his implementation of [Gr86].

References

[BHS93] B. D. Bojanov, H. A. Hakopian, and A. A. Sahakian, “Spline functions and multivari-

ate interpolations”, Kluwer Academic Publishers, , 1993.

[B94] C. de Boor, Gauss elimination by segments and multivariate polynomial interpolation, in “Approximation and Computation: A Festschrift in Honor of Walter Gautschi” (R. V. M. Zahar Ed.), pp. 1–22, ISNM 119, Birkhäuser Verlag, Basel-Boston-Berlin, 1994.

[BC97] L. Bos and J.-P. Calvi, Kergin interpolants at the roots of unity approximate C^{2} functions, J. Anal. Math. 72 (1997), 203–221.

[BW01] L. Bos and S. Waldron, On the structure of Kergin interpolation for points in general position, To appear, 2001.

[BS94] S. C. Brenner and L. R. Scott, “The mathematical theory of finite element methods”,

Springer-Verlag, New York, 1994.

[CGMS83] A. S. Cavaretta, T. N. T. Goodman, C. A. Micchelli, and A. Sharma, Multivariate interpolation and the Radon transform, part III: Lagrange representation, Canad. Math. Soc. Conf. Proceedings 3 (1983), 37–50.

[C78] P. G. Ciarlet, “The finite element method for elliptic problems”, North Holland, ,

1978.

[DM83] W. A. Dahmen and C. A. Micchelli, On the linear independence of multivariate B-

splines II: complete configurations, Math. Comp. 41(163) (1983), 143–163.

[G83] T. N. T. Goodman, Interpolation in minimum semi-norm, and multivariate B-splines, J. Approx. Theory 37 (1983), 212–223.

[Gr86] T. A. Grandine, The stable evaluation of multivariate simplex splines, Math. Comp. 50 (1988), 197–205.

[H81] H. Hakopian, Les differences divisées de plusieurs variables et les interpolations multidimensionnelles de types Lagrangien et Hermitien, C. R. Acad. Sci. Paris Ser. I 292 (1981), 453-456.

[H84] Hakop Hakopian, On multivariate spline functions, B-spline bases and polynomial interpolation II, Studia Math. 79 (1984), 91–102.

[K80] P. Kergin, A natural interpolation of C^{k} functions, J. Approx. Theory 29 (1980), 278–293.

[L92] R. A. Lorentz, “Multivariate Birkhoff Interpolation”, Lecture Notes in Mathematics, No. 1516, Springer Verlag, Heidelberg, 1992.

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

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

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

variate Lagrange interpolation schemes, SIAM J. Math. Anal. 28(1) (1997), 233–258.