<<

. 7
( 83 .)



>>

But this implies that the homogeneous leading term of the “fundamental” polynomials det Gk coin-
k
cides, after this particular choice of gi; j , with
1 @Fj
g= det : i; j = 1; : : : ; d ;
deg f · · · deg f @i
1 d

which is independent of k now; in other words, there exist polynomials gk ; k = 1; : : : ; N , such that
ˆ
deg gk ¡ deg g and det Gk = g + gk . Moreover, g is a homogeneous polynomial of degree at most
ˆ ˆ
deg f + · · · + deg f ’ d. Now, let p be any polynomial, then
1 d
N N N
det Gj (·) p(zj ) p(zj )
Kp = p(zj ) =g + g:
ˆ (8)
det Gj (zj ) j
det Gj (zj ) det Gj (zj )
j=1 j=1 j=1

Combining (8) with (6) then yields the existence of polynomials q1 ; : : : ; qd such that
N N d
p(zj ) p(zj )
p=g + g+
ˆ qj fj
det Gj (zj ) j
det Gj (zj )
j=1 j=1 j=1

and comparing homogeneous terms of degree deg g Kronecker realized that either, for any p such
that deg p ¡ deg g,
N
p(zj )
=0 (9)
det Gj (zj )
j=1

or there exist homogeneous polynomials h1 ; : : : ; hd such that
d
g= hj det F : (10)
j
j=1

The latter case, Eq. (10), says (in algebraic terminology) that there is a syzygy among the leading
terms of the polynomials F ; j = 1; : : : ; d, and is equivalent to the fact that N ¡ deg f · · · deg f ,
j 1 d
while (9) describes and even characterizes the complete intersection case that N = deg f · · · deg f .
1 d
In his paper, Kronecker also mentions that the condition (10) has been overlooked in [21]. Jacobi
dealt there with the common zeros of two bivariate polynomials and derived explicit representations
for the functional
N
f(zj )
[z1 ; : : : ; zN ]f:= ; (11)
det Gj (zj )
j=1
M. Gasca, T. Sauer / Journal of Computational and Applied Mathematics 122 (2000) 23“35 27


which behaves very much like a divided di erence, since it is a combination of point evaluations
d
which, provided that (9) hold true, annihilates deg g’1 .
In addition, Kronecker refers to a paper [6] which he says treats the case of symmetric functions,
probably elementary symmetric polynomials. Unfortunately, this paper is unavailable to us so far.


3. Bivariate tables, the natural approach

Only very few research papers on multivariate polynomial interpolation were published during the
ÿrst part of this century. In the classical book Interpolation [45], where one section (Section 19) is
devoted to this topic, the author only refers to two related papers, recent at that time (1927), namely
[27,28]. The latter one [28], turned out to be inaccessible to us, unfortunately, but it is not di cult
to guess that it might have pursued a tensor product approach, because this is the unique point of
view of [45] (see also [31]).
The formulas given in [27] are Newton formulas for tensor product interpolation in two variables,
and the author, Narumi, claims (correctly) that they can be extended to “many variables”. Since it is
a tensor product approach, the interpolation points are of the form (xi ; yj ), 06i6m; 06j6n, with
xi ; yj arbitrarily distributed on the axes OX and OY , respectively. Bivariate divided di erences for
these sets of points are obtained in [27], by recurrence, separately for each variable. With the usual
notations, the interpolation formula from [27] reads as
j’1
m n i’1
p(x; y) = f[x0 ; : : : ; xi ; y0 ; : : : ; yj ] (x ’ xh ) (y ’ xk ); (12)
i=0 j=0 h=0 k=0

where empty products have the value 1. Remainder formulas based on the mean value theorem are
also derived recursively from the corresponding univariate error formulas in [27]. For f su ciently
smooth there exist values ; ; Á; Á such that
@m+1 f( ; y) m (x ’ xh ) @n+1 f(x; Á) n (y ’ yk )
h=0 k=0
R(x; y) = +
@x (m + 1)! @y (n + 1)!
m+1 n+1


@m+n+2 f( ; Á ) m (x ’ xh ) n (y ’ yk )
h=0 k=0
’ : (13)
@x m+1 @y n+1 (m + 1)! (n + 1)!
The special case of equidistant points on both axes is particularly considered in [27], and since
the most popular formulas at that time were based on ÿnite di erences with equally spaced argu-
ments, Narumi shows how to extend Gauss, Bessel and Stirling univariate interpolation formulas for
equidistant points to the bivariate case by tensor product. He also applies the formulas he obtained
to approximate the values of bivariate functions, but he also mentions that some of his formulas had
been already used in [49].
In [45], the Newton formula (12) is obtained in the same way, with the corresponding remainder
formula (13). Moreover, Ste ensen considers a more general case, namely when for each i; 06i6m,
the interpolation points are of the form y0 ; : : : ; yni , with 06ni 6n. Now with a similar argument the
interpolating polynomial becomes
j’1
ni
m i’1
p(x; y) = f[x0 ; : : : ; xi ; y0 ; : : : ; yj ] (x ’ xh ) (y ’ xk ) (14)
i=0 j=0 h=0 k=0
28 M. Gasca, T. Sauer / Journal of Computational and Applied Mathematics 122 (2000) 23“35


with a slightly more complicated remainder formula. The most interesting particular cases occur
when ni = n, which is the Cartesian product considered above, and when ni = m ’ i. This triangular
case (triangular not because of the geometrical distribution of the interpolation points, but of the
indices (i; j)), gives rise to the interpolating polynomial
j’1
m m’i i’1
p(x; y) = f[x0 ; : : : ; xi ; y0 ; : : : ; yj ] (x ’ xh ) (y ’ xk ); (15)
i=0 j=0 h=0 k=0

that is
j’1
i’1
p(x; y) = f[x0 ; : : : ; xi ; y0 ; : : : ; yj ] (x ’ xh ) (y ’ xk ): (16)
06i+j6m h=0 k=0

Ste ensen refers for this formula to Biermann™s lecture notes [4] from 1905, and actually it seems
that Biermann has been the ÿrst who considered polynomial interpolation on the triangular grid in
a paper [3] from 1903 (cf. [44]) in the context of cubature.
Since the triangular case corresponds to looking at the “lower triangle” of the tensor product
situation only, this case can be resolved by tensor product methods. In particular, the respective
error formula can be written as
m+1 i’1 m’i
@m+1 f( i ; Ái ) h=0 (x ’ xh ) k=0 (y’ yk )
R(x; y) = : (17)
@xi @ym+1’i i! (m ’ i + 1)!
i=0

In the case of Cartesian product Ste ensen also provides the Lagrange formula for (12), which
can be obviously obtained by tensor product of univariate formulas.
Remainder formulas based on intermediate points ( i ; Ái ) can be written in many di erent forms.
For them we refer to Stancu™s paper [44] which also contains a brief historical introduction where the
author refers, among others, to [3,15,27,40,41]. Multivariate remainder formulas with Peano (spline)
kernel representation, however, have not been derived until very recently in [42] and, in particular,
in [43] which treats the triangular situation.


4. Salzer™s papers: from bivariate tables to general sets

In 1944, Salzer [33] considered the interpolation problem at points of the form (x1 + s1 h1 ; : : : ; x n +
sn hn ) where
(i) (x1 ; : : : ; x n ) is a given point in Rn ,
(ii) h1 ; : : : ; hn are given real numbers,
(iii) s1 ; : : : ; sn are nonnegative integers summing up to m.
This is the multivariate extension of the triangular case (16) for equally spaced arguments, where
ÿnite di erences can be used. Often, di erent names are used for the classical Newton interpolation
formula in the case of equally spaced arguments using forward di erences: Newton“Gregory, Harriot“
Briggs, also known by Mercator and Leibnitz, etc. See [18] for a nice discussion of this issue. In
[33], Salzer takes the natural multivariate extension of this formula considering the polynomial
q(t1 ; : : : ; tn ) := p(x1 + t1 h1 ; : : : ; x n + tn hn ) of total degree not greater than m in the variables t1 ; : : : ; tn ,
which interpolates a function f(x1 +t1 h1 ; : : : ; x n +tn hn ) at the points corresponding to ti =si , i=1; : : : ; n,
M. Gasca, T. Sauer / Journal of Computational and Applied Mathematics 122 (2000) 23“35 29


where the si are all nonnegative integers such that 06s1 + · · · + sn 6m. The formula, which is called
in [33] a multiple Gregory“Newton formula, is rewritten there in terms of the values of the function
f at the interpolation points, i.e., in the form
t1 tn m ’ t1 ’ · · · ’ tn
q(t1 ; : : : ; tn ) = ··· f(x1 + s1 h1 ; : : : ; x n + sn hn ): (18)
s1 sn m ’ s1 ’ · · · ’ sn
s1 +···+sn 6m

Note that (18) is the Lagrange formula for this interpolation problem. Indeed, each function
t1 tn m ’ t1 ’ · · · ’ tn
··· (19)
s1 sn m ’ s1 ’ · · · ’ sn
is a polynomial in t1 ; : : : ; tn of total degree m which vanishes at all points (t1 ; : : : ; tn ) with ti non-
negative integers 06t1 + · · · + tn 6m, except at the point (s1 ; : : : ; sn ), where it takes the value 1. In
particular, for n = 1 we get the well-known univariate Lagrange polynomials
t’i
t m’t
˜s (t) = =
s m’s s’i
06i6m;
i=s


for s = 0; : : : ; m.
Salzer used these results in [34] to compute tables for the polynomials (18) and, some years later
in [35], he studied in a similar form how to get the Lagrange formula for the more general case of
formula (16), even starting with this formula. He obtained the multivariate Lagrange polynomials
by a rather complicated expression involving the univariate ones.
It should be noted that several books related to computations and numerical methods published
around this time include parts on multivariate interpolation to some extent, surprisingly, more than
most of the recent textbooks in Numerical Analysis. We have already mentioned Ste ensen™s book
[45], but we should also mention Whittaker and Robinson [51, pp. 371“374], Mikeladze [26, Chapter
XVII] and especially Kunz [23, pp. 248“274], but also Isaacson and Keller [20, pp. 294 “299] and
Berezin and Zhidkov [2, pp. 156 “194], although in any of them not really much more than in [45]
is told.
In [36,37], Salzer introduced a concept of bivariate divided di erences abandoning the idea of
iteration for each variable x and y taken separately. Apparently, this was the ÿrst time (in spite
of the similarity with (11)), that bivariate divided di erences were explicitly deÿned for irregularly
distributed sets of points. Divided di erences with repeated arguments are also considered in [37] by
coalescence of the ones with di erent arguments. Since [36] was just a ÿrst attempt of [37], we only
explain the latter one. Salzer considers the set of monomials {xi yj }, with i; j nonnegative integers,
ordered in a graded lexical term order, that is,
(i; j) ¡ (h; k) ” i + j ¡ h + k or i + j = h + k; i ¿ h: (20)
Hence, the monomials are listed as
{1; x; y; x2 ; xy; y2 ; x3 ; : : :}: (21)
For any set of n + 1 points (xi ; yi ), Salzer deÿnes the associated divided di erence
n
[01 : : : n]f:= Ak f(xk ; yk ); (22)
k=0
30 M. Gasca, T. Sauer / Journal of Computational and Applied Mathematics 122 (2000) 23“35


choosing the coe cients Ak in such a form that (22) vanishes when f is any of the ÿrst n monomials
of list (21) and takes the value 1 when f is the (n + 1)st monomial of that list. In other words, the
coe cients Ak are the solution of the linear system
n
ij
xi yj any of the ÿrst n monomials of (21);
Ak xk yk = 0;
k=0
n
ij
xi yj the (n + 1)th monomial of (21):
Ak xk yk = 1; (23)
k=0

These generalized divided di erences share some of the properties of the univariate ones but not
all. Moreover, they have some limitations, for example, they exist only if the determinant of the
coe cients in (23) is di erent from zero, and one has no control of that property in advance. On
the other hand, observe that for example the simple divided di erence with two arguments (x0 ; y0 )
and (x; y), which is
f(x; y) ’ f(x0 ; y0 )
;
x ’ x0
gives, when applied to f(x; y) = xy, the rational function
xy ’ x0 y0
x ’ x0
and not a polynomial of lower degree. In fact, Salzer™s divided di erences did not have great success.
Several other deÿnitions of multivariate divided di erences had appeared since then, trying to keep
as many as possible of the good properties of univariate divided di erences, cf. [16].


5. Reduction of a problem to other simpler ones

Around the 1950s an important change of paradigm happened in multivariate polynomial inter-
polation, as several people began to investigate more general distributions of points, and not only
(special) subsets of Cartesian products. So, when studying cubature formulae [32], Radon observed
the following in 1948: if a bivariate interpolation problem with respect to a set T ‚ R2 of ( k+2 )
2
2
interpolation points is unisolvent in k , and U is a set of k + 2 points on an arbitrary straight
line ˜ ‚ R2 such that ˜ © T = …, then the interpolation problem with respect to T ∪ U is unisolvent
2
in k+1 . Radon made use of this observation to build up point sets which give rise to unisolvent
interpolation problems for m recursively by degree. Clearly, these interpolation points immediately
yield interpolatory cubature formulae.
The well-known BÃ zout theorem, cf. [50], states that two planar algebraic curves of degree m and
e

<<

. 7
( 83 .)



>>