You are on page 1of 38

# D.

Levy
3 Interpolation
3.1 What is Interpolation?
Imagine that there is an unknown function f(x) for which someone supplies you with
its (exact) values at (n+1) distinct points x
0
< x
1
< · · · < x
n
, i.e., f(x
0
), . . . , f(x
n
) are
given. The interpolation problem is to construct a function Q(x) that passes through
these points, i.e., to ﬁnd a function Q(x) such that the interpolation requirements
Q(x
j
) = f(x
j
), 0 j n, (3.1)
are satisﬁed (see Figure 3.1). One easy way of obtaining such a function, is to connect the
given points with straight lines. While this is a legitimate solution of the interpolation
problem, usually (though not always) we are interested in a diﬀerent kind of a solution,
e.g., a smoother function. We therefore always specify a certain class of functions from
which we would like to ﬁnd one that solves the interpolation problem. For example,
we may look for a function Q(x) that is a polynomial, Q(x). Alternatively, the function
Q(x) can be a trigonometric function or a piecewise-smooth polynomial, and so on.
x
0
x
1
x
2
f(x
0
)
f(x
1
)
f(x
2
)
f(x)
Q(x)
Figure 3.1: The function f(x), the interpolation points x
0
, x
1
, x
2
, and the interpolating
polynomial Q(x)
As a simple example let’s consider values of a function that are prescribed at two
points: (x
0
, f(x
0
)) and (x
1
, f(x
1
)). There are inﬁnitely many functions that pass through
these two points. However, if we limit ourselves to polynomials of degree less than or
equal to one, there is only one such function that passes through these two points: the
1
3.2 The Interpolation Problem D. Levy
line that connects them. A line, in general, is a polynomial of degree one, but if the
two given values are equal, f(x
0
) = f(x
1
), the line that connects them is the constant
Q
0
(x) ≡ f(x
0
), which is a polynomial of degree zero. This is why we say that there is a
unique polynomial of degree 1 that connects these two points (and not “a polynomial
of degree 1”).
The points x
0
, . . . , x
n
are called the interpolation points. The property of “passing
through these points” is referred to as interpolating the data. The function that
interpolates the data is an interpolant or an interpolating polynomial (or whatever
function is being used).
There are cases were the interpolation problem has no solution, e.g., if we look for a
linear polynomial that interpolates three points that do not lie on a straight line. When
a solution exists, it can be unique (a linear polynomial and two points), or the problem
can have more than one solution (a quadratic polynomial and two points). What we are
going to study in this section is precisely how to distinguish between these cases. We
are also going to present diﬀerent approaches to constructing the interpolant.
Other than agreeing at the interpolation points, the interpolant Q(x) and the under-
lying function f(x) are generally diﬀerent. The interpolation error is a measure on
how diﬀerent these two functions are. We will study ways of estimating the interpolation
error. We will also discuss strategies on how to minimize this error.
It is important to note that it is possible to formulate the interpolation problem
without referring to (or even assuming the existence of) any underlying function f(x).
For example, you may have a list of interpolation points x
0
, . . . , x
n
, and data that is
experimentally collected at these points, y
0
, y
1
, . . . , y
n
, which you would like to interpo-
late. The solution to this interpolation problem is identical to the one where the values
are taken from an underlying function.
3.2 The Interpolation Problem
We begin our study with the problem of polynomial interpolation: Given n + 1
distinct points x
0
, . . . , x
n
, we seek a polynomial Q
n
(x) of the lowest degree such that
the following interpolation conditions are satisﬁed:
Q
n
(x
j
) = f(x
j
), j = 0, . . . , n. (3.2)
Note that we do not assume any ordering between the points x
0
, . . . , x
n
, as such an order
will make no diﬀerence. If we do not limit the degree of the interpolation polynomial
it is easy to see that there any inﬁnitely many polynomials that interpolate the data.
However, limiting the degree of Q
n
(x) to be deg(Q
n
(x)) n, singles out precisely one
interpolant that will do the job. For example, if n = 1, there are inﬁnitely many
polynomials that interpolate (x
0
, f(x
0
)) and (x
1
, f(x
1
)). However, there is only one
polynomial Q
n
(x) with deg(Q
n
(x)) 1 that does the job. This result is formally stated
in the following theorem:
2
D. Levy 3.2 The Interpolation Problem
Theorem 3.1 If x
0
, . . . , x
n
∈ R are distinct, then for any f(x
0
), . . . f(x
n
) there exists a
unique polynomial Q
n
(x) of degree n such that the interpolation conditions (3.2) are
satisﬁed.
Proof. We start with the existence part and prove the result by induction. For n = 0,
Q
0
= f(x
0
). Suppose that Q
n−1
is a polynomial of degree n −1, and suppose also
that
Q
n−1
(x
j
) = f(x
j
), 0 j n −1.
Let us now construct from Q
n−1
(x) a new polynomial, Q
n
(x), in the following way:
Q
n
(x) = Q
n−1
(x) + c(x −x
0
) · . . . · (x −x
n−1
). (3.3)
The constant c in (3.3) is yet to be determined. Clearly, the construction of Q
n
(x)
implies that deg(Q
n
(x)) n. (Since we might end up with c = 0, Q
n
(x) could actually
be of degree that is less than n.) In addition, the polynomial Q
n
(x) satisﬁes the
interpolation requirements Q
n
(x
j
) = f(x
j
) for 0 j n −1. All that remains is to
determine the constant c in such a way that the last interpolation condition,
Q
n
(x
n
) = f(x
n
), is satisﬁed, i.e.,
Q
n
(x
n
) = Q
n−1
(x
n
) + c(x
n
−x
0
) · . . . · · · (x
n
−x
n−1
). (3.4)
The condition (3.4) implies that c should be deﬁned as
c =
f(x
n
) −Q
n−1
(x
n
)
n−1

j=0
(x
n
−x
j
)
, (3.5)
and we are done with the proof of existence.
As for uniqueness, suppose that there are two polynomials Q
n
(x), P
n
(x) of degree n
that satisfy the interpolation conditions (3.2). Deﬁne a polynomial H
n
(x) as the
diﬀerence
H
n
(x) = Q
n
(x) −P
n
(x).
The degree of H
n
(x) is at most n which means that it can have at most n zeros (unless
it is identically zero). However, since both Q
n
(x) and P
n
(x) satisfy all the
interpolation requirements (3.2), we have
H
n
(x
j
) = (Q
n
−P
n
)(x
j
) = 0, 0 j n,
which means that H
n
(x) has n + 1 distinct zeros. This contradiction can be resolved
only if H
n
(x) is the zero polynomial, i.e.,
P
n
(x) ≡ Q
n
(x),
and uniqueness is established.

3
3.3 Newton’s Form of the Interpolation Polynomial D. Levy
3.3 Newton’s Form of the Interpolation Polynomial
One good thing about the proof of Theorem 3.1 is that it is constructive. In other
words, we can use the proof to write down a formula for the interpolation polynomial.
We follow the procedure given by (3.4) for reconstructing the interpolation polynomial.
We do it in the following way:
• Let
Q
0
(x) = a
0
,
where a
0
= f(x
0
).
• Let
Q
1
(x) = a
0
+ a
1
(x −x
0
).
Following (3.5) we have
a
1
=
f(x
1
) −Q
0
(x
1
)
x
1
−x
0
=
f(x
1
) −f(x
0
)
x
1
−x
0
.
We note that Q
1
(x) is nothing but the straight line connecting the two points
(x
0
, f(x
0
)) and (x
1
, f(x
1
)).
In general, let
Q
n
(x) = a
0
+ a
1
(x −x
0
) + . . . + a
n
(x −x
0
) · . . . · (x −x
n−1
) (3.6)
= a
0
+
n

j=1
a
j
j−1

k=0
(x −x
k
).
The coeﬃcients a
j
in (3.6) are given by
a
0
= f(x
0
),
a
j
=
f(x
j
) −Q
j−1
(x
j
)

j−1
k=0
(x
j
−x
k
)
, 1 j n.
(3.7)
We refer to the interpolation polynomial when written in the form (3.6)–(3.7) as the
Newton form of the interpolation polynomial. As we shall see below, there are
various ways of writing the interpolation polynomial. The uniqueness of the interpola-
tion polynomial as guaranteed by Theorem 3.1 implies that we will only be rewriting
the same polynomial in diﬀerent ways.
Example 3.2
The Newton form of the polynomial that interpolates (x
0
, f(x
0
)) and (x
1
, f(x
1
)) is
Q
1
(x) = f(x
0
) +
f(x
1
) −f(x
0
)
x
1
−x
0
(x −x
0
).
4
D. Levy 3.4 The Interpolation Problem and the Vandermonde Determinant
Example 3.3
The Newton form of the polynomial that interpolates the three points (x
0
, f(x
0
)),
(x
1
, f(x
1
)), and (x
2
, f(x
2
)) is
Q
2
(x) = f(x
0
)+
f(x
1
) −f(x
0
)
x
1
−x
0
(x−x
0
)+
f(x
2
) −
_
f(x
0
) +
f(x
1
)−f(x
0
)
x
1
−x
0
(x
2
−x
0
)
_
(x
2
−x
0
)(x
2
−x
1
)
(x−x
0
)(x−x
1
).
3.4 The Interpolation Problem and the Vandermonde Deter-
minant
An alternative approach to the interpolation problem is to consider directly a polynomial
of the form
Q
n
(x) =
n

k=0
b
k
x
k
, (3.8)
and require that the following interpolation conditions are satisﬁed
Q
n
(x
j
) = f(x
j
), 0 j n. (3.9)
In view of Theorem 3.1 we already know that this problem has a unique solution, so we
should be able to compute the coeﬃcients of the polynomial directly from (3.8). Indeed,
the interpolation conditions, (3.9), imply that the following equations should hold:
b
0
+ b
1
x
j
+ . . . + b
n
x
n
j
= f(x
j
), j = 0, . . . , n. (3.10)
In matrix form, (3.10) can be rewritten as
_
_
_
_
_
1 x
0
. . . x
n
0
1 x
1
. . . x
n
1
.
.
.
.
.
.
.
.
.
1 x
n
. . . x
n
n
_
_
_
_
_
_
_
_
_
_
b
0
b
1
.
.
.
b
n
_
_
_
_
_
=
_
_
_
_
_
f(x
0
)
f(x
1
)
.
.
.
f(x
n
)
_
_
_
_
_
. (3.11)
In order for the system (3.11) to have a unique solution, it has to be nonsingular.
This means, e.g., that the determinant of its coeﬃcients matrix must not vanish, i.e.
¸
¸
¸
¸
¸
¸
¸
¸
¸
1 x
0
. . . x
n
0
1 x
1
. . . x
n
1
.
.
.
.
.
.
.
.
.
1 x
n
. . . x
n
n
¸
¸
¸
¸
¸
¸
¸
¸
¸
= 0. (3.12)
The determinant (3.12), is known as the Vandermonde determinant. In Lemma 3.4
we will show that the Vandermonde determinant equals to the product of terms of the
form x
i
− x
j
for i > j. Since we assume that the points x
0
, . . . , x
n
are distinct, the
determinant in (3.12) is indeed non zero. Hence, the system (3.11) has a solution that
is also unique, which conﬁrms what we already know according to Theorem 3.1.
5
3.4 The Interpolation Problem and the Vandermonde Determinant D. Levy
Lemma 3.4
¸
¸
¸
¸
¸
¸
¸
¸
¸
1 x
0
. . . x
n
0
1 x
1
. . . x
n
1
.
.
.
.
.
.
.
.
.
1 x
n
. . . x
n
n
¸
¸
¸
¸
¸
¸
¸
¸
¸
=

i>j
(x
i
−x
j
). (3.13)
Proof. We will prove (3.13) by induction. First we verify that the result holds in the
2 ×2 case. Indeed,
¸
¸
¸
¸
1 x
0
1 x
1
¸
¸
¸
¸
= x
1
−x
0
.
We now assume that the result holds for n −1 and consider n. We note that the index
n corresponds to a matrix of dimensions (n + 1) ×(n + 1), hence our induction
assumption is that (3.13) holds for any Vandermonde determinant of dimension n ×n.
We subtract the ﬁrst row from all other rows, and expand the determinant along the
ﬁrst column:
¸
¸
¸
¸
¸
¸
¸
¸
¸
1 x
0
. . . x
n
0
1 x
1
. . . x
n
1
.
.
.
.
.
.
.
.
.
1 x
n
. . . x
n
n
¸
¸
¸
¸
¸
¸
¸
¸
¸
=
¸
¸
¸
¸
¸
¸
¸
¸
¸
1 x
0
. . . x
n
0
0 x
1
−x
0
. . . x
n
1
−x
n
0
.
.
.
.
.
.
.
.
.
0 x
n
−x
0
. . . x
n
n
−x
n
0
¸
¸
¸
¸
¸
¸
¸
¸
¸
=
¸
¸
¸
¸
¸
¸
¸
x
1
−x
0
. . . x
n
1
−x
n
0
.
.
.
.
.
.
x
n
−x
0
. . . x
n
n
−x
n
0
¸
¸
¸
¸
¸
¸
¸
For every row k we factor out a term x
k
−x
0
:
¸
¸
¸
¸
¸
¸
¸
x
1
−x
0
. . . x
n
1
−x
n
0
.
.
.
.
.
.
x
n
−x
0
. . . x
n
n
−x
n
0
¸
¸
¸
¸
¸
¸
¸
=
n

k=1
(x
k
−x
0
)
¸
¸
¸
¸
¸
¸
¸
¸
¸
¸
¸
¸
¸
¸
¸
¸
¸
¸
1 x
1
+ x
0
. . .
n−1

i=0
x
n−1−i
1
x
i
0
1 x
2
+ x
0
. . .
n−1

i=0
x
n−1−i
2
x
i
0
.
.
.
.
.
.
.
.
.
1 x
n
+ x
0
. . .
n−1

i=0
x
n−1−i
n
x
i
0
¸
¸
¸
¸
¸
¸
¸
¸
¸
¸
¸
¸
¸
¸
¸
¸
¸
¸
Here, we used the expansion
x
n
1
−x
n
0
= (x
1
−x
0
)(x
n−1
1
+ x
n−2
1
x
0
+ x
n−3
1
x
2
0
+ . . . + x
n−1
0
),
for the ﬁrst row, and similar expansions for all other rows. For every column l, starting
from the second one, subtracting the sum of x
i
0
times column i (summing only over
“previous” columns, i.e., columns i with i < l), we end up with
n

k=1
(x
k
−x
0
)
¸
¸
¸
¸
¸
¸
¸
¸
¸
1 x
1
. . . x
n−1
1
1 x
2
. . . x
n−1
2
.
.
.
.
.
.
.
.
.
1 x
n
. . . x
n−1
n
¸
¸
¸
¸
¸
¸
¸
¸
¸
. (3.14)
6
D. Levy 3.5 The Lagrange Form of the Interpolation Polynomial
Since now we have on the RHS of (3.14) a Vandermonde determinant of dimension n×n,
we can use the induction to conclude with the desired result.

3.5 The Lagrange Form of the Interpolation Polynomial
The form of the interpolation polynomial that we used in (3.8) assumed a linear com-
bination of polynomials of degrees 0, . . . , n, in which the coeﬃcients were unknown. In
this section we take a diﬀerent approach and assume that the interpolation polyno-
mial is given as a linear combination of n + 1 polynomials of degree n. This time, we
set the coeﬃcients as the interpolated values, {f(x
j
)}
n
j=0
, while the unknowns are the
polynomials. We thus let
Q
n
(x) =
n

j=0
f(x
j
)l
n
j
(x), (3.15)
where l
n
j
(x) are n+1 polynomials of degree n. We use two indices in these polynomials:
the subscript j enumerates l
n
j
(x) from 0 to n and the superscript n is used to remind
us that the degree of l
n
j
(x) is n. Note that in this particular case, the polynomials l
n
j
(x)
are precisely of degree n (and not n). However, Q
n
(x), given by (3.15) may have a
lower degree. In either case, the degree of Q
n
(x) is n at the most. We now require that
Q
n
(x) satisﬁes the interpolation conditions
Q
n
(x
i
) = f(x
i
), 0 i n. (3.16)
By substituting x
i
for x in (3.15) we have
Q
n
(x
i
) =
n

j=0
f(x
j
)l
n
j
(x
i
), 0 i n.
In view of (3.16) we may conclude that l
n
j
(x) must satisfy
l
n
j
(x
i
) = δ
ij
, i, j = 0, . . . , n, (3.17)
where δ
ij
is the Kr¨onecker delta, deﬁned as
δ
ij
=
_
1, i = j,
0, i = j.
Each polynomial l
n
j
(x) has n + 1 unknown coeﬃcients. The conditions (3.17) provide
exactly n+1 equations that the polynomials l
n
j
(x) must satisfy and these equations can
be solved in order to determine all l
n
j
(x)’s. Fortunately there is a shortcut. An obvious
way of constructing polynomials l
n
j
(x) of degree n that satisfy (3.17) is the following:
l
n
j
(x) =
(x −x
0
) · . . . · (x −x
j−1
)(x −x
j+1
) · . . . · (x −x
n
)
(x
j
−x
0
) · . . . · (x
j
−x
j−1
)(x
j
−x
j+1
) · . . . · (x
j
−x
n
)
, 0 j n. (3.18)
7
3.5 The Lagrange Form of the Interpolation Polynomial D. Levy
The uniqueness of the interpolating polynomial of degree n given n + 1 distinct
interpolation points implies that the polynomials l
n
j
(x) given by (3.17) are the only
polynomials of degree n that satisfy (3.17).
Note that the denominator in (3.18) does not vanish since we assume that all inter-
polation points are distinct. The Lagrange form of the interpolation polynomial
is the polynomial Q
n
(x) given by (3.15), where the polynomials l
n
j
(x) of degree n are
given by (3.18). A compact form of rewriting (3.18) using the product notation is
l
n
j
(x) =
n

i=0
i=j
(x −x
i
)
n

i=0
i=j
(x
j
−x
i
)
, j = 0, . . . , n. (3.19)
Example 3.5
We are interested in ﬁnding the Lagrange form of the interpolation polynomial that
interpolates two points: (x
0
, f(x
0
)) and (x
1
, f(x
1
)). We know that the unique interpola-
tion polynomial through these two points is the line that connects the two points. Such
a line can be written in many diﬀerent forms. In order to obtain the Lagrange form we
let
l
1
0
(x) =
x −x
1
x
0
−x
1
, l
1
1
(x) =
x −x
0
x
1
−x
0
.
The desired polynomial is therefore given by the familiar formula
Q
1
(x) = f(x
0
)l
1
0
(x) + f(x
1
)l
1
1
(x) = f(x
0
)
x −x
1
x
0
−x
1
+ f(x
1
)
x −x
0
x
1
−x
0
.
Example 3.6
This time we are looking for the Lagrange form of the interpolation polynomial, Q
2
(x),
that interpolates three points: (x
0
, f(x
0
)), (x
1
, f(x
1
)), (x
2
, f(x
2
)). Unfortunately, the
Lagrange form of the interpolation polynomial does not let us use the interpolation
polynomial through the ﬁrst two points, Q
1
(x), as a building block for Q
2
(x). This
means that we have to compute all the polynomials l
n
j
(x) from scratch. We start with
l
2
0
(x) =
(x −x
1
)(x −x
2
)
(x
0
−x
1
)(x
0
−x
2
)
,
l
2
1
(x) =
(x −x
0
)(x −x
2
)
(x
1
−x
0
)(x
1
−x
2
)
,
l
2
2
(x) =
(x −x
0
)(x −x
1
)
(x
2
−x
0
)(x
2
−x
1
)
.
8
D. Levy 3.5 The Lagrange Form of the Interpolation Polynomial
The interpolation polynomial is therefore given by
Q
2
(x) = f(x
0
)l
2
0
(x) + f(x
1
)l
2
1
(x) + f(x
2
)l
2
2
(x)
= f(x
0
)
(x −x
1
)(x −x
2
)
(x
0
−x
1
)(x
0
−x
2
)
+ f(x
1
)
(x −x
0
)(x −x
2
)
(x
1
−x
0
)(x
1
−x
2
)
+ f(x
2
)
(x −x
0
)(x −x
1
)
(x
2
−x
0
)(x
2
−x
1
)
.
It is easy to verify that indeed Q
2
(x
j
) = f(x
j
) for j = 0, 1, 2, as desired.
Remarks.
1. One instance where the Lagrange form of the interpolation polynomial may seem
to be advantageous when compared with the Newton form is when there is a
need to solve several interpolation problems, all given at the same interpolation
points x
0
, . . . x
n
but with diﬀerent values f(x
0
), . . . , f(x
n
). In this case, the
polynomials l
n
j
(x) are identical for all problems since they depend only on the
points but not on the values of the function at these points. Therefore, they have
to be constructed only once.
2. An alternative form for l
n
j
(x) can be obtained in the following way. Deﬁne the
polynomials w
n
(x) of degree n + 1 by
w
n
(x) =
n

i=0
(x −x
i
).
Then it its derivative is
w

n
(x) =
n

j=0
n

i=0
i=j
(x −x
i
). (3.20)
When w

x
(x) is evaluated at an interpolation point, x
j
, there is only one term in
the sum in (3.20) that does not vanish:
w

n
(x
j
) =
n

i=0
i=j
(x
j
−x
i
).
Hence, in view of (3.19), l
n
j
(x) can be rewritten as
l
n
j
(x) =
w
n
(x)
(x −x
j
)w

n
(x
j
)
, 0 j n. (3.21)
3. For future reference we note that the coeﬃcient of x
n
in the interpolation
polynomial Q
n
(x) is
n

j=0
f(x
j
)
n

k=0
k=j
(x
j
−x
k
)
. (3.22)
9
3.6 Divided Diﬀerences D. Levy
For example, the coeﬃcient of x in Q
1
(x) in Example 3.5 is
f(x
0
)
x
0
−x
1
+
f(x
1
)
x
1
−x
0
.
3.6 Divided Diﬀerences
We recall that Newton’s form of the interpolation polynomial is given by (see (3.6)–(3.7))
Q
n
(x) = a
0
+ a
1
(x −x
0
) + . . . + a
n
(x −x
0
) · . . . · (x −x
n−1
),
with a
0
= f(x
0
) and
a
j
=
f(x
j
) −Q
j−1
(x
j
)
j−1

k=0
(x
j
−x
k
)
, 1 j n.
From now on, we will refer to the coeﬃcient, a
j
, as the j
th
-order divided diﬀerence.
The j
th
-order divided diﬀerence, a
j
, is based on the points x
0
, . . . , x
j
and on the values
of the function at these points f(x
0
), . . . , f(x
j
). To emphasize this dependence, we use
the following notation:
a
j
= f[x
0
, . . . , x
j
], 1 j n. (3.23)
We also denote the zeroth-order divided diﬀerence as
a
0
= f[x
0
],
where
f[x
0
] = f(x
0
).
Using the divided diﬀerences notation (3.23), the Newton form of the interpolation
polynomial becomes
Q
n
(x) = f[x
0
] + f[x
0
, x
1
](x −x
0
) + . . . + f[x
0
, . . . x
n
]
n−1

k=0
(x −x
k
). (3.24)
There is a simple recursive way of computing the j
th
-order divided diﬀerence from
divided diﬀerences of lower order, as shown by the following lemma:
Lemma 3.7 The divided diﬀerences satisfy:
f[x
0
, . . . x
n
] =
f[x
1
, . . . x
n
] −f[x
0
, . . . x
n−1
]
x
n
−x
0
. (3.25)
10
D. Levy 3.6 Divided Diﬀerences
Proof. For any k, we denote by Q
k
(x), a polynomial of degree k, that interpolates
f(x) at x
0
, . . . , x
k
, i.e.,
Q
k
(x
j
) = f(x
j
), 0 j k.
We now consider the unique polynomial P(x) of degree n −1 that interpolates f(x)
at x
1
, . . . , x
n
, and claim that
Q
n
(x) = P(x) +
x −x
n
x
n
−x
0
[P(x) −Q
n−1
(x)]. (3.26)
In order to verify this equality, we note that for i = 1, . . . , n −1, P(x
i
) = Q
n−1
(x
i
) so
that
RHS(x
i
) = P(x
i
) = f(x
i
).
At x
n
, RHS(x
n
) = P(x
n
) = f(x
n
). Finally, at x
0
,
RHS(x
0
) = P(x
0
) +
x
0
−x
n
x
n
−x
0
[P(x
0
) −Q
n−1
(x
0
)] = Q
n−1
(x
0
) = f(x
0
).
Hence, the RHS of (3.26) interpolates f(x) at the n + 1 points x
0
, . . . , x
n
, which is also
true for Q
n
(x) due to its deﬁnition. Since the RHS and the LHS in equation (3.26) are
both polynomials of degree n, the uniqueness of the interpolating polynomial (in
this case through n + 1 points) implies the equality in (3.26).
Once we established the equality in (3.26) we can compare the coeﬃcients of the
monomials on both sides of the equation. The coeﬃcient of x
n
on the left-hand-side of
(3.26) is f[x
0
, . . . , x
n
]. The coeﬃcient of x
n−1
in P(x) is f[x
1
, . . . , x
n
] and the
coeﬃcient of x
n−1
in Q
n−1
(x) is f[x
0
, . . . , x
n−1
]. Hence, the coeﬃcient of x
n
on the
right-hand-side of (3.26) is
1
x
n
−x
0
(f[x
1
, . . . , x
n
] −f[x
0
, . . . , x
n−1
]),
which means that
f[x
0
, . . . x
n
] =
f[x
1
, . . . x
n
] −f[x
0
, . . . x
n−1
]
x
n
−x
0
.

Remark. In some books, instead of deﬁning the divided diﬀerence in such a way that
they satisfy (3.25), the divided diﬀerences are deﬁned by the formula
f[x
0
, . . . x
n
] = f[x
1
, . . . x
n
] −f[x
0
, . . . x
n−1
].
If this is the case, all our results on divided diﬀerences should be adjusted accordingly
as to account for the missing factor in the denominator.
11
3.6 Divided Diﬀerences D. Levy
Example 3.8
The second-order divided diﬀerence is
f[x
0
, x
1
, x
2
] =
f[x
1
, x
2
] −f[x
0
, x
1
]
x
2
−x
0
=
f(x
2
)−f(x
1
)
x
2
−x
1

f(x
1
)−f(x
0
)
x
1
−x
0
x
2
−x
0
.
Hence, the unique polynomial that interpolates (x
0
, f(x
0
)), (x
1
, f(x
1
)), and (x
2
, f(x
2
))
is
Q
2
(x) = f[x
0
] + f[x
0
, x
1
](x −x
0
) + f[x
0
, x
1
, x
2
](x −x
0
)(x −x
1
)
= f(x
0
) +
f(x
1
) −f(x
0
)
x
1
−x
0
(x −x
0
) +
f(x
2
)−f(x
1
)
x
2
−x
1

f(x
1
)−f(x
0
)
x
1
−x
0
x
2
−x
0
(x −x
0
)(x −x
1
).
For example, if we want to ﬁnd the polynomial of degree 2 that interpolates (−1, 9),
(0, 5), and (1, 3), we have
f(−1) = 9,
f[−1, 0] =
5 −9
0 −(−1)
= −4, f[0, 1] =
3 −5
1 −0
= −2,
f[−1, 0, 1] =
f[0, 1] −f[−1, 0]
1 −(−1)
=
−2 + 4
2
= 1.
so that
Q
2
(x) = 9 −4(x + 1) + (x + 1)x = 5 −3x + x
2
.
The relations between the divided diﬀerences are schematically portrayed in Table 3.1
(up to third-order). We note that the divided diﬀerences that are being used as the
coeﬃcients in the interpolation polynomial are those that are located in the top of every
column. The recursive structure of the divided diﬀerences implies that it is required to
compute all the low order coeﬃcients in the table in order to get the high-order ones.
One important property of any divided diﬀerence is that it is a symmetric function
of its arguments. This means that if we assume that y
0
, . . . , y
n
is any permutation of
x
0
, . . . , x
n
, then
f[y
0
, . . . , y
n
] = f[x
0
, . . . , x
n
].
This property can be clearly explained by recalling that f[x
0
, . . . , x
n
] plays the role of
the coeﬃcient of x
n
in the polynomial that interpolates f(x) at x
0
, . . . , x
n
. At the same
time, f[y
0
, . . . , y
n
] is the coeﬃcient of x
n
at the polynomial that interpolates f(x) at the
same points. Since the interpolation polynomial is unique for any given data set, the
order of the points does not matter, and hence these two coeﬃcients must be identical.
12
D. Levy 3.7 The Error in Polynomial Interpolation
x
0
f(x
0
)

f[x
0
, x
1
]

x
1
f(x
1
) f[x
0
, x
1
, x
2
]

f[x
1
, x
2
] f[x
0
, x
1
, x
2
, x
3
]

x
2
f(x
2
) f[x
1
, x
2
, x
3
]

f[x
2
, x
3
]

x
3
f(x
3
)
Table 3.1: Divided Diﬀerences
3.7 The Error in Polynomial Interpolation
Our goal in this section is to provide estimates on the “error” we make when interpolating
data that is taken from sampling an underlying function f(x). While the interpolant and
the function agree with each other at the interpolation points, there is, in general, no
reason to expect them to be close to each other elsewhere. Nevertheless, we can estimate
the diﬀerence between them, a diﬀerence which we refer to as the interpolation error.
We let Π
n
denote the space of polynomials of degree n, and let C
n+1
[a, b] denote the
space of functions that have n + 1 continuous derivatives on the interval [a, b].
Theorem 3.9 Let f(x) ∈ C
n+1
[a, b]. Let Q
n
(x) ∈ Π
n
such that it interpolates f(x) at
the n + 1 distinct points x
0
, . . . , x
n
∈ [a, b]. Then ∀x ∈ [a, b], ∃ξ
n
∈ (a, b) such that
f(x) −Q
n
(x) =
1
(n + 1)!
f
(n+1)

n
)
n

j=0
(x −x
j
). (3.27)
Proof. Fix a point x ∈ [a, b]. If x is one of the interpolation points x
0
, . . . , x
n
, then the
left-hand-side and the right-hand-side of (3.27) are both zero, and the result holds
trivially. We therefore assume that x = x
j
0 j n, and let
w(x) =
n

j=0
(x −x
j
).
We now let
F(y) = f(y) −Q
n
(y) −λw(y),
13
3.7 The Error in Polynomial Interpolation D. Levy
where λ is chosen as to guarantee that F(x) = 0, i.e.,
λ =
f(x) −Q
n
(x)
w(x)
.
Since the interpolation points x
0
, . . . , x
n
and x are distinct, w(x) does not vanish and
λ is well deﬁned. We now note that since f ∈ C
n+1
[a, b] and since Q
n
and w are
polynomials, then also F ∈ C
n+1
[a, b]. In addition, F vanishes at n + 2 points:
x
0
, . . . , x
n
and x. According to Rolle’s theorem, F

has at least n + 1 distinct zeros in
(a, b), F

has at least n distinct zeros in (a, b), and similarly, F
(n+1)
has at least one
zero in (a, b), which we denote by ξ
n
. We have
0 = F
(n+1)

n
) = f
(n+1)

n
) −Q
(n+1)
n

n
) −λ(x)w
(n+1)

n
) (3.28)
= f
(n+1)

n
) −
f(x) −Q
n
(x)
w(x)
(n + 1)!
Here, we used the fact that the leading term of w(x) is x
n+1
, which guarantees that its
(n + 1)
th
derivative is
w
(n+1)
(x) = (n + 1)! (3.29)
Reordering the terms in (3.28) we conclude with
f(x) −Q
n
(x) =
1
(n + 1)!
f
(n+1)

n
)w(x).

In addition to the interpretation of the divided diﬀerence of order n as the coeﬃcient
of x
n
in some interpolation polynomial, it can also be characterized in another important
way. Consider, e.g., the ﬁrst-order divided diﬀerence
f[x
0
, x
1
] =
f(x
1
) −f(x
0
)
x
1
−x
0
.
Since the order of the points does not change the value of the divided diﬀerence, we can
assume, without any loss of generality, that x
0
< x
1
. If we assume, in addition, that
f(x) is continuously diﬀerentiable in the interval [x
0
, x
1
], then this divided diﬀerence
equals to the derivative of f(x) at an intermediate point, i.e.,
f[x
0
, x
1
] = f

(ξ), ξ ∈ (x
0
, x
1
).
In other words, the ﬁrst-order divided diﬀerence can be viewed as an approximation
of the ﬁrst derivative of f(x) in the interval. It is important to note that while this
interpretation is based on additional smoothness requirements from f(x) (i.e. its be-
ing diﬀerentiable), the divided diﬀerences are well deﬁned also for non-diﬀerentiable
functions.
This notion can be extended to divided diﬀerences of higher order as stated by the
following lemma.
14
D. Levy 3.7 The Error in Polynomial Interpolation
Lemma 3.10 Let x, x
0
, . . . , x
n−1
be n + 1 distinct points. Let a = min(x, x
0
, . . . , x
n−1
)
and b = max(x, x
0
, . . . , x
n−1
). Assume that f(y) has a continuous derivative of order n
in the interval (a, b). Then
f[x
0
, . . . , x
n−1
, x] =
f
(n)
(ξ)
n!
, (3.30)
where ξ ∈ (a, b).
Proof. Let Q
n
(y) interpolate f(y) at x
0
, . . . , x
n−1
, x. Then according to the
construction of the Newton form of the interpolation polynomial (3.24), we know that
Q
n
(y) = Q
n−1
(y) + f[x
0
, . . . , x
n−1
, x]
n−1

j=0
(y −x
j
).
Since Q
n
(y) interpolated f(y) at x, we have
f(x) = Q
n−1
(x) + f[x
0
, . . . , x
n−1
, x]
n−1

j=0
(x −x
j
).
By Theorem 3.9 we know that the interpolation error is given by
f(x) −Q
n−1
(x) =
1
n!
f
(n)

n−1
)
n−1

j=0
(x −x
j
),
which implies the result (3.30).

Remark. In equation (3.30), we could as well think of the interpolation point x as
any other interpolation point, and name it, e.g., x
n
. In this case, the equation (3.30)
takes the somewhat more natural form of
f[x
0
, . . . , x
n
] =
f
(n)
(ξ)
n!
.
In other words, the n
th
-order divided diﬀerence is an n
th
-derivative of the function f(x)
at an intermediate point, assuming that the function has n continuous derivatives. Sim-
ilarly to the ﬁrst-order divided diﬀerence, we would like to emphasize that the n
th
-order
divided diﬀerence is also well deﬁned in cases where the function is not as smooth as
required in the theorem, though if this is the case, we can no longer consider this divided
diﬀerence to represent a n
th
-order derivative of the function.
15
3.8 Interpolation at the Chebyshev Points D. Levy
3.8 Interpolation at the Chebyshev Points
In the entire discussion up to now, we assumed that the interpolation points are given.
There may be cases where one may have the ﬂexibility of choosing the interpolation
points. If this is the case, it would be reasonable to use this degree of freedom to
minimize the interpolation error.
We recall that if we are interpolating values of a function f(x) that has n continuous
derivatives, the interpolation error is of the form
f(x) −Q
n
(x) =
1
(n + 1)!
f
(n+1)

n
)
n

j=0
(x −x
j
). (3.31)
Here, Q
n
(x) is the interpolating polynomial and ξ
n
is an intermediate point in the
interval of interest (see (3.27)).
It is important to note that the interpolation points inﬂuence two terms on the
right-hand-side of (3.31). The obvious one is the product
n

j=0
(x −x
j
). (3.32)
The second term that depends on the interpolation points is f
(n+1)

n
) since the value
of the intermediate point ξ
n
depends on {x
j
}. Due to the implicit dependence of ξ
n
on the interpolation points, minimizing the interpolation error is not an easy task. We
will return to this “full” problem later on in the context of the minimax approximation
problem. For the time being, we are going to focus on a simpler problem, namely, how
to choose the interpolation points x
0
, . . . , x
n
such that the product (3.32) is minimized.
The solution of this problem is the topic of this section. Once again, we would like to
emphasize that a solution of this problem does not (in general) provide an optimal choice
of interpolation points that minimizes the interpolation error. All that it guarantees is
that the product part of the interpolation error is minimal.
The tool that we are going to use is the Chebyshev polynomials. The solution of
the problem will be to choose the interpolation points as the Chebyshev points. We will
ﬁrst introduce the Chebyshev polynomials and the Chebyshev points and then explain
why interpolating at these points minimizes (3.32).
The Chebyshev polynomials can be deﬁned using the following recursion relation:
_
_
_
T
0
(x) = 1,
T
1
(x) = x,
T
n+1
(x) = 2xT
n
(x) −T
n−1
(x), n 1.
(3.33)
For example, T
2
(x) = 2xT
1
(x)−T
0
(x) = 2x
2
−1, and T
3
(x) = 4x
3
−3x. The polynomials
T
1
(x), T
2
(x) and T
3
(x) are plotted in Figure 3.2.
Instead of writing the recursion formula, (3.33), it is possible to write an explicit
formula for the Chebyshev polynomials:
16
D. Levy 3.8 Interpolation at the Chebyshev Points
! "#\$ "#% "#& "#' " "#' "#& "#% "#\$ !
!
"#\$
"#%
"#&
"#'
"
"#'
"#&
"#%
"#\$
!
x
T
1
(x)
T
2
(x)
T
3
(x)
Figure 3.2: The Chebyshev polynomials T
1
(x), T
2
(x) and T
3
(x)
Lemma 3.11 For x ∈ [−1, 1],
T
n
(x) = cos(ncos
−1
x), n 0. (3.34)
Proof. Standard trigonometric identities imply that
cos(n + 1)θ = cos θ cos nθ −sin θ sin nθ,
cos(n −1)θ = cos θ cos nθ + sin θ sin nθ.
Hence
cos(n + 1)θ = 2 cos θ cos nθ −cos(n −1)θ. (3.35)
We now let θ = cos
−1
x, i.e., x = cos θ, and deﬁne
t
n
(x) = cos(ncos
−1
x) = cos(nθ).
Then by (3.35)
_
_
_
t
0
(x) = 1,
t
1
(x) = x,
t
n+1
(x) = 2xt
n
(x) −t
n−1
(x), n 1.
Hence t
n
(x) = T
n
(x).

17
3.8 Interpolation at the Chebyshev Points D. Levy
What is so special about the Chebyshev polynomials, and what is the connection
between these polynomials and minimizing the interpolation error? We are about to
answer these questions, but before doing so, there is one more issue that we must
clarify.
We deﬁne a monic polynomial as a polynomial for which the coeﬃcient of the
leading term is one, i.e., a polynomial of degree n is monic, if it is of the form
x
n
+ a
n−1
x
n−1
+ . . . + a
1
x + a
0
.
Note that Chebyshev polynomials are not monic: the deﬁnition (3.33) implies that the
Chebyshev polynomial of degree n is of the form
T
n
(x) = 2
n−1
x
n
+ . . .
This means that T
n
(x) divided by 2
n−1
is monic, i.e.,
2
1−n
T
n
(x) = x
n
+ . . .
A general result about monic polynomials is given by the following theorem
Theorem 3.12 If p
n
(x) is a monic polynomial of degree n, then
max
−1x1
|p
n
(x)| 2
1−n
. (3.36)
Proof. We prove (3.36) by contradiction. Suppose that
|p
n
(x)| < 2
1−n
, |x| 1.
Let
q
n
(x) = 2
1−n
T
n
(x),
and let x
j
be the following n + 1 points
x
j
= cos
_

n
_
, 0 j n.
Since
T
n
_
cos

n
_
= (−1)
j
,
we have
(−1)
j
q
n
(x
j
) = 2
1−n
.
18
D. Levy 3.8 Interpolation at the Chebyshev Points
Hence
(−1)
j
p
n
(x
j
) |p
n
(x
j
)| < 2
1−n
= (−1)
j
q
n
(x
j
).
This means that
(−1)
j
(q
n
(x
j
) −p
n
(x
j
)) > 0, 0 j n.
Hence, the polynomial (q
n
−p
n
)(x) oscillates (n+1) times in the interval [−1, 1], which
means that (q
n
−p
n
)(x) has at least n distinct roots in the interval. However, p
n
(x) and
q
n
(x) are both monic polynomials which means that their diﬀerence is a polynomial of
degree n − 1 at most. Such a polynomial cannot have more than n − 1 distinct roots,
which leads to a contradiction. Note that p
n
−q
n
cannot be the zero polynomial because
that will imply that p
n
(x) and q
n
(x) are identical which again is not possible due to the
assumptions on their maximum values.

We are now ready to use Theorem 3.12 to ﬁgure out how to reduce the interpolation
error. We know by Theorem 3.9 that if the interpolation points x
0
, . . . , x
n
∈ [−1, 1],
then there exists ξ
n
∈ (−1, 1) such that the distance between the function whose values
we interpolate, f(x), and the interpolation polynomial, Q
n
(x), is
max
|x|1
|f(x) −Q
n
(x)|
1
(n + 1)!
max
|x|1
|f
(n+1)

n
)| max
|x|1
¸
¸
¸
¸
¸
n

j=0
(x −x
j
)
¸
¸
¸
¸
¸
.
We are interested in minimizing
max
|x|1
¸
¸
¸
¸
¸
n

j=0
(x −x
j
)
¸
¸
¸
¸
¸
.
We note that

n
j=0
(x −x
j
) is a monic polynomial of degree n + 1 and hence by Theo-
rem 3.12
max
|x|1
¸
¸
¸
¸
¸
n

j=0
(x −x
j
)
¸
¸
¸
¸
¸
2
−n
.
The minimal value of 2
−n
can be actually obtained if we set
2
−n
T
n+1
(x) =
n

j=0
(x −x
j
),
which is equivalent to choosing x
j
as the roots of the Chebyshev polynomial T
n+1
(x).
Here, we have used the obvious fact that |T
n
(x)| 1.
19
3.8 Interpolation at the Chebyshev Points D. Levy
What are the roots of the Chebyshev polynomial T
n+1
(x)? By Lemma 3.11
T
n+1
(x) = cos((n + 1) cos
−1
x).
The roots of T
n+1
(x), x
0
, . . . , x
n
, are therefore obtained if
(n + 1) cos
−1
(x
j
) =
_
j +
1
2
_
π, 0 j n,
i.e., the (n + 1) roots of T
n+1
(x) are
x
j
= cos
_
2j + 1
2n + 2
π
_
, 0 j n. (3.37)
The roots of the Chebyshev polynomials are sometimes referred to as the Chebyshev
points. The formula (3.37) for the roots of the Chebyshev polynomial has the following
geometrical interpretation. In order to ﬁnd the roots of T
n
(x), deﬁne α = π/n. Divide
the upper half of the unit circle into n + 1 parts such that the two side angles are α/2
and the other angles are α. The Chebyshev points are then obtained by projecting these
points on the x-axis. This procedure is demonstrated in Figure 3.3 for T
4
(x).
-1 x
0
x
1
0 x
2
x
3
1
The unit circle

8

8

8
π
8
x
Figure 3.3: The roots of the Chebyshev polynomial T
4
(x), x
0
, . . . , x
3
. Note that they
become dense next to the boundary of the interval
The following theorem summarizes the discussion on interpolation at the Chebyshev
points. It also provides an estimate of the error for this case.
20
D. Levy 3.8 Interpolation at the Chebyshev Points
Theorem 3.13 Assume that Q
n
(x) interpolates f(x) at x
0
, . . . , x
n
. Assume also that
these (n + 1) interpolation points are the (n + 1) roots of the Chebyshev polynomial of
degree n + 1, T
n+1
(x), i.e.,
x
j
= cos
_
2j + 1
2n + 2
π
_
, 0 j n.
Then ∀|x| 1,
|f(x) −Q
n
(x)|
1
2
n
(n + 1)!
max
|ξ|1
¸
¸
f
(n+1)
(ξ)
¸
¸
. (3.38)
Example 3.14
Problem: Let f(x) = sin(πx) in the interval [−1, 1]. Find Q
2
(x) which interpolates f(x)
in the Chebyshev points. Estimate the error.
Solution: Since we are asked to ﬁnd an interpolation polynomial of degree 2, we need
3 interpolation points. We are also asked to interpolate at the Chebyshev points, and
hence we ﬁrst need to compute the 3 roots of the Chebyshev polynomial of degree 3,
T
3
(x) = 4x
3
−3x.
The roots of T
3
(x) can be easily found from x(4x
2
−3) = 0, i.e.,
x
0
= −

3
2
, , x
1
= 0, x
2
=

3
2
.
The corresponding values of f(x) at these interpolation points are
f(x
0
) = sin
_

3
2
π
_
≈ −0.4086,
f(x
1
) = 0,
f(x
2
) = sin
_

3
2
π
_
≈ 0.4086.
The ﬁrst-order divided diﬀerences are
f[x
0
, x
1
] =
f(x
1
) −f(x
0
)
x
1
−x
0
≈ 0.4718,
f[x
1
, x
2
] =
f(x
2
) −f(x
1
)
x
2
−x
1
≈ 0.4718,
and the second-order divided diﬀerence is
f[x
0
, x
1
, x
2
] =
f[x
1
, x
2
] −f[x
0
, x
1
]
x
2
−x
0
= 0.
21
3.8 Interpolation at the Chebyshev Points D. Levy
The interpolation polynomial is
Q
2
(x) = f(x
0
) + f[x
0
, x
1
](x −x
0
) + f[x
0
, x
1
, x
2
](x −x
0
)(x −x
1
) ≈ 0.4718x.
The original function f(x) and the interpolant at the Chebyshev points, Q
2
(x), are
plotted in Figure 3.4.
As of the error estimate, ∀|x| 1,
| sin πx −Q
2
(x)|
1
2
2
3!
max
|ξ|1
|(sin πt)
(3)
|
π
3
2
2
3!
1.292
A brief examination of Figure 3.4 reveals that while this error estimate is correct, it is
far from being sharp.
! "#\$ "#% "#& "#' " "#' "#& "#% "#\$ !
!
"#\$
"#%
"#&
"#'
"
"#'
"#&
"#%
"#\$
!
Q
2
(x)
f(x)
x
Figure 3.4: The function f(x) = sin(π(x)) and the interpolation polynomial Q
2
(x) that
interpolates f(x) at the Chebyshev points. See Example 3.14.
Remark. In the more general case where the interpolation interval for the function
f(x) is x ∈ [a, b], it is still possible to use the previous results by following the
following steps: Start by converting the interpolation interval to y ∈ [−1, 1]:
x =
(b −a)y + (a + b)
2
.
This converts the interpolation problem for f(x) on [a, b] into an interpolation problem
for f(x) = g(x(y)) in y ∈ [−1, 1]. The Chebyshev points in the interval y ∈ [−1, 1] are
22
D. Levy 3.9 Hermite Interpolation
the roots of the Chebyshev polynomial T
n+1
(x), i.e.,
y
j
= cos
_
2j + 1
2n + 2
π
_
, 0 j n.
The corresponding n + 1 interpolation points in the interval [a, b] are
x
j
=
(b −a)y
j
+ (a + b)
2
, 0 j n.
In this case, the product term in the interpolation error is
max
x∈[a,b]
¸
¸
¸
¸
¸
n

j=0
(x −x
j
)
¸
¸
¸
¸
¸
=
¸
¸
¸
¸
b −a
2
¸
¸
¸
¸
n+1
max
|y|1
¸
¸
¸
¸
¸
n

j=0
(x −y
j
)
¸
¸
¸
¸
¸
,
and the interpolation error is given by
|f(x) −Q
n
(x)|
1
2
n
(n + 1)!
¸
¸
¸
¸
b −a
2
¸
¸
¸
¸
n+1
max
ξ∈[a,b]
¸
¸
f
(n+1)
(ξ)
¸
¸
. (3.39)
3.9 Hermite Interpolation
We now turn to a slightly diﬀerent interpolation problem in which we assume that
in addition to interpolating the values of the function at certain points, we are also
interested in interpolating its derivatives. Interpolation that involves the derivatives is
called Hermite interpolation. Such an interpolation problem is demonstrated in the
following example:
Example 3.15
Problem: Find a polynomials p(x) such that p(1) = −1, p

(1) = −1, and p(0) = 1.
Solution: Since three conditions have to be satisﬁed, we can use these conditions to
determine three degrees of freedom, which means that it is reasonable to expect that
these conditions uniquely determine a polynomial of degree 2. We therefore let
p(x) = a
0
+ a
1
x + a
2
x
2
.
The conditions of the problem then imply that
_
_
_
a
0
+ a
1
+ a
2
= −1,
a
1
+ 2a
2
= −1,
a
0
= 1.
Hence, there is indeed a unique polynomial that satisﬁes the interpolation conditions
and it is
p(x) = x
2
−3x + 1.
23
3.9 Hermite Interpolation D. Levy
In general, we may have to interpolate high-order derivatives and not only ﬁrst-
order derivatives. Also, we assume that for any point x
j
in which we have to satisfy an
interpolation condition of the form
p
(l)
(x
j
) = f(x
j
),
(with p
(l)
being the l
th
-order derivative of p(x)), we are also given all the values of the
lower-order derivatives up to l as part of the interpolation requirements, i.e.,
p
(i)
(x
j
) = f
(i)
(x
j
), 0 i l.
If this is not the case, it may not be possible to ﬁnd a unique interpolant as demonstrated
in the following example.
Example 3.16
Problem: Find p(x) such that p

(0) = 1 and p

(1) = −1.
Solution: Since we are asked to interpolate two conditions, we may expect them to
uniquely determine a linear function, say
p(x) = a
0
+ a
1
x.
However, both conditions specify the derivative of p(x) at two distinct points to be
of diﬀerent values, which amounts to a contradicting information on the value of a
1
.
Hence, a linear polynomial cannot interpolate the data and we must consider higher-
order polynomials. Unfortunately, a polynomial of order 2 will no longer be unique
because not enough information is given. Note that even if the prescribed values of the
derivatives were identical, we will not have problems with the coeﬃcient of the linear
term a
1
, but we will still not have enough information to determine the constant a
0
.
A simple case that you are probably already familiar with is the Taylor series.
When viewed from the point of view that we advocate in this section, one can consider
the Taylor series as an interpolation problem in which one has to interpolate the value
of the function and its ﬁrst n derivatives at a given point, say x
0
, i.e., the interpolation
conditions are:
p
(j)
(x
0
) = f
(j)
(x
0
), 0 j n.
The unique solution of this problem in terms of a polynomial of degree n is
p(x) = f(x
0
) + f

(x
0
)(x −x
0
) + . . . +
f
(n)
(x
0
)
n!
(x −x
0
)
n
=
n

j=0
f
(j)
(x
0
)
j!
(x −x
0
)
j
,
which is the Taylor series of f(x) expanded about x = x
0
.
24
D. Levy 3.9 Hermite Interpolation
3.9.1 Divided diﬀerences with repetitions
We are now ready to consider the Hermite interpolation problem. The ﬁrst form we
study is the Newton form of the Hermite interpolation polynomial. We start by extend-
ing the deﬁnition of divided diﬀerences in such a way that they can handle derivatives.
We already know that the ﬁrst derivative is connected with the ﬁrst-order divided dif-
ference by
f

(x
0
) = lim
x→x
0
f(x) −f(x
0
)
x −x
0
= lim
x→x
0
f[x, x
0
].
Hence, it is natural to extend the notion of divided diﬀerences by the following deﬁnition.
Deﬁnition 3.17 The ﬁrst-order divided diﬀerence with repetitions is deﬁned as
f[x
0
, x
0
] = f

(x
0
). (3.40)
In a similar way, we can extend the notion of divided diﬀerences to high-order derivatives
as stated in the following lemma (which we leave without a proof).
Lemma 3.18 Let x
0
x
1
. . . x
n
. Then the divided diﬀerences satisfy
f[x
0
, . . . x
n
] =
_
¸
_
¸
_
f[x
1
, . . . , x
n
] −f[x
0
, . . . , x
n−1
]
x
n
−x
0
, x
n
= x
0
,
f
(n)
(x
0
)
n!
, x
n
= x
0
.
(3.41)
We now consider the following Hermite interpolation problem: The interpolation
points are x
0
, . . . , x
l
(which we assume are ordered from small to large). At each inter-
polation point x
j
, we have to satisfy the interpolation conditions:
p
(i)
(x
j
) = f
(i)
(x
j
), 0 i m
j
.
Here, m
j
denotes the number of derivatives that we have to interpolate for each point
x
j
(with the standard notation that zero derivatives refers to the value of the function
only). In general, the number of derivatives that we have to interpolate may change
from point to point. The extended notion of divided diﬀerences allows us to write the
solution to this problem in the following way:
We let n denote the total number of points including their multiplicities (that cor-
respond to the number of derivatives we have to interpolate at each point), i.e.,
n = m
1
+ m
2
+ . . . + m
l
.
We then list all the points including their multiplicities (that correspond to the number
of derivatives we have to interpolate). To simplify the notations we identify these points
with a new ordered list of points y
i
:
{y
0
, . . . , y
n−1
} = {x
0
, . . . , x
0
. ¸¸ .
m
1
, x
1
, . . . , x
1
. ¸¸ .
m
2
, . . . , x
l
, . . . , x
l
. ¸¸ .
m
l
}.
25
3.9 Hermite Interpolation D. Levy
The interpolation polynomial p
n−1
(x) is given by
p
n−1
(x) = f[y
0
] +
n−1

j=1
f[y
0
, . . . , y
j
]
j−1

k=0
(x −y
k
). (3.42)
Whenever a point repeats in f[y
0
, . . . , y
j
], we interpret this divided diﬀerence in terms
of the extended deﬁnition (3.41). In practice, there is no need to shift the notations to
y’s and we work directly with the original points. We demonstrate this interpolation
procedure in the following example.
Example 3.19
Problem: Find an interpolation polynomial p(x) that satisﬁes
_
_
_
p(x
0
) = f(x
0
),
p(x
1
) = f(x
1
),
p

(x
1
) = f

(x
1
).
Solution: The interpolation polynomial p(x) is
p(x) = f(x
0
) + f[x
0
, x
1
](x −x
0
) + f[x
0
, x
1
, x
1
](x −x
0
)(x −x
1
).
The divided diﬀerences:
f[x
0
, x
1
] =
f(x
1
) −f(x
0
)
x
1
−x
0
.
f[x
0
, x
1
, x
1
] =
f[x
1
, x
1
] −f[x
1
, x
0
]
x
1
−x
0
=
f

(x
1
) −
f(x
1
)−f(x
0
)
x
1
−x
0
x
1
−x
0
.
Hence
p(x) = f(x
0
)+
f(x
1
) −f(x
0
)
x
1
−x
0
(x−x
0
)+
(x
1
−x
0
)f

(x
1
) −[f(x
1
) −f(x
0
)]
(x
1
−x
0
)
2
(x−x
0
)(x−x
1
).
3.9.2 The Lagrange form of the Hermite interpolant
In this section we are interested in writing the Lagrange form of the Hermite interpolant
in the special case in which the nodes are x
0
, . . . , x
n
and the interpolation conditions
are
p(x
i
) = f(x
i
), p

(x
i
) = f

(x
i
), 0 i n. (3.43)
We look for an interpolant of the form
p(x) =
n

i=0
f(x
i
)A
i
(x) +
n

i=0
f

(x
i
)B
i
(x). (3.44)
26
D. Levy 3.9 Hermite Interpolation
In order to satisfy the interpolation conditions (3.43), the polynomials p(x) in (3.44)
must satisfy the 2n + 2 conditions:
_
_
_
A
i
(x
j
) = δ
ij
, B
i
(x
j
) = 0,
i, j = 0, . . . , n.
A

i
(x
j
) = 0, B

i
(x
j
) = δ
ij
,
(3.45)
We thus expect to have a unique polynomial p(x) that satisﬁes the constraints (3.45)
assuming that we limit its degree to be 2n + 1.
It is convenient to start the construction with the functions we have used in the
Lagrange form of the standard interpolation problem (Section 3.5). We already know
that
l
i
(x) =
n

j=0
j=i
x −x
j
x
i
−x
j
,
satisfy l
i
(x
j
) = δ
ij
. In addition, for i = j,
l
2
i
(x
j
) = 0, (l
2
i
(x
j
))

= 0.
The degree of l
i
(x) is n, which means that the degree of l
2
i
(x) is 2n. We will thus
assume that the unknown polynomials A
i
(x) and B
i
(x) in (3.45) can be written as
_
A
i
(x) = r
i
(x)l
2
i
(x),
B
i
(x) = s
i
(x)l
2
i
(x).
The functions r
i
(x) and s
i
(x) are both assumed to be linear, which implies that deg(A
i
) =
deg(B
i
) = 2n + 1, as desired. Now, according to (3.45)
δ
ij
= A
i
(x
j
) = r
i
(x
j
)l
2
i
(x
j
) = r
i
(x
j

ij
.
Hence
r
i
(x
i
) = 1. (3.46)
Also,
0 = A

i
(x
j
) = r

i
(x
j
)[l
i
(x
j
)]
2
+ 2r
i
(x
j
)l
i
(x
J
)l

i
(x
j
) = r

i
(x
j

ij
+ 2r
i
(x
j

ij
l

i
(x
j
),
and thus
r

i
(x
i
) + 2l

i
(x
i
) = 0. (3.47)
Assuming that r
i
(x) is linear, r
i
(x) = ax + b, equations (3.46),(3.47), imply that
a = −2l

i
(x
i
), b = 1 + 2l

i
(x
i
)x
i
.
27
3.9 Hermite Interpolation D. Levy
Therefore
A
i
(x) = [1 + 2l

i
(x
i
)(x
i
−x)]l
2
i
(x).
As of B
i
(x) in (3.44), the conditions (3.45) imply that
0 = B
i
(x
j
) = s
i
(x
j
)l
2
i
(x
j
) =⇒ s
i
(x
i
) = 0, (3.48)
and
δ
ij
= B

i
(x
j
) = s

i
(x
j
)l
2
i
(x
j
) + 2s
i
(x
j
)(l
2
i
(x
j
))

=⇒ s

i
(x
i
) = 1. (3.49)
Combining (3.48) and (3.49), we obtain
s
i
(x) = x −x
i
,
so that
B
i
(x) = (x −x
i
)l
2
i
(x).
To summarize, the Lagrange form of the Hermite interpolation polynomial is given by
p(x) =
n

i=0
f(x
i
)[1 + 2l

i
(x
i
)(x
i
−x)]l
2
i
(x) +
n

i=0
f

(x
i
)(x −x
i
)l
2
i
(x). (3.50)
The error in the Hermite interpolation (3.50) is given by the following theorem.
Theorem 3.20 Let x
0
, . . . , x
n
be distinct nodes in [a, b] and f ∈ C
2n+2
[a, b]. If p ∈
Π
2n+1
, such that ∀0 i n,
p(x
i
) = f(x
i
), p

(x
i
) = f

(x
i
),
then ∀x ∈ [a, b], there exists ξ ∈ (a, b) such that
f(x) −p(x) =
f
(2n+2)
(ξ)
(2n + 2)!
n

i=0
(x −x
i
)
2
. (3.51)
Proof. The proof follows the same techniques we used in proving Theorem 3.9. If x is
one of the interpolation points, the result trivially holds. We thus ﬁx x as a
non-interpolation point and deﬁne
w(y) =
n

i=0
(y −x
i
)
2
.
We also have
φ(y) = f(y) −p(y) −λw(y),
28
D. Levy 3.10 Spline Interpolation
and select λ such that φ(x) = 0, i.e.,
λ =
f(x) −p(x)
w(x)
.
φ has (at least) n + 2 zeros in [a, b]: (x, x
0
, . . . , x
n
). By Rolle’s theorem, we know that
φ

has (at least) n + 1 zeros that are diﬀerent than (x, x
0
, . . . , x
n
). Also, φ

vanishes at
x
0
, . . . , x
n
, which means that φ

has at least 2n + 2 zeros in [a, b].
Similarly, Rolle’s theorem implies that φ

has at least 2n + 1 zeros in (a, b), and by
induction, φ
(2n+2)
has at least one zero in (a, b), say ξ.
Hence
0 = φ
(2n+2)
(ξ) = f
(2n+2)
(ξ) −p
(2n+2)
(ξ) −λw
(2n+2)
(ξ).
Since the leading term in w(y) is x
2n+2
, w
(2n+2)
(ξ) = (2n + 2)!. Also, since
p(x) ∈ Π
2n+1
, p
(2n+2)
(ξ) = 0. We recall that x was an arbitrary (non-interpolation)
point and hence we have
f(x) −p(x) =
f
(2n+2)
(ξ)
(2n + 2)!
n

i=0
(x −x
i
)
2
.

Example 3.21
Assume that we would like to ﬁnd the Hermite interpolation polynomial that satisﬁes:
p(x
0
) = y
0
, p

(x
0
) = d
0
, p(x
1
) = y
1
, p

(x
1
) = d
1
.
In this case n = 1, and
l
0
(x) =
x −x
1
x
0
−x
1
, l

0
(x) =
1
x
0
−x
1
, l
1
(x) =
x −x
0
x
1
−x
0
, l

1
(x) =
1
x
1
−x
0
.
According to (3.50), the desired polynomial is given by (check!)
p(x) = y
0
_
1 +
2
x
0
−x
1
(x
0
−x)
_ _
x −x
1
x
0
−x
1
_
2
+ y
1
_
1 +
2
x
1
−x
0
(x
1
−x)
_ _
x −x
0
x
1
−x
0
_
2
+d
0
(x −x
0
)
_
x −x
1
x
0
−x
1
_
2
+ d
1
(x −x
1
)
_
x −x
0
x
1
−x
0
_
2
.
3.10 Spline Interpolation
So far, the only type of interpolation we were dealing with was polynomial interpolation.
In this section we discuss a diﬀerent type of interpolation: piecewise-polynomial interpo-
lation. A simple example of such interpolants will be the function we get by connecting
data with straight lines (see Figure 3.5). Of course, we would like to generate functions
29
3.10 Spline Interpolation D. Levy
(x
0
, f(x
0
))
(x
1
, f(x
1
))
(x
2
, f(x
2
))
(x
3
, f(x
3
))
(x
4
, f(x
4
))
x
Figure 3.5: A piecewise-linear spline. In every subinterval the function is linear. Overall
it is continuous where the regularity is lost at the knots
that are somewhat smoother than piecewise-linear functions, and still interpolate the
data. The functions we will discuss in this section are splines.
You may still wonder why are we interested in such functions at all? It is easy to
motivate this discussion by looking at Figure 3.6. In this ﬁgure we demonstrate what a
high-order interpolant looks like. Even though the data that we interpolate has only one
extrema in the domain, we have no control over the oscillatory nature of the high-order
interpolating polynomial. In general, high-order polynomials are oscillatory, which rules
them as non-practical for many applications. That is why we focus our attention in this
section on splines.
Splines, should be thought of as polynomials on subintervals that are connected in
a “smooth way”. We will be more rigorous when we deﬁne precisely what we mean by
smooth. First, we pick n+1 points which we refer to as the knots: t
0
< t
1
< · · · < t
n
. A
spline of degree k having knots t
0
, . . . , t
n
is a function s(x) that satisﬁes the following
two properties:
1. On [t
i−1
, t
i
) s(x) is a polynomial of degree k, i.e., s(x) is a polynomial on every
subinterval that is deﬁned by the knots.
2. Smoothness: s(x) has a continuous (k −1)
th
derivative on the interval [t
0
, t
n
].
A spline of degree 0 is a piecewise-constant function (see Figure 3.7). A spline of
30
D. Levy 3.10 Spline Interpolation
! " # \$ % & % \$ # " !
&'!
&
&'!
%
%'!
\$
1
1+x
2
Q
10
(x)
x
Figure 3.6: An interpolant “goes bad”. In this example we interpolate 11 equally spaced
samples of f(x) =
1
1+x
2
with a polynomial of degree 10, Q
10
(x)
(t
0
, f(t
0
))
(t
1
, f(t
1
))
(t
2
, f(t
2
))
(t
3
, f(t
3
))
(t
4
, f(t
4
))
x
Figure 3.7: A zeroth-order (piecewise-constant) spline. The knots are at the interpola-
tion points. Since the spline is of degree zero, the function is not even continuous
31
3.10 Spline Interpolation D. Levy
degree 1 is a piecewise-linear function that can be explicitly written as
s(x) =
_
¸
¸
¸
_
¸
¸
¸
_
s
0
(x) = a
0
x + b
0
, x ∈ [t
0
, t
1
),
s
1
(x) = a
1
x + b
1
, x ∈ [t
1
, t
2
),
.
.
.
.
.
.
s
n−1
(x) = a
n−1
x + b
n−1
, x ∈ [t
n−1
, t
n
],
(see Figure 3.5 where the knots {t
i
} and the interpolation points {x
i
} are assumed to
be identical). It is now obvious why the points t
0
, . . . , t
n
are called knots: these are
the points that connect the diﬀerent polynomials with each other. To qualify as an
interpolating function, s(x) will have to satisfy interpolation conditions that we will
discuss below. We would like to comment already at this point that knots should not be
confused with the interpolation points. Sometimes it is convenient to choose the knots
to coincide with the interpolation points but this is only optional, and other choices can
be made.
3.10.1 Cubic splines
A special case (which is the most common spline function that is used in practice) is
the cubic spline. A cubic spline is a spline for which the function is a polynomial of
degree 3 on every subinterval, and a function with two continuous derivatives overall
(see Figure 3.8).
Let’s denote such a function by s(x), i.e.,
s(x) =
_
¸
¸
¸
_
¸
¸
¸
_
s
0
(x), x ∈ [t
0
, t
1
),
s
1
(x), x ∈ [t
1
, t
2
),
.
.
.
.
.
.
s
n−1
(x), x ∈ [t
n−1
, t
n
],
where ∀i, the degree of s
i
(x) is 3.
We now assume that some data (that s(x) should interpolate) is given at the knots,
i.e.,
s(t
i
) = y
i
, 0 i n. (3.52)
The interpolation conditions (3.52) in addition to requiring that s(x) is continuous,
imply that
s
i−1
(t
i
) = y
i
= s
i
(t
i
), 1 i n −1. (3.53)
We also require the continuity of the ﬁrst and the second derivatives, i.e.,
s

i
(t
i+1
) = s

i+1
(t
i+1
), 0 i n −2, (3.54)
s

i
(t
i+1
) = s

i+1
(t
i+1
), 0 i n −2.
32
D. Levy 3.10 Spline Interpolation
(t
0
, f(t
0
))
(t
1
, f(t
1
))
(t
2
, f(t
2
))
(t
3
, f(t
3
))
(t
4
, f(t
4
))
x
Figure 3.8: A cubic spline. In every subinterval [t
i−1
, t
i
, the function is a polynomial of
degree 2. The polynomials on the diﬀerent subintervals are connected to each other
in such a way that the spline has a second-order continuous derivative. In this example
we use the not-a-knot condition.
33
3.10 Spline Interpolation D. Levy
Before actually computing the spline, let’s check if we have enough equations to
determine a unique solution for the problem. There are n subintervals, and in each
subinterval we have to determine a polynomial of degree 3. Each such polynomial
has 4 coeﬃcients, which leaves us with 4n coeﬃcients to determine. The interpolation
and continuity conditions (3.53) for s
i
(t
i
) and s
i
(t
i+1
) amount to 2n equations. The
continuity of the ﬁrst and the second derivatives (3.54) add 2(n−1) = 2n−2 equations.
Altogether we have 4n − 2 equations but 4n unknowns which leaves us with 2 degrees
of freedom. These indeed are two degrees of freedom that can be determined in various
ways as we shall see below.
We are now ready to compute the spline. We will use the following notation:
h
i
= t
i+1
−t
i
.
We also set
z
i
= s

(t
i
).
Since the second derivative of a cubic function is linear, we observe that s

i
(x) is the line
connecting (t
i
, z
i
) and (t
i+1
, z
i+1
), i.e.,
s

i
(x) =
x −t
i
h
i
z
i+1

x −t
i+1
h
i
z
i
. (3.55)
Integrating (3.55) once, we have
s

i
(x) =
1
2
(x −t
i
)
2
z
i+1
h
i

1
2
(x −t
i+1
)
2
z
i
h
i
+ ˜ c.
Integrating again
s
i
(x) =
z
i+1
6h
i
(x −t
i
)
3
+
z
i
6h
i
(t
i+1
−x)
3
+ C(x −t
i
) + D(t
i+1
−x).
The interpolation condition, s(t
i
) = y
i
, implies that
y
i
=
z
i
6h
i
h
3
i
+ Dh
i
,
i.e.,
D =
y
i
h
i

z
i
h
i
6
.
Similarly, s
i
(t
i+1
) = y
i+1
, implies that
y
i+1
=
z
i+1
6h
i
h
3
i
+ Ch
i
,
34
D. Levy 3.10 Spline Interpolation
i.e.,
C =
y
i+1
h
i

z
i+1
6
h
i
.
This means that we can rewrite s
i
(x) as
s
i
(x) =
z
i+1
6h
i
(x−t
i
)
3
+
z
i
6h
i
(t
i+1
−x)
3
+
_
y
i+1
h
i

z
i+1
6
h
i
_
(x−t
i
)+
_
y
i
h
i

z
i
6
h
i
_
(t
i+1
−x).
All that remains to determine is the second derivatives of s(x), z
0
, . . . , z
n
. We can set
z
1
, . . . , z
n−1
using the continuity conditions on s

(x), i.e., s

i
(t
i
) = s

i−1
(t
i
). We ﬁrst
compute s

i
(x) and s

i−1
(x):
s

i
(x) =
z
i+1
2h
i
(x −t
i
)
2

z
i
2h
i
(t
i+1
−x)
2
+
y
i+1
h
i

z
i+1
6
h
i

y
i
h
i
+
z
i
h
i
6
.
s

i−1
(x) =
z
i
2h
i−1
(x −t
i−1
)
2

z
i−1
2h
i−1
(t
i
−x)
2
+
y
i
h
i−1

z
i
6
h
i−1

y
i−1
h
i−1
+
z
i−1
h
i−1
6
.
So that
s

i
(t
i
) = −
z
i
2h
i
h
2
i
+
y
i+1
h
i

z
i+1
6
h
i

y
i
h
i
+
z
i
h
i
6
= −
h
i
3
z
i

h
i
6
z
i+1

y
i
h
i
+
y
i+1
h
i
,
s

i−1
(t
i
) =
z
i
2h
i−1
h
2
i−1
+
y
i
h
i−1

z
i
6
h
i−1

y
i−1
h
i−1
+
z
i−1
h
i−1
6
=
h
i−1
6
z
i−1
+
h
i−1
3
z
i

y
i−1
h
i−1
+
y
i
h i−1
.
Hence, for 1 i n −1, we obtain the system of equations
h
i−1
6
z
i−1
+
h
i
+ h
i−1
3
z
i
+
h
i
6
z
i+1
=
1
h
i
(y
i+1
−y
i
) −
1
h
i−1
(y
i
−y
i−1
). (3.56)
These are n−1 equations for the n+1 unknowns, z
0
, . . . , z
n
, which means that we have
2 degrees of freedom. Without any additional information about the problem, the only
way to proceed is by making an arbitrary choice. There are several standard ways to
proceed. One option is to set the end values to zero, i.e.,
z
0
= z
n
= 0. (3.57)
This choice of the second derivative at the end points leads to the so-called, natural
cubic spline. We will explain later in what sense this spline is “natural”. In this case,
we end up with the following linear system of equations
_
_
_
_
_
_
_
h
0
+h
1
3
h
1
6
h
1
6
h
1
+h
2
3
h
2
6
.
.
.
.
.
.
.
.
.
h
n−3
6
h
n−3
+h
n−2
3
h
n−2
6
h
n−2
6
h
n−2
+h
n−1
3
_
_
_
_
_
_
_
_
_
_
_
_
_
_
z
1
z
2
.
.
.
z
n−2
z
n−1
_
_
_
_
_
_
_
=
_
_
_
_
_
_
_
y
2
−y
1
h
1

y
1
−y
0
h
0
y
3
−y
2
h
2

y
2
−y
1
h
1
.
.
.
y
n−1
−y
n−2
h
n−2

y
n−2
−y
n−3
h
n−3
yn−y
n−1
h
n−1

y
n−1
−y
n−2
h
n−2
_
_
_
_
_
_
_
35
3.10 Spline Interpolation D. Levy
The coeﬃcients matrix is symmetric, tridiagonal, and diagonally dominant (i.e., |a
ii
| >

n
j=1,j=i
|a
ij
|, ∀i), which means that it can always be (eﬃciently) inverted.
In the special case where the points are equally spaced, i.e., h
i
= h, ∀i, the system
becomes
_
_
_
_
_
_
_
4 1
1 4 1
.
.
.
.
.
.
.
.
.
1 4 1
1 4
_
_
_
_
_
_
_
_
_
_
_
_
_
_
z
1
z
2
.
.
.
z
n−2
z
n−1
_
_
_
_
_
_
_
=
6
h
2
_
_
_
_
_
_
_
y
2
−2y
1
+ y
0
y
3
−2y
2
+ y
1
.
.
.
y
n−1
−2y
n−2
+ y
n−3
y
n
−2y
n−1
+ y
n−2
_
_
_
_
_
_
_
(3.58)
In addition to the natural spline (3.57), there are other standard options:
1. If the values of the derivatives at the endpoints are known, one can specify them
s

(t
0
) = y

0
, s

(t
n
) = y

n
.
2. The not-a-knot condition. Here, we require the third-derivative s
(3)
(x) to be
continuous at the points t
1
, t
n−1
. In this case we end up with a cubic spline with
knots t
0
, t
2
, t
3
, . . . , t
n−2
, t
n
. The points t
1
and t
n−1
no longer function as knots.
The interpolation requirements are still satisﬁed at t
0
, t
1
, . . . , t
n−1
, t
n
. Figure 3.9
shows two diﬀerent cubic splines that interpolate the same initial data. The spline
that is plotted with a solid line is the not-a-knot spline. The spline that is plotted
with a dashed line is obtained by setting the derivatives at both end-points to
zero.
3.10.2 What is natural about the natural spline?
The following theorem states that the natural spline cannot have a larger L
2
-norm of
the second-derivative than the function it interpolates (assuming that that function has
a continuous second-derivative). In fact, we are minimizing the L
2
-norm of the second-
derivative not only with respect to the “original” function which we are interpolating,
but with respect to any function that interpolates the data (and has a continuous second-
derivative). In that sense, we refer to the natural spline as “natural”.
Theorem 3.22 Assume that f

(x) is continuous in [a, b], and let a = t
0
< t
1
< · · · <
t
n
= b. If s(x) is the natural cubic spline interpolating f(x) at the knots {t
i
} then
_
b
a
(s

(x))
2
dx
_
b
a
(f

(x))
2
dx.
Proof. Deﬁne g(x) = f(x) −s(x). Then since s(x) interpolates f(x) at the knots {t
i
}
their diﬀerence vanishes at these points, i.e.,
g(t
i
) = 0, 0 i n.
36
D. Levy 3.10 Spline Interpolation
(t
0
, f(t
0
))
(t
1
, f(t
1
))
(t
2
, f(t
2
))
(t
3
, f(t
3
))
(t
4
, f(t
4
))
x
Figure 3.9: Two cubic splines that interpolate the same data. Solid line: a not-a-knot
spline; Dashed line: the derivative is set to zero at both end-points
Now
_
b
a
(f

)
2
dx =
_
b
a
(s

)
2
dx +
_
b
a
(g

)
2
dx + 2
_
b
a
s

g

dx. (3.59)
We will show that the last term on the right-hand-side of (3.59) is zero, which will
conclude the proof as the other two terms on the right-hand-side of (3.59) are
non-negative. Splitting that term into a sum of integrals on the subintervals and
integrating by parts on every subinterval, we have
_
b
a
s

g

dx =
n

i=1
_
t
i
t
i−1
s

g

dx =
n

i=1
_
(s

g

)
¸
¸
¸
¸
t
i
t
i−1

_
t
i
t
i−1
s

g

dx
_
.
Since we are dealing with the “natural” choice s

(t
0
) = s

(t
n
) = 0, and since s

(x) is
constant on [t
i−1
, t
i
] (say c
i
), we end up with
_
b
a
s

g

dx = −
n

i=1
_
t
i
t
i−1
s

g

dx = −
n

i=1
c
i
_
t
i
t
i−1
g

dx = −
n

i=1
c
i
(g(t
i
)−g(t
i−1
)) = 0.

We note that f

(x) can be viewed as a linear approximation of the curvature
|f

(x)|
(1 + (f

(x))
2
)
3
2
.
37
3.10 Spline Interpolation D. Levy
From that point of view, minimizing
_
b
a
(f

(x))
2
dx, can be viewed as ﬁnding the curve
with a minimal |f

(x)| over an interval.
38
{"found_math": false, "script_math_tex": 0, "script_math_asciimath": 0, "math_annotations": 0, "math_alttext": 0, "mathml": 0, "mathjax_tag": 0, "mathjax_inline_tex": 0, "mathjax_display_tex": 0, "mathjax_asciimath": 0, "img_math": 0, "codecogs_latex": 0, "wp_latex": 0, "mimetex.cgi": 0, "/images/math/codecogs": 0, "mathtex.cgi": 0, "katex": 0, "math-container": 0, "wp-katex-eq": 0, "align": 0, "equation": 0, "x-ck12": 0, "texerror": 0}