. 8
( 83 .)


n, with no common component, intersect each other at exactly mn points in an algebraic closure
of the underlying ÿeld, counting multiplicities. This theorem has many interesting consequences
for bivariate interpolation problems, extensible to higher dimensions. For example, no unisolvent
interpolation problem in n can have more than n + 1 collinear points. Radon™s method in [32] is
a consequence of this type of observations, and some other more recent results of di erent authors
can also be deduced in a similar form, as we shall see later.
M. Gasca, T. Sauer / Journal of Computational and Applied Mathematics 122 (2000) 23“35 31

Another example of a result which shows the more general point of view taken in multivariate
interpolation at that time is due to Thacher Jr. and Milne [47] (see also [48]). Consider two uni-
variate interpolation problems in n’1 , with T1 ; T2 as respective sets of interpolation points, both of
cardinality n. Assume that T1 © T2 has cardinality n ’ 1, hence T = T1 ∪ T2 has cardinality n + 1. The
univariate Aitken“Neville interpolation formula combines the solutions of the two smaller problems
based on T1 and T2 to obtain the solution in n of the interpolation problem with T as the under-
lying set of interpolation points. The main idea is to ÿnd a partition of unity, in this case a ne
polynomials ˜1 ; ˜2 , i.e., ˜1 + ˜2 = 1, such that
˜1 (T2 \T1 ) = ˜2 (T1 \T2 ) = 0
and then combine the solutions p1 ; p2 with respect to T1 ; T2 , into the solution ˜1 p1 + ˜2 p2 with
respect to T . This method was developed in the 1930s independently by Aitken and Neville with
the goal to avoid the explicit use of divided di erences in the computation of univariate Lagrange
polynomial interpolants.
It was exactly this idea which Thatcher and Milne extended to the multivariate case in [47].
Let us sketch their approach in the bivariate case. For example, consider an interpolation problem
with 10 interpolation points, namely, the set T = {(i; j) : 06i + j63}, where i; j are nonnegative
integers, and the interpolation space 3 . The solution pT of this problem is obtained in [47] from
the solutions pTk ∈ 2 ; k = 1; 2; 3, of the 3 interpolation problems with respect to the six-point sets
Tk ‚ T , k = 1; 2; 3, where
T1 = {(i; j) : 06i + j62};

T2 = {(i; j) : (i; j) ∈ T; i ¿ 0};

T3 = {(i; j) : (i; j) ∈ T; j ¿ 0}:
p T = ˜1 p T 1 + ˜2 p T 2 + ˜3 p T 3 ;
where ˜k ; k = 1; 2; 3 are appropriate polynomials of degree 1. In fact, in this case these polynomials
are the barycentric coordinates relative to the simplex (0, 0), (3, 0), (0, 3) and thus a partition
of unity. In [47] the problem is studied in d variables and in that case d + 1 “small” problems,
with respective interpolation sets Tk ; k = 1; : : : ; d, with a simplicial structure (the analogue of the
triangular grid), are used to obtain the solution of the full problem with T = T1 © · · · © Td+1 as
interpolation points.
In 1970, Guenter and Roetman [19], among other observations, made a very interesting remark,
which connects to the Radon=BÃ zout context and deserves to be explained here. Let us consider a set
T of ( d ) points in R , where exactly ( m+d’1 ) of these points lie on a hyperplane H . Then T \H
m+d d
m’1+d m m
consists of ( d ) points. Let us denote by d; H the space of polynomials of d with the variables
restricted to H , which is isomorphic to d’1 . If the interpolation problems deÿned by the sets T \H
m’1 m
and T ©H are unisolvent in the spaces d and d; H , respectively, then the interpolation problem
deÿned by T is unisolvent in d . In other words, the idea is to decompose, whenever possible, a
problem of degree m and d variables into two simpler problems, one of degree m and d’1 variables
and the other one with degree m ’ 1 and d variables.
32 M. Gasca, T. Sauer / Journal of Computational and Applied Mathematics 122 (2000) 23“35

6. The ÿnite element approach

In 1943, Courant [11] suggested a ÿnite di erence method applicable to boundary value problems
arising from variational problems. It is considered one of the motivations of the ÿnite element
method, which emerged from the engineering literature along the 1950s. It is a variational method
of approximation which makes use of the Rayleigh“Ritz“Galerkin technique. The method became
very successful, with hundreds of technical papers published (see, e.g., the monograph [52]), even
before its mathematical basis was completely understood at the end of the 1960s.
Involved in the process of the ÿnite element method there are local polynomial interpolation
problems, generally for polynomials of low degree, thus, with only few interpolation data. The
global solution obtained by solving all the local interpolation problems is a piecewise polynomial of
a certain regularity, depending on the amount and type of interpolation data in the common boundary
between pieces. Some of the interest in multivariate polynomial interpolation along the 1960=1970s
was due to this method. Among the most interesting mathematical papers of that time in Finite
Elements, we can mention [53,5], see also the book [46] by Strang and Fix, but, in our opinion, the
most relevant papers and book from the point of view of multivariate polynomial interpolation are
due to Ciarlet et al., for example [7“9].
In 1972, Nicolaides [29,30] put the classical problem of interpolation on a simplicial grid of ( m+d )
points of Rd , regularly distributed, forming what he called a principal lattice, into the ÿnite element
context. He actually used barycentric coordinates for the Lagrange formula, and moreover gave
the corresponding error representations, see also [7]. However, much of this material can already
be found in [3]. In general, taking into account that these results appeared under di erent titles,
in a di erent context and in journals not accessible everywhere, it is not so surprising any more,
how often the basic facts on the interpolation problem with respect to the simplicial grid had been

7. Hermite problems

The use of partial or directional derivatives as interpolation data in the multivariate case had not
received much attention prior to the ÿnite element method, where they were frequently used. It
seems natural to approach partial derivatives by coalescence, as in univariate Hermite interpolation
problems. However, things are unfortunately much more complicated in several variables. As it was
already pointed out by Salzer and Kimbro [39] in 1958, the Hermite interpolation problem based on
the values of a bivariate function f(x; y) at two distinct points (x1 ; y1 ); (x2 ; y2 ) and on the values of
the partial derivatives @f=@x; @f=@x at each of these two points is not solvable in the space 2 for
any choice of points, although the number of interpolation conditions coincides with the dimension
of the desired interpolation space. Some years later, Ahlin [1] circumvented some of these problems
by using a tensor product approach: k 2 derivatives @p+q f=@xp @yq with 06p; q6k ’ 1 are prescribed
at the n2 points of a Cartesian product. The interpolation space is the one spanned by x yÿ with
06 ; ÿ6nk ’ 1 and a formula for the solution is easily obtained.
We must mention that Salzer came back to bivariate interpolation problems with derivatives in [38]
studying hyperosculatory interpolation over Cartesian grids, that is, interpolation problems where
all partial derivatives of ÿrst and second order and the value of the function are known at the
M. Gasca, T. Sauer / Journal of Computational and Applied Mathematics 122 (2000) 23“35 33

interpolation points. Salzer gave some special conÿgurations of points which yield solvability of this
type of interpolation problem in an appropriate polynomial space and also provided the corresponding
remainder formulae.
Nowadays, Hermite and Hermite“Birkho interpolation problems have been studied much more
systematically, see [16,25] for references.

8. Other approaches

In 1966, Coatmelec [10] studied the approximation of functions of several variables by linear
operators, including interpolation operators. At the beginning of the paper, he only considered in-
terpolation operators based on values of point evaluations of the function, but later he also used
values of derivatives. In this framework he obtained some qualitative and quantitative results on the
approximation order of polynomial interpolation. At the end of [10], Coatmelec also includes some
examples in R2 of points which are distributed irregularly along lines: n + 1 of the points on a line
r0 , n of them on another line r1 , but not on r0 , and so on until 1 point is chosen on a line rn but
not on r0 ∪ · · · ∪ rn’1 . He then points out the unisolvence of the corresponding interpolation problem
in n which is, in fact, again a consequence of BÃ zout™s theorem as in [32].
In 1971, Glaeser [17] considers Lagrange interpolation in several variables from an abstract
algebraic=analytic point of view and acknowledges the inconvenience of working with particular
systems of interpolation points due to the possibility of the nonexistence of a solution, in contrast
to the univariate case. This is due to the nonexistence of polynomial spaces of dimension k ¿ 1 in
more than one variable such that the Lagrange interpolation problem has a unique solution for any
system of k interpolation points. In other words, there are no nontrivial Haar (or Chebychev) spaces
any more for two and more variables, cf. [12] or [24]. In [17], polynomial spaces with dimension
greater than the number of interpolation conditions are considered in order to overcome this problem.
Glaeser investigated these underdetermined systems which he introduced as interpolation schemes
in [17] and also studied the problem of how to particularize the a ne space of all solutions of a
given interpolation problem in order to obtain a unique solution. This selection process is done in
such a way that it controls the variation of the solution when two systems of interpolation points
are very “close” to each other, with the goal to obtain a continuous selection process.


1. We thank Carl de Boor for several references, in particular, for pointing out to us the paper [32]
with the result mentioned at the beginning of Section 4, related to Bezout™s theorem. We are also
grateful the help of Elena Ausejo, from the group of History of Sciences of the University of
Zaragoza, in the search of references.
2. M. Gasca has been partially supported by the Spanish Research Grant PB96 0730.
3. T. Sauer was supported by a DFG Heisenberg fellowship, Grant SA 627=6.
34 M. Gasca, T. Sauer / Journal of Computational and Applied Mathematics 122 (2000) 23“35


[1] A.C. Ahlin, A bivariate generalization of Hermite™s interpolation formula, Math. Comput. 18 (1964) 264“273.
[2] I.S. Berezin, N.P. Zhidkov, Computing Methods, Addison-Wesley, Reading, MA, 1965 (Russian version in 1959).
[3] O. Biermann, Uber naherungsweise Kubaturen, Monatshefte Math. Phys. 14 (1903) 211“225.
[4] O. Biermann, Vorlesungen uber Mathematische Naherungsmethoden, Vieweg, Braunschweig, 1905.
[5] G. Birkho , M.H. Schultz, R.S. Varga, Piecewise Hermite interpolation in one and two variables with applications
to partial di erential equations, Numer. Math. 11 (1968) 232“256.
[6] W. Borchardt, Uber eine Interpolationsformel fur eine Art symmetrischer Funktionen und deren Anwendung, Abh.
d. Preu . Akad. d. Wiss. (1860) 1“20.
[7] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, 1978.
P.G. Ciarlet, P.A. Raviart, General Lagrange and Hermite interpolation in Rn with applications to ÿnite element
methods, Arch. Rational Mech. Anal. 46 (1972) 178“199.
[9] P.G. Ciarlet, C. Wagschal, Multipoint Taylor formulas and applications to the ÿnite element method, Numer. Math.
17 (1971) 84“100.
[10] C. Coatmelec, Approximation et interpolation des fonctions di erentiables de plusieurs variables, Ann. Sci. Ecole
Norm. Sup. 83 (1966) 271“341.
[11] R. Courant, Variational methods for the solution of problems of equilibrium and vibrations, Bull. Amer. Math. Soc.
49 (1943) 1“23.
[12] P.J. Davis, Interpolation and Approximation, Blaisdell, Walthan, MA, 1963 (2nd Edition, Dover, New York, 1975).
[13] Encyklopadie der mathematischen Wissenschaften, Teubner, Leipzig, pp. 1900 “1904.
[14] Encyclopà die des Sciences Mathematiques, Gauthier-Villars, Paris, 1906.
[15] I.A. Ezrohi, General forms of the remainder terms of linear formulas in multidimensional approximate analysis I,
II Mat. Sb. 38 (1956) 389 “ 416 and 43 (1957) 9 “28 (in Russian).
[16] M. Gasca, T. Sauer, Multivariate polynomial interpolation, Adv. Comput. Math., 12 (2000) 377“410.
[17] G. Glaeser, L™interpolation des fonctions di erentiables de plusieurs variables, in: C.T.C. Wall (Ed.), Proceedings of
Liverpool Singularities Symposium II, Lectures Notes in Mathematics, Vol. 209, Springer, Berlin, 1971, pp. 1“29.
[18] H.H. Goldstine, A History of Numerical Analysis from the 16th Through the 19th Century, Springer, Berlin, 1977.
[19] R.B. Guenter, E.L. Roetman, Some observations on interpolation in higher dimensions, Math. Comput. 24 (1970)
[20] E. Isaacson, H.B. Keller, Analysis of Numerical Methods, Wiley, New York, 1966.
[21] C.G.J. Jacobi, Theoremata nova algebraica circa systema duarum aequationum inter duas variabiles propositarum,
Crelle J. Reine Angew. Math. 14 (1835) 281“288.
[22] L. Kronecker, Uber einige Interpolationsformeln fur ganze Funktionen mehrerer Variabeln. Lecture at the academy
of sciences, December 21, 1865, in: H. Hensel (Ed.), L. Kroneckers Werke, Vol. I, Teubner, Stuttgart, 1895,
pp. 133“141. (reprinted by Chelsea, New York, 1968).
[23] K.S. Kunz, Numerical Analysis, McGraw-Hill, New York, 1957.
[24] G.G. Lorentz, Approximation of Funtions, Chelsea, New York, 1966.
[25] R. Lorentz, Multivariate Hermite interpolation by algebaic polynomials: a survey, J. Comput. Appl. Math. (2000)
this volume.
[26] S.E. Mikeladze, Numerical Methods of Mathematical Analysis, Translated from Russian, O ce of Tech. Services,
Department of Commerce, Washington DC, pp. 521“531.
[27] S. Narumi, Some formulas in the theory of interpolation of many independent variables, Tohoku Math. J. 18 (1920)
[28] L. Neder, Interpolationsformeln fur Funktionene mehrerer Argumente, Skandinavisk Aktuarietidskrift (1926) 59.
[29] R.A. Nicolaides, On a class of ÿnite elements generated by Lagrange interpolation, SIAM J. Numer. Anal. 9 (1972)
[30] R.A. Nicolaides, On a class of ÿnite elements generated by Lagrange interpolation II, SIAM J. Numer. Anal.
10 (1973) 182“189.
[31] K. Pearson, On the construction of tables and on interpolation, Vol. 2, Cambridge University Press, Cambridge,
[32] J. Radon, Zur mechanischen Kubatur, Monatshefte Math. Physik 52 (1948) 286“300.
M. Gasca, T. Sauer / Journal of Computational and Applied Mathematics 122 (2000) 23“35 35

[33] H.E. Salzer, Note on interpolation for a function of several variables, Bull. AMS 51 (1945) 279“280.
[34] H.E. Salzer, Table of coe cients for interpolting in functions of two variables, J. Math. Phys. 26 (1948) 294“305.
[35] H.E. Salzer, Note on multivariate interpolation for unequally spaced arguments with an application to double
summation, J. SIAM 5 (1957) 254“262.
[36] H.E. Salzer, Some new divided di erence algorithms for two variables, in: R.E. Langer (Ed.), On Numerical
Approximation, University of Wisconsin Press, Madison, 1959, pp. 61“98.
[37] H.E. Salzer, Divided di erences for functions of two variables for irregularly spaced arguments, Numer. Math.
6 (1964) 68“77.
[38] H.E. Salzer, Formulas for bivariate hyperosculatory interpolation, Math. Comput. 25 (1971) 119“133.
[39] H.E. Salzer, G.M. Kimbro, Tables for Bivariate Osculatory Interpolation over a Cartesian Grid, Convair Astronautics,
[40] A. Sard, Remainders: functions of several variables, Acta Math. 84 (1951) 319“346.
[41] A. Sard, Remainders as integrals of partial derivatives, Proc. Amer. Math. Soc. 3 (1952) 732“741.
[42] T. Sauer, Yuan Xu, On multivariate Lagrange interpolation, Math. Comput. 64 (1995) 1147“1170.
[43] T. Sauer, Yuan Xu, A case study in multivariate Lagrange interpolation, in: S.P. Singh (Ed.), Approximation Theory,
Wavelets and Applications, Kluwer Academic Publishers, Dordrecht, 1995, pp. 443“452.
[44] D.D. Stancu, The remainder of certain linear approximation formulas in two variables, J. SIAM Numer. Anal.
1 (1964) 137“163.
[45] I.F. Ste ensen, Interpolation, Chelsea, New York, 1927 (2nd Edition, 1950).
[46] G. Strang, G.J. Fix, An Analysis of the Finite Element Method, Prentice-Hall, Englewood Cli s, NJ, 1973.
[47] H.C. Thacher Jr., W.E. Milne, Interpolation in several variables, J. SIAM 8 (1960) 33“42.
[48] H.C. Thacher Jr., Derivation of interpolation formulas in several independent variables, Ann. N.Y. Acad. Sci.
86 (1960) 758“775.
[49] T.N. Thiele, Interpolationsrechnung, Teubner, Leipzig, 1909.
[50] R.S. Walker, Algebraic Curves, Springer, Berlin, 1978.
[51] E.T. Whittaker, G. Robinson, Calculus of Observations, 4th Edition, Blackie and Sons, London, 1944.
[52] O.C. Zienkiewicz, The Finite Element Method in Structural and Continuum Mechanics, McGraw-Hill, London, 1967.
[53] M. Zlamal, On the ÿnite element method, Numer. Math. 12 (1968) 394“409.
Journal of Computational and Applied Mathematics 122 (2000) 37“50

Elimination techniques: from extrapolation to totally positive
matrices and CAGD
M. Gasca — , G. Muhlbach
Department of Applied Mathematics, University of Zaragoza, 5009 Zaragoza, Spain
Received 20 May 1999; received in revised form 22 September 1999

In this survey, we will show some connections between several mathematical problems such as extrapolation, linear
systems, totally positive matrices and computer-aided geometric design, with elimination techniques as the common tool
to deal with all of them. c 2000 Elsevier Science B.V. All rights reserved.


. 8
( 83 .)