Barycentric lagrange interpolation

SIAM REVIEW Vol. 46, No. 3, pp. 501–517

c 2004 Society for Industrial and Applied Mathematics

Barycentric Lagrange Interpolation∗
Jean-Paul Berrut† Lloyd N. Trefethen‡
Dedicated to the memory of Peter Henrici (1923–1987)
Abstract. Barycentric interpolation is a variant of Lagrange polynomial interpolation that is fast and stable. It deserves to be known as the standard method ofpolynomial interpolation. Key words. barycentric formula, interpolation AMS subject classifications. 65D05, 65D25 DOI. 10.1137/S0036144502417715

1. Introduction. “Lagrangian interpolation is praised for analytic utility and beauty but deplored for numerical practice.” This heading, from the extended table of contents of one of the most enjoyable textbooks of numerical analysis [1], expresses awidespread view. In the present work we shall show that, on the contrary, the Lagrange approach is in most cases the method of choice for dealing with polynomial interpolants. The key is that the Lagrange polynomial must be manipulated through the formulas of barycentric interpolation. Barycentric interpolation is not new, but most students, most mathematical scientists, and even many numerical analysts donot know about it. This simple and powerful idea deserves a place at the heart of introductory courses and textbooks in numerical analysis.1 As always with polynomial interpolation, unless the degree of the polynomial is low, it is usually not a good idea to use uniformly spaced interpolation points. Instead one should use Chebyshev points or other systems of points clustered at the boundary ofthe interval of approximation. 2. Lagrange and Newton Interpolation. Let n + 1 distinct interpolation points (nodes) xj , j = 0, . . . , n, be given, together with corresponding numbers fj , which may or may not be samples of a function f . Unless stated otherwise, we assume that the nodes are real, although most of our results and comments generalize to the
∗ Received by the editors November 8,2002; accepted for publication (in revised form) December 24, 2003; published electronically July 30, 2004. http://www.siam.org/journals/sirev/46-3/41771.html † D´partement de Math´matiques, Universit´ de Fribourg, 1700 Fribourg/P´rolles, Switzerland e e e e (Jean-Paul.Berrut@unifr.ch). The work of this author was supported by the Swiss National Science Foundation, grant 21–59272.99. ‡ ComputingLaboratory, Oxford University, Oxford, UK (LNT@comlab.ox.ac.uk). 1 We are speaking here of a method of interpolating data in one space dimension by a polynomial. The “barycentric coordinates” popular in computational geometry for representing data on triangles and tetrahedra are related, but different.

501

502

JEAN-PAUL BERRUT AND LLOYD N. TREFETHEN

complex plane. Let Πn denote thevector space of all polynomials of degree at most n. The classical problem addressed here is that of finding the polynomial p ∈ Πn that interpolates f at the points xj , i.e., p(xj ) = fj , j = 0, . . . , n.

The problem is well-posed; i.e., it has a unique solution that depends continuously on the data. Moreover, as explained in virtually every introductory numerical analysis text, the solution canbe written in Lagrange form [44]:
n

(2.1)

p(x) =
j=0

fj j (x),

j (x) =

n k=0, k=j (x − xk ) . n k=0, k=j (xj − xk )

The Lagrange polynomial (2.2)
j (xk )

j

corresponding to the node xj has the property 1, j = k, 0, otherwise, j, k = 0, . . . , n.

=

At this point, many authors specialize Lagrange’s formula (2.1) for small numbers n of nodes, before asserting thatcertain shortcomings make it a bad choice for practical computations. Among the shortcomings sometimes claimed are these: 1. Each evaluation of p(x) requires O(n2 ) additions and multiplications. 2. Adding a new data pair (xn+1 , fn+1 ) requires a new computation from scratch. 3. The computation is numerically unstable. From here it is commonly concluded that the Lagrange form of p is mainly a...