Chapter 3

Interpolation and Curve Fitting

Engineering Numerical Methods Chapter 3: Interpolation and curve fitting

3.1 Introduction
The Figure below illustrates a set of tabular data in the form of a set of [x, f(x)]
pairs. The function f(x) is known at a finite set (actually eight) of discrete values
of x. The value of the function can be determined at any of the eight values of x
simply by a table lookup. However, a problem arises when the value of the
function is needed at any value of x between the discrete values in the table. The
actual function is not known and cannot be determined form the tabular values.
However, the actual function can be approximated by some known function, and
the value of the approximating function can be determined at any desired value
of x. This process, which is called interpolation, is the subject of this Chapter.
The discrete data of the figure below are actually values of the
function f ( x) = 1 / x .
In many applications, the values of the discrete data at the specific points
are not all that is needed. Values of the function at points other than the known
discrete points may be needed (i.e. interpolation). The derivative of the function
may be required (i.e. differentiation). The integral of the function may be of
interest (i.e. integration). These processes are performed by fitting an
approximating function to the set of discrete data and performing the desired
process on the approximating function.

a) interpolation b) differentiation c) integration

62
Engineering Numerical Methods Chapter 3: Interpolation and curve fitting

Many types of approximating functions exist. In fact, any analytical function can
be used as an approximating function. Three of the more common
approximating functions are:

1. Polynomials.
2. Trigonometric functions.
3. Exponential functions.

Approximating functions should have the following properties:

1. The approximating function should be easy to determine.

2. It should be easy to evaluate.
3. It should be easy to differentiate.
4. It should be easy to integrate.

Polynomials satisfy all four of these properties. Consequently, polynomial

approximating functions are used here to fit sets of discrete data for
interpolation, differentiation, and integration.
There are tow fundamentally different ways to fit a polynomial to a set of
discrete data:

1. Exact fits.
2. Approximate fits.

An exact fit yields a polynomial that passes exactly through all of the discrete
points, as illustrated in the figure below. This type of fit is useful for small sets
of smooth data. An approximate fit yields a polynomial that passes through the
set of data in the best manner possible, without being required to pass exactly
through any of the data points. Approximate fits are useful for large sets of
smooth data and small or large sets of rough data.

63
Engineering Numerical Methods Chapter 3: Interpolation and curve fitting

3.2 Direct Fit Polynomials

First let's consider a completely general procedure for fitting a polynomial to a

set of equally spaced or unequally spaced data. Given n+1 sets of data
xo , f ( xo ), x1 , f ( x1 ), ....... , xn , f ( xn ), which will be written as ( xo , f o ) ,
( x1 , f1 ) , ...... , ( xn , f n ) , determine the unique nth-degree polynomial Pn(x) that
passes exactly through the n+1 points:

Pn ( x) = ao + a1 x + a2 x 2 + .........+ an x n ....................................... (1)

For simplicity of notation, let f ( xi ) = f i . Substituting each data point into

Eq. (1) yields n+1 equations:

f o = ao + a1 xo + a2 xo + .........+ an xo
2 n
....................................... (2. 0)
f1 = ao + a1 x1 + a2 x1 + .........+ an x1
2 n
....................................... (2.1)
..........................................................
f n = ao + a1 xn + a2 xn + .........+ an xn
2 n
....................................... (2. n)

There are n+1 linear equations containing the n+1 coefficients a0 to an .

Equation (2) can be solved for ao to an by Gauss elimination. The resulting
polynomial is the unique nth-degree polynomial that passes exactly through the
n+1 data points. The direct fit polynomial procedure work for both equally
spaced data and unequally spaced data.

Example: Direct fit polynomials.

Solution: To illustrate interpolation by a direct fit polynomial, consider the

simple function y = f ( x) = 1 / x , and construct the following set of six
significant figure data:

x f(x)
3.35 0.298507
3.40 0.294118
3.50 0.285714
3.60 0.277778

64
Engineering Numerical Methods Chapter 3: Interpolation and curve fitting

Let's interpolate for y at x = 3.44 using linear, quadratic, and cubic interpolation.
The exact value is

1
y (3.44) = f (3.44) = = 0.290698...
3.44

Let's illustrate the procedure in detail for a quadratic polynomial:

P2 ( x) = a + b x + c x 2

To center data around x = 3.44, the first three points are used. Applying P2(x) at
each of these data points gives the following three equations:

0.298507 = a + b (3.35) + c (3.35) 2

0.294118 = a + b (3.40) + c (3.40) 2
0.285714 = a + b (3.50) + c (3.50) 2

Solving theses three equations for a, b, and c by Gauss elimination without

scaling or pivoting yields

P2 ( x) = 0.876561− 0.256080x + 0.0249333x 2

Substituting x = 3.44 in the polynomial gives

P2 (3.44) = 0.876561− 0.256080(3.44) + 0.0249333(3.44) 2 = 0.290697

The error is Error (3.44) = P2 (3.44) − f (3.44) = 0.290697 − 0.290698 = −0.000001 .

For a linear polynomial, use x = 3.4 and 3.5 to center that data around
x = 3.44. The resulting linear polynomial is

P1 ( x) = 0.579854 − 0.0840400 x

Substituting x = 3.44 in the polynomial gives P1 (3.44) = 0.290756

65
Engineering Numerical Methods Chapter 3: Interpolation and curve fitting

For cubic polynomial, all four points must be used. The resulting cubic
polynomial is

P3 ( x) = 1.121066− 0.470839x + 0.0878000x 2 − 0.00613333x 3

Substituting x = 3.44 in the polynomial gives P3 (3.44) = 0.290698 .

The results are summarized below, where the results of linear, quadratic, and
cubic interpolation and the errors, Error(3.44) = pn (3.44) − 0.290698, are
tabulated. The advantages of higher-degree interpolation are obvious.

P(3.44) = 0.290756 Linear interpolation Error = 0.000058

=0.290697 Quadratic interpolation =-0.000001
=0.290698 Cubic interpolation =0.000000

________________________________________________________________

The main advantage of direct fit polynomials is that the explicit form of the
approximating function is obtained, and the interpolation at several values of x
can be accomplished simply by evaluating the polynomial at each value of x.
The work required to obtain the polynomial does not have to be redone for each
value of x. A second advantage is that the data can be unequally spaced.
The main disadvantage of direct fit polynomials is that each time the degree
of the polynomial is changed, all of the work required to fit the new polynomial
must be redone.

66
Engineering Numerical Methods Chapter 3: Interpolation and curve fitting

3.2.1 Direct Multivariate Polynomial Approximation

Many problems arise in engineering and science when the dependant variable is
a function of two or more independent variables, for example, z = f(x.y) is a two
variable or bivariate, function. Such functions in general called multivariate
functions. When multivariate function is given by tabular data, multivariate
approximation is required for interpolation, differentiation, and integration.

Consider the bivariate function, z = f(x.y), and the set of tabular data presented
in the following table. The tabular data can be fit be a multivariate polynomial
of the form

z = f ( x, y) = a + bx + cy + dxy + ex 2 + fy 2 + gx2 y + hxy2 + ix3 + jy 3 + .......

The number of data points must equal the number of coefficients in the
polynomial. A linear bivariate polynomial in x and y is obtained by including the
first four terms in the equation. A quadratic bivariate polynomial in x and y is
obtained by including the first eight terms in the equation. The number of terms
in approximating polynomial increases rapidly as the degree of approximations
increases. This leads to ill-conditioned systems of linear equations for
determining the coefficients. Consequently, multivariate high-degree
approximation must be used with caution.

Example: Direct multivariate linear interpolation.

Solution: Consider the following values, use direct multivariate linear

interpolation to calculate z(x,y) = z(1100, 1225),

x
y 800 1000 1200
1150 1380.4 1500.2 1614.5
1200 1377.7 1499.0 1613.6
1250 1375.2 1497.1 1612.6

The form of the approximating polynomial is

z = f ( x, y ) = a + bx + cy + dxy

67
Engineering Numerical Methods Chapter 3: Interpolation and curve fitting

Substituting the four data points that bracket x = 1100 and y = 1225 into the
polynomial gives

1499.0 = a + (1000) b + (1200) c + (1000) (1200) d

1497.1 = a + (1000) b + (1250) c + (1000) (1250) d
1613.6 = a + (1200) b + (1200) c + (1200) (1200) d
1612.6 = a + (1200) b + (1250) c + (1200) (1250) d

Solving for a, b, c, and d by Gauss elimination yields

z = f ( x, y) = 1079.6 + 0.4650x − 0.1280y + 0.0900 10−3 xy

Substitute x = 1100 and y = 1225 in the approximating polynomial gives

z(1100, 1225)= 1555.5. The error in this result is Error = 1555.5 − 1556.0 = −0.5 .
The advantage of this approach is that we can use the same approximating
polynomial to evaluate other values of z(x,y),if required, without reevaluating
the polynomial coefficients.

3.3 Least Squares Approximation

3.3.1 Introduction
An approximate fit yields a polynomial that passes through the set of points in
the best possible manner without being required to pass exactly through any of
the points. Several definitions of best possible manner exist. Consider the set of
discrete points, xi , Y ( xi ) = ( xi , Yi ) , and the approximate polynomial y(x)
chosen to represent the set of discrete points, as illustrated in the figure below.
The discrete points do not fall on the approximating polynomial. The deviations
(i.e. distances) of the points from the approximating function must be minimized
in some manner.

68
 Engineering Numerical Methods Chapter 3: Interpolation and curve fitting

If the values of the independent variable xi are considered exact, then all the
deviation is assigned to the dependant variable Yi, and the deviation ei is the
vertical distance between Yi and yi=f(xi). Thus,

ei = Yi − yi

It is certainly possible that the values of Yi are quite accurate, but the
corresponding values of xi are in error. In that case, the deviation would be
measured by the horizontal distance illustrated in the above figure. If xi and Yi
both have uncertainties in their values, then the perpendicular distance between
a point and the approximating function would be the deviation. The usual
approach in the approximate fitting of tabular data is to assume that the
deviation is the vertical distance between a point and the approximating
function, as specified by the above equation.

Several best criteria are illustrated in the next figure for a straight line
approximation. From the figure we can see that the least squares criteria, in
which the sum of the squares of the deviations is minimized. The least squares
procedure yields a good compromise criterion for the best fit approximation.

Best fit criteria. (a) Minimize e

i . (b) Minimize e i . (c) Minimax. (d) Least squares.

69
Engineering Numerical Methods Chapter 3: Interpolation and curve fitting

3.3.2 The Straight Line Approximation

The simplest polynomial is a linear polynomial, the straight lint. Least squares
straight line approximations are an extremely useful and common approximate
fit, which is determined as follows. Given N data points, (xi ,Yi), fit the best
straight line through the set data. The approximating function is

y = a + b x ............................................................ (1)

At each value of xi, Eq. (1) gives

yi = a + b xi (i = 1, ............, N )

The deviation ei at each value of xi is

ei = Yi − yi (i = 1, ............, N )

The sum of the squares of the deviations defines the function S(a,b):

N N
S (a, b) =  (ei ) = (Yi − a − bx i ) 2
2

i =1 i =1

The function S(a,b) is a minimum when S / a = S / b = 0 . Thus,

S N
=  2 (Yi − a − bxi ) (−1) = 0
a i =1
......................... (2)
S N
=  2 (Yi − a − bxi ) (− xi ) = 0
b i =1

Dividing Eq. (2) by 2 and rearranging yields

N N
a N + b  xi =  Yi
i =1 i =1
N N N
a  xi + b  xi =  xi Yi
2

i =1 i =1 i =1

Which is called the normal equations of the least squares fit. They can be
solved for a and b by Gauss elimination.

70
Engineering Numerical Methods Chapter 3: Interpolation and curve fitting

Example: Least squares straight line approximation.

Solution: Consider the constant pressure specific heat for air at low
temperatures presented in the following table, where T is the temperature and Cp
is the specific heat. Determine a least squares straight line approximation for the
set of data:
T 300 400 500 600 700 800 900 1000
Cp 1.0045 1.0134 1.0296 1.0507 1.0743 1.0984 1.1212 1.1410
C p = a + bT

For this problem the equation becomes

8 8
8 a + b  Ti =  C p ,i
i =1 i =1
8 8 N
a  Ti + b  Ti =  Ti C p ,i
2

i =1 i =1 i =1

Evaluating the summations and substituting into the above equations gives
8 a + 5200 b = 8.5331
5200 a + 3,800,000 b = 5632.74

Solving for a and b be Gauss elimination without scaling or pivoting yields

C p = 0.933194 + 0.205298  10 −3 T

Substituting the initial values of T into this equation gives the results presented
in the next table and figure, which presents the exact data, the least squares
straight line approximation, and the percent error. The straight line is not a very
good approximation of the data.

71
Engineering Numerical Methods Chapter 3: Interpolation and curve fitting

3.3.3 High-Degree Polynomial Approximation

Given the N data points, (xi ,Yi), fit the best nth-degree polynomial through the
set of data. Consider the nth-degree polynomial:

y = ao + a1 x + a2 x 2 + .........+ an x n

The sum of the squares of the deviation is given by

N N
S (a0 , a1 ,..., an ) =  (ei ) 2 = (Yi − ao − a1 xi − .... − an xi ) 2
n

i =1 i =1

The function S (a0 , a1 ,..., an ) is a minimum when

S N
=  2 (Yi − ao − a1 xi − .... − a n xi ) (−1) = 0
n

 ao i =1
..............................................................................
S N
=  2 (Yi − ao − a1 xi − .... − a n xi ) (− xi ) = 0
n n

 a n i =1

Dividing by 2 and rearranging yields the normal equations:

N N N
ao N + a1  xi + .... + a n  xi =  Yi
n

i =1 i =1 i =1

.............................................................
N N N N
ao  xi + a1  xi + .... + a n  xi =  xi Yi
n n +1 2n n

i =1 i =1 i =1 i =1

And this can be solved for ao to an by Gauss elimination.

Example: Least squares quadratic polynomial approximation.

Solution: Determine a least squares quadratic polynomial approximation for this

set of data:

x 1000 1500 2000 2500 3000

Y 1.1410 1.2095 1.2520 1.2782 1.2955

72
Engineering Numerical Methods Chapter 3: Interpolation and curve fitting

y = a + b x + c x2
For this problem the polynomial becomes

a N + b  xi + c  xi =  Yi
2

a  xi + b  xi + c  xi =  xi Yi
2 3

a  xi +b  xi + c  xi =  xi Yi
2 3 4 2

Or may be represented in matrix form

 N

x x i i
 a    Yi 
2

3    
  xi x x i  b  =   xi Yi 
2
i
 x 2 x x 3 4
   xi 2 Yi 
 i i i  c 

 5 10  10 3 22.5  10 6  a   6.1762 
    3 
 10  10 22.5  10 55  10  b  =  12.5413  10 
3 6 9

22.5  10 6 55  10 9 12 
142.125  10   c  288.5186  10 6 


Solving for a, b, and c by Gauss elimination yields

y = 0.965460+ 0.211197 10−3 x − 0.0339143 10−6 x 2

Substituting the initial values of x into the approximating polynomial gives the
results presented in the next table and figure. The quadratic polynomial is
a reasonable approximation of the discrete data.

Y y

73
Engineering Numerical Methods Chapter 3: Interpolation and curve fitting

3.3.4 Multivariate Polynomial Approximation

Many problems arises in engineering and science where the dependant variable
is a function of two or more independent variable, for example, z = f(x,y) is a
two-variable, or bivariate, function.
Given the N data points, (xi, yi, Zi), fit the best linear bivariate polynomial
through the set of data. Consider the linear polynomial:

z = a + bx + cy

The sum of the squares of the deviations is given by

S (a, b, c) =  (ei ) 2 =  ( Z i − a − bx i − cyi ) 2

The function S(a,b,c) is a minimum when

S
=  2 ( Z i − a − bx i − cyi ) (−1) = 0
a
S
=  2 ( Z i − a − bx i − cyi ) (− xi ) = 0
b
S
=  2 ( Z i − a − bx i − cyi ) (− yi ) = 0
b

Dividing by 2 and rearranging yields the normal equations:

N N N
a N + b  xi + c  y i =  Z i
i =1 i =1 i =1
N N N N
a  x i + b  xi + c  x i y i =  x i Z i
2

i =1 i =1 i =1 i =1
N N N N
a  y i + b  xi y i + c  y i =  y i Z i
2

i =1 i =1 i =1 i =1

Which can be solved for a, b, and c by Gauss elimination.

74
Engineering Numerical Methods Chapter 3: Interpolation and curve fitting

A linear fit to a set of bivariate data may be inadequate. Consider the quadratic
bivariate polynomial:

z = a + b x + c y + d x 2 + e y 2 + f xy

The sum of the squares of the deviations is given by

S (a, b, c, d , e, f ) =  (ei ) 2 =  ( Z i − a − bx i − cyi − dxi − eyi − f xi yi ) 2

2 2

The function S(a,b,c,d,e,f) is a minimum when

S
=  2 ( Z i − a − bx i − cyi − dxi − eyi − f xi yi ) (−1) = 0
2 2

a
...................................................................................
S
=  2 ( Z i − a − bx i − cyi − dxi − eyi − f xi yi ) (− xi yi ) = 0
2 2

f

Dividing by 2 and rearranging yields the normal equations:

a N + b  xi + c  y i + d  xi + e  y i + f  xi y i =  Z i
2 2

a  xi + b  xi + c  xi yi + d  xi + e xi yi + f  xi yi =  xi Z i
2 3 2 2

a  y i + b  xi y i + c  y i + d  xi y i + e  y i + f  xi y i =  y i Z i
2 2 3 2

a  xi + b  xi + c  xi yi + d  xi + e xi yi + f  xi yi =  xi Z i
2 3 2 4 2 2 3 2

a  yi + b  xi yi + c  yi + d  xi yi + e yi + f  xi yi =  yi Z i
2 2 3 2 2 4 3 2

a  xi yi + b  xi yi + c  xi yi + d  xi yi + e xi yi + f  xi yi =  xi yi Z i
2 2 3 3 2 2

Which can be solved for a to f by Gauss elimination.

75
Engineering Numerical Methods Chapter 3: Interpolation and curve fitting

Example: Least squares quadratic bivariate polynomial approximation.

Solution: Consider the following values, use Least squares quadratic bivariate
polynomial to calculate z(x,y) = z(1100, 1225),

x
y 800 1000 1200
1150 1380.4 1500.2 1614.5
1200 1377.7 1499.0 1613.6
1250 1375.2 1497.1 1612.6

The form of the approximating polynomial is

z = a + b x + c y + d x 2 + e y 2 + f xy

 N x y x y x y  a 
    i 
2 2
Z 

i i i i i i

  xi x x y x x y x y     i i 
2 3 2 2
i i i i i i i i  b x Z
 2 
    yi Z i 
  i2 x y y x y y x y
2 2 3
y i i  c
 = 
i i i i i i

  xi x x y x x y x y i  d    xi Z i 
3 2 4 2 2 3 2
i i i i i i i
 3 
 e    yi 2 Z i 
  yi x y y x y y x y
2 2 3 2 2 4
i i 
i i i i i i
   
 xy  i i i 
2
 i i x y x y x y x y x y
i
2
i i i
2
i
3
i i i
3
i
2

i   
f   x y Z

Evaluating the summation and substituting in the matrix yields

 9 E0 10.800 E 3 9.000 E 3 12.975 E 6 9.240 E 6 10.800 E 6   a 

10.800 E 3 12.975 E 6 10.800 E 6 15.606 E 9 11.088 E 9 12.975 E 9   b 

 9.000 E 3 10.800 E 6 9.240 E 6 12.975 E 9 9.720 E 9 11.088 E 9   c 
  
12.975 E 6 15.606 E 9 12.975 E 9 18.792 E12 13.321 E12 15.606 E12  d 
 9.240 E 6 11.088 E 9 9.720 E 9 13.321 E12 10.450 E12 11.664 E12  e 
  
10.800 E 6 12.975 E 9 11.088 E 9 15.606 E12 11.664 E12 13.321 E12   f 
13.4703 E 3
16.1638 E 6
 
13.6118 E 6
= 
19.4185 E 9
14.1122 E 9
 
16.3337 E 9

76
 Engineering Numerical Methods Chapter 3: Interpolation and curve fitting

A problem arises for high-degree polynomials. The coefficients in the matrix are
varying over a range of several orders of magnitude, which gives rise to
ill-conditioned system. Normalizing each equation helps the situation. Double
precision calculations are required.

Each raw in the matrix should be normalized by the exponential term in the first
coefficient of each raw. Solving the normalized equations by Gauss eliminations
yields

z( x, y) = 914.033 + 0.645500x − 0.020500 y − 0.0000775x 2 − 0.000040y 2 + 0.0000825 xy

Evaluating z(1100, 1225) = 1556.3.

The error is, Error = 1556.3 - 1556.0 = 0.3, which is smaller than the error
incurred by interpolation using direct multivariate linear approximation.

77