DEBRE TABOR UNIVERSITY

Gafat Institute of Technology

DEPARTMENT OF MECHANICAL ENGINEERING

COURSE MODULE FOR THE COURSE

NUMERICAL METHODS (MEng2053)

BY: Endalkachew Addis

DEBRE TABOR, ETHIOPIA

2016 E.C

P.O. BOX: 272

NM Module 2016 E.C

CHAPTER ONE

1 INTRODUCTION TO NUMERICAL METHODS AND BASIC CONCEPTS

The chapter gives the use of Numerical methods and introduces the fundamental concepts of
Numerical methods such as Relative error, absolute error, approximate error, round off error and
truncation error. Also it introduces about significant digits, Taylor’s theorem etc.

Objective

After completing this chapter students will be able to

 Describe the use of Numerical

Methods (K, L2)
 Describe various errors (K, L2)
 Describe some preliminaries (K, L2)

Introduction
Mathematical models are an integral part in solving engineering problems. Many times, these
mathematical models are derived from engineering and science principles, while at other times
the models may be obtained from experimental data.
Mathematical models generally result in need of using mathematical procedures that
include but are not limited to
(A) differentiation,
(B) nonlinear equations,
(C) simultaneous linear equations,
(D) curve fitting by interpolation or regression,
(E) integration, and
(F) differential equations.

These mathematical procedures may be suitable to be solved exactly as you must have
experienced in the series of calculus or applied mathematics courses you have taken, but in most
cases, the procedures need to be solved approximately using numerical methods. Numerical or
Computational methods is a subject concerned with developing methods for computing, usually
approximately, the solution of mathematical expressions which can’t be solved analytically.

DTU-MED Compiled BY Endalkachw A. |Numerical 1

Methods
NM Module 2016 E.C

True or Absolute error:

True error denoted by Et is the difference between the true value (also called the exact value)
and the approximate value.
Relative true error:
Relative true error is denoted by t and is defined as the ratio between the true error and the true
value.
Approximate error:
Approximate error is denoted by E a and is defined as the difference between the present
approximation and previous approximation.

Relative approximate error:

Relative approximate error is denoted by a and is defined as the ratio between the approximate
error and the present approximation
Round off error:
1
A computer can only represent a number approximately. For example, a number like may be
3
represented as 0.333333 on a PC. Then the round off error in this case is
1
 0.333333  0.00000033 . Then there are other numbers that cannot be represented exactly.
3
For example,  and 2 are numbers that need to be approximated in computer calculations.

Truncation error:
Truncation error is defined as the error caused by truncating a mathematical procedure. For
example, the Maclaurin series for e x is given as

x2 x3
ex  1 x    ....................
2! 3!

This series has an infinite number of terms but when using this series to calculate e x , only a
finite number of terms can be used. For example, if one uses three terms to calculate e x , then
x2
ex  1 x  .
2!
the truncation error for such an approximation is

DTU-MED Compiled BY Endalkachw A. |Numerical 2

Methods
NM Module 2016 E.C
 x2 
Truncation error = e x  1  x  ,
 2! 

Significant Figures
The significant figures (sig figs) of a number are the digits from the first nonzero digit on the
left to either last digit (zero or non) on the right if there is a decimal, or last nonzero digit of a
number if there is no decimal.
• 2300 or 2.3×103 has 2 sig figs
• 2300. or 2.300×103 has 4 sig figs
• 2300.0 or 2.3000×103 has 5 sig figs
• 23,040 or 2.304×104 has 4 sig figs
• 0.035 or 3.5×10-2 has 2 sig figs
• 0.03500 or 3.500×10-2 has 4 sig figs

Taylor Series:
From Maclaurin series, we have
x2 x4 x6
cos( x)  1     (1)
2! 4! 6!
x3 x5 x7
sin( x)  x     (2)
3! 5! 7!
x 2 x3
ex  1 x    (3)
2! 3!
All the above expressions are actually a special case of Taylor series called the Maclaurin series.
Why are these applications of Taylor’s theorem important for numerical methods? Expressions
such as given in Equations (1), (2) and (3) give you a way to find the approximate values of these
functions by using the basic arithmetic operations of addition, subtraction, division, and
multiplication.

Questions

Question1: Define various errors?

Question2: Explain the working of numerical methods and define significant digits?
References
1. M.K. Jain, S.R.K Iyengar and R.K Jain “ Numerical Methods for Scientific and
engineering computation” Sixth edition, New age international publishers. P.P: 1-15.

2. Chapra, S.C and Canale, R.P., Numerical methods for engineers, McGraw – Hill
international edition, 6th edition, 2010.

DTU-MED Compiled BY Endalkachw A. |Numerical 3

Methods
NM Module 2016 E.C

CHAPTER TWO

ROOTS OF NON LINEAR EQUATIONS

This chapter discusses about various methods such as Bisection method, False position method,
Secant method, Newton- Raphson method their advantages and drawbacks for solving non linear
equations.

Learning objectives
After completing this course the students will be able to:
 Explain Bisection method and apply to solve engineering problems (K, L3)
 Derive False position method Apply False position method to
solve engineering problems (K, L3)
 Derive Secant method and apply to solve engineering problems (K, L3)
 Derive Newton- Raphson method and apply to solve engineering problems (K, L3)

Introduction
Most of Engineering Problems are complex and difficult to solve analytically. If the modelling
equation is linear and quadratic equation it is possible to solve analytical and using the general
formula of quadratic equation. But if our modelling equation is Non-linear equation we have no
general formula to solve our modelling problem. There are two common numerical method two
solve non-linear equation the first one is closed – domain (Bracketing) the other one is Open
domain method.

DTU-MED Compiled BY Endalkachw A. |Numerical 4

Methods
NM Module 2016 E.C

2.2. Bracketing Methods

2.2.1. Bisection Method of Solving a Nonlinear Equation

What is the bisection method and what is it based on?

One of the first numerical methods developed to find the root of a nonlinear equation f ( x)  0
was the bisection method (also called binary-search method). The method is based on the
following theorem.
Theorem
 An equation f ( x)  0 , where f (x) is a real continuous function, has at least one root
between x and xu if f ( x ) f ( xu )  0 (See Figure 1).
 Note that if f ( x ) f ( xu )  0 , there may or may not be any root between x and xu
(Figures 2 and 3).
 If f ( x ) f ( xu )  0 , then there may be more than one root between x and xu (Figure 4).
So the theorem only guarantees one root between x and xu .
Bisection method
Since the method is based on finding the root between two points, the method falls under
the category of bracketing methods.

Since the root is bracketed between two points, x and xu , one can find the mid-point,
xm between x and xu . This gives us two new intervals

1. x and xm , and
2. xm and xu .

Advantages of bisection method

a) The bisection method is always convergent. Since the method brackets the root, the method
is guaranteed to converge.
b) As iterations are conducted, the interval gets halved. So one can guarantee the error in the
solution of the equation.

Drawbacks of bisection method

a) The convergence of the bisection method is slow as it is simply based on halving the
interval.
b) If one of the initial guesses is closer to the root, it will take larger number of iterations
to reach the root.

DTU-MED Compiled BY Endalkachw A. |Numerical 5

Methods
NM Module 2016 E.C

2.2.2. FALSE-POSITION METHOD

Introduction
In chapter 03.03, the bisection method was described as one of the simple bracketing methods
for solving a nonlinear equation of the general form
f ( x)  0 (1)

f x 

f xU 

Exact root

xL
O xr xU x

f xL 

Figure 2.1 False-Position Method

The above nonlinear equation can be stated as finding the value of xsuch that Equation (1) is
satisfied. In the bisection method, we identify proper values of x L (lower bound value) and xU
(upper bound value) for the current bracket, such that

f ( xL ) f ( xU )  0 . (2)

The next predicted/improved root xm can be computed as the midpoint between x L and xU as

x L  xU
xm  (3)
2

DTU-MED Compiled BY Endalkachw A. |Numerical 6

Methods
NM Module 2016 E.C
The new upper and lower bounds are then established, and the procedure is repeated until the
convergence is achieved (such that the new lower and upper bounds are sufficiently close to each
other).
However, in the example shown in Figure 1, the bisection method may not be efficient because it
does not take into consideration that f ( x L ) is much closer to the zero of the function f (x) as
compared to f ( xU ) . In other words, the next predicted root xr would be closer to x L (in the
example as shown in Figure 1), than the mid-point between x L and xU . The false-position
method takes advantage of this observation mathematically by drawing a secant from the
function value at x L to the function value at xU , and estimates the root as where it crosses the x-
axis.

False-Position Method
Based on two similar triangles, shown in Figure 1, one gets
0  f ( x L ) 0  f ( xU )
 (4)
xm  x L xm  xU

From Equation (4), one obtains

xm  xL  f xU   xm  xU  f xL 
xU f xL   xL f xU   xm f xL   f xU 

The above equation can be solved to obtain the next predicted root xm as

xU f  x L   x L f xU 
xm  (5)
f  x L   f  xU 

The above equation, through simple algebraic manipulations, can also be expressed as
f  xU 
x m  xU  (6)
 f  x L   f  xU  
 
 x L  xU 
or
f x L 
xm  x L 
 f  xU   f  x L  
 
 xU  x L 

DTU-MED Compiled BY Endalkachw A. |Numerical 7

Methods
NM Module 2016 E.C

2.3.1. Newton-Raphson Method of Solving a Nonlinear Equation

Introduction
Methods such as the bisection method and the false position method of finding roots of a
nonlinear equation f ( x)  0 require bracketing of the root by two guesses. Such methods are
called bracketing methods. These methods are always convergent since they are based on
reducing the interval between the two guesses so as to zero in on the root of the equation.
In the Newton-Raphson method, the root is not bracketed. In fact, only one initial guess
of the root is needed to get the iterative process started to find the root of an equation. The
method hence falls in the category of open methods. Convergence in open methods is not
guaranteed but if the method does converge, it does so much faster than the bracketing methods.

Derivation
The Newton-Raphson method is based on the principle that if the initial guess of the root of
f ( x)  0 is at xi , then if one draws the tangent to the curve at f ( xi ) , the point xi 1 where the
tangent crosses the x -axis is an improved estimate of the root (Figure 1).

Using the definition of the slope of a function, at x  xi

f  xi   0
f xi  = tan θ =
xi  xi 1

which gives
f xi 
xi 1 = xi  (1)
f xi 

Equation (1) is called the Newton-Raphson formula for solving nonlinear equations of the form
f x   0 . So starting with an initial guess, xi , one can find the next guess, xi 1 , by using
Equation (1). One can repeat this process until one finds the root within a desirable tolerance.

Drawbacks of the Newton-Raphson Method

1. Divergence at inflection points

2. Division by zero
3. Oscillations near local maximum and minimum
4. Root jumping

DTU-MED Compiled BY Endalkachw A. |Numerical 8

Methods
NM Module 2016 E.C
2.3.2. Secant Method of Solving Nonlinear Equations

The Newton-Raphson method of solving a nonlinear equation f ( x)  0 is given by the iterative

formula

f ( xi )
xi 1 = xi  (1)
f ( xi )

One of the drawbacks of the Newton-Raphson method is that you have to evaluate the derivative
of the function. With availability of MATLAB etc, this process has become more convenient.
However, it still can be a laborious process, and even intractable if the function is derived as part
of a numerical scheme. To overcome these drawbacks, the derivative of the function, f (x) is
approximated as

f ( xi )  f ( xi 1 )
f ( xi )  (2)
xi  xi 1

Substituting Equation (2) in Equation (1) gives

f ( xi )( xi  xi 1 )
xi 1  xi  (3)
f ( xi )  f ( xi 1 )

The above equation is called the secant method. This method now requires two initial guesses,
but unlike the bisection method, the two initial guesses do not need to bracket the root of the
equation. The secant method is an open method and may or may not converge. However, when
secant method converges, it will typically converge faster than the bisection method. However,
since the derivative is approximated as given by Equation (2), it typically converges slower than
the Newton-Raphson method.

The secant method can also be derived from geometry,as shown in Figure 1. Taking two initial
guesses, xi 1 and xi , one draws a straight line between f ( xi ) and f ( xi 1 ) passing through the x -
axis at xi 1 . ABE and DCE are similar triangles.

Hence
AB DC

AE DE
f ( xi ) f ( xi 1 )

xi  xi 1 xi 1  xi 1

DTU-MED Compiled BY Endalkachw A. |Numerical 9

Methods
NM Module 2016 E.C

On rearranging, the secant method is given as

f ( xi )( xi  xi 1 )
xi 1  xi 
f ( xi )  f ( xi 1 )

Questions

Question 1: Derive the Bisection method, False position method, Secant method, Newton-
Raphson methods.
Question 2: Bisection method and Newton - Raphson methods advantages and drawbacks for
solving non linear equations.
Question 3: Find the roots of the equation x 3  0.165x 2  3.993 104  0
by using the following methods?
i) Bisection method ii) False position method
iii) Secant method iv) Newton - Raphson method
Question 4: Find the root of equation x + sinx -1 = 0 using Newton - Raphson method.
Question 5: Find the real root of equation ex – 4x2 = 0 using Bisection method.

References
1. M.K. Jain, S.R.K Iyengar and R.K Jain “ Numerical Methods for Scientific and
engineering computation” Sixth edition, New age international publishers. P.P: 18 - 50.
2. Chapra, S.C and Canale, R.P., Numerical methods for engineers, McGraw – Hill
international edition, 6th edition, 2010.

DTU-MED Compiled BY Endalkachw A. |Numerical 10

Methods
NM Module 2016 E.C

CHAPTER THREE

SYSTEM OF LINEAR EQUATIONS

The chapter gives formation of matrices from linear equations and discusses about solution of
different matrices. Various methods such as Gauss elimination, Gauss – Seidel method and LU
Decomposition methods are applied to solve engineering problems.

Learning objectives

After completing this course the students will be able to:

 Illustrate formation of matrices from system of linear equations and consistent and
Inconsistent systems (K, L2)
 Illustrate Gauss elimination method and apply to solve some problems (K, L3)
 Describe Gauss - Seidel method and apply to solve some problems (K, L3)
 Describe LU Decomposition method and apply to solve some problems (K, L3)

Introduction

3.1.1 Matrix algebra is used for solving systems of equations

Matrix algebra is used to solve a system of simultaneous linear equations. In fact, for many
mathematical procedures such as the solution to a set of nonlinear equations, interpolation,
integration, and differential equations, the solutions reduce to a set of simultaneous linear
equations
A general set of m linear equations and n unknowns,

a11x1  a12 x2    a1n xn  c1

a21x1  a22 x2    a2n xn  c2

……………………………………

…………………………………….

am1 x1  am2 x2  ........ amnxn  cm

can be rewritten in the matrix form as

DTU-MED Compiled BY Endalkachw A. |Numerical 11

Methods
NM Module 2016 E.C
 a11 a12 . . a1n   x1   c1 
a a 22 . . a 2 n  x  c 
 21  2  2 
        
     
       
a m1 am2 . . a mn   x n  c m 

Denoting the matrices by  A ,  X  , and C  , the system of equation is

AX   C , where  A is called the coefficient matrix, C  is called the right hand side vector and
X  is called the solution vector.
Sometimes AX   C  systems of equations are written in the augmented form. That is

a a12 ...... a1n  c1 

 11 
 a 21 a 22 ...... a 2 n  c 2 
A C      
 
   
a m1 am2 ...... a mn  c 
n

3.1.2 A system of equations can be consistent or inconsistent

A system of equations AX   C  is consistent if there is a solution, and it is inconsistent if
there is no solution. However, a consistent system of equations does not mean a unique solution,
that is, a consistent system of equations may have a unique solution or infinite solutions (Figure
1).

[A][X]= [B]

Consistent System Inconsistent System

Unique Solution Infinite Solutions

Figure 3.1. Consistent and inconsistent system of equations flow chart.

DTU-MED Compiled BY Endalkachw A. |Numerical 12

Methods
NM Module 2016 E.C

3.2.1) Gauss Elimination Method

One of the most popular techniques for solving simultaneous linear equations is the Gaussian
elimination method. The approach is designed to solve a general set of n equations and n
unknowns

a11x1  a12x2  a13x3  ...  a1n xn  b1

a21x1  a22x2  a23x3  ...  a2n xn  b2

. .
. .
. .

an1x1  an 2 x2  an3 x3  ...  ann xn  bn

Gaussian elimination consists of two steps

1. Forward Elimination of Unknowns: In this step, the unknown is eliminated in each
equation starting with the first equation. This way, the equations are reduced to one
equation and one unknown in each equation.
2. Back Substitution: In this step, starting from the last equation, each of the unknowns is
found.

Forward Elimination of Unknowns:

In the first step of forward elimination, the first unknown, x1 is eliminated from all rows below
the first row. The first equation is selected as the pivot equation to eliminate x1 . So, to eliminate
x1 in the second equation, one divides the first equation by a11 (hence called the pivot element)
and then multiplies it by a 21 . This is the same as multiplying the first equation by a 21 / a11 to
give

a 21 a a
a 21 x1  a12 x2  ...  21 a1n xn  21 b1
a11 a11 a11

Now, this equation can be subtracted from the second equation to give

 a   a  a
 a 22  21 a12  x 2  ...   a 2 n  21 a1n  x n  b2  21 b1
 a11   a11  a11

or
 x2  ...  a2 n xn  b2
a22

where

DTU-MED Compiled BY Endalkachw A. |Numerical 13

Methods
NM Module 2016 E.C
a21
  a22 
a22 a12
a11

a21
a2 n  a2 n  a1n
a11

This procedure of eliminating x1 , is now repeated for the third equation to the n th equation to
reduce the set of equations as

a11x1  a12 x2  a13 x3  ...  a1n xn  b1

 x2  a23
a22  x3  ...  a2 n xn  b2

 x2  a33
a32  x3  ...  a3 n xn  b3

. . .
. . .

an 2 x2  an 3 x3  ...  ann

 xn  bn

This is the end of the first step of forward elimination. Now for the second step of forward
elimination, we start with the second equation as the pivot equation and a22 as the pivot element.
 (the pivot
So, to eliminate x 2 in the third equation, one divides the second equation by a22
 . This is the same as multiplying the second equation by
element) and then multiply it by a32
 / a22
a32  and subtracting it from the third equation. This makes the coefficient of x 2 zero in the
third equation. The same procedure is now repeated for the fourth equation till the n th equation
to give

a11x1  a12 x2  a13 x3  ...  a1n xn  b1

 x2  a23
a22  x3  ...  a2 n xn  b2

 x3  ...  a3n xn  b3

a33

. . .
. . .

an3 x3  ...  ann

 xn  bn

The next steps of forward elimination are conducted by using the third equation as a pivot
equation and so on. That is, there will be a total of n  1 steps of forward elimination. At the
end of n  1 steps of forward elimination, we get a set of equations that look like

DTU-MED Compiled BY Endalkachw A. |Numerical 14

Methods
NM Module 2016 E.C

a11 x1  a12 x2  a13 x3  ...  a1n xn  b1

 x2  a23
a22  x3  ...  a2 n xn  b2

 x3  ... a3n xn  b3

a33

. . .
. . .
n 1
ann xn  bnn 1

Back Substitution:
Now the equations are solved starting from the last equation as it has only one unknown.

bn( n 1)
xn  ( n 1)
a nn

Then the second last equation, that is the (n  1) th equation, has two unknowns: x n and xn 1 , but
x n is already known. This reduces the (n  1) th equation also to one unknown. Back
substitution hence can be represented for all equations by the formula

bii 1   aiji 1 x j

n

j i 1
xi  for i  n 1, n  2,,1
aiii 1

and

bn( n 1)
xn  ( n 1)
a nn

3.2.2) Gauss-Seidel Method

In certain cases, such as when a system of equations is large, iterative methods of solving
equations are more advantageous. Elimination methods, such as Gaussian elimination, are prone
to large round-off errors for a large set of equations. Iterative methods, such as the Gauss-Seidel
method, give the user control of the round-off error. Also, if the physics of the problem are well
known, initial guesses needed in iterative methods can be made more judiciously leading to
faster convergence.
What is the algorithm for the Gauss-Seidel method? Given a general set of n equations and n
unknowns, we have

DTU-MED Compiled BY Endalkachw A. |Numerical 15

Methods
NM Module 2016 E.C
a11x1  a12 x2  a13 x3  ...  a1n xn  c1

a21x1  a22 x2  a23 x3  ...  a2n xn  c2

. . .
. . .

an1 x1  an 2 x2  an3 x3  ...  ann xn  cn

If the diagonal elements are non-zero, each equation is rewritten for the corresponding unknown,
that is, the first equation is rewritten with x1 on the left hand side, the second equation is
rewritten with x 2 on the left hand side and so on as follows.

c2  a21 x1  a23 x3   a2 n xn
x2 
a22


cn 1  an 1,1 x1  an 1, 2 x2   an 1,n  2 xn  2  an 1,n xn
xn 1 
an 1,n 1
cn  an1 x1  an 2 x2    an ,n 1 xn 1
xn 
ann

These equations can be rewritten in a summation form as

n
c1   a1 j x j
j 1
j 1
x1 
a11
n
c2   a2 j x j
j 1
j 2
x2 
a 22

.
.
n
c n 1  a
j 1
n 1, j xj
j  n 1
x n 1 
a n 1,n 1

DTU-MED Compiled BY Endalkachw A. |Numerical 16

Methods
NM Module 2016 E.C
n
cn   a nj x j
j 1
j n
xn 
a nn

Hence for any row i ,

n
ci   aij x j
j 1
j i
xi  , i  1,2,, n.
aii

Now to find xi ’s, one assumes an initial guess for the xi ’s and then uses the rewritten equations
to calculate the new estimates. Remember, one always uses the most recent estimates to
calculate the next estimates, xi . At the end of each iteration, one calculates the absolute relative
approximate error for each xi as

x inew  x iold
a  100
i
x inew

where xinew is the recently obtained value of xi , and xiold is the previous value of xi .

When the absolute relative approximate error for each xi is less than the pre-specified
tolerance, the iterations are stopped.

3.2.3 LU Decompostion Method

It is possible to show that any square matrix A can be expressed as a product of a lower
triangular matrix L and an upper triangular matrix U:
A= LU
The process of computing L and U for a given A is known as LU decomposition or LU
factorization. LU decomposition is not unique (the combinations of L and U for a prescribed A
are endless), unless certain constraints are placed on L or U.
After decomposing A, it is easy to solve the equations Ax =B. We first rewrite the equations as
LUx = B. Upon using the notation Ux = D, the equations become
LD = B
which can be solved for D by forward substitution. Then
Ux = D
will yield x by the back substitution process.
The advantage of LU decomposition over the Gauss elimination method is that once A is
decomposed, we can solve Ax = B for as many constant vectors B as we please. The cost of each
additional solution is relatively small, since the forward and back substitution operations are
much less time consuming than the decomposition process.

DTU-MED Compiled BY Endalkachw A. |Numerical 17

Methods
NM Module 2016 E.C

Questions

Question1: Write the formation of matrices from the system of linear eqautions?
Question2: Explain the step by step procedure involved Gauss elimination method and
Gauss – Seidel methods.
Question3: Solve the following system of linear equations by using Gauss elimination
Method and LU Decomposition methods
4x – 2y + z = 15, -3x –y + 4z = 8 and x – y + 3z = 13
Question4: Solve by Gauss – Siedel method, the following system of equations
20x + y -2z = 17, 3x + 20y – z = -18, 2x -3y + 20 z = 25

References:
1. M.K. Jain, S.R.K Iyengar and R.K Jain “Numerical Methods for Scientific and
engineering computation” Sixth edition, New age international publishers. P.P: 104 -
152.
2. Chapra, S.C and Canale, R.P., Numerical methods for engineers, McGraw – Hill
international edition, 6th edition, 2010.

DTU-MED Compiled BY Endalkachw A. |Numerical 18

Methods
NM Module 2016 E.C

CHAPTER FOUR

CURVE FITTING

This chapter deals with Lagrange and Newton’s divided difference interpolation methods to
solve curve fitting problems and least squares method to fit regression linear and polynomial
models.

Objectives

After completing this course the students will be able to:

 Apply Lagrange and Newton’s divided difference interpolation methods
to solve problems (K, L3)
 Apply least squares regression analysis to solve problems (K, L3)

4.1.1 What is interpolation?

Many times, data is given only at discrete points such as x0 , y0 , x1 , y1  , ......, xn 1, yn 1  ,
xn , yn  . So, how then does one find the value of y at any other value of x ? Well, a continuous
function f x  may be used to represent the n  1 data values with f x  passing through the
n  1 points (Figure 1). Then one can find the value of y at any other value of x . This is called
interpolation.
Polynomial interpolation involves finding a polynomial of order n that passes through
the n  1 points. One of the methods of interpolation is called the direct method. Other methods
include Newton’s divided difference polynomial method and the Lagrangian interpolation
method. We will discuss the direct method and Newton’s divided difference polynomial method
in this chapter.

4.1.1 Direct Method

The direct method of interpolation is based on the following premise. Given n  1 data points,
fit a polynomial of order n as given below

y  a0  a1 x  ............... a n x n (1)

through the data, where a0 , a1 ,........., an are n  1 real constants. Since n  1 values of y are
given at n  1 values of x , one can write n  1 equations. Then the n  1 constants,
a0 , a1 ,........., an can be found by solving the n  1 simultaneous linear equations. To find the
value of y at a given value of x , simply substitute the value of x in Equation 1.

DTU-MED Compiled BY Endalkachw A. |Numerical 19

Methods
NM Module 2016 E.C
But, it is not necessary to use all the data points. How does one then choose the order of
the polynomial and what data points to use? This concept and the direct method of interpolation
are best illustrated using examples.
4.2 Interpolation methods
a) Lagrange method

DTU-MED Compiled BY Endalkachw A. |Numerical 20

Methods
NM Module 2016 E.C

DTU-MED Compiled BY Endalkachw A. |Numerical 21

Methods
NM Module 2016 E.C

b) Newton’s divided difference method

To illustrate this method, linear and quadratic interpolation is presented first. Then, the general
form of Newton’s divided difference polynomial method is presented.
Linear Interpolation

Given ( x0 , y0 ) and ( x1 , y1 ), fit a linear interpolant through the data. Noting y  f (x) and
y1  f ( x1 ) , assume the linear interpolant f1 ( x) is given by (Figure 2)

f1 ( x)  b0  b1 ( x  x0 )

Since at x  x0 ,

f1 ( x0 )  f ( x0 )  b0  b1 ( x0  x0 )  b0

and at x  x1 ,

f1 ( x1 )  f ( x1 )  b0  b1 ( x1  x0 )

 f ( x0 )  b1 ( x1  x0 )

giving

f ( x1 )  f ( x0 )
b1 
x1  x0

So

b0  f ( x0 )

f ( x1 )  f ( x0 )
b1 
x1  x0

giving the linear interpolant as

DTU-MED Compiled BY Endalkachw A. |Numerical 22
Methods
NM Module 2016 E.C
f1 ( x)  b0  b1 ( x  x0 )

f ( x1 )  f ( x0 )
f 1 ( x)  f ( x0 )  ( x  x0 )
x1  x0

x1 , y1 

f1  x 

x0 , y0 
x

Figure 4.1 Linear interpolation.

Quadratic Interpolation
Given ( x0 , y0 ), ( x1 , y1 ), and ( x 2 , y 2 ), fit a quadratic interpolant through the data. Noting y  f (x),
y0  f ( x0 ), y1  f ( x1 ), and y 2  f ( x 2 ), assume the quadratic interpolant f 2 ( x) is given by

f 2 ( x)  b0  b1 ( x  x0 )  b2 ( x  x0 )( x  x1 )

At x  x0 ,

f 2 ( x0 )  f ( x0 )  b0  b1 ( x0  x0 )  b2 ( x0  x0 )( x0  x1 )

 b0

b0  f ( x0 )

At x  x1

f 2 ( x1 )  f ( x1 )  b0  b1 ( x1  x0 )  b2 ( x1  x0 )( x1  x1 )

DTU-MED Compiled BY Endalkachw A. |Numerical 23

Methods
NM Module 2016 E.C
f ( x1 )  f ( x0 )  b1 ( x1  x0 )

giving

f ( x1 )  f ( x0 )
b1 
x1  x0

At x  x2

f 2 ( x2 )  f ( x2 )  b0  b1 ( x2  x0 )  b2 ( x2  x0 )( x2  x1 )

f ( x1 )  f ( x0 )
f ( x 2 )  f ( x0 )  ( x2  x0 )  b2 ( x2  x0 )( x2  x1 )
x1  x0

Giving

f ( x 2 )  f ( x1 ) f ( x1 )  f ( x0 )

x 2  x1 x1  x0
b2 
x 2  x0

Hence the quadratic interpolant is given by

f 2 ( x)  b0  b1 ( x  x0 )  b2 ( x  x0 )( x  x1 )

f ( x2 )  f ( x1 ) f ( x1 )  f ( x0 )

f ( x1 )  f ( x0 ) x2  x1 x1  x0
 f ( x0 )  ( x  x0 )  ( x  x0 )( x  x1 )
x1  x0 x2  x0

x1 , y1 
 x2 , y 2 

f 2 x 

x0 , y0 
x

Figure 4.2 Quadratic interpolation

DTU-MED Compiled BY Endalkachw A. |Numerical 24

Methods
NM Module 2016 E.C

General Form of Newton’s Divided Difference Polynomial

In the two previous cases, we found linear and quadratic interpolants for Newton’s divided difference
method. Let us revisit the quadratic polynomial interpolant formula

f 2 ( x)  b0  b1 ( x  x0 )  b2 ( x  x0 )( x  x1 )

where

b0  f ( x0 )

f ( x1 )  f ( x0 )
b1 
x1  x0

f ( x 2 )  f ( x1 ) f ( x1 )  f ( x0 )

x 2  x1 x1  x0
b2 
x 2  x0

Note that b0 , b1 , and b2 are finite divided differences. b0 , b1 , and b2 are the first, second, and third
finite divided differences, respectively. We denote the first divided difference by

f [ x0 ]  f ( x0 )

the second divided difference by

f ( x1 )  f ( x0 )
f [ x1 , x0 ] 
x1  x0

and the third divided difference by

f [ x2 , x1 ]  f [ x1 , x0 ]
f [ x2 , x1 , x0 ] 
x 2  x0

f ( x 2 )  f ( x1 ) f ( x1 )  f ( x0 )

x 2  x1 x1  x0

x 2  x0

where f [ x0 ], f [ x1 , x0 ], and f [ x2 , x1 , x0 ] are called bracketed functions of their variables enclosed in

square brackets.

Rewriting,

f 2 ( x)  f [ x0 ]  f [ x1 , x0 ]( x  x0 )  f [ x2 , x1 , x0 ]( x  x0 )( x  x1 )

This leads us to writing the general form of the Newton’s divided difference polynomial for n  1 data
points, x0 , y0 , x1 , y1 ,......,xn1 , yn1 , xn , yn  , as

DTU-MED Compiled BY Endalkachw A. |Numerical 25

Methods
NM Module 2016 E.C
f n ( x)  b0  b1 ( x  x0 )  ....  bn ( x  x0 )( x  x1 )...(x  xn1 )

where

b0  f [ x0 ]

b1  f [ x1 , x0 ]

b2  f [ x2 , x1 , x0 ]

bn1  f [ xn1 , xn2 ,....,x0 ]

bn  f [ xn , xn1 ,....,x0 ]

where the definition of the mth divided difference is

bm  f [ xm ,........,x0 ]

f [ xm ,........,x1 ]  f [ xm 1 ,........,x0 ]

xm  x0

From the above definition, it can be seen that the divided differences are calculated recursively.

For an example of a third order polynomial, given ( x0 , y0 ), ( x1 , y1 ), ( x 2 , y 2 ), and ( x3 , y3 ),

f 3 ( x)  f [ x0 ]  f [ x1 , x0 ]( x  x0 )  f [ x 2 , x1 , x0 ]( x  x0 )( x  x1 )
 f [ x3 , x 2 , x1 , x0 ]( x  x0 )( x  x1 )( x  x 2 )

b0
b1
b2
x0 f x0 
f x1 , x0 
b3

x1 f x1  f x2 , x1 , x0 

f x2 , x1  f x3 , x2 , x1, x0 

x2 f x2  f x3 , x2 , x1 
f x3 , x2 

x3 f x3 

DTU-MED Compiled BY Endalkachw A. |Numerical 26

Methods
NM Module 2016 E.C
Figure 4.3 Table of divided differences for a cubic polynomial.

Least square regression

Where substantial error is associated with data, polynomial interpolation is inappropriate and
may yield unsatisfactory results when used to predict intermediate values. A more appropriate
strategy for such cases is to derive an approximating function that fit the shape or general trend
of the data without necessarily matching the individual points. According to this method one
has to determine the best approximating function which minimizes the sum of the square of the
errors, i.e,

Linear Regression
Linear regression is the most popular regression model. In this model, we wish to predict
response to n data points ( x1 , y1 ), ( x2 , y2 ),......,( xn , yn ) by a regression model given by

y  a0  a1 x (1)

where a 0 and a1 are the constants of the regression model.

A measure of goodness of fit, that is, how well a0  a1 x predicts the response variable y
is the magnitude of the residual  i at each of the n data points.

Ei  yi  (a0  a1 xi ) (2)

Ideally, if all the residuals  i are zero, one may have found an equation in which all the
points lie on the model. Thus, minimization of the residual is an objective of obtaining
regression coefficients.
The most popular method to minimize the residual is the least squares methods, where
the estimates of the constants of the models are chosen such that the sum of the squared residuals
n

E
2
is minimized, that is minimize i .
i 1

Why minimize the sum of the square of the residuals? Why not, for instance, minimize the
sum of the residual errors or the sum of the absolute values of the residuals? Alternatively,
constants of the model can be chosen such that the average residual is zero without making
individual residuals small. Will any of these criteria yield unbiased parameters with the smallest
variance? All of these questions will be answered below.
To explain this data by a straight line regression model,
y  a0  a1 x (3)

Let us use the least squares criterion where we minimize

DTU-MED Compiled BY Endalkachw A. |Numerical 27

Methods
NM Module 2016 E.C

n n 2

S r   Ei    yi  a0  a1 xi 
2
(9)
i 1 i 1

S r is called the sum of the square of the residuals.

To find a 0 and a1 , we minimize S r with respect to a 0 and a1 .

S r n
 2  yi  a0  a1 xi  1  0 (10)
a0 i 1

S r n
 2  yi  a0  a1 xi  xi   0 (11)
a1 i 1

giving
n n n
  yi   a0   a1 xi  0 (12)
i 1 i 1 i 1

n n n
  y i x i   a 0 x i   a1 x i2  0 (13)
i 1 i 1 i 1

n
Noting that a
i 1
0  a 0  a 0  . . .  a 0  na0

n n
na0  a1  xi  y i (14)
i 1 i 1

n n n
a 0  x i  a1  x i2   x i y i (15)
i 1 i 1 i 1

DTU-MED Compiled BY Endalkachw A. |Numerical 28

Methods
NM Module 2016 E.C

( xn , yn )
y
xi , yi 
Ei  yi  a0  a1 xi

 x2 , y 2 

x3 , y3 

y  a0  a1 x
x1 , y1 
x

Figure 4.4 Linear regression of y vs. x data showing residuals and square of residual at a
typical point, xi .

Solving the above Equations (14) and (15) gives

n n n
n x i y i  x i  y i
i 1 i 1 i 1
a1  2
(16)
 n 
n
n x   x i  2
i
i 1  i 1 
n n n n

 xi2  y i   xi  xi y i
i 1 i 1 i 1 i 1
a0  2
(17)
n
  n
n x i2   x i 
i 1  i 1 

Redefining

_ x
i 1
i
x (20)
n

DTU-MED Compiled BY Endalkachw A. |Numerical 29

Methods
NM Module 2016 E.C
n

_ y
i 1
i
y (21)
n
we can rewrite

_ _
a 0  y  a1 x

Questions

Question 1: Explain about Interpolation?

Question 2: Find a polynomial P(x) that interpolates the points {(-1,2), (0,3), (2,-1), (5,1)} using
the Lagrange method. Hence find out P(2.5).

Question 3:Derive linear and quadratic interpolation equations by using Newton’s divided
difference method?

Question 4: Using Newton’s divided differences, find the cubic polynomial for the following
data.

x: 0 1 2 3

y: 1 2 1 10

Question 5: Fit the curve of the form y = a + bx by the method of least squares from the
following data.
x 0 5 10 15 20 25
y 12 15 17 22 24 30

References:
1. M.K. Jain, S.R.K Iyengar and R.K Jain “Numerical Methods for Scientific and
engineering computation” Sixth edition, New age international publishers. P.P: 210 -
290.
2. Chapra, S.C and Canale, R.P., Numerical methods for engineers, McGraw – Hill
international edition, 6th edition, 2010.

DTU-MED Compiled BY Endalkachw A. |Numerical 30

Methods
NM Module 2016 E.C

CHAPTER FIVE
NUMERICAL DIFFERENTIATION

5.1) Introduction
In the case of differentiation, we first write the interpolating formula on the interval (x0,xn) and
differentiate the polynomial term by term to get an approximated polynomial to the derivative of
the function. When the tabular points are equidistant, we can use the finite divided differences
such as Central, backward and forward differences formula.

DTU-MED Compiled BY Endalkachw A. |Numerical 31

Methods
NM Module 2016 E.C

DTU-MED Compiled BY Endalkachw A. |Numerical 32

Methods
NM Module 2016 E.C

DTU-MED Compiled BY Endalkachw A. |Numerical 33

Methods
NM Module 2016 E.C

Questions

Question 1: Given the evenly spaced data points find first and second derivatives?

X: 0 0.1 0.2 0.3 0.4

f(x): 0.00 0.0819 0.1341 0.1646 0.1797

References

1. M.K. Jain, S.R.K Iyengar and R.K Jain “Numerical Methods for Scientific and
engineering computation” Sixth edition, New age international publishers. P.P: 348 -
370.
2. Chapra, S.C and Canale, R.P., Numerical methods for engineers, McGraw – Hill
international edition, 6th edition, 2010.

DTU-MED Compiled BY Endalkachw A. |Numerical 34

Methods
NM Module 2016 E.C

CHAPTER SIX

NUMERICAL INTEGRATION

The chapter deals with Trapezoidal rule, Simpson’s 1/3rd & 3/8th rules and Gauss quadrature
rules - single and multi segmental rules to solve numerical integration problems.

Objectives
After completing this course the students will be able to:

 Apply Trapezoidal rule of integration to solve problems (K, L3)

 Apply Simpson’s 1/3rd & 3/8th rules of integration to solve problems (K, L3)
 Apply Gauss Quadrature rule of integration to solve problems (K, L3)

DTU-MED Compiled BY Endalkachw A. |Numerical 35

Methods
NM Module 2016 E.C

DTU-MED Compiled BY Endalkachw A. |Numerical 36

Methods
NM Module 2016 E.C

DTU-MED Compiled BY Endalkachw A. |Numerical 37

Methods
NM Module 2016 E.C

DTU-MED Compiled BY Endalkachw A. |Numerical 38

Methods
NM Module 2016 E.C

DTU-MED Compiled BY Endalkachw A. |Numerical 39

Methods
NM Module 2016 E.C

DTU-MED Compiled BY Endalkachw A. |Numerical 40

Methods
NM Module 2016 E.C

DTU-MED Compiled BY Endalkachw A. |Numerical 41

Methods
NM Module 2016 E.C

Simpson’s 3/8 Rule

The third order Lagrange polynomial can be fit to four points and integrated,

to yield

where h=(b-a)/3. This equation is called Simpson’s 3/8 rule since h is multiplied by 3/8.

Questions

Question 1: Derive Single segment and multi segment Trapezoidal rule?

Question 2: Derive Single segment and multi segment Simpson’s 1/3rd rule?

Question 3: Evaluate ∫ dx/(1 + x2) with Trapezoidal rule from limits 0 to 1 with h = 1/6

Question 4: Evaluate ∫ dx/x with Simpson’s 1/3 rd rule from limits 0 to 1 with h =1/4.

Question 5: Evaluate ∫ eSinx with Simpson’s 3/8 th rule from limits 0 to π/2.

References

3. M.K. Jain, S.R.K Iyengar and R.K Jain “Numerical Methods for Scientific and
engineering computation” Sixth edition, New age international publishers. P.P: 348 -
370.
4. Chapra, S.C and Canale, R.P., Numerical methods for engineers, McGraw – Hill
international edition, 6th edition, 2010.

DTU-MED Compiled BY Endalkachw A. |Numerical 42

Methods
NM Module 2016 E.C

CHAPTER SEVEN
NUMERICAL SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS

These chapter deals with Euler’s method, Modified Euler’s method, Runge – Kutta 4th order
method and Finite difference method to solve ordinary differential equations.

Learning objectives

After completing this chapter the students will be able to:

 Apply Euler’s method to solve problems (K, L3)
 Apply Modified Euler’s method to solve problems (K, L3)
 Apply Runge – Kutta 4th order method to solve problems (K, L3)
 Apply Finite difference method to solve problems (K, L3)

DTU-MED Compiled BY Endalkachw A. |Numerical 43

Methods
NM Module 2016 E.C

DTU-MED Compiled BY Endalkachw A. |Numerical 44

Methods
NM Module 2016 E.C

DTU-MED Compiled BY Endalkachw A. |Numerical 45

Methods
NM Module 2016 E.C

DTU-MED Compiled BY Endalkachw A. |Numerical 46

Methods
NM Module 2016 E.C

DTU-MED Compiled BY Endalkachw A. |Numerical 47

Methods
NM Module 2016 E.C

DTU-MED Compiled BY Endalkachw A. |Numerical 48

Methods
NM Module 2016 E.C

DTU-MED Compiled BY Endalkachw A. |Numerical 49

Methods
NM Module 2016 E.C

Questions

Question1: Derive Euler’s method?

Question2: Derive modified Euler’s method?
Question3: Derive Runge – Kutta 4th order method?
Question4: Solve y1 + 2x + y = 0 , y(0) = -1 by using Euler method find y(0.4),
choose h = 0.1
Question5: Solve (1 + x) y1 + y2 = 0, y(0) = 1 by using modified Euler method by choosing
h = 0.1, find y (0.1) and y(0.2).
Question6: Find y(0.1), y(0.2) and y(0.3) from y1 = x + y2, y(0) = 1 by using Runge – Kutta
method of fourth order.

References

i. M.K. Jain, S.R.K Iyengar and R.K Jain “Numerical Methods for Scientific and
engineering computation” Sixth edition, New age international publishers. P.P:
425 - 585.
ii. Chapra, S.C and Canale, R.P., Numerical methods for engineers, McGraw – Hill
international edition, 6th edition, 2010.

*************** The End ******************

DTU-MED Compiled BY Endalkachw A. |Numerical 50

Methods