I need help with …

Copyright ©2007 by the Society for Industrial and Applied Mathematics
This electronic version is for personal use and may not be duplicated or distributed.
Chapter 3
Chebyshev
Expansions
The best is the cheapest.
—Benjamin Franklin
3.1
Introduction
In Chapter 2, approximations were considered consisting of expansions around a specific
value of the variable (finite or infinite); both convergent and divergent series were described.
These are the preferred approaches when values around these points (either in R or C) are
needed.
In this chapter, approximations in real intervals are considered. The idea is to approximate a function f(x) by a polynomial p(x) that gives a uniform and accurate description
in an interval [a, b].
Let us denote by Pn the set of polynomials of degree at most n and let g be a bounded
function defined on [a, b]. Then the uniform norm ||g|| on [a, b] is given by
||g|| = max |g(x)|.
x∈[a,b]
(3.1)
For approximating a continuous function f on an interval [a, b], it is reasonable to
consider that the best option consists in finding the minimax approximation, defined as
follows.
Definition 3.1. q ∈ Pn is the best (or minimax) polynomial approximation to f on [a, b] if
||f − q|| ≤ ||f − p|| ∀p ∈ Pn .
(3.2)
Minimax polynomial approximations exist and are unique (see [152]) when f is
continuous, although they are not easy to compute in general. Instead, it is a more effective
approach to consider near-minimax approximations, based on Chebyshev polynomials.
51
From “Numerical Methods for Special Functions” by Amparo Gil, Javier Segura, and Nico Temme
Copyright ©2007 by the Society for Industrial and Applied Mathematics
This electronic version is for personal use and may not be duplicated or distributed.
52
Chapter 3. Chebyshev Expansions
Chebyshev polynomials form a special class of polynomials especially suited for
approximating other functions. They are widely used in many areas of numerical analysis:
uniform approximation, least-squares approximation, numerical solution of ordinary and
partial differential equations (the so-called spectral or pseudospectral methods), and so on.
In this chapter we describe the approximation of continuous functions by Chebyshev
interpolation and Chebyshev series and how to compute efficiently such approximations.
For the case of functions which are solutions of linear ordinary differential equations with
polynomial coefficients (a typical case for special functions), the problem of computing
Chebyshev series is efficiently solved by means of Clenshaw’s method, which is also presented in this chapter.
Before this, we give a very concise overview of well-known results in interpolation theory, followed by a brief summary of important properties satisfied by Chebyshev
polynomials.
3.2
Basic results on interpolation
Consider a real function f that is continuous on the real interval [a, b]. When values of this
function are known at a finite number of points xi , one can consider the approximation by
a polynomial Pn such that f(xi ) = Pn (xi ). The next theorem gives an explicit expression
for the lowest degree polynomial (the Lagrange interpolation polynomial) satisfying these
interpolation conditions.
Theorem 3.2 (Lagrange interpolation). Given a function f that is defined at n + 1 points
x0 < x1 < · · · < xn ∈ [a, b], there exists a unique polynomial of degree smaller than or equal to n such that Pn (xi ) = f(xi ), i = 0, . . . , n. (3.3) This polynomial is given by Pn (x) = n f(xi )Li (x), (3.4) i=0 where Li (x) is defined by !n πn+1 (x) j=0,j=i (x − xj ) = !n , (3.5) Li (x) = (x − xi )πn+1 (xi ) j=0,j=i (xi − xj ) ! πn+1 (x) being the nodal polynomial, πn+1 (x) = nj=0 (x − xj ). Additionally, if f is continuous on [a, b] and n + 1 times differentiable in (a, b), then for any x ∈ [a, b] there exists a value ζx ∈ (a, b), depending on x, such that Rn (x) = f(x) − Pn (x) = f n+1 (ζx ) πn+1 (x). (n + 1)! (3.6) Proof. The proof of this theorem can be found elsewhere [45, 48]. Li are called the fundamental Lagrange interpolation polynomials. From "Numerical Methods for Special Functions" by Amparo Gil, Javier Segura, and Nico Temme Copyright ©2007 by the Society for Industrial and Applied Mathematics This electronic version is for personal use and may not be duplicated or distributed. 3.2. Basic results on interpolation 53 The first part of the theorem is immediate and Pn satisfies the interpolation conditions, because the polynomials Li are such that Li (xj ) = δij . The formula for the remainder can be proved from repeated application of Rolle’s theorem (see, for instance, [48, Thm. 3.3.1]). For the particular case of Lagrange interpolation over n nodes, a simple expression for the interpolating polynomial can be given in terms of forward differences when the nodes are equally spaced, that is, xi+1 − xi = h, i = 0, . . . , n − 1. In this case, the interpolating polynomials of Theorem 3.2 can be written as n s i f0 , (3.7) Pn (x) = i i=0 where s= x − x0 , h i−1 s 1" (s − j), = i! j=0 i fj = fj+1 − fj , fj = f(xj ), (3.8) 2 fj = (fj+1 − fj ) = fj+2 − 2fj+1 + fj , . . . . This result is easy to prove by noticing that fs = ( + I)s f0 , s = 0, 1, . . . , n, and by expanding the binomial of commuting operators and I (I being the identity, Ifi = fi ). The formula for the remainder in Theorem 3.2 resembles that for the Taylor formula of degree n (Lagrange form), except that the nodal polynomials in the latter case contain only one node, x0 , which is repeated n + 1 times (in the sense that the power (x − x0 )n+1 appears). This interpretation in terms of repeated nodes can be generalized; both the Taylor formula and the Lagrange interpolation formula can be seen as particular cases of a more general interpolation formula, which is Hermite interpolation. Theorem 3.3 (Hermite interpolation). Let f be n times differentiable with continuity in [a, b] and n + 1 times differentiable in (a, b). Let x0 < x1 < · · · < xk ∈ [a, b], and let ni ∈ N such that n0 + n1 + · · · + nk = n − k. Then, there exists a unique polynomial Pn of degree not larger than n such that Pn(j) (xi ) = f (j) (xi ), j = 0, . . . , ni , i = 0, . . . , k. (3.9) Furthermore, given x ∈ [a, b], there exists a value ζx ∈ (a, b) such that f(x) = Pn (x) + f n+1 (ζx ) πn+1 (x), (n + 1)! (3.10) where πn+1 (x) is the nodal polynomial πn+1 (x) = (x − x0 )n0 +1 · · · (x − xk )nk +1 (3.11) in which each node xi is repeated ni + 1 times. Proof. For the proof we refer to [45]. An explicit expression for the interpolating polynomial is, however, not so easy as for Lagrange’s case. A convenient formalism is that of Newton’s divided difference formula, also for Lagrange interpolation (see [45] for further details). From "Numerical Methods for Special Functions" by Amparo Gil, Javier Segura, and Nico Temme Copyright ©2007 by the Society for Industrial and Applied Mathematics This electronic version is for personal use and may not be duplicated or distributed. 54 Chapter 3. Chebyshev Expansions For the case of a single interpolation node x0 which is repeated n times, the corresponding interpolating polynomial is just the Taylor polynomial of degree n at x0 . It is very common that successive derivatives of special functions are known at a certain point x = x0 (Taylor’s theorem, (2.1)), but it is not common that derivatives are known at several points. Therefore, in practical evaluation of special functions, Hermite interpolation different from the Taylor case is seldom used. Lagrange interpolation is, however, a very frequently used method of approximation and, in addition, will be behind the quadrature methods to be discussed in Chapter 5. For interpolating a function in a number of nodes, we need, however, to know the values which the function takes at these points. Therefore, in general we will need to rely on an alternative (high-accuracy) method of evaluation. However, for functions which are solutions of a differential equation, Clenshaw’s method (see §3.6.1) provides a way to compute expansions in terms of Chebyshev polynomials. Such infinite expansions are related to a particular and useful type of Lagrange interpolation that we discuss in detail in §3.6.1 and introduce in the next section. 3.2.1 The Runge phenomenon and the Chebyshev nodes Given a function f which is continuous on [a, b], we may try to approximate the function by a Lagrange interpolating polynomial. We could naively think that as more nodes are considered, the approximation will always be more accurate, but this is not always true. The main question to be addressed is whether the polynomials Pn that interpolate a continuous function f in n + 1 equally spaced points are such that lim ||f − Pn || = lim ||Rn || = 0, n→∞ n→∞ (3.12) where, if f is sufficiently differentiable, the error can be estimated through (3.6). A pathological example for which the Lagrange interpolation does not converge is provided by f(x) = |x| in the interval [−1, 1], for which equidistant interpolation diverges for 0 < |x| < 1 (see [189, Thm. 4.7]), as has been proved by Bernstein. A less pathological example, studied by Runge, showing the Runge phenomenon, gives a clear warning on the problems of equally spaced nodes. Considering the problem of interpolation of f(x) = 1 1 + x2 (3.13) on [−5, 5], Runge observed that limn→∞ ||f − Pn || = ∞, but that convergence takes place in a smaller interval [−a, a] with a 3.63. This bad behavior in Runge’s example is due to the values of the nodal polynomial πn+1 (x), which tends to present very strong oscillations near the endpoints of the interval (see Figure 3.1). From "Numerical Methods for Special Functions" by Amparo Gil, Javier Segura, and Nico Temme Copyright ©2007 by the Society for Industrial and Applied Mathematics This electronic version is for personal use and may not be duplicated or distributed. 3.2. Basic results on interpolation 55 2 6 10 1.5 5 10 1 0 0.5 5 -10 0 6 -10 -4 -2 0 2 -4 4 -2 0 2 4 x x Figure 3.1. Left: the function f(x) = 1/(1 + x2 ) is plotted in [−5, 5] together with the polynomial of degree 10 which interpolates f at x = 0, ±1, ±2, ±3 ± 4, ±5. Right: the nodal polynomial π(x) = x(x − 1)(x − 4)(x − 9)(x − 16)(x − 25). The uniformity of the error in the interval of interpolation can be considerably improved by choosing the interpolation nodes xi in a different way. Without loss of generality, we will restrict our study to interpolation on the interval [−1, 1]; the problem of interpolating f with nodes xi in the finite interval [a, b] is equivalent to the problem of interpolating g(t) = f(x(t)), where a+b b−a x(t) = + t (3.14) 2 2 with nodes ti ∈ [−1, 1]. Theorem 3.4 explains how to choose the nodes in [−1, 1] in order to minimize uniformly the error due to the nodal polynomial and to quantify this error. The nodes are given by the zeros of a Chebyshev polynomial. Theorem 3.4. Let xk ! = cos((k + 1/2)π/(n + 1)), k = 0, 1, . . . , n. Then the monic polynomial T̂n+1 (x) = nk=0 (x − xk ) is the polynomial of degree n + 1 with the smallest possible uniform norm (3.1) in [−1, 1] in the sense that ||T̂n+1 || ≤ ||qn+1 || (3.15) for any other monic polynomial qn+1 of degree n + 1. Furthermore, ||T̂n+1 || = 2−n . (3.16) The selection of these nodes will not guarantee convergence as the number of nodes tends to infinity, because it also depends on how the derivatives of the function f behave, but certainly enlarges the range of functions for which convergence takes place and eliminates the problem for the example provided by Runge. Indeed, taking as nodes xk = 5 cos ((k + 1/2)π/11) , k = 0, 1, . . . , 10, (3.17) instead of the 11 equispaced points, the behavior is much better, as illustrated in Figure 3.2. From "Numerical Methods for Special Functions" by Amparo Gil, Javier Segura, and Nico Temme Copyright ©2007 by the Society for Industrial and Applied Mathematics This electronic version is for personal use and may not be duplicated or distributed. 56 Chapter 3. Chebyshev Expansions 1 2 0.8 1.5 0.6 1 0.4 0.5 0.2 0 0 -4 -2 0 2 4 -4 -2 x 0 2 4 x Figure 3.2. Left: the function f(x) = 1/(1 + x2 ) is plotted together with the interpolation polynomial for the 11 Chebyshev points (see (3.17)). Right: the interpolation errors for equispaced points and Chebyshev points is shown. Before proving Theorem 3.4 and further results, we summarize the basic properties of Chebyshev polynomials, the zeros of which are the nodes in Theorem 3.4. 3.3 Chebyshev polynomials: Basic properties Let us first consider a definition and some properties of the Chebyshev polynomials of the first kind. Definition 3.5 (Chebyshev polynomial of the first kind Tn (x)). The Chebyshev polynomial of the first kind of order n is defined as follows: Tn (x) = cos n cos−1 (x) , x ∈ [−1, 1], n = 0, 1, 2, . . . . (3.18) From this definition the following property is evident: Tn (cos θ) = cos (nθ), 3.3.1 θ ∈ [0, π], n = 0, 1, 2, . . . . (3.19) Properties of the Chebyshev polynomials Tn (x) The polynomials Tn (x), n ≥ 1, satisfy the following properties, which follow straightforwardly from (3.19). (i) The Chebyshev polynomials Tn (x) satisfy the following three-term recurrence relation: Tn+1 (x) = 2xTn (x) − Tn−1 (x), n = 1, 2, 3, . . . , (3.20) with starting values T0 (x) = 1, T1 (x) = x. From "Numerical Methods for Special Functions" by Amparo Gil, Javier Segura, and Nico Temme Copyright ©2007 by the Society for Industrial and Applied Mathematics This electronic version is for personal use and may not be duplicated or distributed. 3.3. Chebyshev polynomials: Basic properties 57 1 0.8 0.6 0.4 T0(x) T (x) 1 T (x) 2 T (x) 3 T (x) 4 T (x) 0.2 0 −0.2 5 −0.4 −0.6 −0.8 −1 −1 −0.8 −0.6 −0.4 −0.2 0 x 0.2 0.4 0.6 0.8 1 Figure 3.3. Chebyshev polynomials of the first kind Tn (x), n = 0, 1, 2, 3, 4, 5. Explicit expressions for the first six Chebyshev polynomials are T0 (x) = 1, T1 (x) = x, T2 (x) = 2x2 − 1, T3 (x) = 4x3 − 3x, (3.21) T4 (x) = 8x4 − 8x2 + 1, T5 (x) = 16x5 − 20x3 + 5x. The graphs of these Chebyshev polynomials are plotted in Figure 3.3. (ii) The leading coefficient (of xn ) in Tn (x) is 2n−1 and Tn (−x) = (−1)n Tn (x). (iii) Tn (x) has n zeros which lie in the interval (−1, 1). They are given by 2k + 1 xk = cos π , k = 0, 1, . . . , n − 1. 2n (3.22) Tn (x) has n + 1 extrema in the interval [−1, 1] and they are given by xk = cos kπ , n k = 0, 1, . . . , n. (3.23) At these points, the values of the polynomials are Tn (xk ) = (−1)k . With these properties, it is easy to prove Theorem 3.4, which can also be expressed in the following way. From "Numerical Methods for Special Functions" by Amparo Gil, Javier Segura, and Nico Temme Copyright ©2007 by the Society for Industrial and Applied Mathematics This electronic version is for personal use and may not be duplicated or distributed. 58 Chapter 3. Chebyshev Expansions Theorem 3.6. The polynomial T̂n (x) = 21−n Tn (x) is the minimax approximation on [−1, 1] to the zero function by a monic polynomial of degree n and ||T̂n || = 21−n . (3.24) Proof. Let us suppose that there exists a monic polynomial pn of degree n such that |pn (x)| ≤ 21−n for all x ∈ [−1, 1], and we will arrive at a contradiction. Let xk , k = 0, . . . , n, be the abscissas of the extreme values of the Chebyshev polynomial of degree n. Because of property (ii) of this section we have pn (x1 ) > 21−n Tn (x1 ),
pn (x0 ) < 21−n Tn (x0 ), pn (x2 ) > 21−n Tn (x2 ), . . . .
Therefore, the polynomial
Q(x) = pn (x) − 21−n Tn (x)
changes sign between each two consecutive extrema of Tn (x). Thus, it changes sign n
times. But this is not possible because Q(x) is a polynomial of degree smaller than n (it is
a subtraction of two monic polynomials of degree n).
Remark 1. The monic Chebyshev polynomial T̂n (x) is not the minimax approximation
in Pn (Definition 3.1) of the zero function. The minimax approximation in Pn of the zero
function is the zero polynomial.
Further properties
Next we summarize additional properties of the Chebyshev polynomials of the first kind
that will be useful later. For further properties and proofs of these results see, for instance,
[148, Chaps. 1–2].
(a) Relations with derivatives.




T0 (x) = T1 (x),




T1 (x) = 41 T2 (x),







(x)
Tn+1 (x) Tn−1

Tn (x) = 1
−
2
n+1
n−1 ,
(3.25)
n ≥ 2,
(1 − x2 )Tn (x) = n [xTn (x) − Tn+1 (x)] = n [Tn−1 (x) − xTn (x)] .
(3.26)
(b) Multiplication relation.
2Tr (x)Tq (x) = Tr+q (x) + T|r−q| (x),
(3.27)
with the particular case q = 1,
2xTr (x) = Tr+1 (x) + T|r−1| (x).
(3.28)
From “Numerical Methods for Special Functions” by Amparo Gil, Javier Segura, and Nico Temme
Copyright ©2007 by the Society for Industrial and Applied Mathematics
This electronic version is for personal use and may not be duplicated or distributed.
3.3. Chebyshev polynomials: Basic properties
59
(c) Orthogonality relation.
1
−1
Tr (x)Ts (x)(1 − x2 )−1/2 dx = Nr δrs ,
(3.29)
with N0 = π and Nr = 12 π if r = 0.
(d) Discrete orthogonality relation.
1. With the zeros of Tn+1 (x) as nodes: Let n > 0, r, s ≤ n, and let xj =
cos((j + 1/2)π/(n + 1)). Then
n

Tr (xj )Ts (xj ) = Kr δrs ,
(3.30)
j=0
where K0 = n + 1 and Kr = 12 (n + 1) when 1 ≤ r ≤ n.
2. With the extrema of Tn (x) as nodes: Let n > 0, r, s ≤ n, and xj = cos(πj/n).
Then
n

Tr (xj )Ts (xj ) = Kr δrs ,
(3.31)
j=0
where K0 = Kn = n and Kr = 12 n when 1 ≤ r ≤ n − 1.
The double prime indicates that the terms with suffixes j = 0 and j = n are to be
halved.
(e) Polynomial representation.
The expression of Tn (x) in terms of powers of x is given by (see [38, 201])
[n/2]
Tn (x) =

dk(n) xn−2k ,
(3.32)
k=0
where
dk(n) = (−1)k 2n−2k−1

n
n−k
,
n−k
k
2k < n, (3.33) and dk(2k) = (−1)k , k ≥ 0. (3.34) (f) Power representation. The power xn can be expressed in terms of Chebyshev polynomials as follows: x =2 n 1−n [n/2] n k=0 k Tn−2k (x), (3.35) where the prime indicates that the term for k = 0 is to be halved. From "Numerical Methods for Special Functions" by Amparo Gil, Javier Segura, and Nico Temme Copyright ©2007 by the Society for Industrial and Applied Mathematics This electronic version is for personal use and may not be duplicated or distributed. 60 Chapter 3. Chebyshev Expansions The first three properties are immediately obtained from the definition of Chebyshev polynomials. Property (c) means that the set of Chebyshev polynomials {Tn (x)} is an orthogonal set with respect to the weight function w(x) = (1 − x2 )−1/2 in the interval (−1, 1). This concept is developed in Chapter 5, and it is shown that this orthogonality implies the first discrete orthogonality of property (d); see (5.86). This property, as well as the second discrete orthogonality, can also be easily proved using trigonometry (see [148, Chap. 4]). See also [148, Chap. 2] for a proof of the last two properties. Shifted Chebyshev polynomials Shifted Chebyshev polynomials are also of interest when the range of the independent variable is [0, 1] instead of [−1, 1]. The shifted Chebyshev polynomials of the first kind are defined as Tn∗ (x) = Tn (2x − 1), 0 ≤ x ≤ 1. (3.36) Similarly, one can also build shifted polynomials for a generic interval [a, b]. Explicit expressions for the first six shifted Chebyshev polynomials are T0∗ (x) = 1, T1∗ (x) = 2x − 1, T2∗ (x) = 8x2 − 8x + 1, T3∗ (x) = 32x3 − 48x2 + 18x − 1, (3.37) T4∗ (x) = 128x4 − 256x3 + 160x2 − 32x + 1, T5∗ (x) = 512x5 − 1280x4 + 1120x3 − 400x2 + 50x − 1. 3.3.2 Chebyshev polynomials of the second, third, and fourth kinds Chebyshev polynomials of the first kind are a particular case of Jacobi polynomials Pn(α,β) (x) (up to a normalization factor). Jacobi polynomials, which can be defined through the Gauss hypergeometric function (see §2.3) as −n, n + α + β + 1 1 − x n+α ; F Pn(α,β) (x) = , (3.38) 2 1 α+1 n 2 are orthogonal polynomials on the interval [−1, 1] with respect to the weight function w(x) = (1 − x)α (1 + x)β , α, β > −1, that is,
1
−1
Pr(α,β) (x)Ps(α,β) (x)w(x) dx = Mr δrs .
(3.39)
In particular, for the case α = β = −1/2 we recover the orthogonality relation (3.29).
Furthermore,

−n, n 1 − x
Tn (x) = 2 F1
.
(3.40)
;
1/2
2
From “Numerical Methods for Special Functions” by Amparo Gil, Javier Segura, and Nico Temme
Copyright ©2007 by the Society for Industrial and Applied Mathematics
This electronic version is for personal use and may not be duplicated or distributed.
3.3. Chebyshev polynomials: Basic properties
61
As we have seen, an important property satisfied by the polynomials Tn (x) is that,
with the change x = cos θ, the zeros and extrema are equally spaced in the θ variable. The
zeros of Tn (x) (see (3.18)) satisfy
θk − θk−1 = | cos−1 (xk ) − cos−1 (xk−1 )| = π/n,
(3.41)
and similarly for the extrema.
This is not the only case of Jacobi polynomials with equally spaced zeros (in the θ
variable), but it is the only case with both zeros and extrema equispaced. Indeed, considering
the Liouville–Green transformation (see §2.2.4) with the change of variable x = cos θ, we
can prove that

α+1/2
β+1/2
1
1
u(α,β)
cos
(θ)
=
sin
θ
θ
Pn(α,β) (cos θ),
n
2
2
0 ≤ θ ≤ π,
(3.42)
satisfies the differential equation
d 2 u(α,β)
(θ)
n
+ (θ)u(α,β)
(θ) = 0,
n
dθ 2 

1
1
− α2
− β2
4
4
1
2
.
(θ) = 4 (2n + α + β + 1) +
+
cos2 21 θ
sin2 21 θ
(3.43)
From this, we observe that for the values |α| = |β| = 12 , and only for these values, (θ) is
constant and therefore the solutions are trigonometric functions
u(α,β)
= C(α,β) cos(θw(α,β)
+ φ(α,β) ),
n
n
w(α,β)
= n + (α + β + 1)/2
n
(3.44)
with C(α,β) and φ(α,β) values not depending on θ. The solutions u(α,β)
, |α| = |β| = 12 have
n
therefore equidistant zeros and extrema. The distance between zeros is
θk − θk−1 =
π
.
n + (α + β + 1)/2
(3.45)
Jacobi polynomials have the same zeros as the solutions u(α,β)
(except that θ = 0, π may
n
also be zeros for the latter). Therefore, Jacobi polynomials have equidistant zeros for
|α| = |β| = 12 . However, due to the sine and cosine factors in (3.42), the extrema of Jacobi
polynomials are only equispaced when α = β = − 12 .
The four types of Chebyshev polynomials are the only classical orthogonal (hypergeometric) polynomials for which the elementary change of variables x = cos θ makes all
zeros equidistant. Furthermore, these are the only possible cases for which equidistance
takes place, not only in the θ variable but also under more general changes of variable (also
including confluent cases) [52, 53].
Chebyshev polynomials are proportional to the Jacobi polynomials with equispaced
θ zeros. From (3.42) such Chebyshev polynomials can be written as
cos(θw(α,β)
+ φ(α,β) )
n
α+1/2
β+1/2 .
sin 12 θ
cos 12 θ
Tnα,β (θ) = C(α,β)
(3.46)
From “Numerical Methods for Special Functions” by Amparo Gil, Javier Segura, and Nico Temme
Copyright ©2007 by the Society for Industrial and Applied Mathematics
This electronic version is for personal use and may not be duplicated or distributed.
62
Chapter 3. Chebyshev Expansions
C(α,β) can be arbitrarily chosen and it is customary to take C(α,β) = 1, except when
α = β = 12 , in which case C(α,β) = 12 . On the other hand, for each selection of α and β
(with |α| = |β| = 12 ) there is only one possible selection of φ(α,β)) in [0, π) which gives
a polynomial solution. This phase is easily selected by requiring that Tnα,β (θ) be finite as
θ → 0, π. With the standard normalization considered, the four families of polynomials
Tnα,β (θ) (proportional to Pn(α,β) ) can be written as
(−1/2,−1/2)
Tn
(1/2,1/2)
Tn
(θ) = cos(nθ) = Tn (x),
(θ) = sin((n + 1)θ) = Un (x),
sin θ
(−1/2,1/2)
(θ) =
Tn
(1/2,−1/2)
Tn
(θ) =
cos((n + 12 )θ)
cos( 12 θ)
sin((n + 12 )θ)
sin( 12 θ)
= Vn (x),
(3.47)
= Wn (x).
These are the Chebyshev polynomials of first (T ), second (U), third (V ), and fourth (W )
kinds. The third- and fourth-kind polynomials are trivially related because Pn(α,β) (x) =
(−1)n Pn(β,α) (−x).
Particularly useful for some applications are Chebyshev polynomials of the second
(1/2,1/2)
(θ(x)))
kind. The zeros of Un (x) plus the nodes x = −1, 1 (that is, the x zeros of un
are the nodes of the Clenshaw–Curtis quadrature rule (see §9.6.2). All Chebyshev polynomials satisfy three-term recurrence relations, as is the case for any family of orthogonal
polynomials; in particular, the Chebyshev polynomials of the second kind satisfy the same
recurrence as the polynomials of the first kind. See [2] or [148] for further properties.
3.4
Chebyshev interpolation
Because the scaled Chebyshev polynomial T̂n+1 (x) = 2−n Tn+1 (x) is the monic polynomial
of degree n + 1 with the smallest maximum absolute value in [−1, 1] (Theorem 3.6), the
selection of its n zeros for Lagrange interpolation leads to interpolating polynomials for
which the Runge phenomenon is absent.
By considering the estimation for the Lagrange interpolation error (3.6) under the
condition of Theorem 3.2, taking as interpolation nodes the zeros of Tn+1 (x),

π
1
, k = 0, . . . , n,
(3.48)
k+
xk = cos
2 n+1
and considering the minimax property of the nodal polynomial T̂n+1 (x) (Theorem 3.6), the
following error bound can be obtained:
|Rn (x)| =
1
|f (n+1) (ζx )|
|f (n+1) (ζx )|
|T̂n+1 (x)| ≤ 2−n
≤ n
||f (n+1) ||,
(n + 1)!
(n + 1)!
2 (n + 1)!
(3.49)
From “Numerical Methods for Special Functions” by Amparo Gil, Javier Segura, and Nico Temme
Copyright ©2007 by the Society for Industrial and Applied Mathematics
This electronic version is for personal use and may not be duplicated or distributed.
3.4. Chebyshev interpolation
63
where ||f (n+1) || = maxx∈[−1,1] |f (n+1) (x)|. By considering a linear change of variables
(3.14), an analogous result can be given for Chebyshev interpolation in an interval [a, b].
Interpolation with Chebyshev nodes is not as good as the best approximation (Definition 3.1), but usually it is the best practical possibility for interpolation and certainly much
better than equispaced interpolation. The best polynomial approximation is characterized
by the Chebyshev equioscillation theorem.
Theorem 3.7 (Chebyshev equioscillation theorem). For any continuous function f in
[a, b], a unique minimax polynomial approximation in Pn (the space of the polynomials of
degree n at most) exists and is uniquely characterized by the alternating or equioscillation
property that there are at least n + 2 points at which f(x) − Pn (x) attains its maximum
absolute value, with alternating signs.
Proof. Proofs of this theorem can be found, for instance, in [48, 189].
Because the function f(x) − Pn (x) alternates signs between each two consecutive
extrema, it has at least n + 1 zeros; therefore Pn is a Lagrange interpolating polynomial,
interpolating f at n + 1 points in [a, b]. The specific location of these points depends on
the particular function f , which makes the computation of best approximations difficult in
general.
Chebyshev interpolation by a polynomial in Pn , interpolating the function f at the
n + 1 zeros of Tn+1 (x), can be a reasonable approximation and can be computed in an
effective and stable way. Given the properties of the error for Chebyshev interpolation on
[−1, 1] and the uniformity in the deviation of Chebyshev polynomials with respect to zero
(Theorem 3.6), one can expect that Chebyshev interpolation gives a fair approximation to the
minimax approximation when the variation of f is soft. In addition, the Runge phenomenon
does not occur.
Uniform convergence (in the sense of (3.12)) does not necessarily hold but, in fact,
there is no system of preassigned nodes that can guarantee uniform convergence for any
continuous function f (see [189, Thm. 4.3]). The sequence of best uniform approximations pn for a given continuous function f does uniformly converge. For the Chebyshev
interpolation we need to consider some additional “level of continuity” in the form of the
modulus of continuity.
Definition 3.8. Let f be a function defined in an interval [a, b]. We define the modulus of
continuity as
ω(δ) =
sup
x1 ,x2 ∈[a,b]
|x1 −x2 | 1,
(3.76)
From “Numerical Methods for Special Functions” by Amparo Gil, Javier Segura, and Nico Temme
Copyright ©2007 by the Society for Industrial and Applied Mathematics
This electronic version is for personal use and may not be duplicated or distributed.
3.6. Computing the coefficients of a Chebyshev expansion
69
0.04
Abs. error
0.03
0.02
0.01
0
-1
-0.5
0
0.5
1
x
Figure 3.4. Error after truncating the series in (3.73) after the term k = 5.
then |f(x) − Sn (x)| = O(r −n ) as n → ∞ for all x ∈ [−1, 1].
The ellipse Er has semiaxis of length (r + 1/r)/2 on the real z-axis and of length
(r − 1/r)/2 on the imaginary axis.
For entire functions f we can take any number r in this theorem, and in fact the rate
of convergence can be of order O(1/n!). For example, we have the generating function for
the modified Bessel coefficients In (z) given by
ezx = I0 (z)T0 (x) + 2
∞

In (z)Tn (x),
−1 ≤ x ≤ 1,
(3.77)
n=1
where z can be any complex number. The Bessel functions behave like In (z) = O((z/2)n /n!)
as n → ∞ with z fixed, and the error |ezx − Sn (x)| has a similar behavior. The absolute
error when the series is truncated after the n = 5 term is shown in Figure 3.5.
3.6
Computing the coefficients of a Chebyshev expansion
In general, the Chebyshev coefficients of the Chebyshev expansion of a function f can be
approximately obtained by the numerical computation of the integral of (3.67). To improve
the speed of computation, fast Fourier cosine transform algorithms for evaluating the sums
in (3.68) can be considered. For numerical aspects of the fast Fourier transform we refer
the reader to [226].
From “Numerical Methods for Special Functions” by Amparo Gil, Javier Segura, and Nico Temme
Copyright ©2007 by the Society for Industrial and Applied Mathematics
This electronic version is for personal use and may not be duplicated or distributed.
70
Chapter 3. Chebyshev Expansions
0.004
Abs. Err.
0.003
0.002
0.001
0
-1
-0.5
0
0.5
1
x
Figure 3.5. Error after truncating the series for e2x in (3.77) after the term n = 5.
Compare with Figure 3.4.
In the particular case when the function f is a solution of an ordinary linear differential equation with polynomial coefficients, Clenshaw [37] proposed an alternative method,
which we will now discuss.
3.6.1
Clenshaw’s method for solutions of linear differential
equations with polynomial coefficients
The method works as follows.
Let us assume that f satisfies a linear differential equation in the variable x with
polynomial coefficients pk (x),
m

pk (x)f (k) (x) = h(x),
(3.78)
k=0
and where the coefficients of the Chebyshev expansion of the function h are known. In
general, conditions on the solution f will be given at x = 0 or x = ±1.
Let us express formally the sth derivative of f as follows:
f (s) (x) = 12 c0(s) + c1(s) T1 (x) + c2(s) T2 (x) + · · · .
(3.79)
Then the following expression can be obtained for the coefficients:
(s+1)
(s+1)
− cr+1
,
2rcr(s) = cr−1
r ≥ 1.
(3.80)
From “Numerical Methods for Special Functions” by Amparo Gil, Javier Segura, and Nico Temme
Copyright ©2007 by the Society for Industrial and Applied Mathematics
This electronic version is for personal use and may not be duplicated or distributed.
3.6. Computing the coefficients of a Chebyshev expansion
71
To see how to arrive to this equation, let us start with
(1)
Tn−1 (x)
f (x) = 12 c0(1) + c1(1) T1 (x) + c2(1) T2 (x) + · · · + cn−1
(3.81)
(1)
+ cn(1) Tn (x) + cn+1
Tn+1 (x) + · · ·
and integrate this expression. Using the relations in (3.25), we obtain
1
1
1
f(x) = c0 + c0(1) T1 (x) + c1(1) T2 (x) + · · ·
2
4
2

1 (1) Tn (x) Tn−2 (x)
+ cn−1
−
n
n−2
2

1 (1) Tn+1 (x) Tn−1 (x)
+ cn
−
n+1
n−1
2

1 (1) Tn+2 (x) Tn (x)
+ cn+1
−
+ ···.
2
n+2
n
(3.82)
Comparing the coefficients of the Chebyshev polynomials in this expression and the
Chebyshev expansion of f , we arrive at (3.80) for s = 1. Observe that a relation for c0
is not obtained in this way. Substituting in (3.82) given values of f at, say, x = 0 gives a
relation between c0 and an infinite number of coefficients cn(1) .
A next element in Clenshaw’s method is using (3.28) to handle the powers of x
occurring in the differential equation satisfied by f . Denoting the coefficients of Tr (x) in
the expansion of g(x) by Cr (g) when r > 0 and twice this coefficient when r = 0, and using
(3.28), we infer that

(s)
(s)
Cr xf (s) = 12 cr+1
.
(3.83)
+ c|r−1|
This expression can be generalized as follows:
p

1 p (s)
c
.
Cr xp f (s) = p
2 j=0 j |r−p+2j|
(3.84)
When the expansion (3.79) is substituted into the differential equation (3.78) together
with (3.80), (3.84), and the associated boundary conditions, it is possible to obtain an infinite
set of linear equations for the coefficients cr(s) . Two strategies can be used for solving these
equations.
Recurrence method. The equations can be solved by recurrence for r = N − 1, N −
2, . . . , 0, where N is an arbitrary (large) positive integer, by assuming that cr(s) = 0
(s)
. This is done as follows.
for r > N and by assigning arbitrary values to cN
(s)
Consider r = N in (3.80) and compute cN−1
, s = 1, . . . , m. Then, considering
(3.84) and the differential equation (3.78), select r appropriately in order to compute
(0)
cN−1
= cN−1 . We repeat the process by considering r = N − 1 in (3.80) and
(s)
computing cN−2
, etc. Obviously and unfortunately, the computed coefficients cr will
not satisfy, in general, the boundary conditions, and we will have to take care of these
in each particular case.
From “Numerical Methods for Special Functions” by Amparo Gil, Javier Segura, and Nico Temme
Copyright ©2007 by the Society for Industrial and Applied Mathematics
This electronic version is for personal use and may not be duplicated or distributed.
72
Chapter 3. Chebyshev Expansions
Iterative method. The starting point in this case is an initial guess for cr which satisfies
the boundary conditions. Using these values, we use (3.80) to obtain the values of
cr(s) , s = 1, . . . , m, and then the relation (3.84) and the differential equation (3.78) to
compute corrected values of cr .
The method based on recursions is, quite often, more rapidly convergent than the
iterative method; therefore, and in general, the iterative method could be useful for correcting
the rounding errors arising in the application of the method based on recursions.
Example 3.14 (Clenshaw’s method for the J -Bessel function). Let us consider, as a simple example (due to Clenshaw), the computation of the Bessel function J0 (t) in the range
0 ≤ t ≤ 4. This corresponds to solving the differential equation for J0 (4x), that is,
xy + y + 16xy = 0
(3.85)
in the range 0 ≤ x ≤ 1 with conditions y(0) = 1, y (0) = 0. This is equivalent to solving
the differential equation in [−1, 1], because J0 (x) = J0 (−x), x ∈ R.
Because J0 (4x) is an even function of x, the Tr (x) of odd order do not appear in
its Chebyshev expansion. By substituting the Chebyshev expansion into the differential
equation, we obtain
Cr (xy ) + Cr (y ) + 16Cr (xy) = 0,
r = 1, 3, 5, . . . ,
and considering (3.84),

1
cr−1 + cr+1
+ cr + 8 (cr−1 + cr+1 ) = 0,
2
r = 1, 3, 5, . . . .
(3.86)
(3.87)
This equation can be simplified. First, we see that by replacing r → r − 1 and r → r + 1
in (3.87) and subtracting both expressions, we get

1

cr−2 + cr − cr − cr+2
+ cr−1 − cr+1
2
(3.88)
+ 8 (cr−2 + cr − cr − cr+2 ) = 0, r = 2, 4, 6, . . . .
It is convenient to eliminate the terms with the second derivatives. This can be done by
using (3.80). In this way,

r cr−1
(3.89)
+ 8 (cr−2 − cr+2 ) = 0, r = 2, 4, 6, . . . .
+ cr+1
Now, expressions (3.80) and (3.89) can be used alternatively in the recurrence process, as
follows:



cr−1
= cr+1
+ 2rcr
r = N, N − 2, N − 4, . . . , 2.
(3.90)


cr−2 = cr+2 − 18 r cr−1
+ cr+1
As an illustration, let us take as first trial coefficient c̃20 = 1 and all higher order
coefficients zero. Applying the recurrences (3.90) and considering the calculation with 15
significant digits, we obtain the values of the trial coefficients given in Table 3.1.
Using the coefficients in Table 3.1, the trial solution of (3.85) at x = 0 is given by
ỹ(0) = 12 c̃0 − c̃2 + c̃4 − c̃6 + c̃8 − · · · = 8050924923505.5,
(3.91)
From “Numerical Methods for Special Functions” by Amparo Gil, Javier Segura, and Nico Temme
Copyright ©2007 by the Society for Industrial and Applied Mathematics
This electronic version is for personal use and may not be duplicated or distributed.
3.6. Computing the coefficients of a Chebyshev expansion
73
Table 3.1. Computed coefficients in the recurrence processes (3.90). We take as
starting values c̃20 = 1 and 15 significant digits in the calculations.
r

c̃r+1
c̃r
0
807138731281 −8316500240280
2 −5355660492900 13106141731320
4
2004549104041 −2930251101008
6 −267715177744
282331031920
8
18609052225
−15413803680
10
−797949504
545186400
12
23280625
−13548600
14
−492804
249912
16
7921
−3560
18
−100
40
20
1
Table 3.2. Computed coefficients of the Chebyshev expansion of the solution of (3.85).
r
cr
0
0.1002541619689529 10−0
2 −0.6652230077644372 10−0
4
0.2489837034982793 10−0
6 −0.3325272317002710 10−1
8
0.2311417930462743 10−2
10 −0.9911277419446611 10−4
12
0.2891670860329331 10−5
14 −0.6121085523493186 10−7
16
0.9838621121498511 10−9
18 −0.1242093311639757 10−10
20
0.1242093311639757 10−12
and the final values for the coefficients cr of the solution y(x) of (3.85) will be obtained by
dividing the trial coefficients by ỹ(0). This gives the requested values shown in Table 3.2.
The value of y(1) will then be given by
y(1) = 12 c0 + c2 + c4 + c6 + c8 + · · · = −0.3971498098638699,
the relative error being 0.57 10−13 when compared with the value of J0 (4).
(3.92)

From “Numerical Methods for Special Functions” by Amparo Gil, Javier Segura, and Nico Temme
Copyright ©2007 by the Society for Industrial and Applied Mathematics
This electronic version is for personal use and may not be duplicated or distributed.
74
Chapter 3. Chebyshev Expansions
Remark 2. Several questions arise in this successful method. The recursion given in (3.90)
is rather simple, and we can find its exact solution; cf. the expansion of the J0 in (3.139).
In Chapter 4 we explain that in this case the backward recursion scheme for computing the
Bessel coefficients is stable. In more complicated recursion schemes this information is not
available. The scheme may be of large order and may have several solutions of which the
asymptotic behavior is unknown. So, in general, we don’t know if Clenshaw’s method for
differential equations computes the solution that we want, and if for the wanted solution the
scheme is stable in the backward direction.
Example 3.15 (Clenshaw’s method for the Abramowitz function). Another but not so
easy example of the application of Clenshaw’s method is provided by MacLeod [146]
for the computation of the Abramowitz functions [1],
∞
Jn (x) =
t n e−t −x/t dt,
2
n integer.
(3.93)
0
Chebyshev expansions for J1 (x) for x ≥ 0 can be obtained by considering the following
two cases depending on the range of the argument x.
If 0 ≤ x ≤ a,
J1 (x) = f1 (x) −
√
πxg1 (x) − x2 h1 (x) log x,
(3.94)
where f1 , g1 , and h1 satisfy the system of equations
xg1 + 3g1 + 2g1 = 0,

x2 h
1 + 6xh1 + 6h1 + 2xh1 = 0,
(3.95)
xf1 + 2f1 = 3×2 h1 + 9xh1 + 2h1 ,
with appropriate initial
f1 , g1 , and h1 are expanded in
conditions at x = 0. The2 functions
2
a series of the form ∞
c
T
(t),
where
t
=
(2x
/a
)
−
1.
k
k
k=0
If x > a,
π ν −ν
e q1 (ν)
(3.96)
3 3
with ν = 3 (x/2)2/3 . The function q1 (ν) can be expanded in a Chebyshev series of
the variable
a 2/3
2B
t=
,
(3.97)
− 1, B = 3
ν
2
and q1 satisfies the differential equation
J1 (x) ∼
4ν3 q1 − 12ν3 q1 + (12ν3 − 5ν)q1 + (5ν + 5)q1 = 0,
(3.98)
where the derivatives
are taken with respect to ν. The function q1 is expanded in a
series of the form ∞
k=0 ck Tk (t), where t is given in (3.97).
The transition point a is selected in such a way that a and B are exactly represented.
Also, the number of terms needed for the evaluation of the Chebyshev expansions for a
prescribed accuracy is taken into account.
The differential equations (3.95) and (3.98) are solved by using Clenshaw’s
method.

From “Numerical Methods for Special Functions” by Amparo Gil, Javier Segura, and Nico Temme
Copyright ©2007 by the Society for Industrial and Applied Mathematics
This electronic version is for personal use and may not be duplicated or distributed.
3.7. Evaluation of a Chebyshev sum
3.7
75
Evaluation of a Chebyshev sum
Frequently one has to evaluate a partial sum of a Chebyshev expansion, that is, a finite series
of the form
SN (x) = 12 c0 +
N

ck Tk (x).
(3.99)
k=1
Assuming we have already computed the coefficients ck , k = 0, . . . , N, of the expansion, it would be nice to avoid the explicit computation of the Chebyshev polynomials
appearing in (3.99), although they easily follow from the relations (3.18) and (3.19). A first
possibility for the computation of this sum is to rewrite the Chebyshev polynomials Tk (x)
in terms of powers of x and then use the Horner scheme for the evaluation of the resulting
polynomial expression. However, one has to be careful when doing this because for some
expansions there is a considerable loss of accuracy due to cancellation effects.
3.7.1
Clenshaw’s method for the evaluation of a Chebyshev sum
An alternative and efficient method for evaluating this sum is due to Clenshaw [36]. This
scheme of computation, which can also be used for computing partial sums involving other
types of polynomials, corresponds to the following algorithm.
Algorithm 3.1. Clenshaw’s method for a Chebyshev sum.
Input: x; c0 , c1 , . . . , cN .

Output: .
SN (x) = N
k=0 ck Tk (x).
•
bN+1 = 0; bN = cN .
•
DO r = N − 1, N − 2, . . . , 1:
br = 2xbr+1 − br+2 + cr .
•
.
SN (x) = xb1 − b2 + c0 .
Let us explain how we arrived at this algorithm. For simplicity, let us first consider
the evaluation of
N

.
ck Tk (x) = 12 c0 + SN (x).
(3.100)
SN (x) =
k=0
This expression can be written in vector form as follows:


T0 (x)
 T1 (x) 


.
SN (x) = cT t = (c0 , c1 , . . . , cN ) 
.
..


.
(3.101)
TN (x)
From “Numerical Methods for Special Functions” by Amparo Gil, Javier Segura, and Nico Temme
Copyright ©2007 by the Society for Industrial and Applied Mathematics
This electronic version is for personal use and may not be duplicated or distributed.
76
Chapter 3. Chebyshev Expansions
On the other hand, the three-term recurrence relation satisfied by the Chebyshev polynomials
(3.20) can also be written in matrix form,
 



1
1
T0 (x)
 −2x
  T1 (x)   −x 
1
 



 1
  T2 (x)   0 
−2x
1
 



(3.102)

  T3 (x)  =  0  ,
1
−2x 1
 









.
.
.
.
.
..
..
..
..
  .. 


TN (x)
0
1 −2x 1
or
At = d,
(3.103)
where A is the (N + 1) × (N + 1) matrix of the coefficients of the recurrence relation and
d is the right-hand side vector of (3.102).
Let us now consider a vector bT = (b0 , b1 , . . . , bN ) such that
Then,
bT A = cT .
(3.104)
.
Sn = cT t = bT At = bT d = b0 − b1 x.
(3.105)
For SN , we have
1
1
1
SN = .
SN − c0 = (b0 − b1 x) − (b0 − 2xb1 + b2 ) = (b0 − b2 ).
2
2
2
(3.106)
The coefficients br can be computed using a recurrence relation if (3.104) is interpreted
as the corresponding matrix equation for the recurrence relation (and considering bN+1 =
bN+2 = 0). In this way,
br − 2xbr+1 + br+2 = cr ,
r = 0, 1, . . . , N.
(3.107)
The three-term recurrence relation is computed in the backward direction, starting from
r = N. With this, we arrive at
.
SN = xb1 − b2 + c0 ,
(3.108)
c0
.
2
(3.109)
SN = xb1 − b2 +
Error analysis
The expressions provided by (3.105) or (3.106), together with the use of (3.107), are simple
and avoid the explicit computation of the Chebyshev polynomials Tn (x) (with the exception
of T0 (x) = 1 and T1 (x) = x). However, these relations will be really useful if one can
be sure that the influence of error propagation when using (3.107) is small. Let us try to
quantify this influence by following the error analysis due to Elliott [62].
From “Numerical Methods for Special Functions” by Amparo Gil, Javier Segura, and Nico Temme
Copyright ©2007 by the Society for Industrial and Applied Mathematics
This electronic version is for personal use and may not be duplicated or distributed.
3.7. Evaluation of a Chebyshev sum
77
Let us denote by
/
Q = Q + δQ
(3.110)
an exact quantity Q computed approximately (δQ represents the absolute error in the computation).
From (3.107), we obtain

(3.111)
b̂n = cˆn + 2x̂b̂n+1 − b̂n+2 + rn ,
where rn is the roundoff error arising from rounding the quantity inside the brackets. We
can rewrite this expression as

b̂n = cˆn + 2xb̂n+1 − b̂n+2 + ηn + rn ,
(3.112)
where, neglecting the terms of order higher than 1 for the errors,
ηn = 2(δx)b̂n+1 ≈ 2(δx)bn+1 .
(3.113)
From (3.107), it is clear that δbn satisfies a recurrence relation of the form
Yn − 2xYn+1 + Yn+2 = δcn + ηn + rn ,
(3.114)
which is the same recurrence relation (with a different right-hand side) as that satisfied by
bn . It follows that
N

(δcn + ηn + rn )Tn (x).
(3.115)
δb0 − δb1 x =
n=0
On the other hand, because of (3.105), the computed .
SN will be given by

/
.
S N = b̂0 − b̂1 x̂ + s,
(3.116)
where s is the roundoff error arising from computing the expression inside the brackets.
Hence,
/
.
(3.117)
S N = (b0 − xb1 ) + (δb0 − x(δb1 )) − b1 (δx) + s,
and using (3.115) it follows that
δ.
SN = (δb0 − xδb1 ) − b1 δx + s =
N

(δcn + ηn + rn ) Tn (x) − b1 δx + s.
(3.118)
n=0
Let us rewrite this expression as
δ.
SN =
N

n=0
(δcn + rn ) Tn (x) + 2δx
N

bn+1 Tn (x) − b1 δx + s.
(3.119)
n=0
From “Numerical Methods for Special Functions” by Amparo Gil, Javier Segura, and Nico Temme
Copyright ©2007 by the Society for Industrial and Applied Mathematics
This electronic version is for personal use and may not be duplicated or distributed.
78
Chapter 3. Chebyshev Expansions
At this point, we can use the fact that the bn coefficients can be written in terms of the
Chebyshev polynomials of the second kind Un (x) as follows:
bn =
N

(3.120)
ck Uk−n (x).
k=n
We see that the term
N

N
n=0 bn+1 Tn (x) in (3.119) can be expressed as
bn+1 Tn (x) =
n=0
=
N

n=1
N

N
N

ck Uk−n (x) Tn−1 (x)
bn Tn−1 (x) =
n=1
k=n
k

ck
Uk−n (x)Tn−1 (x)
(3.121)
n=1
k=1
N

= 12
ck (k + 1)Uk−1 (x),
k=1
where, in the last step, we have used
k

sin(k − n + 1)θ cos(n − 1)θ = 12 (k + 1) sin kθ.
(3.122)
n=1
Substituting (3.121) in (3.119) it follows that
δ.
SN =
N

n=0
(δcn + rn ) Tn (x) + δx
N

ncn Un−1 (x) + s.
(3.123)
n=1
Since |Tn (x)| ≤ 1, |Un−1 (x)| ≤ n, and assuming that the local errors |δcn |, |rn |, |δx|,
|s| are quantities that are smaller than a given > 0, we have

N

2
δ.

SN ≤ (2N + 3) +
n |cn | .
(3.124)
n=1
In the case of a Chebyshev series where the coefficients are slowly convergent, the
second term on the right-hand side of (3.124) can provide a significant contribution to the
error.
On the other hand, when x is close to ±1, there is a risk of a growth of rounding errors,
and, in this case, a modification of Clenshaw’s method [164] seems to be more appropriate.
We describe this modification in the following algorithm.
Algorithm 3.2. Modified Clenshaw’s method for a Chebyshev sum.
Input: x; c0 , c1 , . . . , cN .

Output: .
SN (x) = N
k=0 ck Tk (x).
•
IF (x ≈ ±1) THEN
From “Numerical Methods for Special Functions” by Amparo Gil, Javier Segura, and Nico Temme
Copyright ©2007 by the Society for Industrial and Applied Mathematics
This electronic version is for personal use and may not be duplicated or distributed.
3.7. Evaluation of a Chebyshev sum
•
•
•
•
•
79
bN = cN ; dN = bN .
IF (x ≈ 1) THEN
DO r = N − 1, N − 2, . . . , 1:
dr = 2(x − 1)br+1 + dr+1 + cr .
br = dr + br+1 .
ELSEIF (x ≈ −1) THEN
DO r = N − 1, N − 2, . . . , 1:
dr = 2(x + 1)br+1 − dr+1 + cr .
br = dr − br+1 .
ENDIF
ELSE
Use Algorithm 3.1.
•
.
SN (x) = xb1 − b2 + c0 .
Oliver [165] has given a detailed analysis of Clenshaw’s method for evaluating a
Chebyshev sum, also by comparing it with other polynomial evaluation schemes for the
evaluation of (3.100); in addition, error bounds are derived. Let us consider two expressions
for (3.100),
.
SN (x) =
N

cn Tn (x),
(3.125)
dn x n ,
(3.126)
n=0
.
SN (x) =
N

n=0
and let ŜN (x) be the actually computed quantity, assuming that errors are introduced at each
stage of the computation process of .
SN (x) using (3.125) (considering Clenshaw’s method)
or (3.126) (considering Horner’s scheme). Then,
N

ρn (x)|cn |
ŜN (x) − SN (x) ≤
(3.127)
n=0
or
N

σn (x)|dn |,
ŜN (x) − SN (x) ≤
(3.128)
n=0
depending on the choice of the polynomial expression; is the accuracy parameter and
ρn and σn are error amplification factors. Reference [165] analyzed the variation of these
factors with x. Two conclusions of this study were that
• the accuracy of the methods of Clenshaw and Horner are sensitive to values of x and
the errors tend to reach their extreme at the endpoints of the interval;
• when a polynomial has coefficients of constant sign or strictly alternating sign, converting into the Chebyshev form does not improve upon the accuracy of the evaluation.
From “Numerical Methods for Special Functions” by Amparo Gil, Javier Segura, and Nico Temme
Copyright ©2007 by the Society for Industrial and Applied Mathematics
This electronic version is for personal use and may not be duplicated or distributed.
80
3.8
Chapter 3. Chebyshev Expansions
Economization of power series
Chebyshev polynomials also play a key role in the so-called economization of power series.
Suppose we have at our disposal a convergent Maclaurin series expansion for the evaluation
of a function f(x) in the interval [−1, 1]. Then, a plausible approximation to f(x) may be the
polynomial pn (x) of degree n, which is obtained by truncating the power series after n + 1
terms. It may be possible, however, to obtain a “better” nth-degree polynomial approximation. This is the idea of economization: it involves finding an alternative representation
for the function containing n + 1 parameters that possesses the same functional form as the
initial approximant. This alternative representation also incorporates information present
in the higher orders of the original power series to minimize the maximum error of the new
approximant over the range of x.
Let SN+1 denote
SN+1 =
N+1

ai x i ,
(3.129)
i=0
the original Maclaurin series for f(x) truncated at order N + 1. Then, one can obtain an
“economic” representation by subtracting from SN+1 a suitable polynomial PN+1 of the
same order such that the leading orders cancel. That is,
CN = SN+1 − PN+1 =
N

ai xi ,
(3.130)
i=0
where ai denotes the resulting expansion coefficient of xi . Obviously, the idea is to choose
PN+1 in such a way that the maximum error of the new Nth order series representation is
significantly reduced. Then, an optimal candidate is
PN+1 = aN+1
TN+1
.
2N
(3.131)
The maximum error of this new Nth order polynomial CN is nearly the same as the maximum
error of the (N + 1)th order polynomial SN+1 and considerably less than SN .
Of course, this procedure may be adapted to ranges different from [−1, 1] by using
Chebyshev polynomials adjusted to the required range. For example, the Chebyshev polynomials Tk (x/c) (or shifted Chebyshev polynomials Tk∗ (x/c)) should be used in the range
[−c, c] (or [0, c]).
3.9
Example: Computation of Airy functions of real
variable
The computation of the Airy functions of real variables Ai(x), Bi(x) and their derivatives
[186] for real values of the argument x is a useful example of application of Chebyshev
From “Numerical Methods for Special Functions” by Amparo Gil, Javier Segura, and Nico Temme
Copyright ©2007 by the Society for Industrial and Applied Mathematics
This electronic version is for personal use and may not be duplicated or distributed.
3.9. Example: Computation of Airy functions of real variable
81
expansions for computing special functions. As is usual, the real line is divided into a
number of intervals and we consider expansions on intervals containing the points −∞
or +∞. An important aspect is selecting the quantity that has to be expanded in terms
of Chebyshev polynomials. The Airy functions have oscillatory behavior at −∞, and
exponential behavior at +∞. It is important to expand quantities that are slowly varying in
the interval of interest.
When the argument of the Airy function is large, the asymptotic expansions given in
(10.4.59)–(10.4.64), (10.4.66), and (10.4.67) of [2] can be considered. The coefficients ck
and dk used in these expansions are given by
c0 = 1,
(3k + 12 )
ck =
54k k! (k + 12 )
6k + 1
dk = −
ck ,
6k − 1
d0 = 1,
k = 0, 1, 2, . . . ,
,
(3.132)
k = 1, 2, 3, . . . .
Asymptotic expansions including the point −∞
We have the representations
*
)
Ai(−z) = π−1/2 z−1/4 sin(ζ + 14 π) f(z) − cos(ζ + 14 π) g(z) ,
)
*
Ai (−z) = −π−1/2 z1/4 cos(ζ + 14 π) p(z) + sin(ζ + 14 π) q(z) ,
Bi(−z) = π
−1/2 −1/4
)
cos(ζ + 14 π) f(z) + sin(ζ + 14 π) g(z)
z
(3.133)
*
,
*
)
Bi (−z) = π−1/2 z1/4 sin(ζ + 14 π) p(z) − cos(ζ + 14 π) q(z) ,
where ζ = 23 z3/2 . The asymptotic expansions for the functions f(z), g(z), p(z), q(z) are
f(z) ∼
∞

k
(−1) c2k ζ
−2k
,
g(z) ∼
k=0
p(z) ∼
∞

∞

(−1)k c2k+1 ζ −2k−1 ,
k=0
k
(−1) d2k ζ
−2k
,
k=0
q(z) ∼
∞

(3.134)
k
(−1) d2k+1 ζ
−2k−1
,
k=0
as z → ∞, |ph z| < 23 π. Asymptotic expansions including the point +∞ Now we use the representations Ai (z) = 12 π−1/2 z−1/4 e−ζ f̃ (z), Ai (z) = − 12 π−1/2 z1/4 e−ζ p̃(z), Bi (z) = 12 π−1/2 z−1/4 eζ g̃(z), Bi (z) = 12 π−1/2 z1/4 eζ q̃(z), (3.135) From "Numerical Methods for Special Functions" by Amparo Gil, Javier Segura, and Nico Temme Copyright ©2007 by the Society for Industrial and Applied Mathematics This electronic version is for personal use and may not be duplicated or distributed. 82 Chapter 3. Chebyshev Expansions where the asymptotic expansions for the functions f̃ (z), g̃(z), p̃(z), q̃(z) are f̃ (z) ∼ g̃(z) ∼ ∞ (−1)k ck ζ −k , p̃(z) ∼ ∞ k=0 k=0 ∞ ∞ −k ck ζ , q̃(z) ∼ k=0 (−1)k dk ζ −k , (3.136) −k dk ζ , k=0 as z → ∞, with |ph z| < π (for f̃ (z) and p̃(z)) and |ph z| < 13 π (for g̃(z) and q̃(z)). Chebyshev expansions including the point −∞ The functions f(z), g(z), p(z), q(z) are the slowly varying quantities in the representations in (3.133), and these functions are computed approximately by using Chebyshev expansions. We write z = x. A possible selection is the x-interval [7, +∞) and for the shifted Chebyshev polynomials we take the argument t = (7/x)3 (the third power also arises in the expansions in (3.134)). We have to obtain the coefficients in the approximation f(x) ≈ p(x) ≈ m1 ar Tr∗ (t), m2 g(x) ≈ 1ζ r=0 r=0 m3 m4 cr Tr∗ (t), q(x) ≈ 1ζ r=0 br Tr∗ (t), (3.137) dr Tr∗ (t). r=0 Chebyshev expansions including the point +∞ Now we consider Chebyshev expansions for the functions f̃ (x), g̃(x), p̃(x), q̃(x), and we have f̃ (x) ≈ p̃(x) ≈ n1 ãr Tr∗ (t̃), g̃(x) ≈ n2 r=0 r=0 n3 n4 r=0 c̃r Tr∗ (t̃), q̃(x) ≈ (−1)r ãr Tr∗ (t̃), (3.138) (−1) c̃r Tr∗ (t̃), r r=0 where t̃ = (7/x)3/2 . The number of terms in the Chebyshev expansions mj , nj , j = 1, 2, 3, 4, are determined for a prescribed accuracy of the functions. The Chebyshev coefficients of functions defined by convergent power series can be computed by rearrangement of the corresponding power series expansions using (3.35). For power series that represent asymptotic expansions this method is not available. In the present case we have used the Maple program chebyshev with default accuracy of 10−10 . The first three coefficients of the Chebyshev expansions (3.137) and (3.138) are given in Table 3.3. For ãr and c̃r , see Table 3.4 for more—and more precise—values. From "Numerical Methods for Special Functions" by Amparo Gil, Javier Segura, and Nico Temme Copyright ©2007 by the Society for Industrial and Applied Mathematics This electronic version is for personal use and may not be duplicated or distributed. 3.10. Chebyshev expansions with coefficients in terms of special functions 83 Table 3.3. The first coefficients of the Chebyshev expansions (3.137) and (3.138). 3.10 Coef. r=0 r=1 r=2 ar br cr dr ãr c̃r 1.0001227513 0.0695710154 0.9998550220 −0.0973635868 0.9972733954 1.0038355798 0.0001230753 0.0001272386 −0.0001453297 −0.0001420772 −0.0026989587 0.0038027377 0.0000003270 0.0000006769 −0.0000003548 −0.0000007225 0.0000271274 0.0000322597 Chebyshev expansions with coefficients in terms of special functions As we have seen in §3.6.1, the coefficients in a Chebyshev expansion can be obtained from recurrence relations when the function satisfies a linear differential equation with polynomial coefficients. All special functions of hypergeometric type satisfy such a differential equation, and in §3.6.1 we have explained how the method works for the Bessel function J0 (4x) in the range −1 ≤ x ≤ 1. However, for this particular function we can obtain expansions in which the coefficients can be expressed in terms of known special functions, which in fact are again Bessel functions. We have (see [143, p. 37]) J0 (ax) = ∞ n (−1)n Jn2 (a/2) T2n (x), n=0 J1 (ax) = 2 ∞ (3.139) n (−1) Jn (a/2)Jn+1 (a/2) T2n+1 (x), n=0 where −1 ≤ x ≤ 1 and 0 = 1, n = 2 if n > 0. The parameter a can be any complex
number. Similar expansions are available for J-Bessel functions of any complex order, in
which the coefficients are 1 F2 -hypergeometric functions, and explicit recursion relations
are available for computing the coefficients. For general integer order, the coefficients are
products of two J-Bessel functions, as in (3.139). See again [143].
Another example is the expansion for the error function,
2 2
ea x erf (ax) =
∞

√ 1 a2
πe 2
In+ 1 12 a2 T2n+1 (x),
2
−1 ≤ x ≤ 1,
(3.140)
n=0
in which the modified Bessel function is used. Again, a can be any complex number.
The expansions in (3.139) and (3.140) can be viewed as expansions near the origin.
Other expansions are available that can be viewed as expansions at infinity, and these may
be considered as alternatives for asymptotic expansions of special functions. For example,
for the confluent hypergeometric U-functions we have the convergent expansion in terms
of shifted Chebyshev polynomials (see (3.36)):
(ωz)a U(a, c, ωz) =
∞

Cn (z)Tn∗ (1/ω), z = 0, |ph z| < 32 π, 1 ≤ ω ≤ ∞. (3.141) n=0 From "Numerical Methods for Special Functions" by Amparo Gil, Javier Segura, and Nico Temme Copyright ©2007 by the Society for Industrial and Applied Mathematics This electronic version is for personal use and may not be duplicated or distributed. 84 Chapter 3. Chebyshev Expansions Furthermore, a, 1 + a − c = 0, −1, −2, . . . . When equalities hold for these values of a and c, the Kummer U-function reduces to a Laguerre polynomial. This follows from U(a, c, z) = z1−c U(1 + a − c, 2 − c, z) (3.142) and U(−n, α + 1, z) = (−1)n n!Lαn (z), n = 0, 1, 2, . . . . (3.143) The expansion (3.141) is given in [143, p. 25]. The coefficients can be represented in terms of generalized hypergeometric functions, in fact, Meijer G-functions, and they can be computed from the recurrence relation 2Cn (z) = 2(n + 1)A1 Cn+1 (z) + A2 Cn+2 (z) + A3 Cn+3 (z), n (3.144) where b = a + 1 − c, 0 = 12 , n = 1 (n ≥ 1), and A1 = 1 − 2z (2n + 3)(n + a + 1)(n + b + 1) − , 2(n + 2)(n + a)(n + b) (n + a)(n + b) A2 = 1 − 2(n + 1)(2n + 3 − z) , (n + a)(n + b) A3 = − (3.145) (n + 1)(n + 3 − a)(n + 3 − b) . (n + 2)(n + a)(n + b) For applying the backward recursion algorithm it is important to know that ∞ (−1)n Cn (z) = 1, |ph z| < 32 π. (3.146) n=0 This follows from lim (ωz)a U(a, c, ωz) = 1 ω→∞ Tn∗ (0) = (−1)n . and (3.147) The standard backward recursion scheme (see Chapter 4) for computing the coefficients Cn (z) works only for |ph z| < π, and for ph z = ±π a modification seems to be possible; see [143, p. 26]. Although the expansion in (3.141) converges for all z = 0 in the indicated sector, it is better to avoid small values of the argument of the U-function. Luke gives an estimate of the coefficients Cn (z) of which the dominant factor that determines the speed of convergence is given by 2 1 Cn (z) = O n2(2a−c−1)/3 e−3n 3 z 3 , n → ∞, (3.148) and we see that large values of z1/3 improve the convergence. For example, if we denote ζ = ωz and we want to use the expansion for |ζ| ≥ R(> 0), we should choose z and ω
(ω ≥ 1) such that z and ζ have the same phase, say θ. We cannot choose z with modulus larger
than R, and an appropriate choice is z = Reiθ . Then the expansion gives an approximation
of the U-function on the half-ray with |ζ| ≥ R with phase θ, and the coefficients Cn (z) can
be used for all ζ on this half-ray. For a single evaluation we can take ω = 1.
From “Numerical Methods for Special Functions” by Amparo Gil, Javier Segura, and Nico Temme
Copyright ©2007 by the Society for Industrial and Applied Mathematics
This electronic version is for personal use and may not be duplicated or distributed.
3.10. Chebyshev expansions with coefficients in terms of special functions
85
Table 3.4. Coefficients of the Chebyshev expansion (3.153).
n
Cn (z)
Dn (z)
0
1
2
3
4
5
6
7
8
9
0.99727 33955 01425
−0.00269 89587 07030
0.00002 71274 84648
−0.00000 05043 54523
0.00000 00134 68935
−0.00000 00004 63150
0.00000 00000 19298
−0.00000 00000 00938
0.00000 00000 00052
−0.00000 00000 00003
1.00383 55796 57251
0.00380 27374 06686
−0.00003 22598 78104
0.00000 05671 25559
−0.00000 00147 27362
0.00000 00004 97977
−0.00000 00000 20517
0.00000 00000 00989
−0.00000 00000 00054
0.00000 00000 00003
The expansion in (3.141) can be used for all special cases of the Kummer U-function,
that is, for Bessel functions (Hankel functions and K-modified Bessel function), for the
incomplete gamma function (a, z), with special cases the complementary error function
and exponential integrals.
Example 3.16 (Airy function). For the Airy function Ai(x) we have the relations
√ 1 1
1
ξ 6 U( 16 , 13 , ξ) = 2 πx 4 e 2 ξ Ai(x),
√
1
1 1
ξ − 6 U(− 16 , − 13 , ξ) = −2 πx− 4 e 2 ξ Ai (x),
(3.149)
3
where ξ = 43 x 2 . For the expansions of the functions f̃ (x) and p̃(x) in (3.138) we take
ω = (x/7)3/2 and z = 43 73/2 = 24.69 . . . . To generate the coefficients with this value of z
we determine the smallest value of n for which the exponential factor in (3.148) is smaller
than 10−15 . This gives n = 8.
Next we generate for both U-functions in (3.149) the coefficients Cn (z) by using
(3.144) in the backward direction, with starting values
.
C19 (z) = 1,
.
C20 (z) = 0,
.
C21 (z) = 0.
(3.150)
We also compute for normalization
S=
18

Cn (z) = −0.902363242772764 1025 ,
(−1)n .
(3.151)
n=0
where the numerical value is for the Ai case. Finally we compute
Cn (z) = .
Cn (z)/S,
n = 0, 1, 2, . . . , 9.
(3.152)
From “Numerical Methods for Special Functions” by Amparo Gil, Javier Segura, and Nico Temme
Copyright ©2007 by the Society for Industrial and Applied Mathematics
This electronic version is for personal use and may not be duplicated or distributed.
86
Chapter 3. Chebyshev Expansions
This gives the coefficients Cn (z) of the expansions
9

√ 1 2 3/2
Cn Tn∗ (7/x)3/2 ,
2 πx 4 e 3 x Ai(x) ≈
n=0
√
1
4
−2 πx e
2 3/2
3x

Ai (x) ≈
9

Dn Tn∗

(7/x)
3/2

(3.153)
,
n=0
x ≥ 7, of which the coefficients are given in Table 3.4.
The values of the coefficients in the first three rows of Table 3.4 correspond approximately with the same as those of the coefficients ãr and c̃r in Table 3.3.

From “Numerical Methods for Special Functions” by Amparo Gil, Javier Segura, and Nico Temme
c 2004 Society for Industrial and Applied Mathematics

SIAM REVIEW
Vol. 46, No. 3, pp. 501–517
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 of polynomial interpolation.
Key words. barycentric formula, interpolation
AMS subject classiﬁcations. 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
a widespread 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 do not 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
of the 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
† Département de Mathématiques, Université de Fribourg, 1700 Fribourg/Pérolles, Switzerland
(Jean-Paul.Berrut@unifr.ch). The work of this author was supported by the Swiss National Science
Foundation, grant 21–59272.99.
‡ Computing Laboratory, 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 diﬀerent.
501
502
JEAN-PAUL BERRUT AND LLOYD N. TREFETHEN
complex plane. Let Πn denote the vector space of all polynomials of degree at most
n. The classical problem addressed here is that of ﬁnding 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 can be written in Lagrange form [44]:
(2.1)
p(x) =
n

n
fj j (x),
k=0, k=j (x − xk )
j (x) = n
k=0, k=j (xj − xk )
j=0
.
The Lagrange polynomial j corresponding to the node xj has the property

(2.2)
j (xk ) =
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 that certain 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 theoretical tool for proving theorems. For computations, it is generally recommended that
one should instead use Newton’s formula, which requires only O(n) ﬂops for each
evaluation of p once some numbers, which are independent of the evaluation point x,
have been computed.2
Newton’s approach consists of two steps. First, compute the Newton tableau of
divided diﬀerences
f [x0 ]
f [x0 , x1 ]
f [x1 ]
f [x0 , x1 , x2 ]
f [x1 , x2 ]
f [x2 ]
f [x2 , x3 ]
f [x3 ]
..
.
f [x0 , x1 , x2 , x3 ]
f [x1 , x2 , x3 ]
..
.
..
.
..
.
f [xn−3 , . . . , xn ]
f [xn−2 , xn−1 , xn ]
..
.
…
f [x0 , x1 , x2 , . . ., xn ]
f [xn−1 , xn ]
f [xn ]
2 We deﬁne a ﬂop—ﬂoating point operation—as a multiplication or division plus an addition or
subtraction. This is quite a simpliﬁcation of the reality of computer hardware, since divisions are
much slower than multiplications on many processors, but this matter is not of primary importance
for our discussion since we are mainly interested in the larger diﬀerences between O(n) and O(n2 )
or O(n2 ) and O(n3 ). As usual, O(np ) means ≤ cnp for some constant c.
BARYCENTRIC LAGRANGE INTERPOLATION
503
by the recursive formula
(2.3)
f [xj , xj+1 , . . . , xk−1 , xk ] =
f [xj+1 , . . . , xk ] − f [xj , . . . , xk−1 ]
xk − xj
with initial condition f [xj ] = fj . This requires about n2 subtractions and n2 /2
divisions. Second, for each x, evaluate p(x) by the Newton interpolation formula
p(x) = f [x0 ] + f [x0 , x1 ](x − x0 ) + f [x0 , x1 , x2 ](x − x0 )(x − x1 ) + · · ·
(2.4)
+ f [x0 , x1 , . . . , xn ](x − x0 )(x − x1 ) · · · (x − xn−1 ).
This requires a mere n ﬂops when carried out by nested multiplication, much less
than the O(n2 ) required for direct application of (2.1).
3. An Improved Lagrange Formula. The ﬁrst purpose of the present work is
to emphasize that the Lagrange formula (2.1) can be rewritten in such a way that it
too can be evaluated and updated in O(n) operations, just like its Newton counterpart (2.4). To this end, it suﬃces to note that the numerator of j in (2.1) can be
written as the quantity
(3.1)
(x) = (x − x0 )(x − x1 ) · · · (x − xn )
divided by x − xj . ( is denoted by Ω or ω in many texts.) If we deﬁne the barycentric
weights by3
(3.2)
wj =
1
,
k=j (xj − xk )
j = 0, . . . , n,
that is, wj = 1/ (xj ) [37, p. 243], we can thus write j as
j (x) = (x)
wj
.
x − xj
Now we note that all the terms of the sum in (2.1) contain the factor (x), which does
not depend on j. This factor may therefore be brought in front of the sum to yield
(3.3)
p(x) = (x)
n

wj
fj .
x − xj
j=0
That’s it! Now Lagrange interpolation too is a formula requiring O(n2 ) ﬂops for
calculating some quantities independent of x, the numbers wj , followed by O(n) ﬂops
for evaluating p once these numbers are known. Rutishauser [51] called (3.3) the “ﬁrst
form of the barycentric interpolation formula.”
What about updating? A look at (3.2) shows that incorporating a new node xn+1
entails two calculations:
• Divide each wj , j = 0, . . . , n, by xj − xn+1 (one ﬂop for each point), for a
cost of n + 1 ﬂops.
• Compute wn+1 with formula (3.2), for another n + 1 ﬂops.
3 In practice, in view of the remarks of the previous footnote, one would normally choose to
implement this formula by n multiplications followed by a single division.
504
JEAN-PAUL BERRUT AND LLOYD N. TREFETHEN
The Lagrange interpolant can thus also be updated with O(n) ﬂops! By applying this
update recursively while minimizing the number of divisions, we get the following
(n)
algorithm for computing wj = wj , j = 0, . . . , n:
(0)
w0 = 1
for j = 1 to n do
for k = 0 to j − 1 do
(j)
(j−1)
wk = (xk − xj )wk
end
(j)
j−1
wj = k=0 (xj − xk )
end
for j = 0 to n do
(j)
(j)
wj = 1/wj
end
This algorithm performs the same operations as (3.2), only in a diﬀerent order.4
A remarkable advantage of Lagrange over Newton interpolation, rarely mentioned
in the literature, is that the quantities that have to be computed in O(n2 ) operations
do not depend on the data fj . This feature permits the interpolation of as many
functions as desired in O(n) operations each once the weights wj are known, whereas
Newton interpolation requires the recomputation of the divided diﬀerence tableau for
each new function.
Another advantage of the Lagrange formula is that it does not depend on the order
in which the nodes are arranged. In the Newton formula, the divided diﬀerences do
have such a dependence. This seems odd aesthetically, and it has a computational
consequence: for larger values of n, many orderings lead to numerical instability. For
stability, it is necessary to select the points in a van der Corput or Leja sequence or in
another related sequence based on a certain equidistribution property [18, 23, 61, 67].
This is not to say that Newton interpolation has no advantages. One is that it
leads to elegant methods for incorporating information on derivatives f (p) (xj ) when
solving so-called Hermite interpolation problems. Another may be its ability to accommodate vector or matrix interpolation from scalar data [60]; using a representation
like (3.3) would require costly matrix inversion. There are also some theoretical instances in which Newton interpolation is preferable, e.g., the construction of multistep
formulas for the solution of ordinary diﬀerential equations [33, 45].
4. The Barycentric Formula. Equation (3.3) is not the end of the story: it can
be modiﬁed to an even more elegant formula, the one that is often used in practice.
Suppose we interpolate, besides the data fj , the constant function 1, whose interpolant is of course itself. Inserting into (3.3), we get
(4.1)
1 =
n

j=0
j (x) = (x)
n

wj
.
x
−
xj
j=0
Dividing (3.3) by this expression and cancelling the common factor (x), we obtain
4 This contrasts with a somewhat faster but unstable [28, p. 96] algorithm that takes advantage
of the relation
j
k=0
(j)
wk
= 0 [67].
BARYCENTRIC LAGRANGE INTERPOLATION
505
the barycentric formula for p:
n

(4.2)
p(x) =
wj
fj
x
−
xj
j=0
n

wj
x
−
xj
j=0
BARYCENTRIC
FORMULA
where wj is still deﬁned by (3.2). Rutishauser [51] called (4.2) the “second (true) form
of the barycentric formula.” See p. 229 of [39] for comments on this terminology.
We see that the barycentric formula is a Lagrange formula, but one with a special
and beautiful symmetry. The weights wj appear in the denominator exactly as in the
numerator, except without the data factors fj . This means that any common factor
in all the weights wj may be cancelled without aﬀecting the value of pn , and in the
discussion to follow we will use this freedom.
Like (3.3), (4.2) can also take advantage of the updating of the weights wj in
O(n) ﬂops to incorporate a new data pair (xn+1 , fn+1 ).
5. Chebyshev and Other Point Distributions. For certain special sets of nodes
xj , one can give explicit formulas for the barycentric weights wj . The obvious place
to start is equidistant nodes with spacing h = 2/n on the
interval [−1, 1]. Here the

weights can be directly calculated to be wj = (−1)n−j nj (hn n!) [55], which after
cancelling the factors independent of j yields

n
(5.1)
.
wj = (−1)j
j
For the interval [a, b] we would multiply the original formula for wj by 2n (b − a)−n ,
but this constant factor too can be dropped, so we end up with (5.1) again, regardless
of a and b.
A glance at (5.1) reveals that if n is large, the weights wj for equispaced barycentric interpolation vary by exponentially large factors, of order approximately 2n . This
sounds dangerous, and it is. The eﬀect will be that even small data near the center of
the interval are associated with large oscillations in the interpolant, on the order of 2n
times bigger, near the edge of the interval [40, 65]. This so-called Runge phenomenon
is not a problem with the barycentric formula, but is intrinsic in the underlying interpolation problem. Among other things it implies that polynomial interpolation in
equally spaced points is highly ill-conditioned: small changes in the data may cause
huge changes in the interpolant.
For polynomial interpolation to be a well-conditioned process, unless n is rather
small, one must dispense with equally spaced points [59, Thm. 6.21.3]. As is well
known in approximation theory, the right approach is to use point sets that are clustered at the endpoints of the interval with an asymptotic density proportional to
(1 − x2 )−1/2 as n → ∞. Remarkably, this is the same asymptotic density one gets if
[−1, 1] is interpreted as a conducting wire and the points xj are interpreted as point
charges that repel one another with an inverse-linear force and that are allowed to
move along the wire to ﬁnd their equilibrium conﬁguration [64, Chap. 5]. And it
is precisely the same asymptotic density that is required to make the weights wj of
comparable scale in the sense that although they may not all be exactly equal, they
do not vary by factors exponentially large in n.
506
JEAN-PAUL BERRUT AND LLOYD N. TREFETHEN
The simplest examples of clustered point sets are the families of Chebyshev points,
obtained by projecting equally spaced points on the unit circle down to the unit
interval [−1, 1]. Four standard varieties of such points have been deﬁned, and for each,
there is an explicit formula for in (3.1) that may be easily diﬀerentiated [9, 35, 47].
From the identity
(5.2)
wj =
1
(xj )
mentioned earlier, we accordingly obtain explicit formulas for the weights wj .
The Chebyshev points of the ﬁrst kind are given by
xj = cos
(2j + 1)π
,
2n + 2
j = 0, . . . , n.
In this case after cancelling factors independent of j we ﬁnd [40, p. 249]
wj = (−1)j sin
(5.3)
(2j + 1)π
.
2n + 2
Note that these numbers vary by factors O(n), not exponentially, reﬂecting the good
distribution of the points. The Chebyshev points of the second kind are given by
xj = cos
jπ
,
n
j = 0, . . . , n.
Here we ﬁnd [52]

(5.4)
j
wj = (−1) δj ,
δj =
1/2, j = 0 or j = n,
1,
otherwise;
all but two of the weights are exactly equal. Formulas for the Chebyshev points of
the third and fourth kinds can be found in [9].
For all of these sets of Chebyshev points, if the interval [−1, 1] is linearly transformed to [a, b], the weights as deﬁned by (3.2) all get multiplied by 2n (b − a)−n .
However, as this factor cancels out in the barycentric formula, there is again no need
to include it (and it is safer not to—see section 7).
One sees that, with equidistant or Chebyshev points, no expensive computations
are needed to get the weights wj , and thus only O(n) operations are required for
evaluating pn . No other interpolation method seems to achieve this, for as mentioned
in section 3, Newton interpolation always requires O(n2 ) operations for the divided
diﬀerences.
Here is a Matlab code segment that samples the function f (x) = |x| + x/2 − x2
in 1001 Chebyshev points of the second kind and evaluates the barycentric interpolant
in 5000 points. One would rarely use so many points in practice; we pick such large
values just to illustrate the eﬀectiveness of the formula. (At the end of section 7 we
shall modify this program to avoid failure when x = xi .)
n = 1000;
fun = inline(’abs(x)+.5*x-x.ˆ2’);
x = cos(pi*(0:n)’/n);
f = fun(x);
c = [1/2; ones(n-1,1); 1/2].*(-1).ˆ((0:n)’);
507
BARYCENTRIC LAGRANGE INTERPOLATION
n=20
0.5
0
−0.5
−1
−0.5
0
0.5
1
−0.5
0
0.5
1
n=100
0.5
0
−0.5
−1
Fig. 5.1
Barycentric interpolation of the function f (x) = |x| + x/2 − x2 in 21 and 101 Chebyshev
points of the second kind on [−1, 1]. The dots mark the interpolated values fj .
xx = linspace(-1,1,5000)’;
numer = zeros(size(xx));
denom = zeros(size(xx));
for j = 1:n+1
xdiff = xx-x(j);
temp = c(j)./xdiff;
numer = numer + temp*f(j);
denom = denom + temp;
end
ff = numer./denom;
plot(x,f,’.’,xx,ff,’-’)
On one of our workstations this code runs in about 1 second, producing a plot like
those shown in Figure 5.1 but for the case n = 1000. For further numerical examples
of barycentric interpolation, see [6] and [40].
Other point sets with the asymptotic distribution (1 − x2 )−1/2 also lead to wellconditioned polynomial approximation, notably the Legendre points—zeros or extrema
of the Legendre polynomials. However, explicit formulas for the weights wj are not
known for Legendre points.
6. Convergence Rates for Smooth Functions. If a smooth function f deﬁned
on the interval [−1, 1] is interpolated by polynomials in Chebyshev points or in any
other system of points with the asymptotic density (1 − x2 )−1/2 , the rate of convergence of the interpolants to f as n → ∞ is remarkably fast. In particular, suppose
that f is a function that can be analytically continued to a function f (z) that is
analytic (holomorphic) in a neighborhood of [−1, 1] in the complex plane. Then the
508
JEAN-PAUL BERRUT AND LLOYD N. TREFETHEN
0
10
−4
error
10
−8
10
−12
10
−16
10
0
50
100
150
200
n
Fig. 6.1
Convergence of the error (6.1) in polynomial interpolation of two smooth functions f
in Chebyshev points of the second kind in [−1, 1]. The lower curve corresponds to f (x) =
exp(x)/ cos(x), and the upper curve to f (x) = (1+16×2 )−1 . In both cases we observe steady
convergence down to the level of rounding errors. The wiggles in the second function arise
from the missing central node for odd n.
interpolants pn satisfy an error estimate
(6.1)
max |f (x) − pn (x)| ≤ CK −n
x∈[−1,1]
for some constants C and K > 1 [24, p. 173]. The bigger the region of analyticity,
the bigger we may take K to be. To be precise, if f is analytic on and inside an
ellipse in the complex plane with foci ±1 and axis lengths 2L and 2, then we may
take K = L + [64, Chap. 5].
Figure 6.1 illustrates this fast convergence for two smooth functions. The convergence rate for f (x) = exp(x)/ cos(x) is determined by the poles of f (z) at z = ±π/2;
we have
π
+ π 2 − 1 ≈ 2.7822.
K =
2
The convergence rate for f (x) = (1 + 16×2 )−1 is determined by the poles of this
function at z = ±i/4; we have
K =
1
+
4
17
≈ 1.2808.
16
These facts about convergence pertain to polynomial interpolation in general, not
barycentric interpolation per se, but the latter has a special feature that might be
well suited to taking advantage of such convergence rates under certain circumstances.
Suppose we apply the barycentric formula to Chebyshev points associated with values
n = 2, 4, 8, 16, . . . . Then it is possible to compute the interpolants recursively, reusing
the previous function values fj each time n is doubled [39]. If f is analytic in a
neighborhood of [−1, 1], each doubling of n results approximately in a squaring of the
error; one has quadratic convergence of the overall process.
A function f may be less smooth than the class we have considered (e.g., twice
continuously diﬀerentiable, or inﬁnitely diﬀerentiable); or it may be smoother (e.g.,
BARYCENTRIC LAGRANGE INTERPOLATION
509
analytic in the entire complex plane). A great deal is known about exactly how
the convergence rates of polynomial interpolants depend on the precise degree of
smoothness of f . See Chapter 5 of [64] for an introduction and [58] for a more
advanced treatment.
7. Numerical Stability. Polynomial interpolation is an area in which the errors
introduced by computer arithmetic are often a problem. One of the major advantages
of the barycentric formula is its good behavior in this respect. But this is a subject
where confusion is widespread, and one must be sure to distinguish several diﬀerent
phenomena.
One phenomenon to be clear about is that for improperly distributed sets of nodes,
as discussed in section 5, the underlying interpolation problem is ill-conditioned: small
changes in the data may cause large changes in the interpolant, typically manifested
as large oscillations near the endpoints. Unless the nodes converge to the (1 − x2 )−1/2
distribution, the condition numbers grow exponentially as n → ∞. In such circumstances polynomial interpolation is usually not a good method for applications,
regardless of whether one uses a Lagrange, Newton, or any other formulation. Even
though the interpolants may converge to f in theory as n → ∞ when f is suﬃciently
smooth, rounding errors will destroy the convergence in practice, since the rounding
errors are not smooth and thus get multiplied by exponentially large factors.
Let us suppose, then, that we are working with a set of points that are clustered
in the above sense, such as Chebyshev or Legendre points. The barycentric formula
is now stable provided that two matters described below are attended to. For years
this conclusion has been “folklore” in certain circles, and a rigorous analysis will be
provided in a forthcoming paper by N. J. Higham that was motivated by reading
an early draft of the present article [41]. (In a word, Higham shows that (3.3) is
unconditionally stable and that (4.2) is stable too, provided one has a clustered set
of interpolation points.)
The ﬁrst matter is the question of underﬂow and overﬂow in the computation of
the weights in cases where we are working from the general formula (3.2) rather than
from well-behaved explicit expressions like (5.3) and (5.4). In most cases diﬃculties
can be avoided as follows. Suppose we are working on an interval [a, b] of length 4C.
(In the ﬁeld of complex analysis, C is called t…

Turn in your highest-quality paper
Get a qualified writer to help you with

“ I need help with … ”

Get high-quality paper

Guarantee! All work is written by expert writers!

Still stressed from student homework?

Get quality assistance from academic writers!

Order now