Piecewise polynomial interpolation

For certain x-values x1 x2 xn we are given the function values yi = f (xi ). In some cases below we will also assume
that we are additionally given some derivatives si = f 0 (xi ). We want to find an interpolating function p(x) which satisfies all
the given data and is hopefully close to the function f (x).
We could use a single interpolating polynomial p(x). But this is usually a bad idea: for a large value of n we will obtain
large oscillations.
We should only use an interpolating polynomial if we know that this will not be a problem and several of the following
conditions hold
the derivatives f (k) do not grow very fast (e.g., f (x) = sin x)
the points x1 , . . . , xn are close together, and we evaluate p(x) at a point x inside of the interval [x1 , xn ]
for equidistant nodes, we only evaluate p(x) for x near the center of the interval [x1 , . . . , xn ]
if we want to evaluate p(x) over a whole interval [a, b] we should choose x1 , . . . , xn as Chebyshev nodes for this interval.
In all other cases it is much better to use a piecewise polynomial: We break the interval [a, b] into smaller subintervals, and
use polynomial interpolation with low degree polynomials on each subinterval. Typically we choose polynomial degree of
about 3. This is a good compromise between small errors and control of oscillations.

Piecewise linear interpolation

We are given x-values x1 , . . . , xn and y-values yi = f (xi ) for i = 1, . . . , n. Let
yi+1 yi
Li (x) := yi +
(x xi ).
xi+1 xi

(1)

We then define p(x) as the piecewise linear function with

for x [xi , xi+1 ] :

p(x) = Li (x)

We then have from the error formula for polynomial interpolation with 2 points that with hi := xi+1 xi
for x [xi , xi+1 ] :

f 00 (t)
(x xi )(x xi+1 )
2!

 h2
| f (x) p(x)| 21 max  f 00 (t) i
4
t[xi ,xi+1 ]
f (x) p(x) =

since the function |(x xi )(x xi+1 )| has its maximum

hi
2

h2i in the midpoint of the interval [xi , xi+1 ].

We see that the interpolation error satisfies | f (x) p(x)| Ch2i where C = 18 maxt[x1 ,xn ] | f 00 (t)|. If we choose equidistant
points with hi = (b a)/(n 1) we have | f (x) p(x)| C(b a)2 /n2 , i.e., doubling the number of points reduces the error
bound by a factor of 4.
However, if the function f (x) has different behavior on different parts of the interval we can get better results by choosing
the points x1 , . . . , xn accordingly: If | f 00 (x)| is small in a certain region we can use a wider spacing hi ; if | f 00 (x)| is large in
another reason we should place the nodes more closely, so that hi is small there. In this way we can achieve a small overall
error


 00 
2
1


|
max f (x) p(x)| 8 max
hi max f (t)
x[x1 ,xn ]

i=1,...,n1

t[xi ,xi+1 ]

with a small number of nodes. We say the choice of the nodes x1 , . . . , xn is adapted to the behavior of the function f .
One advantage of piecewise linear interpolation is that the behavior of p resembles the behavior of f :
whereever the function f is increasing/decreasing, we have that the function p is increasing/decreasing
However, we have drawbacks:
the function p(x) is not smooth: it has kinks (jumps of p0 (x)) at the nodes x2 , . . . , xn1 in general
the error | f (x) p(x)| Ch2i for x [xi , xi+1 ] only decreases fairly slowly with decreasing spacing hi . We would rather
have a higher power like Ch4i .

Piecewise cubic Hermite interpolation

Both of these drawbacks can be fixed by using a piecewise cubic polynomial p(x).
We assume that we are given
x1 , . . . , xn
y1 , . . . , yn where yi = f (xi )
s1 , . . . , sn where si = f 0 (xi )
In this case we can construct on each interval [xi , xi+1 ] a cubic Hermite polynomial pi (x) with
p0 (xi ) = si ,

pi (xi ) = yi ,

p(xi+1 ) = yi+1 ,

p0 (xi+1 ) = si+1 .

We can construct this polynomial pi (x) using the divided difference method (using the nodes xi , xi , xi+1 , xi+1 ). We will now
explain a different way to express this polynomial pi (x) since this will be useful for cubic splines later. The idea is to write
pi (x) as the linear function Li (x) plus to extra terms. The function Li (x) already has the correct function values yi , yi+1 at the
endpoints, so the two extra terms have to be chosen such that we get the given derivatives si , si+1 at the endpoints. We let
pi (x) = Li (x) + A (x xi )(x xi+1 )2 + B (x xi )2 (x xi+1 )
with certain values A, B which we will figure out now. Note that the two functions (x xi )(x xi+1 )2 and (x xi )2 (x xi+1 )
have function value zero at both endpoints, and the following values of its derivatives:
q(x)
(x xi )(x xi+1 )2
(x xi )2 (x xi+1 )

q0 (xi )
2
hi
0

q0 (xi+1 )
0
h2i

q00 (xi )
4hi
2hi

q00 (xi+1 )
2hi
4hi

Note that
Li0 (x) =

yi+1 yi
=: si
xi+1 xi

and hence
!

p0 (xi ) = si + Ai h2i = si
and similarly Bi =

si+1 si
h2i

Ai =

si si
h2i

so that we have an explicit formula for the cubic Hermite polynomial pi (x):
pi (x) = Li (x) +

si si
si+1 si
(x xi )(x xi+1 )2 +
(x xi )2 (x xi+1 ).
2
2
hi
hi

(2)

We define the piecewise cubic Hermite polynomial p(x) by

for x [xi , xi+1 ] :

p(x) = pi (x)

(3)

Then we obtain from the error formula for polynomial interpolation with 4 points xi , xi , xi+1 , xi+1 that
for x [xi , xi+1 ] :

f (4) (t)
(x xi )2 (x xi+1 )2
4!

 h4
1
| f (x) p(x)| 24
max  f 00 (t) i
16
t[xi ,xi+1 ]
f (x) p(x) =

hi
2

h2i in the midpoint of the interval [xi , xi+1 ].



1
We see that the interpolation error satisfies | f (x) p(x)| Ch4i where C = 2416
maxt[x1 ,xn ]  f (4) (t). If we choose equidistant
points with hi = (b a)/(n 1) we have | f (x) p(x)| C(b a)4 /n4 , i.e., doubling the number of points reduces the error
bound by a factor of 16.
since the function |(x xi )(x xi+1 )| has its maximum

However, if the function f (x) has different

 (4) behavior on different parts of the interval we can get better results
 (4) by
 choosing



the points x1 , . . . , xn accordingly: If f (x) is small in a certain region we can use a wider spacing hi ; if f (x) is large in

another reason we should place the nodes more closely, so that hi is small there.
error

1
h4i max
max | f (x) p(x)| 2416 max
i=1,...,n1

x[x1 ,xn ]

t[xi ,xi+1 ]

In this way we can achieve a small overall



 (4) 
 f (t)

with a small number of nodes.

Again, a major advantage of using piecewise polynomials is that we can pick a nonuniform spacing of the nodes adapted to
the behavior of the function f .
The cubic Hermite spline has the following drawbacks:
We need all the derivatives si = f 0 (xi ) for i = 1, . . . , n. In many cases these values are not available.
We have that p0 (x) is continuous, but p00 (x) has jumps at the points x2 , . . . , xn1 in general. We would like to have a
smoother function p(x).

Complete cubic spline

The complete cubic spline fixes these two problems. We now assume that we are given
x1 , . . . , xn
y1 , . . . , yn where yi = f (xi )
s1 = f 0 (x1 ) and sn = f 0 (xn ),
i.e., we only need the derivatives at the endpoints (see the section Not-a-knot spline below if these are not available). If
these values are given, we can pick arbitrary numbers s2 , . . . , sn1 and obtain with (3), (2) a piecewise cubic function p(x)
which interpolates all the given data values.
How should we pick n 2 the numbers s2 , . . . , sn1 to obtain a nice function p(x)? We can actually use this freedom to
achieve a function p(x) where p00 (x) is continuous at the points x2 , . . . , xn1 : We want to pick the n 2 numbers x2 , . . . , xn1
such that the n 2 equations
p00i1 (xi ) = p00i (xi )
i = 2, . . . , n 2
(4)
are satisfied. Using (2) we obtain with Li00 (x) = 0 and the values of q00 (xi ), q00 (xi+1 ) from the table above
si si
si+1 si
(4hi ) +
(2hi )
2
hi
h2i
si si
si+1 si
p00i (xi+1 ) =
(2hi ) +
(4hi )
2
hi
h2i
p00i (xi ) =

Therefore (4) gives (substituting i 1 for i in the second equation above) for i = 2, . . . , n 1
si si
si+1 si
si1 si1
si si1
(4hi ) +
(2hi ) =
(2hi1 ) +
(4hi1 )
2
2
2
hi
hi
hi1
h2i1


2
4
4
2
6
6
si1 +
+
si + si+1 =
si1 + si
hi1
hi1 hi
hi
hi1
hi
Note that the value s1 in the last equation and the value sn in the last equation are given, and should therefore be moved to
the right hand side. Hence we obtain the tridiagonal linear system (after diving each equation by 2)



2

3 hs11 + hs22 hs11

2
1



s2
h1 + h2
h2

s3
s2
1
2
2
1

+
3

s3
h2
h3
h2
h2
h3
h3

.
..
..
..

. =
.
(5)

.
.
.
.
.




1
2
2
1

sn2
3 hn3
+ hn2

hn3
hn3 + hn2
hn2
n3

n2

1
2
2
s
n1
+
s

s
n2
n1
n
hn2
hn2
hn1
3 hn2 + hn1 hn1
This gives the following algorithm for finding the cubic spline interpolation:

 for i = 1, . . . , n 1: let hi := xi+1 xi , si :=

yi+1 yi
hi

define the matrix A on the left hand side of (5) and the vector b on the right hand side of (5)

s2

solve the tridiagonal linear system A ... = b using Gaussian elimination without pivoting
sn1
For a given point x [x1 , xn ] we evaluate the cubic spline as follows:
find the interval [xi , xi+1 ] containing x
evaluate pi (x)
using (2), (1)
In Matlab we can find the complete cubic spline as follows: yt = spline([x1 , . . . , xn ],[s1 , y1 , . . . , yn , sn ],xt)
Here xt is a vector of points where we want to evaluate the spline, and yt is the corresponding vector of function values.

Optimal energy property for complete cubic spline

It turns out that a complete cubic spline gives a smooth function p(x) without large oscillations. In fact, the complete
cubic spline is the optimal interpolating curve in a certain sense.
Historically, people constructing ships used thin flexible rulers made of wood (called splines) to find smooth curves
passing through given points (xi , yi ). For a thin piece of wood of length L one needs a certain energy to bend it into a curve
with curvature (s) along the arc length s [0, L]:
Z L

C (s)2 ds

E=
s=0

Here C is a stiffness constant. If one tries to pass a thin piece of wood through a number of points and allows it to relax it
will assume the shape with lowest possible energy E.
If we describe the curve by a function y = p(x) we have for small slopes p0 (x) that (s) p00 (x) and
E E0 :=

Z xn

C p00 (x)2 dx.

x=x1

It turns out that the complete cubic spline is the smoothest possible interpolating function in the following sense:
Among all functions p(x) (not only piecewise polynomials) satisfying
p0 (x1 ) = s1 ,

p(x1 ) = y1 , . . . , p(xn ) = yn ,
Z xn

the complete cubic spline has the lowest possible energy

p0 (xn ) = sn

p00 (x)2 dx.

x=x1

Not-a-knot cubic spline

Now assume that we are not given any derivatives values. We are given only x1 , . . . , xn and the function values y1 , . . . , yn . In
this case the best way to proceed is as follows: First drop x2 and xn1 and consider only the x-values x1 , x3 , x4 , . . . , xn3 , xn2 , xn
with the corresponding y-values. If we pick arbitrary values s1 , sn we can find the interpolating cubic spline function p(x) as
explained above. The function p(x) is a cubic function on the interval [x1 , x3 ] given by (2) (with x1 , x3 in place of xi , xi+1 )
(x2 is not a knot). Similarly p(x) is a cubic function on the interval [xn2 , xn ] given by (2) (with xn2 , xn in place of xi , xi+1 )
(xn1 is not a knot).
In order to determine s1 , sn we need two additional equations: We get them from the points (x2 , y2 ) and (xn1 , yn1 ) and
require
p(x2 ) = y2 ,
p(xn1 ) = yn1
The first equation depends on s1 , s3 . The last equation depends on sn2 , sn . We therefore obtain a tridiagonal linear system
for the unknowns s1 , s3 , s4 , . . . , sn3 , sn2 , sn . We solve this linear system using Gaussian elimination without pivoting and
obtain a cubic spline function called the not-a-not cubic spline.
In Matlab we can find the not-a-knot cubic spline as follows: yt = spline([x1 , . . . , xn ],[y1 , . . . , yn ],xt)
Here xt is a vector of points where we want to evaluate the spline, and yt is the corresponding vector of function values.