zar97/09/03

Numerical Methods For Engineers

page 4-21

4.2 Interpolation
The shortest way is not necessarily the safest way.- zar97
Interpolation - a process of determining an approximating value between precise data
points. Two situations may arise
(1) Given some function f(x) and an interval I, find the approximation to f(x) for any x in I.
(2) Given a set of n data points (xi , yi), i =1,2,... n, estimate the value of f(x) at some
point x not among the xis.
Three methods to be investigated: (a) Newtons interpolating polynomial, (b) Lagrange
interpolating polynomial, and (c) Spline interpolation (fitting data in a piecewise fashion).

The most common method used for this purpose is polynomial interpolation. The
Polynomial of nth degree has a general form as

y a0 a1x a2 x 2 . . . an x n

(4.2-1)

For (n+1) data points, there is only one polynomial of order n that may pass through all the
points. For example,

However, there are several mathematical formats in which the polynomial can be
expressed: Newton and Lagrange polynomials are two alternatives suited for computer
implementation.

zar97/09/03

Numerical Methods For Engineers

page 4-22

4.2.1 Newtons Divided-Difference Interpolating Polynomials

a. First-order interpolating polynomial : Linear interpolation - connecting two data points
with a straight line.

Using similar triangle, we can derive

p1 ( x) f ( x0 )

f ( x1 ) f ( x0 )
( x x0 )
( x1 x0 )

(4.2-2)

f
is a finite-divided-difference approximation of the first derivative. In general,
x
the smaller the interval between the data points, the better the approximation.
Note

Example 4.2-1:

For f(x) = ln x, estimate ln(2) using linear interpolation in the interval

(i) between x= 1 and x=5 and (ii) between x=1 and x= 3.

Solution:
x
1
5

f(x)
0
1.609638

x
1
3

f(x)
0
1.098612

Using linear interpolation formula,

p1 ( 2 ) 0

1609438
.
0
( 2 1 )= 0.402360
51

p1 ( 2 ) 0

1.098612 0
( 2 1 ) = 0.549306
31

Percent true error, respectively (i) t = 42.0 % (ii) t = 20.8 %

zar97/09/03

Worksheet 4.2-

Numerical Methods For Engineers

page 4-22.1

Estimate p(4) using linear interpolation in the interval between x= 2 and x=6
and compute the percent true error knowing the exact value f(4) = 0.6021.
x
2
6

f(x)
0.3010
0.7782

Solution:

Worksheet 4.2-

Estimate p(4) using quadratic for the following set of data also compute the
percent true error.
x
1
2
6

Solution:

f(x)
0
0.3010
0.7782

zar97/09/03

Worksheet 4.2-

Numerical Methods For Engineers

Estimate p(4) using cubic interpolation for the following set of data also
compute the percent true error.
x
1
2
6
10

Solution:

page 4-26.1

f(x)
0
0.3010
0.7782
1

zar97/09/03

Numerical Methods For Engineers

page 4-23

b. Second-order interpolating polynomial : Quadratic interpolation - connecting three

data points with a curvature.
A second-order polynomial for this purpose can be conveniently expressed as

p2 ( x) b0 b1 ( x x 0 ) b2 ( x xo )( x x1 )

(4.2-3)

Note eq.(4.2-3) and (4.2-1) are equivalent. Multiplying the terms in eq.(4.2-3), we can show
that
a o b0 b1 x 0 b2 x 0 x1
a1 b1 b2 x 0 b2 x1
a 3 b2

The coefficients in eq.(4.2-3) can be solved in a simple manner:

at

x = x0 ,
b0 = f(x0)

(4.2-4)

x = x1 ,
b1

f ( x1 ) f ( x 0 )
x1 x 0

(4.2-5)

x = x2 ,
b2

f ( x2 ) f ( x1 )
x2 x1

f ( x1 ) f ( x0 )
x1 x0

x2 x0

(4.2-6)

zar97/09/03

Example 4.2-2:

Numerical Methods For Engineers

page 4-24

Fit a second-order polynomial to three points used in example 4.2-1 and

use the polynomial to evaluate ln(2).

Solution:
x
1
3
5

f(x)
0
1.098612
1.609638

Eq.(4.2-4) yields
b0 = 0
Eq.(4.2-5) gives
b1

1.098612 0
0549306
.
31

and eq.(4.2-6) gives

1.609638 1.098612
0.549306
5

3
b2
0.073448
51
The quadratic polynomial is then

p2(x) = 0.549306(x-1) - 0.073448(x-1)(x-3)

and at x = 2,
p(2) = 0.622754
which has a percent true error of t = 10.2 %. Thus, the quadratic interpolation improves
the result as compared to the linear interpolation.

zar97/09/03

Numerical Methods For Engineers

page 4-25

c. General Form of Newtons Interpolating polynomial : - connecting (n+1) data points

with nth-order polynomial.
The nth-order polynomial is

pn ( x) b0 b1 ( x x0 ) b2 ( x x0 )( x x1 ) . . . bn ( x x0 )( x x1 ) . . . ( x xn 1 ) (4.2-7)
The (n+1) data points, x0, x1, x2, ... , xn can be used to determine the coefficients b 0, b1, b2,
... , and bn in the following manner
b0 = f (x0)
b1 = f [x1, x0]
b2 = f [x2, x1, x0]

(4.2-8)
(4.2-9)
(4.2-10)

.
.
.

bn = f [xn, xn-1, xn-2, ... , x1, x0]

where the bracketed function evaluations are finite divided differences (FDD).
For example,
f ( xi ) f ( x j )
1st FDD:
f [ xi , x j ]
xi x j
f [ xi , x j ] f [ x j , x k ]

2nd FDD:

f [ xi , x j , x k ]

and nth FDD:

f [ x n , x n 1 ,..., x1 , x 0 ]

xi x k
f [ x n , xn 1 ,... , x1 ] f [ x n 1 , x n 2 ,... , x 0 ]
xn x0

(4.2-11)

(4.2-12)

(4.2-13)

(4.2-14)

th

Then the n -order Newtons Interpolating polynomial is

pn ( x) f ( x 0 ) f [ x1 , x 0 ]( x x 0 ) f [ x 2 , x1 , x 0 ]( x x 0 )( x x1 ) . . .
f [ x n , xn 1 ,... , x1 , x 0 ]( x x 0 )( x x1 )( x x 2 )... ( x xn 1 )

Remark:

(4.2-15)

1. Data points may not be equally spaced or x values be in ascending order.

2. Eq.(4.2-12) through (4.2-14) are recursive i.e. higher-order differences are
composed of lower-order differences which is very useful in computer
implementation.

zar97/09/03

Numerical Methods For Engineers

page 4-26

Determine a third-order Newtons interpolating polynomial to problem in

example 4.2-2 with x3= 4 and use the result to evaluate ln(2).

Example 4.2-3:

Solution:
i
0
1
2
3

xi
1
3
5
4

f(xi)
0
1.098612
1.609638
1.386294

st

1 FDD

nd

2 FDD

rd

3 FDD

With n = 3,

p3 ( x) b0 b1 ( x x 0 ) b2 ( x x 0 )( x x1 ) b3 ( x x 0 )( x x1 )( x x3 )
The 1st FDD [eq.(4.2-12)]
1.098612 0
0549306
.
31
1609638
.
1.098612
f [ x 2 , x1 ]
0.255513
53
1386294
.
1.609638
f [ x3 , x2 ]
0.223344
45
f [ x1 , x 0 ]

The 2nd FDD [eq.(4.2-12)]

0.255513 0.549306
0.073448
51
0.223344 0.255513
f [ x 3 , x 2 , x1 ]
0.032169
43
f [ x 2 , x1 , x 0 ]

and eq.(4.2-6) gives

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

0.032169 ( 0.073448)
0.013760
4 1

Then
p3(x) = 0.549306(x-1) - 0.073448(x-1)(x-3) + 0.013760(x-1)(x-3)(x-5)
p(2) = 0.622754 + 0.04128 = 0.664034
which has a percent true error of t = 4.2 %. Thus, the cubic interpolation polynomial
improves the result as compared to the linear and quadratic interpolating polynomials.

zar97/09/03

Numerical Methods For Engineers

page 4-27

Errors of Newtons Interpolating Polynomials

Recall that the truncation error for the Taylor Series could be expressed as

Rn

f n 1 ( )
( x i 1 x i ) n 1
(n 1)!

where xi xi 1

An equivalent expression for the error for the nth-order interpolating polynomial may be
written as

f n 1 ( )
Rn
( x x 0 )( x x1 ). . .( x x n )
(n 1)!

where x 0 x n

(4.2-16)

This expression requires that the function is known and differentiable which is not usually
the case. An alternative form using a finite divided difference to approximate the (n+1)th
derivative is recommended.

Rn f [ x , xn , xn 1 , . . ., x0 ]( x x0 )( x x1 ) . . .( x xn )
The above error expression may be then be used to estimate the error if an additional data
point f(xn+1) is available, as in the form

Rn f [ xn 1 , x n , x n1 ,. . ., x 0 ]( x x 0 )( x x1 ). . .( x x n )

(4.2-17)

zar97/09/03

Numerical Methods For Engineers

page 4-28

Estimate the error in both a quadratic and a cubic interpolating polynomial in

the previous examples.

Example 4.2-4:

Solution:

xi

f(xi)

0
1
2
3

1
3
5
4

0
1.098612
1.609638
1.386294

1st FDD

2nd FDD

3rd FDD

4th FDD

zar97/09/03

Numerical Methods For Engineers

page 4-29

4.2.2 Lagrange Interpolating Polynomials

The Lagrange interpolating polynomial is a reformulation of the Newton polynomial that
avoids the calculation of finite divided differences. It can be represented as

pn ( x ) Li ( x ) f ( xi )

(4.2-18)

i 0

where
n

Li ( x )
j 0
j i

x xj
xi x j

(4.2-19)

with denotes the product of. For example,

(i) Linear version (n=1)

p1 ( x )

x x1
x x0
f ( x0 )
f ( x1 )
x0 x1
x1 x0

(ii) the second-order version (n=2)

x x1 x x 2
x x0 x x2
x x0 x x1
p2 ( x )

f ( x0 )

f ( x1 )

f ( x2 )
x0 x1 x0 x2
x1 x0 x1 x2
x2 x0 x2 x1

zar97/09/03

Example 4.2-5:

Numerical Methods For Engineers

page 4-30

For f(x) = ln x, estimate ln(2) using Lagrange interpolating polynomials of the

first and second order
i
0
1
2

x
1
3
6

f(x)
0
1.609638
1.791760

Solution:
First order,

p1 ( 2)

23
2 1
( 0)
(1.098612) 0.549306
1 3
31

Second order,

2 3 2 6
2 1 2 6
2 1 2 3
p2 ( 2)
.
)

(0)

(1098612

(1.791760) 0.612957
1 3 1 6
3 1 3 6
6 1 6 3

zar97/09/03

Numerical Methods For Engineers

page 4-31

Error Estimate of Lagranges Interpolating Polynomials

An error associated with the Lagrange interpolating polynomial can be estimated in a
similar manner as with the Newtons method and is given as
n

Rn f [ xn 1 , xn , xn 1 , . . ., x0 ] x xi

(4.2-20)

i 0

However, since this estimate contains the finite divided differences, it is rarely used with
the Lagrange interpolation.
Read section 12.3 and section 12.4 (p.384-387).

Remark:
1. Interpolation using both the Newton and the Lagrange polynomials may be used with
equally or arbitrarily spaced data.. However, the method with equally spaced data is
preferable since it may later be used to derive numerical integration formulas.
2. Most systems in engineering are notoriously ill-conditioned, especially when we attempt
to determine the coefficients of the polynomial in the conventional form. Therefore, lowerorder polynomials are recommended and results must be checked carefully.

zar97/09/03

Numerical Methods For Engineers

page 4-32

4.2.3 Spline Interpolation

So far we have used the nth-order polynomials to interpolate between (n+1) data points.
However, there are cases when the polynomials, especially of higher degree, lead to
erroneous results. Thus, an alternative approach connecting lower order polynomials to
subsets of data points is used. Such connecting polynomials are called spline functions.
Spline functions have special characteristics that the connections between adjacent curves
are visually smooth. They also provide superior approximation when there is abrupt change
in the function.
I. Linear Splines - connecting two points with a straight line successively. The first-order
splines for a set of ordered data points (x0, x1, x2,. . ., xn) can be defined as a set of (n-1)
linear functions,
f(x) = f(x0) + m0 (x - x0 )

x0 x x1

f(x) = f(x1) + m1 (x - x0 )
.
.
.
f(x) = f(xn-1) + mn-1 (x - xn-1 )

x1 x x2

(4.2-21)

xn-1 x xn

where

mi

f ( xi 1 ) f ( xi )
xi 1 xi

(4.2-22)

These equations may be used to evaluate the function at any point in the corresponding
interval. Note that the approach is similar to linear interpolation.

zar97/09/03

Example 4.2-6 :

Numerical Methods For Engineers

page 4-33

Fit the data with first-order splines and evaluate the function at x = 2.7.
x
f(x)

1
1

2
4

3
9

4
16

Solution:
The interval 2 x 3 contains the point x = 2.7, therefore, the slope

94
5
3 2

and the resulting linear spline,

f(x) = f(2) + 5 (x - 2)

2 x 3

f(2) = 4 + 5 (2.7 - 2) = 7.5

Note: The first-order splines are neither continuous in the first derivative nor smooth at the
connections (called a knot). This is the main disadvantage of the linear spline.

zar97/09/03

Numerical Methods For Engineers

page 4-34

II. Quadratic Splines. To ensure that the mth derivatives are continuous at the knots, we
have to use a spline of at least degree (m+1). Therefore, the quadratic splines will have
continuous first derivatives at the knots.
In this section, we wish to derive a second-order polynomial for each interval between data
points. The general form of the polynomial for each interval can be represented as

f i ( x ) ai x 2 bi x ci

(4.2-23)

For (n+1) data points (x0, x1, x2, . . . , xn), there are n intervals. In each interval, there are 3
unknown constants to be determined, consequently, 3n unknown constants (the as, bs,
and cs) to evaluate. They require 3n simultaneous equations or conditions to be solved.

Figure 4-5 : Quadratic splines over three intervals.

zar97/09/03

Numerical Methods For Engineers

page 4-35

These conditions are as follows

1. The value of the functions must be equal at the interior nodes, that is

f ( xi ) ai xi2 bi xi ci ai 1xi2 bi 1xi ci 1

for

i 1 to n 1.

(4.2-24)

These equations provide 2n - 2 conditions since only interior nodes are used.

2. The first and the last functions must pass through the end points, that is

f ( x0 ) a1 x02 b1 x0 c1

and

f ( xn ) an xn2 bn xn cn

(4.2-25)

This provides two additional conditions for a total of 2n conditions.

3. The first derivatives at the interior points must be equal, that is

2ai xi bi 2ai 1xi bi 1

for

i 1 to n 1.

(4.2-26)

This provides another n - 1 conditions for a total of 3n - 1.

4. Assume that the second derivative is zero at the first point, that is
a1 = 0.

(4.2-27)

Therefore, we now have 3n conditions which can be used to solve for 3n unknown
constants.
Note: the last condition restricts the function to be linear for the first two points.

zar97/09/03

Example 4.2-7 :

Numerical Methods For Engineers

page 4-36

Fit the data with the quadratic splines and evaluate the function at x = 2.7.

x
1
2
3
4

f(x)
1
4
9
16

Solution:
Since we have 4 data points, we need 3 intervals and 9 conditions to solve for the
quadratic splines.
1. Functions are equal at the interior nodes.

@ x = 2,

4a1 + 2b1 + c1 = 4a2 + 2b2 + c2 = 4

@ x = 3,

9a2 + 3b2 + c2 = 9a3 + 3b3 + c3 = 9

2. First and last functions pass through the end points.

@ x = 1,

a 1 + b 1 + c1 = 1

@ x = 4,

16a3 + 4b3 + c3 = 16

3. The first derivatives are equal at the interior nodes.

@ x = 2,

4a1 + b1 = 4a2 + b2

@ x = 3,

6a2 + b2 = 6a3 + b3

4. The second derivative at the first point is zero.

a1 = 0

zar97/09/03

Numerical Methods For Engineers

page 4-37

We now only have 8 unknowns since a1 = 0. We may express the preceding conditions in
a matrix form as

2
0

0
1

0
1

0 b1 4
4 2 1 0 0 0 c1 4

9 3 1 0 0 0 a2 9

0 0 0 9 3 1 b2 9

0 0 0 0 0 0 c2 1

0 0 0 16 4 1 a 3 16

4 1 0 0 0 0 b3 0

6 1 0 6 1 0 c3 0

0
0
0
1
0
0
0

Solving for the unknown constants yield

a1 = 0

b1 = 3

c1 = - 2

a2 = 2

b2 = - 5

c2 = 6

a3 = 0

b3 = 7

c3 = - 12

The quadratic splines for each interval are then,

f1 (x) = 3x - 2

1 x 2

f2 (x) = 2x2 - 5x + 6

2 x 3

f3 (x) = 7x - 12

3 x 4

The interpolation at x = 2.7,

f2 (2.7) = 2 (2.7)2 - 5 (2.7) + 6 = 7.08

and the exact value is 7.29.

zar97/09/03

Numerical Methods For Engineers

page 4-38

III. Cubic Splines. To ensure that the mth derivatives are continuous at the nodes, we have
to use a spline of at least degree (m+1). Therefore, the cubic splines will have continuous
first and second derivatives at the knots.
In this section, we wish to derive a third-order polynomial for each interval between data
points. The general form of the cubic spline for each interval can be represented as

f i (x ) ai x 3 bi x 2 ci x di

(4.2-28)

For (n+1) data points (x0, x1, x2, . . . , xn), there are n intervals. In each interval, there are 4
unknown constants to be determined, consequently, 4 n unknown constants (the as, bs,
cs and ds) to evaluate. They require 4n simultaneous equations or conditions.

Figure 4-6 : Cubic splines over three intervals.

zar97/09/03

Numerical Methods For Engineers

page 4-39

These conditions are as follows

1. The value of the functions are equal at the interior nodes, that is

f ( xi ) ai xi3 bi xi2 ci xi di ai 1xi3 bi 1xi2 ci 1 xi di 1 for i 1 to n 1.

(4.2-29)

These equations provide 2n - 2 conditions since only interior nodes are used.
2. The first and the last functions must pass through the end points, that is

f ( x0 ) a1 x03 b1x02 c1 x0 d1 and

f ( xn ) an xn3 bn xn2 cn xn di

(4.2-30)

This provides two additional conditions for a total of 2n conditions.

3. The first derivatives at the interior nodes must be equal, that is

3ai xi2 2bi xi ci 3ai 1 xi2 2bi 1xi ci 1

for

i 1 to n 1.

(4.2-31)

This provides another n - 1 conditions for a total of 3n - 1.

4. The second derivative at the interior nodes must be equal, that is

6ai xi 2bi 6ai 1 xi 2bi 1

for

i 1 to n 1.

(4.2-32)

This provides an additional n - 1 conditions for a total of 4n - 2.

5. Assume that the third derivatives at the first and last nodes are zero, that is

a1 = 0

and

an = 0.

We now have 4n conditions which can be used to solve for 4n unknown constants.
Note: the last condition restricts the function to be quadratic for the first dan last two points
and there are several assumptions could be made.

zar97/09/03

Example 4.2-8 :

Numerical Methods For Engineers

page 4-40

Fit the data with the quadratic splines and evaluate the function at x = 2.7.

x
1
2
3
4

f(x)
1
4
9
16

Solution:
Since we have 4 data points, we need 3 intervals and 12 conditions to solve for the cubic
splines.
1. Functions are equal at the interior nodes.

@ x = 2,

8a1 + 4b1 + 2c1 + d1 = 8a2 + 4b2 + 2c2 + d2 = 4

@ x = 3,

27a2 + 9b2 + 3c2 + d2 = 27a3 + 9b3 + 3c3 + d3 = 9

2. First and last functions pass through the end points.

@ x = 1,

a 1 + b 1 + c1 + d 1 = 1

@ x = 4,

64a3 + 16b3 + 4c3 + d3 = 16

3. The first derivatives are equal at the interior nodes.

@ x = 2,

12a1 + 4b1 + c1 = 12a2 + 4b2 + c2

@ x = 3,

27a2 + 6b2 + c2 = 27a3 + 6b3 + c3

4. The second derivative at the interior nodes are equal.

@ x = 2,

12a1 + 2b1 = 12a2 + 2b2

@ x = 3,

18a2 + 2b2 = 18a3 + 2b3

zar97/09/03

Numerical Methods For Engineers

page 4-41

5. The third derivative is zero at the end nodes.

a1 = 0

and

a3 = 0.

We now only have 10 unknowns since a1 = 0 and a3 = 0. We may express the preceding
conditions in a matrix form as

4
0

0
1

0
4

0
0

5 1

0 0

0 0

4 2

0 0 27 9 3

0 0

0 0

1 1

0 0

0 0

0 0

16

1 0

0 0 12 4

0 0 27 6 1
2 0

0 0 12 2

0 0 18 2 0

0 b1 4
0 0 c1 4

0 0 d1 9

3 1 a 2 9
0 0 b2 1

4 1 c2 16
1 0 d 2 0

1 0 b3 0
0 0 c3 0

0 0 d 3 0

Solving for the unknown constants yields

a1 = 0

b1 = 0.625

c1 = 1.125

d1 = -0.75

a2 = 0.25

b2 = - 1.125

c2 = 5.875

d2 = - 5.25

a3 = 0

b3 = 1.125

c3 = - 0.875

d3 = 1.5

The quadratic splines for each interval are then,

f1 (x) = 0.625x2 + 1.125x - 0.75

1 x 2

f2 (x) = 0.25x3 - 1.125x2 + 5.875x - 5.25

2 x 3

f3 (x) = 1.125x2 - 0.875x + 1.5

3 x 4

The interpolation at x = 2.7,

f2 (2.7) = 0.25 (2.7)3 - 1.125 (2.7)2 + 5.875 (2.7) - 5.25 = 7.332.