CHAPTER 3

Interpolation
Exercises 3.1
1 (a) The Lagrange interpolating polynomial through the points (0, 1), (2, 3), (3, 0) is
P(x) = 1
(x 2)(x 3)
(0 2)(0 3)
+ 3
(x 0)(x 3)
(2 0)(2 3)
+ 0
(x 0)(x 2)
(3 0)(3 2)
1 (b) The Lagrange interpolating polynomial is
P(x) = 1
(x + 1)(x 3)(x 5)
(2 + 1)(2 3)(2 5)
+ 1
(x + 1)(x 2)(x 5)
(3 + 1)(3 2)(3 5)
+ 2
(x + 1)(x 2)(x 3)
(5 + 1)(5 2)(5 3)
1 (c) The Lagrange interpolating polynomial is
P(x) = 2
(x 2)(x 4)
(0 2)(0 4)
+ 1
(x 0)(x 4)
(2 0)(2 4)
+ 4
(x 0)(x 2)
(4 0)(4 2)
2 (a) P(x) = 1 + x
4
3
x(x 2)
2 (b) P(x) =
1
3
(x + 1)
1
12
(x + 1)(x 2) +
1
24
(x + 1)(x 2)(x 3)
2 (c) P(x) = 2 +
3
2
x
3 (a) Compute the Newtons divided differences of the data points:
1 3
1
1 1 1
2 0
2 3 1
4
3 7
The interpolating polynomial is given by P(x) = 3 (x + 1) + (x + 1)(x 1).
3 (b) By Theorem 3.2, there is exactly one polynomial of degree 3 that passes through all four
points. This is the degree 2 polynomial in (a). No degree 3 polynomial exists.
3 (c) There are innitely many interpolating polynomials of degree 6, each of form
P(x) = 3 (x + 1) + (x + 1)(x 1) + (x + 1)(x 1)(x 2)(x 3)Q(x)
52 CHAPTER 3 INTERPOLATION
for an arbitrary degree 2 polynomial Q(x).
4 (a) P(x) = x +
2
3
x(x 1)(x 2)
4 (b) For example, P
1
(x) = P(x)+x(x1)(x2)(x3), P
2
(x) = P(x)+2x(x1)(x2)(x3)
4 (c) No, since P(4) = 20 from part (a).
5 (a) Compute the Newtons divided differences of the data points:
2 8
2
0 4 0
2 0
1 2 0
2
3 2
The polynomial P(x) = 8 2(x + 2) passes through the four data points.
5 (b) According to Theorem 3.2, there are no other degree 3 polynomials through the four data
points, but P
4
(x) = 8 2(x + 2) + c(x + 2)x(x 1)(x 3) interpolates for all c.
6 P(x) = 1+2(x1) (x1)(x2) +
1
2
(x1)(x2)(x3) +x(x1)(x2)(x3)(x4)
7 By denition, P(x) passes through the points (1, 0), ..., (10, 0), (12, 44). In Lagrange form, the
interpolating polynomial can be written
P(x) = 44
(x 1) (x 10)
(12 1) (12 10)
Substituting x = 0 gives
P(x) = 44
(0 1) (0 10)
(12 1) (12 10)
= 44
1 2 10
11 10 2
= 4
8 1118
9 (a) For the data points (1, 0), (2, 0), (3, 0), (4, 0), (5, 0), (6, 0), (7, 10), the Lagrange interpolat-
ing polynomial is
L(x) = 10
(x 1)(x 2) (x 6)
(7 1)(7 2) (7 6)
9 (b) The degree 6 polynomial asked for contains the seven data points in (a). According to
Theorem 3.2, there is only one polynomial that does so, the one computed above. Therefore
it sufces to check whether the L(x) from part (a) passes through (8, 70). Since
L(8) = 10
7 6 2
6 5 1
= 70
Section 3.1 53
the polynomial L(x) from part (a) also works for part (b).
10 40
11 Since the four points have distinct x-coordinates, Theorem 3.2 implies that only one poly-
nomial of degree 3 passes through them. This polynomial is obviously the parabola
y = ax
2
+ bx + c. Therefore no cubic polynomials pass through the four points.
12 No, a set of 5 points can be contained by only one polynomial of degree 4 or less.
13 (a) The Newtons divided difference triangle follows:
1800 280
0.06
1850 283 0.0010
0.16 0.000016
1900 291 0.0042
0.79
2000 370
The interpolating polynomial is
P(x) = 280 + 0.06(x 1800) + 0.001(x 1800)(x 1850)
+ 0.000016(x 1800)(x 1850)(x 1900)
Estimate the 1950 carbon dioxide concentration by
P(1950) = 280 + 0.06(150) + 0.001(150)(100) + 0.000016(150)(100)(50) = 316
13 (b) Estimate the 2050 concentration by
P(2050) = 280 + 0.06(250) + 0.001(250)(200) + 0.000016(250)(200)(150) = 465
14 (a) 48
14 (b) 49.66
15 (a) A reasonably simple approach to evaluating the Lagrange interpolating polynomial is to
rst form the n sums x x
1
, ..., x x
n
and n(n 1)/2 sums x
i
x
j
. Then each term of the
polynomial requires n2 multiplications in the numerator, n2 in the denominator, plus one
multiplication and one division (which we count as a multiplication), for a total of 2n 2 for
each of the n terms. Therefore there are n(2n2) multiplications and n+n(n1)/2+n1 =
1
2
n
2
+
3
2
n 1 additions.
15 (b) The nested Newton divided difference polynomial needs 2 adds and 1 multiplication in
each nesting, for a total of n 1 multiplications and 2n 2 additions.
Computer Problems 3.1
1 (a) The newtdd MATLAB code can be used to construct the coefcients of the degree one
interpolating polynomial. The commands
54 CHAPTER 3 INTERPOLATION
>> x=[1970 1990];y=[3707475887 5281653820];
>> c=newtdd(x,y,2);
>> nest(1,c,1980,x)
return 4494564853.5 as the 1980 population linear estimate, compared to the actual 4452584502,
a difference of only 400 million.
1 (b) Follow the steps in (a), but use the three data points and construct a degree two interpolating
polynomial. The resulting 1980 estimate is 4454831983.7, differing by only 2 million or so
from the actual population.
1 (c) Following the steps in (a) for all four data points results in an estimate of 4472888287.8,
about 20 million more than the actual population.
3 The MATLAB function could have the form:
function y0=polyinterp(x,y,x0)
len=length(x);
c=newtdd(x,y,len);
y0=nest(len-1,c,x0,x);
4 The fundamental domain for cosine can be chosen to be [0, /2].
5 (a) tan(

2
x) =
sin(

2
x)
cos(

2
x)
=
sin

2
cos x cos

2
sin x
cos

2
cos x + sin

2
sin x
=
cos x
sin x
.
5 (b) The tan x function is periodic with period , or in other words, tan(x + ) = tan x.
Therefore if we can compute tan x on the domain [

2
,

2
), all other x can be referred to this
interval by subtracting integer multiples of .
Likewise, the identity tan(

2
x) = 1/ tan x shows that for x satisfying

4
< x <

2
,
tan x = 1/ tan(

2
x) where

2
x lies in [0,

4
]. Furthermore, tan x for x in [

2
, 0] is the
negative of tan |x|.
Summarizing, given any real number, rst subtract enough integer multiples of to locate
a corresponding x in [

2
,

2
]. If x < 0, replace x by |x| and remember to attach a minus sign
at the nal step. If the positive x is less than

4
, compute tan from the fundamental domain,
and if x >

4
, set tan x = 1/ tan(

2
x).
5 (c) The following MATLAB code moves the input x to the fundamental domain [0, /4] and
carries out the evaluation of the interpolating polynomial.
function y=tan1(x)
b=pi
*
(0:3)/12;yb=tan(b);
c=newtdd(b,yb,4);
s=1;d=1;
x=mod(x+pi/2,pi)-pi/2;
if x < 0
s=-1;x=abs(x);
end
if x > pi/4
d=-1; x=pi/2-x;
Section 3.2 55
end
y=s
*
nest(3,c,x,b)d;
nest(3,c,x,b)
5 (d) The maximum error in the interval [0, /4] is approximately 0.003.
Exercises 3.2
1 (a) The Newton divided difference triangle is
0 0
2

2
1
4

0
The degree 2 interpolating polynomial is P
2
(x) =
2

x
4

2
x(x

2
).
1 (b) P
2
(

4
) =
2

4

4

4
(

4
) =
1
2
+
1
4
=
3
4
1 (c) The interpolation error formula for Theorem 3.3 applies to give

sin

4
P
2

4
3
4
6
=

3
128
0.242.
1 (d) The error is

sin

4
P
2

2
2

3
4

0.043,
less than the error bound in (c).
2 (a) P
2
(x) = (ln 2)(x 1)
1
6
(ln 2)(x 1)(x 2)
2 (b) P
2
(3) = 5(ln 2)/3 1.155245
2 (c) 2/3
2 (d) 0.0566
3 (a) According to Theorem 3.3,
|f(x) P
9
(x)| =

(x 0)(x
1
9
) (x 1)
10!

2
10
e
2c
where c is between 0 and 1. At x =
1
2
, the error is

1
2

P
9

1
2

9
18
7
18
5
18
3
18
1
18
1
18
3
18
5
18
7
18
9
18
10!
2
10
1 7.06 10
11
.
56 CHAPTER 3 INTERPOLATION
3 (b) Since f(
1
2
) = 1/e, we can use P
9
(
1
2
) as an approximation for 1/e. According to (a), since
7.06 10
11
0.5 10
9
, the approximation will have at least 9 correct decimal places.
4 (a) 3 5 7 9/5
7
0.0121
4 (b) 5 3 3 5/5
7
0.00288
5 By Theorem 3.3, the interpolation error at x is
|f(x) P(x)| |
(x 0.1)(x 0.2) (x 0.6)
6!
||f
(6)
(c)|
when c is between the smallest and largest of 0.1, 0.2, . . . , 0.6, and x. We are given no infor-
mation about f
(6)
(c). All things being equal, the ratio of the errors at x = 0.55 and x = 0.35
is
|(0.55 0.1)(0.55 0.2) (0.55 0.6)|
|(0.35 0.1)(0.35 0.2) (0.35 0.6)|
=
(0.45)(0.35)(0.25)(0.15)(0.05)(0.05)
(0.25)(0.15)(0.05)(0.05)(0.15)(0.25)
= 4.2.
6 0.00000714, making the assumption that f
(6)
(x) is similar in size to f
(8)
(x) on the interval
[0, 1].
Computer Problems 3.2
1 (a) Using the MATLAB program newtdd to generate the coefcients of the interpolating
polynomial results in
P
4
(x) = 1.433329 + 1.989870(x 0.6) + 3.258900(x 0.6)(x 0.7)
+ 3.680667(x 0.6)(x 0.7)(x 0.8)
+ 4.0000417(x 0.6)(x 0.7)(x 0.8)(x 0.9),
or in nested form,
P
4
(x) = 1.433329 + (x 0.6)(1.989870 + (x 0.7)(3.258900
+ (x 0.8)(3.680667 + (x 0.9)4.0000417))).
1 (b) Substitute into P
4
(x) to get P
4
(0.82) = 1.958910 and P
4
(0.98) = 2.612848.
1 (c) According to Theorem 3.3, the interpolation error satises
|e
x
2
P
4
(x)| =
|(x 0.6)(x 0.7)(x 0.8)(x 0.9)(x 1)|
5!
|f
(v)
(c)|
where
f
(v)
(x) = (32x
5
+ 160x
3
+ 120x)e
x
2
.
Section 3.2 57
For c 1, an upper bound for f
(v)
(c) is
|f
(v)
(c)| (32 + 160 + 120)e
1
= 312e 848.1.
Substituting into the error formula gives the error bounds
|e
(0.82)
2
P
4
(0.82)| 0.0000537
and
|e
(0.98)
2
P
4
(0.98)| 0.000217.
The actual errors, using the results of part (b), are
|e
(0.82)
2
P
4
(0.82)| 0.0000234
and
|e
(0.98)
2
P
4
(0.98)| 0.000107.
1 (d)
0.0002
0.0002
0.0001
0.0001
0.5 1.0
(a)
1 2
1
1
2
3
(b)
2 Write a for loop to subtract the output of sin1 from the MATLAB sin function for a grid of
x inputs, and plot the result:
6 3 3 6
0.003
0.003
58 CHAPTER 3 INTERPOLATION
3 The degree 9 interpolating polynomial from Newtons divided differences, written in nested
form, is
P
9
(x) = 67.052 + (x 1994)(0.9556 + (x 1995)(0.4195
+ (x 1996)(0.6883 + (x 1997)(0.03575 + (x 1998)(0.002175
+ (x 1999)(0.0123056 + (x 2000)(0.007915
+ (x 2001)(0.0028646 + (x 2002)(0.0007352)))))))))
and is plotted below.
1994 1999 2004
50
100
Substituting x = 2010 gives P
9
(2010) 1.952 10
12
barrels, which is mathematically
correct but nonsensical, and not consistent with the data shown. This is a typical example of
the Runge phenomenon, which prevents the degree-9 interpolating polynomial from being a
usable model for extrapolating this data
4 The degree 3 interpolating polynomial from Newtons divided differences, written in nested
form, is
P
3
(x) = 67.052 + (x 1994)(0.9556 + (x 1995)(0.4195 + (x 1996)(0.6883)))
and is plotted below.
1994 1999 2004
50
100
Section 3.3 59
Substituting x = 1998 gives the estimate 74.25810
6
for the 1998 oil production, a slight
overestimate of the actual value 73.400 10
6
. The Runge phenomenon has disappeared due
to the lower degree of the interpolating polynomial.
Exercises 3.3
1 (a) Using the formula x = cos(2i 1)/(2n) for i = 1, . . . , n, the interpolation nodes are
cos

12
, cos
3
12
, cos
5
12
, cos
7
12
, cos
9
12
, cos
11
12
.
1 (b) The general formula for Chebyshev nodes on [a, b] is (a+b)/2+((ba)/2) cos(2i 1)/n
for i = 1, . . . , n. The 4 nodes on [2, 2] are therefore 2 cos(2i 1)/8 for i = 1, . . . , 4 or
2 cos

8
, 2 cos
3
8
, 2 cos
5
8
, 2 cos
7
8
.
1 (c) Using the general formula from (b), the 6 nodes in [4, 12] are
8 + 4 cos

12
, 8 + 4 cos
3
12
, 8 + 4 cos
5
12
, 8 + 4 cos
7
12
, 8 + 4 cos
9
12
, 8 + 4 cos
11
12
.
1 (d) Using the general formula from (b), the 5 nodes in [0.3, 0.7] are
0.2 + 0.5 cos

10
, 0.2 + 0.5 cos
3
10
, 0.2 + 0.5 cos
5
10
, 0.2 + 0.5 cos
7
10
, 0.2 + 0.5 cos
9
10
.
2 (a) 1/32
2 (b) 2
2 (c) 128
2 (d) 1/512
3 To nd a degree 5 interpolating polynomial, we need n = 6 nodes. By Theorem 3.3, the
interpolation error is
e
x
Q
5
(x) =
(x x
1
) (x x
6
)
6!
f
(vi)
(c)
where 1 c 1. Using (3.20),
|e
x
Q
5
(x)|
|f
(vi)
(c)|
2
5
6!
=
e
c
2
5
6!

e
2
5
6!
1.1798 10
4
.
Since 0.11798 10
3
< 0.5 10
3
, Q
5
(x) and e
x
will agree to 3 decimal places for 1
x 1.
60 CHAPTER 3 INTERPOLATION
4 7.55 10
9
, 7 correct digits
5 To nd a degree 3 interpolating polynomial P
3
(x), n = 4 nodes are needed. By Theorem 3.3
and (3.20),
| sin x P
3
(x)| =
(x x
1
)(x x
2
)(x x
3
)(x x
4
)
4!
|f
(iv)
(c)|
1
2
3
4!
0.0052.
6 5.358 10
5
, 3 correct digits
7 To get 6 correct digits we want to control the interpolation error to less than 0.5 10
6
. Using
Theorem 3.3 and (3.20), the error for the degree d interpolation polynomial is
| ln x P
d
(x)| =
|(x x
1
) (x x
d+1
)|
(d + 1)!
|f
(d+1)
(c)|
(
e1
2
)
d+1
2
d
(d + 1)!
d!c
d1

(e 1)
d+1
2
2d+1
(d + 1)
where 1 c e. This expression is less than 0.5 10
6
when d = 14.
8 T
n
(0) = (1)
n/2
for n even, 0 for n odd.
9 (a) Using the recursion relation (3.19) with x = 1, we nd T
n+1
(1) = 2T
n
(1)
T
n1
(1). Since T
0
(1) = 1 and T
1
(1) = 1, this implies T
2
(1) = 1, T
3
(1) = 1,
T
4
(1) = 1, and in general T
n
(1) = (1)
n
. Therefore T
999
(1) = 1.
9 (b) T
1000
(1) = 1, from the discussion in (a).
9 (c) Using (3.19), T
n+1
(0) = T
n1
(0). We know T
0
(0) = 1 and T
1
(0) = 0, so T
n
(0) = 0 for
n odd, T
n
(0) = 1 if 4 divides evenly into n, and T
n
(0) = 1 if n is even and not a multiple of
4. Therefore T
999
(0) = 0.
9 (d) T
1000
(0) = 1, from the discussion in (c).
9 (e) Using (3.19), T
n+1
(
1
2
) = T
n
(
1
2
) T
n1
(
1
2
), and T
0
(
1
2
) = 1, T
1
(
1
2
) =
1
2
.
Note that T
2
(
1
2
) =
1
2
, T
3
(
1
2
) = 1, T
4
(
1
2
) =
1
2
, T
5
(
1
2
) =
1
2
establishes the pattern:
T
n
(
1
2
) 1 if 3 divides n and T
n
(
1
2
) =
1
2
otherwise. Thus T
999
(
1
2
) = 1.
9 (f) T
1000
(
1
2
) =
1
2
, from the discussion in (e).
Computer Problems 3.3
1 The base points in Program 3.3 must be changed to be the Chebyshev nodes, using for example
the MATLAB code
>> b=pi/4+(pi/4)
*
cos((1:2:7)
*
pi/8);
Then evaluating the degree 3 interpolating polynomial at a grid of points in the interval [2, 2]
gives the following plot, where the interpolating polynomial is the thick curve and y = sin x
is the thin curve.
Section 3.3 61
2 1 1 2
1
1
3 According to Theorem 3.3 and formula (3.20), the interpolation error of the degree n 1
interpolation polynomial on the fundamental domain [1, e] is at most

e1
2

|f
(n)
(c)|
2
n1
n!
where |f
(n)
(x)| = (n1)!x
n
| and 1 c e. The maximum value of |f(n)(c)| for 1 c e
is (n 1)!, so the interpolation error is bounded above by (e 1)
n
/2
2n1
n. Checking shows
that this expression is less than 0.5 10
10
for n = 26. Program 3.4 can be adapted to nd
the degree 25 interpolating polynomial at 26 Chebyshev nodes in [1, e].
The part of Program 3.4 that moves any input between 10
4
and 10
4
into the fundamental
domain must also be changed. Since e
10
> 10
4
, the exponent k satisfying e
k
x < e
k+1
can
be found using MATLAB code
k = 10;
while exp(k)>x
k = k-1;
end
Now that xe
k
is in [1, e], it can be substituted into the interpolation polynomial, and the result
is added to k to get the approximate natural logarithm.
4 The functions f(x) = e
|x|
(thin curve), the interpolating polynomial at evenly-spaced points
(dashed curve), and the interpolating polynomial at the Chebyshev nodes (thick curve) are
shown for (a) n = 10 and (b) n = 20. The evenly-spaced data points are plotted in (a), and
the Chebyshev points in (b).
62 CHAPTER 3 INTERPOLATION
1 1
1
2
3
(a)
1 1
1
2
3
(b)
Note the Runge phenomenon, especially evident in (b). The empirical interpolation errors
are plotted in (c) n = 10 and (d) n = 20. The dashed curve is the interpolation error for
evenly-spaced interpolation, and the solid curve is the Chebyshev interpolation error.
error
10
0
10
2
10
4
10
6
1 0 1
(c)
error
10
0
10
2
10
4
10
6
1 0 1
(d)
5 The MATLAB commands
>> b=2
*
(0:(n-1))/(n-1)-1;
or
>> b=cos((1:2:2
*
n-1)
*
pi/(2
*
n));
together with
>> yb=exp(-b.2);
>> c=newtdd(b,yb,n);
dene the coefcients of the degree n 1 interpolating polynomial for evenly-spaced or
Chebyshev nodes, respectively. Setting (a) n = 10 or (b) n = 20 gives the interpolation
curves below. For both n = 10 and 20, the evenly-spaced and Chebyshev interpolants are
Section 3.4 63
indistinguishable in the plot. In (a), the 10 evenly-spaced interpolation points are plotted and
in (b), the Chebyshev nodes are shown.
1 1
1
(a)
1 1
1
(b)
The reason the curves are indistinguishable is clear when the empirical interpolation errors
are plotted, in (c) n = 10 and (d) n = 20. The dashed curve is the interpolation error for
evenly-spaced interpolation, and the solid curve is the Chebyshev interpolation error. There is
a small hint of Runge phenomenon at the ends of the interval.
10
0
10
8
10
4
1 0 1
(c)
10
0
10
4
10
8
10
12
10
16
1 0 1
(d)
Exercises 3.4
1 (a) Calculating S

1
(x) = 3x
2
+1 and S

2
(x) = 3(x 1)
2
+6(x 1) +3, note that S

1
(1) = 4
and S

2
(1) = 3, so that S(x) is not a cubic spline.
1 (b) The derivatives are
S

1
(x) = 6x
2
+ 2x + 4
S

1
(x) = 12x + 2
S

2
(x) = 3(x 1)
2
+ 4(x 1) + 12
S

2
(x) = 6(x 1) + 14
64 CHAPTER 3 INTERPOLATION
Calculating S
1
(1) = 12, S

1
(1) = 12, S

1
(1) = 14, and S
2
(1) = 12, S

2
(1) = 12, S

2
(1) = 14,
we conclude that S(x) is a cubic spline.
2 (a) cubic spline
2 (b) not-a-knot only
3 (a) The second derivatives S

1
(x) =
9
2
x and S

2
(x) =
9
2
(x 1) + 2c must agree at x = 1,
implying that
9
2
= 2c, or c =
9
4
. Since S

1
(0) = 0 and S

2
(2) =
9
2
+
9
2
= 0, the spline is
natural. The spline fails the not-a-knot condition because S

1
(1) =
9
2
=
9
2
= S

2
(1), and is
not parabolically terminated because of the degree 3 terms.
3 (b) The second derivatives S

1
(x) = 8 and S

2
(x) = 2c must agree at x = 1, thus c = 4. The
spline is parabolically terminated because of the absence of degree 3 terms, and is not-a-knot
because S

1
(1) = 0 = S

2
(1). It fails to be natural because S

1
(0) = 8 = 0.
3 (c) The rst derivatives S

1
(x) =
3
2
+ 7x 3x
2
and S

2
(x) = c + (x 1) 3(x 1)
2
must
agree at x = 1, so c =
5
2
. Since S

1
(0) = 7 6(0) = 7 = 0, the spline is not natural. It is
not parabolically terminated because of the nonzero degree 3 terms. The spline is not-a-knot
because S

1
(1) = S

2
(1) = 6 and S

2
(2) = S

3
(2) = 6.
4 k
1
= 29/6, k
2
= 3/2, k
3
= 7/6; not-a-knot
5 By Theorem 3.9, there is exactly one natural cubic spline through any three points with distinct
x coordinates. The cubic spline S
1
(x) = S
2
(x) = x on [0, 1] and [1, 2], respectively, satises
the conditions.
6 S(x) = 1. Spline is also not-a-knot and natural.
7 (a) There are n = 3 given data points. First calculate
1
= x
2
x
1
= 10 = 1,
2
= x
3
x
2
=
2 1 = 1,
1
= y
2
y
1
= 1 0 = 1,
2
= y
3
y
2
= 4 1 = 3. Equation (3.30) has the
form

1 0 0
1 4 1
0 0 1

c
1
c
2
c
3

0
6
0

with solutions c
1
= 0, c
2
=
3
2
, c
3
= 0. Equation (3.28) yields
d
1
=
c
2
c
1
3
1
=
3
2
0
3(1)
=
1
2
d
2
=
c
3
c
2
3
2
=
0
3
2
3(1)
=
1
2
and equation (3.29) gives
b
1
=

1

1
3
(2c
1
+ c
2
) =
1
2
Section 3.4 65
b
2
=

2

2
3
(2c
2
+ c
3
) = 2
The cubic spline is therefore
S
1
(x) = 0 +
1
2
x + 0x
2
+
1
2
x
3
on [0, 1]
S
2
(x) = 1 + 2(x 1) +
3
2
(x 1)
2

1
2
(x 1)
3
on [1, 2].
7 (b) There are n = 3 data points. Calculate
1
= 2,
2
= 1,
1
= 0,
2
= 3. System (3.30) is

1 0 0
2 6 1
0 0 1

c
1
c
2
c
3

0
9
0

with solutions c
1
= 0, c
2
=
3
2
, c
3
= 0. Equation (3.28) implies d
1
=
1
4
, d
2
=
1
2
and equation
(3.29) implies b
1
= 1, b
2
= 2. The cubic spline is therefore
S
1
(x) = 1 (x + 1) +
1
4
(x + 1)
3
on [1, 1]
S
2
(x) = 1 + 2(x 1) +
3
2
(x 1)
2

1
2
(x 1)
3
on [1, 2].
8 (a)

1 +
5
3
x
1
6
x
3
on [0, 2]
3
1
3
(x 2) (x 2)
2
+
1
3
(x 2)
3
on [2, 3]
8 (b)

x
3
on [0, 1]
1 + 3(x 1) + 3(x 1)
2
(x 1)
3
on [1, 2]
9 Setting S
1
(1) = S
2
(1) implies 4 +b
1
= 1, or b
1
= 3. Setting S

1
(1) = S

2
(1) implies b
1
+3 =
b
2
, or b
2
= 0. Therefore S

1
(0) = b
1
+3(0)
2
= 3, and S

2
(3) = b
2
+6(31)6(31)
2
= 12.
10 True
11 (a) By Theorem 3.10, there is exactly one parabolically terminated cubic spline through the
three points. The parabola passing through the three points, used on both intervals [0, 1] and
[1, 2], will be the unique spline. By Newton divided differences:
0 2
2
1 0 2
2
2 2
66 CHAPTER 3 INTERPOLATION
The parabola is S(x) = 2 2x + 2x(x 1) = 2x
2
4x + 2.
11 (b) We know from Section 3.1 that there are innitely many degree 3 polynomials through
the three points. Each can be used as a not-a-knot spline. The cubic polynomials are S(x) =
2x
2
4x + 2 + cx(x 1)(x 2) for arbitrary c = 0.
12 One. Cubic interpolating polynomial P(x) = 3
1
3
(x 1)(x 3)
1
24
(x 1)(x 3)(x 4)
is a not-a-knot spline.
13 If a cubic spline is natural and parabolically terminated, then S
1
(x) = a
1
+ b
1
(x x
1
) +
c
1
(x x
1
)
2
and S

1
(x
1
) = 0. Differentiating gives S

1
(x) = 2c
1
, implying c
1
= 0, and S
1
(x)
is linear. The same argument applies to show that S
n1
(x) is linear.
14 No. The leftmost three data points must be collinear, and so must the rightmost three data
points.
15 We have S
1
(x) = 1, S
2
(x) = 1+b
2
x+c
2
x
2
+d
2
x
3
. Differentiating, S

2
(x) = b
2
+2c
2
x+3d
2
x
2
and S

2
(x) = 2c
2
+ 6d
2
x. The spline conditions imply 0 = S

1
(0) = S

2
(0) = b
2
and 0 =
S

1
(0) = S

2
(0) = 2c
2
; in other words, the conditions are satised provided that 0 = b
2
= c
2
,
for arbitrary d
2
. Thus S
2
(x) = 1 + d
2
x
3
works for any d
2
.
16 Clamped.
17 Since there are innitely many parabolas through n = 2 data points with distinct x coordinates,
there are innitely many parabolically terminated cubic splines. Thus existence holds but
uniqueness fails.
18 n = 2: There are innitely many parabolas through two arbitrary points with x
1
= x
2
; each
is a not-a-knot cubic spline. n = 3: Innitely many exist. Let P
2
(x) be the parabola through
the three points. Then S(x) = P
2
(x) +c(x x
1
)(x x
2
)(x x
3
) is a not-a-knot cubic spline
for each c.
Computer Problems 3.4
1 (a) Programs 3.5 and 3.6 can be used to compute the spline coefcients and plot the spline,
respectively. The command
>> c=splinecoeff([0 1 2 3],[3 5 4 1])
returns the spline coefcients

S
1
(x) = 3 +
8
3
x
2
3
x
3
on [0, 1]
S
2
(x) = 5 +
2
3
(x 1) 2(x 1)
2
+
1
3
(x 1)
3
on [1, 2]
S
3
(x) = 4
7
3
(x 2) (x 2)
2
+
1
3
(x 2)
3
on [2, 3]
.
The plot from
Section 3.4 67
>> c=splineplot([0 1 2 3],[3 5 4 1])
is shown as (a) below.
1 (b) Programs 3.5 and 3.6 can be used as in (a) to compute the spline coefcients and plot the
spline. The result is

S
1
(x) = 3 + 2.5629(x + 1) 0.5629(x + 1)
3
on [1, 0]
S
2
(x) = 5 + 0.8742x 1.6887x
2
+ 0.3176x
3
on [0, 3]
S
3
(x) = 1 0.6824(x 3) + 1.1698(x 3)
2
0.4874(x 3)
3
on [3, 4]
S
4
(x) = 1 + 0.1950(x 4) 0.2925(x 4)
2
+ 0.0975(x 4)
3
on [4, 5]
.
The plot of the natural cubic spline is shown in (b) below.
1 1 2 3 4 5
1
2
3
4
5
(a)
1 1 2 3 4 5
1
2
3
4
5
(b)
2 (a) The not-a-knot condition can be implemented by including the appropriate two lines in
Program 3.5. The spline is

S
1
(x) = 3 +
23
6
x 2x
2
+
1
6
x
3
on [0, 1]
S
2
(x) = 5 +
1
3
(x 1)
3
2
(x 1)
2
+
1
6
(x 1)
3
on [1, 2]
S
3
(x) = 4
13
6
(x 2) (x 2)
2
+
1
6
(x 2)
3
on [2, 3]
.
and is plotted below by Program 3.6.
2 (b) Use Program 3.5 with the not-a-knot condition as in (a). The spline is

S
1
(x) = 3 + 3.8867(x + 1) 2.1500(x + 1)
2
+ 0.2633(x + 1)
3
on [1, 0]
S
2
(x) = 5 + 0.3767x 1.3600x
2
+ 0.2633x
3
on [0, 3]
S
3
(x) = 1 0.6733(x 3) + 1.0100(x 3)
2
0.3367(x 3)
3
on [3, 4]
S
4
(x) = 1 + 0.3367(x 4) 0.3367(x 4)
3
on [4, 5]
.
and is plotted below as (b).
68 CHAPTER 3 INTERPOLATION
1 1 2 3 4 5
1
2
3
4
5
(a)
1 1 2 3 4 5
1
2
3
4
5
(b)
3 The spline is natural, and Program 3.5 calculates

S
1
(x) = 1 + 2.6607x 0.6607x
3
on [0, 1]
S
2
(x) = 3 + 0.6786(x 1) 1.9821(x 1)
2
+ 1.3036(x 1)
3
on [1, 2]
S
3
(x) = 3 + 0.6250(x 2) + 1.9286(x 2)
2
1.5536(x 2)
3
on [2, 3]
S
4
(x) = 4 0.1786(x 3) 2.7321(x 3)
2
+ 0.9107(x 3)
3
on [3, 4]
.
1 1 2 3 4 5
1
2
3
4
5
4 This is a curvature-adjusted cubic spline, computed by Program 3.5 after enabling the curvature
adjusted condition, and setting the second derivatives at the endpoints through the variables
v1=3 and vn=2. The spline is

S
1
(x) = 1 + 1.8006x + 1.5000x
2
1.3006x
3
on [0, 1]
S
2
(x) = 3 + 0.8988(x 1) 2.4018(x 1)
2
+ 1.5030(x 1)
3
on [1, 2]
S
3
(x) = 3 + 0.6042(x 2) + 2.1071(x 2)
2
1.7113(x 2)
3
on [2, 3]
S
4
(x) = 4 0.3155(x 3) 3.0268(x 3)
2
+ 1.3423(x 3)
3
on [3, 4]
.
Section 3.4 69
1 1 2 3 4 5
1
2
3
4
5
5 The cubic spline is clamped, and can be found by Program 3.5 with the clamped conditions
enabled and variables v1 = 0 and vn = 1.

S
1
(x) = 1 + 4.6786x
2
2.6786x
3
on [0, 1]
S
2
(x) = 3 + 1.3214(x 1) 3.3571(x 1)
2
+ 2.0357(x 1)
3
on [1, 2]
S
3
(x) = 3 + 0.7143(x 2) + 2.7500(x 2)
2
2.4643(x 2)
3
on [2, 3]
S
4
(x) = 4 1.1786(x 3) 4.6429(x 3)
2
+ 3.8214(x 3)
3
on [3, 4]
.
1 1 2 3 4 5
1
2
3
4
5
6 The cubic spline is clamped as in Computer Problem 5, and can be found by Program 3.5 with
the clamped conditions enabled and variables v1 = -2 and vn = 1.

S
1
(x) = 1 2x +
57
7
x
2

29
7
x
3
on [0, 1]
S
2
(x) = 3 +
13
7
(x 1)
30
7
(x 1)
2
+
17
7
(x 1)
3
on [1, 2]
S
3
(x) = 3 +
4
7
(x 2) + 3(x 2)
2

18
7
(x 2)
3
on [2, 3]
S
4
(x) = 4
8
7
(x 3)
33
7
(x 3)
2
+
27
7
(x 3)
3
on [3, 4]
.
70 CHAPTER 3 INTERPOLATION
1 1 2 3 4 5
1
2
3
4
5
7 Program 3.5 can be used to nd the cubic spline, using input vectors
>> x=pi
*
(0:4)/8;
>> y=cos(x);
The clamps v1 and vn should be chosen to match the derivatives at the endpoints, namely
S

(0) = sin 0 = 0 and S

(/2) = sin /2 = 1, respectively. The cubic spline is plotted

S
1
(x) = 1 0.5065x
2
+ 0.0327x
3
on [0,

8
]
S
2
(x) = 0.9239 0.3826(x

8
) 0.4679(x

8
)
2
+ 0.0931(x

8
)
3
on [

8
,

4
]
S
3
(x) = 0.7071 0.7070(x

4
) 0.3582(x

4
)
2
+ 0.1396(x

4
)
3
on [

4
,
3
8
]
S
4
(x) = 0.3827 0.9237(x
3
4
) 0.1937(x
3
4
)
2
+ 0.1639(x
3
4
)
3
on [
3
8
,

2
]
.
1
1
8 As in Computer Problem 7, Program 3.5 can be used to compute the clamped cubic spline
matching rst derivatives S

(0) = 1 and S

(/2) = 0. The result is

S
1
(x) = 0 + x 0.0006x
2
0.1639x
3
on [0,

8
]
S
2
(x) = 0.3827 + 0.9237(x

8
) 0.1937(x

8
)
2
0.1396(x

8
)
3
on [

8
,

4
]
S
3
(x) = 0.7071 + 0.7070(x

4
) 0.3582(x

4
)
2
0.0931(x

4
)
3
on [

4
,
3
8
]
S
4
(x) = 0.9237 + 0.3826(x
3
8
) 0.4679(x
3
8
)
2
0.0327(x
3
8
)
3
on [
3
8
,

2
]
.
Section 3.4 71
1
1
9 The clamped cubic spline can be calculated from Program 3.5 after setting the clamps S

(1) =
f

(1) = 1 and S

(3) = f

(3) = 1/3. The spline is

S
1
(x) = (x 1) 0.4638(x 1)
2
+ 0.1713(x 1)
3
on [1,
3
2
]
S
2
(x) = 0.4055 + 0.6647(x
3
2
) 0.2068(x
3
2
)
2
+ 0.0563(x
3
2
)
3
on [
3
2
, 2]
S
3
(x) = 0.6931 + 0.5001(x 2) 0.1224(x 2)
2
+ 0.0295(x 2)
3
on [2,
5
2
]
S
4
(x) = 0.9163 + 0.3998(x
5
2
) 0.0782(x
5
2
)
2
+ 0.0155(x
5
2
)
3
on [
5
2
, 3]
.
The maximum difference between the spline and the function f(x) = ln x can be determined
by a slight modication of Program 3.5. The maximum interpolation error on [1, 3] is approx-
imately 0.0005464.
10 Using an n = 48 point spline reduces the maximum interpolation error to 0.496 10
7
.
11 (a) The natural cubic spline through the four data points is

S
1
(x) = 3039585530 + 64621942(x 1960) + 21671(x 1960)
3
S
2
(x) = 3707475887 + 71123224(x 1970) + 650128(x 1970)
2
13542(x 1970)
3
S
3
(x) = 5281653820 + 80877678(x 1990) 162405(x 1990)
2
+ 5414(x 1990)
3
on [1960, 1970], [1970, 1990] and [1990, 2000] respectively. Substituting x = 1980 gives the
estimate 4470178717.
11 (b) Using the data points, we can estimate the derivatives at the left and right endpoints
as (3707475887 3039585530)/10 = 66789036 and (6079603571 5281653820)/10 =
79794975, respectively. The clamped cubic spline through the four points is

S
1
(x) = 3039585530 + 66789036(x 1960) 362212(x 1960)
2
+ 36221(x 1960)
3
S
2
(x) = 3707475887 + 70411154(x 1970) + 724424(x 1970)
2
15477(x 1970)
3
S
3
(x) = 5281653820 + 80815906(x 1990) 204186(x 1990)
2
+ 10209(x 1990)
3
on [1960, 1970], [1970, 1990] and [1990, 2000] respectively. Substituting x = 1980 gives the
estimate 4468552974. The natural cubic spline estimate misses the actual 1980 population by
17.6 million, while the clamped spline misses by 16 million, and is slightly closer.
72 CHAPTER 3 INTERPOLATION
12 (a) The natural cubic spline through the data is
S
1
(x) = 280 + 0.0613(x 1800) 5.2174 10
7
(x 1800)
3
S
2
(x) = 283 + 0.0574(x 1850) 7.826 10
5
(x 1850)
2
+ 4.2607 10
5
(x 1850)
3
S
3
(x) = 291 + 0.3691(x 1900) + 0.006313(x 1900)
2
2.1043 10
5
(x 1900)
3
on [1800, 1850], [1850, 1900] and [1900, 2000] respectively, which evaluates to 322.6 ppm at
x = 1950.
12 (b) The parabolically-terminated cubic spline through the data is
S
1
(x) = 280 + 0.0469(x 1800) + 2.6154 10
4
(x 1800)
2
S
2
(x) = 283 + 0.0731(x 1850) + 2.6154 10
4
(x 1850)
2
+ 2.9538 10
5
(x 1850)
3
S
3
(x) = 291 + 0.3208(x 1900) + 0.004692(x 1900)
2
on [1800, 1850], [1850, 1900] and [1900, 2000] respectively, which evaluates to 318.8 ppm at
x = 1950. The nearly-identical splines from (a) and (b) are shown below.
12 (c) For four data points, there is no difference between the not-a-knot cubic spline and the
degree three interpolating polynomial.
1800 1850 1900 1950 2000
200
300
400
(a)
1800 1850 1900 1950 2000
200
300
400
(b)
13 Program 3.6 can be adapted to plot the three different cubic splines. They are barely distin-
guishable in the following plot; the solid curve is natural, the dashed curve is not-a-knot, and
the dotted curve is parabolically-terminated.
1994 1999 2004
65
70
75
Section 3.5 73
Exercises 3.5
1 (a) Using the B ezier curve equations on page 164, one calculates b
x
= 3(0) = 0, c
x
=
3(20) 0 = 6, d
x
= 1006 = 5, b
y
= 3(20) = 6, c
y
= 3(02) 6 = 12, d
y
=
0 0 6 + 12 = 6. The B ezier curve is

x(t) = x
1
+ b
x
t + c
x
t
2
+ d
x
t
3
= 6t
2
5t
3
y(t) = y
1
+ b
y
t + c
y
t
2
+ d
y
t
3
= 6t 12t
2
+ 6t
3
.
1 (b) Similar to (a). One calculates b
x
= 3, c
x
= 3, d
x
= 3 and b
y
= 3, c
y
= 3, d
y
= 0. The
B ezier curve is

x(t) = 1 3t 3t
2
+ 3t
3
y(t) = 1 3t + 3t
2
.
1 (c) Similar to (a). One calculates b
x
= 0, c
x
= 3, d
x
= 2 and b
y
= 3, c
y
= 3, d
y
= 0. The
B ezier curve is

x(t) = 1 + 3t
2
2t
3
y(t) = 2 + 3t 3t
2
.
2 (a) (1, 1), (1, 2/3), (3, 1/3), (9, 1)
2 (b) (3, 2), (13/3, 5/3), (16/3, 5/3), (8, 5)
2 (c) (2, 1), (2, 2/3), (7/3, 1/3), (2, 2)
3 We use Example 3.16 to draw three line segments as a three-piece B ezier curve. Each piece is
characterized by four two-dimensional points:
(1, 2) (1, 2) (3, 4) (3, 4)
(3, 4) (3, 4) (5, 1) (5, 1)
(5, 1) (5, 1) (1, 2) (1, 2)
Translating to B ezier curves yields

x(t) = 1 + 2t
2
(3 2t)
y(t) = 2 + 2t
2
(3 2t)

x(t) = 3 + 2t
2
(3 2t)
y(t) = 4 3t
2
(3 2t)

x(t) = 5 4t
2
(3 2t)
y(t) = 1 + t
2
(3 2t)
4

x(t) = 0
y(t) = 15t
2
10t
3

x(t) = 15t
2
10t
3
y(t) = 5

x(t) = 5
y(t) = 5 15t
2
+ 10t
3

x(t) = 5 15t
2
+ 10t
3
y(t) = 0
5 There are innitely many B ezier curves satisfying the requirements. The data points are
(1, 0), (1, y
2
), (1, y
3
), (1, 0) which generate the B ezier curve

x(t) = 1 + 6t
2
4t
3
y(t) = 3y
2
t + (3y
3
6y
2
)t
2
+ (3y
2
3y
3
)t
3
.
74 CHAPTER 3 INTERPOLATION
At t =
1
2
, x(
1
2
) = 0, so we must solve 1 = y(
1
2
) =
3
8
(y
2
+ y
3
). Any combination of y
2
and
y
3
=
8
3
y
2
satises the requirements. The solutions are

x(t) = 1 + 6t
2
4t
3
y(t) = 3y
2
t + (8 9y
2
)t
2
+ (6y
2
8)t
3
.
for arbitrary y
2
. A particularly symmetric B ezier curve is given by setting y
2
= y
3
=
4
3
, or
y(t) = 4t 4t
2
.
6

x(t) =
1
2
t + 2t
2

3
2
t
3
y(t) = 1 4t
2
+ 3t
3
7 (a) For a space curve, the B ezier equations are applied to the x, y and z coordinates. Therefore
b
x
= 3, c
x
= 9, d
x
= 5, b
y
= 0, c
y
= 6, d
y
= 5, and b
z
= 0, c
z
= 3, d
z
= 3. The B ezier
space curve is

x(t) = 1 + 3t 9t
2
+ 5t
3
y(t) = 6t
2
5t
3
z(t) = 3t
2
3t
3
.
7 (b) Similar to (a). Calculate b
x
= 0, c
x
= 6, d
x
= 6, b
y
= 3, c
y
= 9, d
y
= 6, and
b
z
= 3, c
z
= 12, d
z
= 8. The B ezier space curve is

x(t) = 1 6t
2
+ 6t
3
y(t) = 1 + 3t 9t
2
+ 6t
3
z(t) = 2 + 3t 12t
2
+ 8t
3
.
7 (c) Similar to (a). Calculate b
x
= 3, c
x
= 12, d
x
= 10, b
y
= c
y
= d
y
= 0, and b
z
= 0, c
z
=
6, d
z
= 4. The B ezier space curve is

x(t) = 2 + 3t 12t
2
+ 10t
3
y(t) = 1
z(t) = 1 + 6t
2
4t
3
.
8 (a) (1, 1, 1), (1, 2/3, 4/3), (3, 1/3, 11/3), (9, 1, 8)
8 (b) (3, 2, 3), (13/3, 5/3, 10/3), (16/3, 5/3, 4), (8, 5, 4)
8 (c) (2, 1, 0), (2, 2/3, 0), (7/3, 1/3, 0), (2, 2, 2)
9 Differentiating the B ezier formula x(t) = x
1
+b
x
t+c
x
t
2
+d
x
t
3
gives x

(t) = b
x
+2c
x
t+3d
x
t
2
.
Substituting yields x

(0) = b
x
= 3(x
2
x
1
) and
x

(1) = b
x
+ 2c
x
+ 3d
x
= b
x
+ 2c
x
+ 3(x
4
x
1
b
x
c
x
)
= 3(x
1
x
1
) 3(x
2
x
2
) + 3(x
4
x
1
)
= 4(x
4
x
3
).
Section 3.5 75
Similarly, y

(0) = 3(y
2
y
1
) and y

(1) = 3(y
4
y
3
). Therefore the direction vector of the
parametric curve at t = 0 is 3(x
2
x
1
, y
2
y
1
), the direction from the rst endpoint to the
rst control point. The direction vector at t = 1 is 3(x
4
x
3
, y
4
y
3
), the direction from the
second control point to the second endpoint. Therefore the B ezier equations are designed to
allow the control points to control the tangent directions at the endpoints.
Computer Problems 3.5
1 A B ezier spline from Exercise 5 that satises the requirements is

x(t) = 1 + 6t
2
4t
3
y(t) = 4t 4t
2
.
A parametric plot of the spline is shown below.
1 1
1
2
2 A Bezier spline from Exercise 6 that satises the requirements is

x(t) =
1
2
t + 2t
2

3
2
t
3
y(t) = 1 4t
2
+ 3t
3
.
A parametric plot of the spline is shown below.
1
1
76 CHAPTER 3 INTERPOLATION
3 Examples of Bezier splines for the four letters are shown. Program 3.7 can be adapted to
accept, instead of input data from the mouse, n 4 arrays xlist and ylist holding rows
of four x and y data points, respectively. For example, in (b) below we dened
>> xlist=[1 1 1 1;1 3 3 1;1 3 3 1];
>> ylist=[1 1 5 5;5 5 3 3;3 3 1 1];
to represent the three splines dened by the four point groups
(1, 1) (1, 1) (1, 5) (1, 5)
(1, 5) (3, 5) (3, 3) (1, 3)
(1, 3) (3, 3) (3, 1) (1, 1)
1 2 3 4 5
1
2
3
4
5
(a)
1 2 3 4 5
1
2
3
4
5
(b)
1 2 3 4 5
1
2
3
4
5
(c)
1 2 3 4 5
1
2
3
4
5
(d)