BACHELOR OF COMPUTER

APPLICATION

BCA-PC(L)-122
Computer Oriented Numerical
Methods

Directorate of Distance Education

Guru Jambheshwar University of
Science& Technology
Hisar - 125001
1
 CONTENTS

1. Computer Arithmetic 1-20

1.0 Learning objective
1.1 Introduction
1.2 Representation of Floating PointNumbers
1.2.1 Fixed PointRepresentation
1.2.2Floating PointRepresentation
1.3 Arithmetic operations with Normalized Floating PointNumbers
1.3.1 Addition
1.3.2 Subtraction
1.3.3 Multiplication
1.3.4 Division
1.4 Error
1.4.1Measurement of errors
1.5 Check YourProgress
1.6 Summary
1.7 Keywords
1.8 Self-AssessmentTest
1.9 Answers to check yourprogress
1.10 References/ SuggestedReadings

2
2. Iterative Methods 21-44
2.0 Objective
2.1 Introduction
2.2 Bisection Method
2.3 Rate of Convergence
2.4 False Position or Regula Falsi Method
2.5 Order of Convergence of False Position or Regula Falsi Method
2.6 Newton Raphson Method
2.7 Convergence of Newton Raphson Method
2.8 Bairstow’s Method
2.9 Check your progress
2.10 Summary
2.11 Keywords
2.12 Self-Assessment Test
2.13 Answer to check your progress
2.14 References/ SuggestedReadings

3. Solution of Simultaneous Linear Equations & Ordinary 45-77

Differential Equations
3.0 Objective
3.1 Introduction
3.2 Gaussian Elimination Method
3.3 Partial and Full Pivoting (Ill – Conditions)
3.4 Gauss – Seidel Iteration Method
3.5 Solution of a Differential equation
3.6 Taylor’s Series Method
3.7 Euler’s Method
3.8 Runge – Kutta Method
3.8.1 First Order Runge – Kutta Method
3.8.2 Second Order Runge – Kutta Method
3.8.3 Fourth Order Runge – Kutta Method

3
 3.9 Predictor – Corrector Methods
3.9.1 Milne’s Method
3.9.2 Adams – Bashforth Method
3.10 Check Your Progress
3.11 Summary
3.12 Keywords
3.13 Answer to Check Your Progress
3.14 Self Assessment Questions
3.15 References/ SuggestedReadings

4. Interpolation-1(Newton) 78-105
4.0 Objective
4.1 Introduction
4.2 Errors in Polynomial Interpolation
4.3 Finite Differences
4.4 Detection of Error by use of Difference Tables
4.5 Differences of a Polynomial
4.6 Newton’s Forward Interpolation Formula
4.7 Newton’s Backward Interpolation Formula
4.8 Check your progress
4.9 Summary
4.10 Keywords
4.11 Self-Assessment Test
4.12 Answer to check your progress
4.13 References/ SuggestedReadings

5. Interpolation-2(Lagranges) 106-124
5.0 Objective
5.1 Introduction
5.2 Central Difference Interpolation Formulae
5.2.1 Gauss’s Forward Interpolation Formula

4
 5.2.1 Gauss’s Backward Interpolation Formula
5.2.3 Stirling’s Formula
5.2.4 Bessel’s Formula
5.3 Interpolation with Unequal Intervals
5.3.1 Lagrange’s Interpolation Formula
5.4 Inverse Interpolation
5.4.1 Lagrange’s Method
5.4.2 Iterative Method
5.5 Check Your Progress
5.6 Summary
5.7 Keywords
5.8 Self Assessment Questions
5.9 Answer to Check Your Progress
5.10 References/ Suggested Readings
6. Numerical Differentiation & Integration 125-153
6.0 Objective
6.1 Introduction
6.2 Numerical Differentiation
6.2.1 Derivatives using Forward Difference Formula
6.2.2 Derivatives using Backward Differences Formula
6.2.3 Derivatives using Central Differences Formulae
6.3 Maxima and Minima of a Tabulated Function
6.4 Numerical Integration
6.4.1 Trapezoidal Rule
6.4.2 Simpson’s one-third Rule
6.4.3 Simpson’s three-eighth Rule
6.5 Errors in Quadrature Formulae
6.6 Check Your Progress
6.7 Summary
6.8 Keywords

5
6.9 Self Assessment Questions
6.10 Answer to Check Your Progress
6.11 References/ SuggestedReading

6
 SUBJECT: Computer Oriented Numerical Methods

COURSE CODE: BCA-PC(L)-122

AUTHOR: Dr. Pradeep K. Bhatia
LESSON NO. 1

Computer Arithmetic
REVISED/UPDATED SLM BY AMIT

STRUCTURE
1.0 Learning objective
1.1 Introduction
1.2 Representation of Floating PointNumbers
1.2.1 Fixed PointRepresentation
1.2.2Floating PointRepresentation

1.3 Arithmetic operations with Normalized Floating PointNumbers

1.3.1 Addition
1.3.2 Subtraction
1.3.3 Multiplication
1.3.4 Division
1.4 Error
1.4.1Measurement of errors
1.5 Check YourProgress
1.6 Summary
1.7 Keywords
1.8 Self-AssessmentTest
1.9 Answers to check yourprogress
1.10 References/ SuggestedReadings

1
 1.0 LEARNING OBJECTIVES

After reading this unit, you should be able to:

1. In computer arithmetic we learn how data is manipulating by using arithmetic operation
2. Data manipulation is learning to produce result of given problem
3. Learn about the error occur during arithmetic problems
4. Detail study of number representation in memory

1.1 INTRODUCTION

This chapter is to emphasis two types of real numbers viz. fixed point
representation,floatingpointrepresentation;withinfloatingpoint (non-normalizedand
normalized) and their representations in computer memory. Rules to perform
arithmetic operations (Addition, Subtraction, Multiplication, and Division) with
normalized floating numbers are also considered. At the last the various types of errors
with measurement that can be introduced during numerical computation are
alsodefined.

1.2 REPRESENTATION OF FLOATING POINT NUMBERS

For easier understanding we assume that computer can store and operate with decimal numbers,
although it does whole work with binary number system. Also only finite number of digits can
be stored in the memory of the computer. We will assume that a computer has a memory in
which each location can store digits with provision for sign (+ or -)

There are two methods for storing of real numbers in the memory of the computer:

1.2.1 Fixed PointRepresentation

1.2.2 Floating PointRepresentation

2
1.2.1 Fixed Point Representation

Memory location storing the number 412456.2465

+ 4 1 2 4 5 6 2 4

Assumed decimal position

Figure 1.1 Fixed point representations in Memory

This representation is called fixed point representation, since the position of decimal point is
fixed after 6 positions from left. In such a representation largest positive number we can store
999999.99 and smallest positive number we can store 000000.01. This range is quite inadequate.

Example 1.1: Following are the examples of fixed point representations in the decimal number
system:

210000

0.0005432

65754.546

234.00345

Example 1.2: Following are the examples of fixed point representations in the binary number
system:

10111

10.11101

111.00011

0.00011

Disadvantages of fixed point representation

3
Inadequate range: Range of numbers that can be represented is restricted by number of
digits or bits used.

1.2.2 Floating Point Representation

Floating point representation overcomes the above mentioned problem and in position to
accommodate a much wider range of numbers than fixed point representation.

In this representation a real number consists of two basic parts:

i) Mantissa part
ii) Exponent part
In such a representation it is possible to float a decimal point within number towards left or right
side.

For example:

53436.256 = 5343.6256 x 10 1

534.36256 x 102

53.436256 x 103

5.3436256 x 104

.53436256 x 105

.054436256 x 106

and so on

= 534362.56 x 10-1

5343625.6 x 10-2

53436256.0 x 10-3

534362560.0 x 10-4

And so on
4
 Floating Point Mantissa Exponent
Number

5343.6256 x 101 5343.6256 1

534.36256 x 102 534.36256 2

53.436256 x 103 53.436256 3

5.3436256 x 104 5.3436256 4

.53436256 x 105 0.53436256 5

Normalized
0.054436256 x 106 0.053436256 6 Floating
Point
……………. ………… …………
Number
534362.56 x 10-1 534362.56 -1

5343625.6 x 10-2 5343625.6 -2

53436256.0 x 10-3 53436256.0 -3

534362560.0 x 10-4 534362560.0 -4

…………. ………… ……

In general floating representation of a number of any base may be written as:

N = ±Mantissa x (Base) ±exponent

Representation of floating point number in computer memory (with four-digit mantissa)

Let us assume we have hypothetical 8-digit computer out of which four digits are used for
mantissa and two digits are used for exponent with a provision of sign of mantissa and sign of
exponent.

5
 Implied decimal point

± ±

Sign of Mantissa Sign of Exponent

Mantissa exponent

Figure 1.2 Floating point representation in memory (4-digit mantissa)

Normalized Floating Representation

It has been noted that a number may have more than one floating point representations. In order
to have unique representation of non-zero numbers a normalized floating point representation is
used.

A floating point representation in decimal number system is normalized floating point iff
mantissa is less than 1 and greater than equal to .1 or 1/10(base of decimal number system).

i.e.

.1≤ |mantissa| <1

A floating point representation in binary number system is normalized floating point iff mantissa
is less than 1 and greater than equal to .5 or 1/2(base of binary number system).

i.e.

.5≤ |mantissa| <1

In general, a floating point representation is called normalized floating point representation iff
mantissa lies in the range:

1/base≤ |mantissa| <1

Representation of normalized floating point number in computer memory with four-digit

mantissa:

6
 ± ±

Exponent
Sign of Implied Mantissa
Sign of
Mantissa Decimal point exponent

Figure 1.3 Normalized floating point representation in memory (4-digit mantissa)

Note: In computer, storage of floating point numbers is taken place in normalized form.

Disadvantages of floating point representation

 All the eight digits cannot be stored, since two digits are required by exponent.

Some specific rules are to be followed when arithmetic operations are performed with such
numbers.

1.3 ARITHMETIC OPERATIONS WITH NORMALIZED

FLOATING POINT
NUMBERS

1.3.1 Addition
For adding two normalized floating point numbers following rules are to be followed:

a) Their exponents are to be made same if they are not same.

b) Add their mantissa to get the mantissa of resultant.
c) Result is written in normalized floating point number.
d) Check the overflow condition.

7
Example 1.3

Add .4567E05 to .3456E05

Sol. Here exponents are equal, we have to add only mantissa and exponent remains
unchanged.

Addend .3456E05
Augend .4567E05
-----------
Sum . 8023E05

Example1.4

Add .3456E05 and .5433E07

Sol. Here exponents are not equal, therefore firstly make exponents same such that mantissa
of number with smaller exponent sifted towards R.H.S. equal to the number of digits smaller
exponent less than with larger exponent i.e. 7-5=2.
.3456E05  .0034E07 Addend
.5456E07 Augend
.5490E07 Sum

Example 1.5

Add .3456E03 and .5433E07

Sol. Here exponents are not equal, therefore firstly make exponents same such that mantissa
of number with smaller exponent sifted towards R.H.S. equal to the number of digits smaller
exponent less than with larger exponent i.e. 7-3=4.

.3456E03  .0000E07 Addend

.5456E07 Augend

.5456E07 Sum

8
Example 1.6

Add .4567E05 to .7456E05

Sol. Here exponents are equal, we have to add only mantissa and exponent remains
unchanged.

.3456E05 Addend

.7567E05 Augend

1.1023E05->.1102E06 Sum (Last digit of mantissa is chopped)

Example 1.7

Add .4567E99 to .7456E99

Sol. Here exponents are equal, we have to add only mantissa and exponent remains
unchanged.

.3456E05 Addend

.7567E05 Augend

1.1023E05->.1102E100 Sum (Last digit of mantissa is chopped)

OVERFLOW

As per exponent part can not store more than two digits, the number is larger than the
largest number that can be stored in a memory location. This condition is called
overflow condition and computer will intimate an error condition.

1.3.2 Subtraction
Rules to subtract a number from other are as follows:

a. Their exponents are to be made same if they are not same.

b. Subtract mantissa of one number from other to get the mantissa of resultant.
c. Result is written in normalized floating point number.
d. Check the underflow condition

9
Example 1.8

Subtract .3456E05 from .4567E05

Sol. Here exponents are equal, we have to subtract mantissa and exponents remain
unchanged.

.4567E05 Minuend

.3456E05 Subtrahend

.18114E05 Difference

Example 1.9

Subtract .3456E05 from .5433E07

Sol. Here exponents are not equal, therefore firstly make exponents same such that mantissa
of number with smaller exponent sifted towards R.H.S. equal to the number of digits smaller
exponent less than with larger exponent i.e. 7-5=2.

.5456E07

.3456E05  .0034E07

.5422E07

Example 1.10

Subtract .3456E03 from .5433E07

Sol. Here exponents are not equal, therefore firstly make exponents same such that mantissa
of number with smaller exponent sifted towards R.H.S. equal to the number of digits smaller
exponent less than with larger exponent i.e. 7-3=4.

.5433E07

.3456E03  .0000E07

10
.5433E07

Example 1.11

Subtract .5345E05 from .5444E05

Sol. Here exponents are equal, we have to subtract only mantissa and exponent remains
unchanged.

.5433E05

.5345E05

.0088E05->.8800E03

Example 1.12

Subtract .5345E-99 from .5433E-99

Sol. Here exponents are equal, we have to subtract only mantissa and exponent remains
unchanged.

.5433E-99

.5345E-99

.0088E-99 -> .8800E-101 (UNDERFLOW)

As per exponent part can not store more than two digits, the number is smaller than the
smallest number that can be stored in a memory location. This condition is called underflow
condition and computer will intimate an error condition.

1.3.3 Multiplication

If two normalized floating point numbers are to be multiplied following rules are followed:

a) Exponents are added to give exponent of the product.

b) Mantissas of two numbers are multiplied to give mantissa of the product.
c) Result is written in normalized form.
d) Check for overflow/underflow condition.
11
 (m1 x 10e1 ) x (m2 x10e2 ) = (m1xm2)x10(e1+e2)

Example 1.13 Find the product of following normalized floating point representation with 4
digit mantissa.

.4454E23 and .3456E-45

Sol.

Product of mantissa

.4454 x .3456 = .1539302

Discarded

Sum of exponents

23-45 = -18

Resultant Product is 0.1539E-18

Example 1.14 Find the product of following normalized floating point representation with 4-
digit mantissa.

.4454E23 and .1456E-45

Sol.

Product of mantissa

.4454 x .1456 = .0648502

Sum of exponent

23-45 = -18

Product is .0648502E-18 -> .648502E-19

Resultant product is 0.6485E-19

Example 1.15 Find the product of following normalized floating point representation with 4-
digit mantissa.
12
.4454E50 and .3456E51

Sol.

Product of mantissa

.4454 x .3456 = .1539302

Sum of exponent

50+51 = 101

Product is .1539E101 (OVERFLOW)

As per exponent part cannot store more than two digits, the number is larger than the largest
number that can be stored in a memory location. This condition is called overflow condition
and computer will intimate an error condition.

Example 1.16 Find the product of following normalized floating point representation with 4
digit mantissa.

.4454E-50 and .3456E-51

Sol.

Product of mantissa

.4454 x .3456 = .1539302

Discarded

Sum of exponent

-50-51 = -101

Product is .1539E-101 (UNDERFLOW)

As per exponent part can not store more than two digits, the number is smaller than the
smallest number that can be stored in a memory location. This condition is called underflow
condition and computer will intimate an error condition.

13
 1.3.4 Division

If two normalized floating point numbers are to be divided following rules are to be
followed:

a. Exponent of second number is subtracted from first number to obtain of the result.
b. Mantissas of first number is divided by second number to obtain mantissa of the result
c. Result is written in normalized form.
d. Check for overflow/underflow condition.
(m1 x 10e1 ) ÷ (m2 x10e2 ) = (m1÷m2)x10(e1- e2)

Example 1.17 Division of .8888E-05 by .2000 E -03

Sol. .8888E-05 ÷ .2000 E -03 = (.8888 ÷ .2222) E-2

= 4.4440E-2

= .4444E-1

1.4 ERRORS IN NUMBER REPRESENTATION

A computer has finite word length and so only a fixed number of digits are stored and used
during computation. This would mean that even in storing an exact decimal number in its
converted form in the computer memory, an error is introduced. This error is machine dependent.
After the computation is over, the result in the machine form is again converted to decimal form
understandable to the users and some more error may be introduced at this stage.

Input Error-1(e1) Error-2(e2) Output

Number Number+

e1+e2
Error introduced Error introduced
when data stored
when
informationretrieve
Figure 1.4 Effect of the errors on the result
d

14
 1.4.1 Measurement of errors

a) Error = True value – Approximate value= Etrue - Ecal

b) Absolute error = | Error| = E true - E ca l

E true  E cal
c) Relative error =
E true

E true  E cal
d) Percentage error = * 100
E true

Note: For numbers close to 1, absolute error and relative error are nearly equal.

For numbers not close to 1 there can be great difference.

Example: If X = 100500 Xcal = 100000

Absolute error = 100500 - 100000 =500

500
Relative error (Rx)= =0.005
10000

Example: If X = 1.0000 Xcal = 0.9898

Absolute error = 1.0000 – 0.9898 = 0.0102

0.0102
Relative error = = 0.0102
1

e) Inherent error

Error arises due to finite representation of numbers.

For example

15
1/3 = 0.333333 ……..

2 = 1.414………

22/7 = 3.141592653589793…………..

It is noticed that every arithmetic operation performed during computation, gives rise
to some error, which once generated may decay or grow in subsequent calculations.
In some cases, error may grow so large as to make the computed result totally
redundant and we call such a procedure numerically unstable. In some case it can be
avoided by changing the calculation procedure, which avoids subtractions of nearly
equal numbers or division by a small number or discarded remaining digits of
mantissa.

Example Compute midpoint of the numbers

A = 4.568 B= 6.762

Using the four-digit arithmetic.

Solution: Method I

AB 4.568  6.762

C= = = .5660x 10
2 2

Method II

B-A 6.762 - 4.568

C= A+ =4.568 + = .5665 x 10
2 2

f) Transaction Error

Transaction error arises due to representation of finite number of terms of an infinite series.

For example, finite representation of series Sin x, Log x, ex etc.

x3 x5 x7 x9
Sin x = x – + - + ………………
3 5 7 9

16
Apart from above type of errors, we face following two types of errors during computation,
we come across with large number of significant digits and it will be necessary to cut number
up to desired significant digits. Meanwhile two types of errors are introduced.

 Round-off Error
 Chopping-off Error
g) Round-off

Round-off a number to n significant digits, discard all digits to the right of the nth digit, and if
this discarded number is:

-less than half a unit in the nth place, leave the nth digit unaltered.

- greater than half a unit in the nth place, increase the nth place digit by unity.

- exactly half a unit in the nth place, increase the nth digit by unity if it is odd, otherwise leave
it is unchanged.

The number thus round-off said to be correct to n significant digits.

h) Chopping-off

In case of chopping-off a number to n significant digits, discard all digits to the right of the
nth digit, and leave the nth digit unaltered.

Note: Chopping-off introduced more error than round-off error.

Example: The numbers given below are rounded-off to five significant digits:

2.45678 to 2.4568

1.45334 to 1.4533

2.45657 to 2.4566

2.45656 to 2.4565

Example: The numbers given below are chopped-off to five significant digits:

2.45678 to 2.4567
17
 1.45334 to 1.4533

2.45657 to 2.4565

2.45656 to 2.4565

1.5 CHECK YOURPROGRESS

1. Convert 1234543into normalizes number.

2. Add .567E03 and .234E04.

3. Subtract .136E04 from 2.34E04.

4. Multiply .12E04 and .24E05.

5. Divide .23456E06 by .23456E02.

6. Calculate the relative error If X = 200500 Xcal = 200000

1.6 SUMMARY

Let’s take a quick review of all that we have learned about Computer Arithmetic.

1. We start from learning of representation of numbers in memory. We learn how to

represent fixed point numbers and floating point numbers. We learn about
disadvantages of fixed and floating point number system.

2. We learn how a number is converted into normalized number. Further we learn about
the mathematic operations (Addition, Subtraction, Multiplications, Division) on
normalized number.

3. We learn in this chapter about the error. We learn how an error is accrued when a
number is converted into normalized number.

18
 1.7 KEYWORDS

Fixed Point Representation: - in this type of rage of numbers that can be represented is
restricted by number of digits or bits used.

Floating Point Representation: - it overcomes the problem associate with fixed point
representation.It consiststwo-part mantissa and Exponent.

Normalized number: - Number whose mantissa part ranging between 1 and .1 are called
normalized number.

Error: -Difference Between actual value (experimental value) to calculated value is called error.

Over flow: - when exponent is greater than two digits is called over flow.

1.8 SELF ASSESSMENTTEST

1. Discuss the errors, if any, introduced by floating point representation of decimal

representation of decimal numbers in computers.
2. Describe the various components of computation errors introduced by the computer.
3. Assuming computer can handle 4 digit mantissa, calculate the absolute and relative errors
in the following operations were p=0.02455 and q=0.001756.
(a) p+q (b) p-q (c) px q (d) p  q
4. Assuming the computer can handle only 4 digits in the mantissa, write an algorithm to
add, subtract, multiply, and divide two numbers using normalized floating point
arithmetic

1.9 ANSWER TO CHECK YOURPROGRESS

1. .1234543E07
2. .2907E04

19
3. .22040E05
4. .0288E09 (.288E08)
5. 1E04 (.1E05)
6. .0025

1.10 REFERENCES/SUGGESTEDREADINGS
1. Computer Oriented Numerical Methods, V. Rajaraman, PHI.
2. Introductory Methods of Numerical Analysis, S.S. Sastry, PHI.
3. Numerical Methods for Scientific and Engineering Computation, M.K.Jain, S.R.Lyenger,
R.K.Jain, Wiley Eastern Limited.
4. Numerical Methods with Programs in C, T. Veerarajan, T. Ramachandarn., Tata McGraw
Hill.

20
 SUBJECT: Computer Oriented Numerical Methods

COURSE CODE: BCA-PC(L)-122

AUTHOR: Prof. Kuldip Bansal
LESSON NO. 2

ITERATIVE METHODS

REVISED/UPDATED SLM BY AMIT

STRUCTURE
2.0 Objective
2.1 Introduction
2.2 Bisection Method
2.3 Rate of Convergence
2.4 False Position or Regula Falsi Method
2.5 Order of Convergence of False Position or Regula Falsi Method
2.6 Newton Raphson Method
2.7 Convergence of Newton Raphson Method
2.8 Bairstow’s Method
2.9 Check your progress
2.10 Summary
2.11 Keywords
2.12 Self-Assessment Test
2.13 Answer to check your progress
2.14 References/ SuggestedReadings

21
 2.0 OBJECTIVE

The objective of this lesson is to develop Iterative methods for finding the roots of the algebraic
and the transcendental equations.

2.1 INTRODUCTION

To find roots of an equation f(x) = 0, up to a desired accuracy, iterative methods are very useful.

2.2 BISECTION METHOD

This method is due to Bolzano and is based on the repeated application of the intermediate value
property. Let the function f(x) be continuous between aand b. For definiteness, let f(a) be
negative &f(b) be positive. Then, the first approximation to the root is

x1 = ½ (a+b)

If f(x1) = 0, then x1 is a root of f(x) = 0. Otherwise, the root lies between a and x1 or x1 and b
according as f(x1) is positive or negative. Then, we bisect the interval as before and continue the
process until the root is found to desired accuracy. If f(x1) is +ve, so that the root lies between
aand x1. Then the second approximation to the root is x2 = ½(a+x1). If f(x2) is –ve, the root lies
1
between x1 and x2. Then the third approximation to the root is x3 = (x1+x2) and so
2
on.Graphically it can be shown as

22

Fig 2.1
At each step the interval is determined in which the root lies. Middle value is next
approximation. The process is carried out till the result upto desired accuracy is obtained.

2.3 RATE OF CONVERGENCE

Let x0, x1, x2, ……… be the values of a root () of an equation at the 0th, 1st, 2nd …….. iterations
while its actual value is 3.5567. The values of this root calculated by three different methods, are
as given below:

Root 1st Method 2nd Method 3rd Method

x0 5 5 5

x1 5.6 3.8527 3.8327

x2 6.4 3.5693 3.56834

x3 8.3 3.55798 3.55743

x4 9.7 3.55687 3.55672

x5 10.6 3.55676

x6 11.9 3.55671

The values in the 1st method do not converge towards the root 3.5567. In the 2nd and 3rd methods,
the values converge to the root after 6th and 4th iterations respectively. Clearly 3rd method
converges faster than the 2nd method. This fastness of convergence in any method is represented
by its rate of convergence.

If e be the error then ei =  - xi = xi+1 – xi.

If ei+1/eiis almost constant, convergence is said to be linear, i.e., slow.

If ei+1/eipis nearly constant, convergence is said to be of order p, i.e., faster.

23
Since in case of Bisection method, the new interval containing the root, is exactly half the length
of the previous one, the interval width is reduced by a factor of ½ at each step. At the end of the
nth step, the new interval will therefore, be of length (b-a)/2n. If on repeating this process n times,
the latest interval is as small as given , then (b – a) /2n

or n  [log(b-a) – log ]/log 2

This gives the number of iterations required for achieving an accuracy . For example, the
minimum number of iterations required for converging to a root in the interval (0, 1) for a given
 are as under:
: 10-2 10-3 10-4

n: 7 10 14

As the error decreases with each step by a factor of ½, (i.e. ex+1/ex = ½), the convergence in the
bisection method is ‘linear’.

Example: Find a root of the equation x3 – 4x – 9 = 0, using the bisection method correct to three
decimal places.

Solution: Let f(x) = x3 – 4x – 9

Since f(2) is –ve and f(3) is +ve, a root lies between 2 and 3.

 First approximation to the root is

x1 = ½(2+3) = 2.5

Thus f(x1) = (2.5)3 – 4(2.5) – 9 = -3.375, i.e., – ve.

 The root lies between x1 and 3. Thus the second approximation to the root is
x2 = ½ (x1+3) = 2.75.

Then f(x2) = (2.75)3 – 4(2.75) – 9 = 0.7969, i.e., + ve.

Therefore, the root lies between x1and x2. Thus the third approximation to the root is

x3 = ½(x1+x2) = 2.625

24
Then f(x3) = (2.625)3- 4(2.625)-9 = -1.4121, i.e, - ve

The root lies between x2 and x3. Thus the fourth approximation to the root is x4 = ½ (x2 + x3) =
2.6875.

Repeating this process, the successive approximations are

x5 = 2.71875, x6 = 2.70313, x7 = 2.71094

x8 = 2.70703, x9 = 2.70508, x10 = 2.70605

x11 = 2.70654,

Hence the root is 2.706.

Example: Find real positive root of the following equation by bisection method:

x3 – 7x +5 = 0

Solution:

Let f(x) = x3 – 7x +5,

f(0) = 5, f(1) = -1

 Root lies between 0 and 1,

ab
Values of a, b, and the signs + of functional values are shown as follows:
2

ab ab
a b f(a) f(b) f( )
2 2

0 1 0.5 + - +

0.5 1 0.75 + - +

0.75 1 0.875 + - -

0.75 0.875 0.8125 + - -

25
0.75 0.8125 0.7812 + - +

0.7812 0.8125 0.7968 + - -

0.7812 0.7968 0.7890 + - -

0.7812 0.7890 0.7851 + - -

0.7812 0.7851 0.7831 + - -

0.7812 0.7831 0.7822 + - +

2.7822 0.7831 0.7826 + - +

Root lies between 0.7826 and 0.7831

 Root is 0.783.

2.4 FALSE POSITION OR REGULA FALSI METHOD

Let f(x) = 0 be the equation to be solved and the graph of y =

f(x) be drawn. If the line joining the two points A  [xi-1, f(xi-
1)] and B  [xi, f(xi)] meets the x-axis at (xi+1, 0), xi+1 is
the approximate value of the root of the equation f(x) = 0.

The equation of the line joining the points A and B is

f ( xi )  f ( xi 1 ) Fig 2.2
y – f(xi –1) =  ( x  xi 1 )
xi  xi 1

Putting y = 0

f ( xi )  f ( xi 1 )
– f(xi –1) =  ( x  xi 1 )
xi  xi 1

26
 xi  xi 1
or x- xi –1 = –  f ( xi 1 )
f ( xi )  f ( xi 1 )

( x i  x i 1 ) f ( x i 1 )
 x = xi –1 –
f ( x i )  f ( x i 1 )
Therefore, the iterative formula is

( x i  x i 1 ) f ( x i 1 )
xi+1 = xi– 1 -
f ( x i )  f ( x i 1 )

x i 1 f ( x i )  x i f ( x i 1 )
or xi + 1 =
f ( x i )  f ( x i 1 )

xi 1 xi
=  [f(xi) – f(xi-1)]
f ( xi 1 ) f ( xi )

2.5 ORDER OF CONVERGENCE OF FALSE POSITION OR

REGULA FALSI METHOD

The iterative formula of Regula Falsi Method is

x i 1 f ( x i )  x i f ( x i 1 )
xi+1= (1)
f ( x i )  f ( x i 1 )

Let  be the root of the equation f(x) = 0 and ei-1, ei, ei+1 be the errors when xi-1, xi, xi+1are the
approximate values of the root .

 ei-1 = xi-1 -  or xi-1 = ei-1+ 

Similarly, xi= ei+, xi+1 = ei+1 +

Substituting the values of xi-1, xi, xi+1 in (1)

27
 (ei 1   ) f (ei   )  (ei   ) f (ei 1   )
ei+1 +  =
f (ei   )  f (ei 1   )

[(ei 1 f (ei   )  ei f (ei 1   )]   [ f (ei   )  f (ei 1   )]

=
f (ei   )  f (ei 1   )

ei 1 f (ei   )  ei f (ei 1   )
= +
f (ei   )  f (ei 1   )

ei 1 f (  ei )  ei f (  ei 1 )
 ei+1 = (2)
f (  ei )  f (  ei 1 )

Now expanding f(+ei) and f(+ei-1) by Taylor’s theorem.

Numerator = ei-1 f(+ei) – eif( +ei-1)

2
e
= ei-1[f() +eif()+ i f " ( )  ........ ]
2!

ei21
–ei[f() +ei-1f  ()+ f " ( )  ........ ]
2!

2
e e  e i e i21
= [ei-1 – ei] f() + i 1 i f " ( )  ........
2!

ei 1ei (ei  ei 1 )
~ f " ( )
2!

[ f() = 0, being the root of f(x) = 0 Terms containing ei3, e3i-

1andhigher degree terms are neglected.]

Denominator = f( +ei) – f(+ei-1)

2
ei
= [ f() +eif () + f " ( )  ........ ]
2!

ei21
– [f() +ei-1f()+ f " ( )  ........ ]
2!
28
 2
e i  e i21
= (ei– ei-1) f() + f " ( )  ........
2!

~(ei– ei-1) f’() [Terms containing ei2, e2i-1and

higher degree terms are neglected.]

 (2) becomes

1
ei 1 ei (ei  ei 1 ) f " (α )
ei+1 = 2!
(ei -ei-1 )f ' (α )

ei 1 ei f " (α )
or ei+1 =  ei 1 ei k , (3)
2 f ' (α )

f " (α )
where k =
2 f ' (α )

If p is the order of convergence

ei epi-1 k or taking ei = epi-1 k (4)

for all i  n, k is a constant.

Eliminating ei-1 from (3) and (4)

1/ p
 e  1
1 k
ei+1 = ei  i  k  ei p (5)
 k'  k '1 / p

Also ei+1 = eipk (6)

Equating the values of ei+1from (5) and (6)

1 k
= eipk
1
ei p (7)
k '1 / p

Choosing k and ksuch that k = k. k1/p= k1+1/p,(7) becomes

ei1+1/p = ei p

29
 1
 1+ =p or p2 – p – 1 = 0
p

1  5
 p=
2

5  1 3.236
Taking +ve sign, p = = = 1.618
2 2

 Order of convergence of Regular Falsi Method is 1.618

Example: Solve x3 – 9x +1 = 0 for the root lying between 2 and 4 by the method of false
position.

Solution: Let f(x) = x3 – 9x+1 = 0

 f(2) = – 9, f(4) = 29

In the iterative formula

( x i  x i 1 ) f ( x i 1 )
xi+1 = xi-1-
f ( x i )  f ( x i 1 )

x0  2, x1  4
Putting i = 1,
f ( x0 )  9 f ( x1 )  29

( x1  x 0 ) f ( x 0 )
x2 = x0 -
f ( x1 )  f ( x 0 )

(4  2)(9) 18
= 2–  2  2.47
29  (9) 38

For second approximation x3,

x1  2.47, x2  4
i = 2,
f ( x1 )  6.063 f ( x 2 )  29

1.53  6.063
x3 = 2.47 + = 2.73
35.063
30
For third approximation x4,

x2  2.73, x3  4
i = 3,
f ( x2 )  3.2 f ( x3 )  29

(1.27)(3.2)
x4 = 2.73 +  2.85
32.2

f(2.85) = –2.07

Putting i= 4, the fourth approximation is

(1.15)(2.07)
x5 = 2.85+  2.92,
31.07

f(2.92) = – 0.37

and for i = 5

(1.08)(3.7)
x6 = 2.92+  2.93
29.37

f(2.93) = 0.21

similarly, for

i=6

(1.07)(0.21)
x7 = 2.93+  2.937
29.21

 Root of f(x) = 0 is 2.94, correct to two significant figures.

31
 2.6 NEWTON-RAPHSON METHOD

Let f(x) = 0 be the equation whose solution is required. If xi be a point near the root, f(x) may be
written as

f(x) = (xi+ x  xi ).
Expanding it by Taylor’s series

( x  xi ) 2 
f(x) = f(xi) + (x-xi) f(xi) + f " ( x i )  ....  0
2!

As a first approximation, (x-xi)2 and higher degree terms are neglected.

 f(xi) + (x-xi) f(xi)= 0

f ( x i )
or x-xi= -
f ' ( xi )

f ( x i )
or x = xi -
f ' ( xi )

Iterative algorithm of Newton-Raphson method is

f ( x i )
xi+1 = xi - when f(xi) 0
f ' ( xi )

Geometrical interpretation of this formula may be given as follows:

Let the graph of y = f(x) be drawn and Pi be any point (xi , yi) on it.

32

Fig 2.3
Equation of the tangent at pi is

y – f(xi) = f(xi). (x – xi)

Putting y = 0, i.e., tangent at Pi meets the x -axis at Mi + 1whose abscissa is given by

f ( x i )
x – xi = –
f ' ( xi )

f ( x i )
or x = xi -
f ' ( xi )

which is nearer to the root .

 Iterative algorithm is

f ( x i )
xi+1=xi - .
f ' ( xi )

Thus in this method, we have replaced the part of the curve between the point Pi and x – axis by
a tangent to the curve at Pi and so on.

Example: Find the real root of the equation xex– 2 = 0 correct to two decimal places, using
Newton –Raphson method.

Solution :Given f(x) = xex –2, we have

f (x) = xex +ex and f(x) = xex + 2ex

Therefore, we obtain

f(0) = –2 and f(1) = e –2 = 0.71828

Hence, the required root lies in the interval(0,1)and is nearer to1. Also f (x) and f(x) do not
vanish in (0,1); f(x) and f(x) will have the same sign at x = 1. Therefore, we take the first
approximation x0= 1, and using Newton-Raphson method, we get

f ( x 0 ) e  2
x1= x0 - = = 0.867879
f ' ( x0 ) 2e

33
and f(x1) = 0.06716

The second approximation is

f ( x1 ) 0.06716
x2=x1 - = 0.867879 - = 0.85278
f ' ( x1 ) 4.44902

and

f(x2) = 7.655 x 10-4

Thus, the required root is 0.853

Example: Find a real root of the equation x3–x – 1 = 0 using Newton Raphson method, correct to
four decimal places.

Solution: Let f(x) = x3 – x – 1, then we note that f(1) = –1, f(2) = 5. Therefore, the root lies in the
interval (1, 2). We also note that

f (x) = 3x2 – 1,f(x) = 6x

and

f(1) = –1, f(1) = 6, f(2) = 5, f(2) = 12

Since f(2) and f(2) are of the same sign, we choose x0 = 2 as the first approximation to the root.
The second approximation is computed using Newton-Raphson method as

f ( x 0 ) 5
x1=x0 - =2- = 1.54545 and f(x1) = 1.14573
f ' ( x0 ) 11

The successive approximations are

1.14573
x2=1.54545 - = 1.35961, f(x2) = 0.15369
6.16525

0.15369
x3=1.35961 - = 1.32579, f(x3) = 4.60959 x 10-3
4.54562

34
 4.60959  10 3 
x4=1.32579 - = 1.32471, f(x4) = -3.39345 x 10-5
4.27316

3.39345  10 5
x5=1.32471+ =1.324718, f(x5) = 1.823 x 10-7
4.26457

Hence, the required root is 1.3247.

2.7 CONVERGENCE OF NEWTON-REPHSON METHOD

To examine the convergence of Newton-Raphson formula

f ( xn )
xn+1 = xn - (1)
f ' ( xn )

We compare it with the general iteration formula xn+1 = (xn), and thus obtain

f ( xn )
(xn) = xn -
f ' ( xn )

or, we write it as

f ( x)
(x) = x -
f ' ( x)

We also know that the iteration method converges if |(x)|<1. Therefore, Newton-Raphson
formula (1) converges, provided

|f(x) f(x)| <|f(x)|2 (2)

in the interval considered. Newton-Raphson formula therefore converges, provided the initial
approximation x0 is chosen sufficiently close to the root and f(x), f(x) and f(x) are continuous
and bounded in any small interval containing the root.

35
Definition:
Let
xn =  + εn, xn+1 =  + εn+1

where  is a root of f(x) = 0. If we can prove that εn+1 = K εp, where K is a constant and  n is the
error involved at the nth step, while finding the root by an iterative method, then the rate of
convergence of the method is p.

We can now establish that Newton-Raphosn method converges quadratically. Let

xn =  + εn, xn+1 =  + εn+1,

where  is a root of f(x) = 0 and εn is the error involved at the nth step, while finding the root by
Newton-Raphson formula (1). Then, Eq. (1) gives

f (   n )
 + εn+1 =  + εn – ,
f ' (   n )

i.e.,

f (   n )  f ' (   n )  f (   n )
εn+1 = εn – = n
f ' (   n ) f ' (   n )

Using Taylor’s expansion, we get

1
εn+1 =  n  f ' ( )   n f ' ' ( )  ....
f ' ( )   n f " ( )  .......

  n2  
–  f ( )   n f ' ( )  f " ( )  ..... 
 2  

Since  is a root, f() = 0. Therefore, the above expression simplifies to

 n2 1
εn+1 = f "( )
2 f ' ( )   n f "( )

36
 1
 n2 f " ( )  f " ( ) 
= 1 n
2 f ' ( )  f ' ( ) 

 n2 f " ( )  f " ( ) 
= 1 n
2 f ' ( )  f ' ( ) 

or

 n2 f " ( )
εn+1 =  o( n3 )
2 f ' ( )

On neglecting terms of order ( n3 ) and higher power, we obtain

εn+1 = K ( n2 ) , where (3)

f " ( )
K= (4)
2 f ' ( )

It shows that Newton-Raphson method has second order convergence or converges quadratically.

Example: Set up Newton’s scheme of iteration for finding the square root of a positive number
N.

Solution: The square root of N can be carried out as a root of the equation x2 – N = 0. Let
f(x) = x2 – N. By Newton’s method, we have

f ( xn )
xn+1 = xn - .
f ' ( xn )

In this problem, f(x) = x2 – N, f(x) = 2x. Therefore,

xn2  N 1  N
xn+1 = xn -   xn   (5)
2 xn 2 xn 

Example: Evaluate 12 , by Newton’s formula.

37
Solution : Since 9 = 3, 16 = 4, we take x0 = (3 + 4)/2 = 3.5. Using equation (5), we have

1 N 1 12 
x1 =  x0     3.5    3.4643
2 x0  2  3.5 

1 12 
x2 =  3 .4643    3 .4641
2 3 .4643 

1 12 
x3 =  3 .4641    3 .4641
2 3 .4641 

Hence, 12 = 3.4641upto four decimal places.

2.8 BAIRSTOW’S METHOD

Lin-Bairstow’s method is often useful in finding quadratic factors of polynomial and finding the
complete roots of a polynomial equation with real coefficients. Let the polynomial equation is
given by
f(x) = A3x3 + A2x2 + A1x + A0 = 0 (1)

Let x2 + Rx + S be a quadratic factor of f(x) and also let x2 + rx + s be an approximate factor.

Usually, first approximations to r and s can be obtained from the last three terms of the given
polynomial. Thus,

A1 A0
r and s (2)
A2 A2

Dividing f(x) by x2 + rx + s, let

f(x) = (x2 + rx +s) (B2x +B1) + Cx +D

= B2x3 +(B2r +B1)x2 + (C +B1r + sB2) x + (B1s+D), (3)

where the constants B1, B2, C and D have to be determined. Equating the coefficients of the like
powers of x in equations (1) and (3), we get

38
 B2 = A3

B1 = A2 - rB2

C = A1 - rB1 – sB2 (4)

D = A0 – sB1

From (4), it is clear that the coefficients B1, B2, C and D are functions of r and s. Since x2 + Rx +
S is a factor of the given polynomial, So
C(R, S) = 0 and D(R, S) = 0 (5)

Taking

R = r + r and S = s + s (6)

Equations (5) can be expanded by Taylor’s series and we obtain

C C
C(R, S) = C(r, s) + r .  s. 0
r s

(7)

D D
D(R, S) = D(r, s) + r .  s. 0
r s

where the derivatives are to be computed at r and s. We solve equation (7) for r and s. Using
these values in (6) will give the next approximation to R and S, respectively. This process can be
repeated until successive values of R and S show no significant change.

Example: Find the quadratic factor of the polynomial

f(x) = x3 – 2x2 + x – 2

Solution:

Here, we have A3 = 1, A2 = –2, A1 = 1 and A0 = –2

1
So, r=– and s=1
2
39
Using Equations (4), we have

B2 = 1; B1 = – 2 – r

C = 1 – r (– 2 –r ) – s = 1 + 2r + r2 – s,

and D = – 2– s (– 2 – r ) = – 2 + 2s +rs.
1 3
Also [C(r, s)](- ½,1) = 1 – 1 + –1=– ;
4 4

1 1
[D(r,s)](- ½,1) = –2 + 2 – =– ;
2 2

 C   C 
   2  2r  1;    1;
 r    1 ,1   s    1 ,1 
 2   2 

 D 
   s  1;
 r    1 ,1
 2 

 D  3
and   2r 
 s    ,1
1 2
 2 

Following equations (7), we get

r – s = ¾ and r +3/2s = ½

13 1
Hence r = and s = 
20 10

Thus, we have

1 13 3
R=– + = = 0.15
2 20 20

1 9
and S=1– = = 0.9
10 10

40
Therefore, the quadratic factor is x2 + 0.15x + 0.9, which is now taken as the approximate
quadratic factor to get the second approximation. So that, for the second approximation r = 0.15
and s = 0.9, then

C = 1 + 2.15 (0.15) – 0.9 = 0.4225,

D = –2 + 2.15 (0.9) = 0.065,

C
= 2 + 2(0.15) = 2.30,
r

C D D
= – 1; = 0.9 and = 2 + r = 2.15.
S r s

Hence, equations (7) give us

2.3 r – s = – 0.4225

and 0.9r + 2.15 s = – 0.065

Solving these equations, we get

r = – 0.1665312

and s = 0.0394783.

Therefore, the second approximations are obtained as

R = 0.15 – 0.16665312 = – 0.0165312

and S = 0.9 + 0.0394783 = 0.9394783

Thus, the second approximation to the quadratic factor is x2– 0.0165312 x + 0.9394783. This
process is repeated until no significant difference in the values of R and S is there, in two
consecutive steps

41
 2.9 CHECK YOUR PROGRESS

1. Find a root of the following equations by using bisection method.

(i) x3 - 2x - 5 = 0 (ii) x – cos x = 0

2. Using Regula Falsi method, compute the real root of the following equations.

(i) xex – 2 = 0 (ii) x3 - 4x - 9 = 0

3. Using Newton Raphson method, evaluate

(i) √32 (ii) 1/3 (iii) 1/15

4. Using Bairstow’s method, obtain the quadratic factor of the polynomial given by

f(x) = x3 – 2x2 + x – 2

2.10 SUMMARY

1. In this chapter we learn about how to find root of an equation. We learn to use different
method to solve equations or to find root of equations.
2. We learn graphical representation of deferent methods.
3. We learn to find rate of convergence of different method (Bisection method, Newton
raphson method, regulafalsi method)

2.11 KEYWORDS

1. Bisection Method:-x3=(x1+x2)/2, Where x1 is root of equation where function is

negative and x2 is root where function value is positive.
42
 f ( x i )
2. Newton Raphson Method:- x = xi -
f ' ( xi )

( x i  x i 1 ) f ( x i 1 )
3. Regula Falsi Method:- x = xi –1 –
f ( x i )  f ( x i 1 )
4. Bairstow’s Method: - Lin-Bairstow’s method is often useful in finding quadratic factors
of polynomial and finding the complete roots of a polynomial equation with real
coefficients.
5. Rate of Convergence: -It is defined how fast a method provide the accurate result of
function or equations.

2.11 SELF ASSESSMENT TEST

1. Find root of X2-20 using bisection method.

2. Find square root of 3 using newton raphson method.
3. Solve X2-12 using regulafasi method.

2.12 ANSWER TO CHECK YOUR PROGRESS

(1) (i) 2.687 (ii) 0.937

(2) (i) 0.853 (ii) 2.7065

(3) (i) 5.6569 (ii) 3.4482 (iii) 0.258

(4) x2 + 1

43
2.13 REFERENCES/ SUGGESTEDREADINGS

1. Computer Oriented Numerical Methods, V. Rajaraman, PHI.

2. Introductory Methods of Numerical Analysis, S.S. Sastry, PHI.
3. Numerical Methods for Scientific and Engineering Computation, M.K.Jain, S.R.Lyenger,
R.K.Jain, Wiley Eastern Limited.
4. Numerical Methods with Programs in C, T. Veerarajan, T. Ramachandarn., Tata McGraw
Hill.

44
 SUBJECT: Computer Oriented Numerical Methods

COURSE CODE: BCA-PC(L)-122

AUTHOR: Prof. KULDIP BANSAL
LESSON NO. 3

SOLUTION OF SIMULTANEOUS LINEAR EQUATIONS

AND ORDINARY DIFFERENTIAL EQUATIONS

REVISED/UPDATED SLM BY AMIT

STRUCTURE

3.0 Objective
3.1 Introduction
3.2 Gaussian Elimination Method
3.3 Partial and Full Pivoting (Ill – Conditions)
3.4 Gauss – Seidel Iteration Method
3.5 Solution of a Differential equation
3.6 Taylor’s Series Method
3.7 Euler’s Method
3.8 Runge – Kutta Method
3.8.1 First Order Runge – Kutta Method
3.8.2 Second Order Runge – Kutta Method
3.8.3 Fourth Order Runge – Kutta Method
3.9 Predictor – Corrector Methods
3.9.1 Milne’s Method
3.9.2 Adams – Bashforth Method
3.10 Check Your Progress
3.11 Summary
3.12 Keywords
3.13 Answer to Check Your Progress

45
3.14 Self Assessment Questions
3.15 References/ SuggestedReadings

46
 3.0 OBJECTIVE

The objective of this lesson is to describe the numerical methods for finding the solution
simultaneous linear equations and ordinary differential equations.

3.1 INTRODUCTION

The general form of a system of m linear equations in n unknowns x1,x2 , x3 , …, x n can be

represented in matrix form as under:

a11 a12 a13  a1n   x 1   b 1 

a    
 11 a 22 a 23  a 2 n   x 2   b2 
 (1)
          
    
a m1 a m 2 a m3  a mn   x n   bm 

Using matrix notation, the above system can be written in compact form as

A (X) = (B) (2)

The solution of the system of equations (2) gives nunknown values x1, x2, …., xn, which satisfy
the system simultaneously. If m>n, we may not be able to find a solution, in principle, which
satisfy all the equations. If m<n, the system usually will have an infinite number of solutions.
However, in this lesson, we shall restrict to the case m = n. In this case, if A 0, then the
system will have a unique solution, while, if A= 0, then there exists no solution.
Various numerical methods are available for finding the solution of the system of equations (2),
and they are classified as direct and iterative methods. In direct methods, we get the solution of
the system after performing all the steps involved in the procedure. The direct method, we
consider is Gaussian elimination.Under iterative methods, the initial approximate solution is

47
assumed to be known and is improved towards the exact solution in an iterative way. We
consider Gauss-Seidel iterative method.

3.2 GAUSSIAN ELIMINATION METHOD

In the Gaussian elimination method, the solution to the system of equations (2) is obtained in two
steps. In the first step, the given system of equations is reduced to an equivalent upper triangular
form using elementary transformations. In the second step, the upper triangular system is solved
using back substitution procedure by which we obtain the solution in the order xn, xn-1, xn-2, …, x2,
x1.
This method is explained by considering a system of n equations in n unknowns in the form as
follows
a11 x1  a12 x 2   a1n x n  b1 
a 21 x1  a 22 x 2   a 2 n x n  b 2 
 (3)
   
a n1 x1  a n 2 x 2   a nn x n  b n 

Step 1: We divide the first equation by a11 and then subtract this equation multiplied by a21,
a31,.....,an1 from 2nd, 3rd,....., nth equation. Then the system (3) reduces to the following form:
x1  a '12 x2    a '1n xn  b'1 
a '22 x2    a '2 n xn  b'2 
 (4)
   
a 'n 2 x2    a 'nn xn  b'n 

Here, we observe that the last (n – 1) equations are independent of x1, that is, x1 is eliminated
from the last (n–1) equations.
This procedure is repeated with the second equation of (4), that is, we divide the second equation
by a'22 and then x2 is eliminated from 3rd, 4th, …, nth equations of (4). The same procedure is
repeated again and again till the given system assumes the following upper triangular form:

48
 c 11 x1  c 12 x 2    c 1n x n  d 1 
c 22 x 2    c 2 n x n  d 2 

   (5)
C nn x n  d n 



Step II: Now, the values of the unknowns are determined by back substitution procedure, in
which we obtain xn from the last equation of (5) and then substituting this value of xnin the
preceding equation, we get the value of xn-1. Continuing this way, we can find the values of all
other unknowns in the order xn, xn-1, ….., x2,x1.
Example: Solve the following system of equations using Gaussian elimination method
2x + 3y – z = 5
4x + 4y – 3z =3
-2x + 3y – z =1
Solution: The given system of equations is solved in two stages.
Step 1: We divide the first equation by 2 and then subtract the resulting equation (multiplied by
4 and –2) from the second and third equations, respectively. Thus, we eliminate x from the 2nd
and 3rd equations. The resulting new system is given by
3 z 5 
x y-  
2 2 2

- 2y - z  - 7  (1)
6y - 2z  6 


Now, we divide the second equation of (1) by –2 and eliminate y from the last equation and the
modified system is given by
3 z 5
x y-  
2 2 2

z 7 
y    (2)
2 2 
- 5z  15


Step II : From the last equation of (2), we get
z =3 (3)
using this value of z, the second equation of (2) gives

49
 7 3
y  2 (4)
2 2
Using these values of y and z in the first equation of (2), we get
5 3
x  3 1
2 2
Hence, the solution is given by x = 1, y = 2, z = 3

3.3 PARTIAL AND FULL PIVOTING (ILL – CONDITIONS)

The Gaussian elimination method fails if any one of the pivot elements becomes zero. In such a
situation, we rewrite the equations in a different order to avoid zero pivots. Changing the order
equations is called pivoting.
Now we introduce the concept of partial pivoting. In this technique, if the pivot aii happens to be
zero, then the ith column elements are searched for the numerically largest element. Let the jth
row (j >i) contains this element, then we interchange the ith equation with the jth equation and
proceed for elimination. This process is continued whenever pivots become zero during
elimination. For example, let us examine the solution of the following simple system:

10-6x1 + x2 = 1
x1 + x2 = 2
Using Gaussian elimination method with and without partial pivoting, assuming that the
solution is required accurate to only four decimal places. The solution by Gaussian
elimination method gives x1 = 0, x2 = 1. If we use partial pivoting, the system takes the
form
x1 + x2 = 2
10-6x1 + x2 = 1
Using Gaussian elimination method, the solution is found to be x1 = 1, x2 = 1, which is a
meaningful and accurate result.
In full pivoting which is also known as complete pivoting, we interchange rows as well as
columns, such that the largest element in the matrix of the system becomes the pivot element. In

50
this process, the position of the unknown variables also get changed. Full pivoting, in fact, is
more complicated than the partial pivoting. Partial pivoting is preferred for hand computation.
Example: Solve the system of equations
x+ y+z=7
3x + 3y+ 4z =24
2x + y + 3z =16
by Gaussion elimination method with partial pivoting.
Solution: In matrix notation, the given system can be written as

1 1 1   x   7 
3 3 4   y    24  (1)
    
 2 1 3  z  16 

To start with, we observe that the pivot element a11= 1 ( 0). However, looking at the first
column, it shows that the numerically largest element is 3 which is in the second row. Hence, we
interchange the first row with the second row and then proceed for elimination. Thus, equation
(1) takes the form

3 3 4   x   24 
1 1 1  y    7  (2)
    
   
2 1 3   z  16 

after partial pivoting.

Step I: Dividing the first row of the system (2) by 3 and then subtracting the resulting row,
multiplied by 1 and 2 from the second and third rows of the system (2), we get

 4  x   8 
1 1
 3    
     
 1    
- y  -1 (3)
 3    
     
 -1 1   z   0 
 3     

The second row in equation (3) cannot be used as the pivot row, as a22 = 0. Interchanging the
second and third rows, we obtain

51
  4  x  8 
1 1 3     
     
 -1 1  y  0  (4)
 3    
     
1
 -   z    1
 3     

which is in the upper triangular form. From the last row of (4), we get

z=3 (5)

The second row of (4) with this value of z gives

-y + 1 = 0 or y = 1 (6)

Using these values of y and z, the first row of (4) gives

x + 1+ 4 = 8 or x =3 (7)

Thus, the required solution is

x = 3, y = 1, z=3

Example: Solve by Gaussian elimination method with partial pivoting, the following system
of equations:
0x1+ 4x2 +2x3 + 8x4 = 24

4x1+ 10x2 +5x3 + 4x4 = 32

4x1+ 5x2 +6.5x3 + 2x4 = 26

9x1+ 4x2 +4x3 + 0x4 = 21

Solution: The given system can be written as

0 4 2 8   x1   24 
4    
 10 5 4   x 2   32 
 (1)
4 5 6.5 2   x 3   26 
    
9 4 4 0  x 4   21 

To start with, we note that the pivot row, that is, the first row has a zero pivot element
(a11 = 0). This row should be interchanged with any row following it, which on becoming a pivot
52
row should not have a zero pivot element. While interchanging rows, it is better to interchange
with a row having largest pivotal element. Thus, we interchange the first and fourth rows, which
is called partial pivoting and we get

9 4 4 0  x1   21 
4    
 10 5 4   x 2   32 
 (2)
4 5 6.5 2   x 3   26 
    
0 4 2 8   x 4   24 

We note that, in partial pivoting, the unknown vector remains unaltered, while the right-hand
side vector gets changed.

Now, carry out Gaussian elimination process in two steps.

Step I: In this case, divide the first row of the system (2) by 9 and then subtracting this resulting
row multiplied by 4 from the second and third rows of equation (2), we get

 4 4   x1   2.333 
1 0   
9 9    
 
0 8.222 3.222 4   x2    22.667  (3)
   
0 3.222 4.722 2  x3  16.667 
 
0 4 2 8   x4   24 


Now we divide the second pivot row by 8.222 and subtract the resultant row multiplied by 3.222
and 4 from the third and fourth rows of equation (3), we obtain

 4 4   x1   2.333 
   
1 9 9
0     
   x    2.757 
0 1 0.392 0.486 
 2  
(4)
0 0 3.459 0.432   x 3   7.784 
     
0 0 0.432 6.054   x 4  12.973 

Finally, we divide the third pivot row by 3.459 and subtract the resultant row multiplied by 0.432
from fourth row of equation (4), thereby getting the triangular form

53
  4 4   x1   2.333 
   
1 9 9
0     
   x    2.757 
0 1 0.392 0.486 
 2  
(5)
0 0 1 0.125   x 3   2.250 
     
0 0 0 5.999   x 4  11.999 

Step II: From the last row of equation (5), we get x4= 2.00. Using this value of x4 into the third
row of equation (5), we obtain

x3 + 0.25 = 2.25 or x3 = 2.00 (6)

Similarly, we get

x2 = 1.00, x1 = 1.00

Thus, the solution of the given system is given by

x1 = 1.0, x2 = 1.0, x3 = 2.0, x4 = 2.0

3.4 GAUSS-SEIDEL ITERATION METHOD

It is a well-known iterative method for solving a system of linear equations.

a11x1 + a12x2+........+ a1nxn= b1,

a21x1 + a22 x2+........+ a2nxn= b2

.............................................................................

an1x1 + an2 x2+........+ annxn= bn

In Gauss-Seidel method, the corresponding elements of xi(r+1) replaces those of xi(r) as soon as
they become available. Hence, it is called the method of successive displacements. In this
method (r+1)th approximation or iteration is computed from

54
 b1 a12 (r) a 
x1(r 1)  - x 2   - 1n x (r)
n 
a11 a11 a11

(r 1) b 2 a21 (r 1) a2n (r) 
x2  - x1   - xn 
a22 a22 a22  (1)
............................................................. 

(r 1) b n an1 (r 1) an, (n -1) (r 1) 
xn  - x1   - x n -1 
ann ann ann 

Thus, the general procedure can be written in the following compact form

b i i -1 a ij (r  1) n a ij
x i(r  1)  - xj   x (r)j
a ii j 1 a ii j i 1 a ii (2)

for all i = 1,2,…….., n and r = 1,2,…

To describe system (1) in the first equation, we substitute the r-th approximation into the right-
( r 1)
hand side and denote the result by x1 . In the second equation, we substitute
( r 1) ( r 1)
( x1 , x 3( r ) ,...., x n( r ) ) and denote the result x2 . In the third equation, we substitute
( r 1)
( x1 , x 2( r 1) , x 4( r 1) ,...., x n( r ) ) and denote the result by x3(r 1) , and so on. This process is continued

till we arrive at the desired result.

Example: Find the solution of the following system of equations.

x1 – 0.25x2 –0.25x3 = 0.5

–0.25x1 +x2 – 0.25x4 = 0.5

–0.25x1+ x3 – 0.25x4 = 0.25

–0.25x2 – 0.25x3 + x4 = 0.25

using Gauss-Seidel method and perform the first four iterations.

Solution: The given system of equations can be rewritten as

x1 = 0.5 +0.25x2 +0.25x3

55
 x2 = 0.5 +0.25x1 +0.25x4 (1)

x3 = 0.25 +0.25x1 +0.25x4

x4 = 0.25 +0.25x2 +0.25x3

Taking x2 = x3 = x4 = 0 on the right- hand side of the first equation of system (1), we get x1(1) =
0.5 Taking x3 = x4 = 0 and the current value of x1, we get

x2(1) = 0.5 + (0.25) (0.5) +0= 0.625

form the second equation of system (1). Further, we take x4 = 0 and the current value of x1, we
obtain

x3(1) = 0.25 + (0.25) (0.5) +0 = 0.375

from the third equation of system (1). Now, using the current values of x2 and x3, the fourth
equation of system (1) gives

x4(1) = 0.25 + (0.25) (0.625) +(0.25) (0.375) = 0.5

The Gauss-Seidel iterations for the given set of equations can be written as

x1(r+1) = 0.5 +0.25x2(r) +0.25x3(r)

x2(r+1) = 0.5 +0.25x1(r+1) +0.25x4(r)

x3(r+1) = 0.25 +0.25x1(r+1) +0.25x4(r)

x4(r+1) = 0.25 +0.25x2(r+1)+0.25x3(r+1)

Now, by Gauss-Seidel procedure, the second and subsequent approximations can be obtained
and the sequence of the first four approximation are given as below:

Iteration Variables
Number (r) x1 x2 x3 x4
1 0.5 0.625 0.375 0.5
2 0.75 0.8125 0.5625 0.5938
3 0.8438 0.8594 0.60941 0.6172
4 0.8672 0.8711 0.6211 0.6231
56
 3.5 SOLUTION OF A DIFFERENTIAL EQUATION

The solution of an ordinary differential equation means finding an explicit expression for y in
terms of a finite number of elementary functions of x.
Let us consider the first order differential equation

dy/dx = f(x,y), given y(x0) = y0 , ---(1)

to study the various numerical methods of solving such equations. In most of these methods, we
replace the differential equation by a difference equation and then solve it. These methods yield
solutions either as a power series in x from which the values of y can be found by direct
substitution, or a set of values of x and y. The method of Taylor series belong to the former class
of solutions. In this method, y in (1) is approximated by a truncated series, each terms of which is
a function of x. The information about the curve at one point is utilized and the solution is not
iterated. As such, these are referred to as single-step methods. The methods of Euler, Runga-
Kutta, Milne, Adams-Bashforth etc. belong to the latter class of solutions. In these methods, the
next point on the curve is evaluated in short steps ahead, by performing iterations till sufficient
accuracy is achieved. As such, these methods are called step-by-step methods.

Euler and Runga-Kutta methods are used for computing y over a limited range of x-values
whereas Milne and Adams methods may be applied for finding y over a wider range of x-values.
Therefore Milne and Adams methods require starting values which are found by Taylor series or
Runga-Kutta methods.

Initial and boundary conditions:

An ordinary differential equations of the nth order is of the form

F(x, y, dy/dx, d2y/dx2, ….,dny / dxn) = 0 (2)

Its general solution contains n arbitrary constants and is of the form

(x,y, c1,c2, ……., cn) = 0 (3)

57
To obtain its particular solution, n conditions must be given so that the constants c1,c2, ……., cn
can be determined. If these conditions are prescribed at one point only (say x0) then the
differential equation together with the conditions constitute an initial value problem of the nth
order. If the conditions are prescribed at two or more points, then the problem is termed as
boundary value problem.

In this lesson, we shall first describe methods for solving initial value problems and then explain
finite difference method for solving boundary value problems.

3.6 TAYLOR’S SERIES METHOD

Consider the first order equation

dy
 f ( x, y ) (1)
dx

Differentiating (1) we have

d 2 y f f dy
  ,
dx 2 x y dx

i.e. , yn  fx  fy f ' (2)

Differentiating this successively, we can get y, yiv etc. Puttingx = x0 and y = 0, the values of
(y)0, (y)0, (y)0 can be obtained. Hence the Taylor’s series

( x  x0 ) 2 ( x  x0 ) 3
y = y0+(x-x0)(y)0+ (y)0 + (y)0 + .... (3)
2! 3!

gives the values of y for every value of x for which (3) converges.

On finding the value y1 for x = x1 from (3), y, y etc. can be evaluated at x = x1 by means of (1),
(2) etc. Then y can be expanded about x =x1. In this way, the solution can be extended beyond
the range of convergence of series (3).

58
Remarks: This is a single step method and works well so long as the successive derivatives can
be calculated easily. It is useful for finding starting values for the application of powerful
methods like Runga-Kutta, Milne and Adams-Bashforth.

Example: Find by Taylor’s series method, the values of y at x = 0.1 and x = 0.2 to five places
of decimals from dy/dx = x2y – 1, y(0) = 1.

Solution:

Here (y0) = 1.

 Differentiating successively and substituting, we get

y = x2y – 1, (y)0 = –1

y = 2xy + x2y, (y)0 = 0

y = 2y + 4xy + x2y, (y)0 = 2

yiv = 6y + 6xy + x2y, (yiv)0 = – 6, etc.

Putting these values in the Taylor’s series, we have

x2 x3 x4
y = 1+ x (-1) + (0)  (2)  (6)  .......
2! 3! 4!

x3 x4
 1 x    .......
3 4

Hence y (0.1) = 0.90033 and y(0.2) = 0.80227.

Example: Apply Taylor’s method to obtain approximate value of y at x = 0.2 for the
differential equation dy/dx = 2y + 3ex, y(0) = 0. Compare the numerical solution obtained with
the exact solution.

Solution:

We have y = 2y + 3ex ; y (0) = 2y (0) + 3e0 = 3.

59
Differentiating successively and substituting x =0, y = 0 we get

y = 2y+3ex , y(0) = 2y (0) +3 =9

y = 2y+3ex, y(0) = 2y(0) +3 =21 yiv = 2y+3ex,

yiv(0) = 2y (0) +3 =45 etc.

Putting these values in the Taylor’s series, we have

x2 x3 x 4 iv
y (x) = y(0) + xy (0) + y" (0)  y'" (0)  y (0)  .......
2! 3! 4!

9 2 21 3 45 4
= 0 +3x  x  x  x  .......
2 6 24

9 2 7 3 15 4
= 3x  x  x  x  .......
2 2 8

Hence y (0.2) = 3 (0.2) + 4.5(0.2)2 +3.5 (0.2)3 + 1.875 (0.4)4 +.....

= 0.8110 (1)

dy
Now  2 y  3e x is a Leibniz’s linear in x. Its I.F. being e-2x , the solution is
dx

ye-2x=  3ex.e-2xdx + c = - 3e-x + c

or y = -3ex + ce2x

Since y = 0 when x = 0, this implies c =3. Thus the exact solution is

y = 3 (e2x – ex)

when x = 0.2, y = 3 (e0.4 – e0.2) = 0.8112 (2)

Comparing (1) and (2) it is clear that (1) approximates to the exact value upto 3 decimal places.

60
 3.7 EULER’S METHOD

Consider the equation

dy
 f ( x, y ), y ( x0 )  y o . (1)
dx

The curve of solution through P(x0, y0) for this differential equation is shown in the following
figure. Further, we want to find the ordinate of any other point Q on this curve

Fig 7.1

Let us divide LM into n sub-intervals each of width h at L1, L2,...... so that h is quite small.

In the interval LL1, we approximate the curve by the tangent at P. If the ordinate through L1
meets this tangent in P1(x0+h,y1), then

y1 = L1P1=LP+R1P1 = y0 + PR1 tan

 dy 
= y0 + h   = y0 + h f(x0,y0)
 dx  P

Let P1Q1 be the curve of solution of (1) through P1 and let its tangent at P1 meet the ordinate
through L2 in P2(x0+2h,y2). Then

61
 y2=y1+h(fx0+h,y1) (2)

Repeating this process n times, we finally reach on an approximation MPn of MQ given by

yn = yn-1 + hf(x0+ n  1 h, yn-1)

or yn+1 = yn+ hf(xn,yn), where xn = x0 + nh.

This is Euler’s method of finding an approximate solution of (1).

Example: Using Euler’s method, find an approximate value of y corresponding to x = 1, given

that dy/dx = x + y and y= 1 when x = 0.
Solution:
We take n = 10 and h = 0.1 which is sufficiently small. The various calculations are arranged as
follows:

x y x + y = dy/dx old y + 0.1 (dy/dx) = new y

0.0 1.00 1.00 1.00+0.1(1.00)=1.10

0.1 1.10 1.20 1.10+0.1(1.20)=1.22

0.2 1.22 1.42 1.22+0.1(1.42)=1.36

0.3 1.36 1.66 1.36+0.1(1.66)=1.53

0.4 1.53 1.93 1.53+0.1(1.93)=1.72

0.5 1.72 2.22 1.72+0.1(2.22)=1.94

0.6 1.94 2.54 1.94+0.1(2.54)=2.19

0.7 2.19 2.89 2.19+0.1(2.89)=2.48

0.8 2.48 3.29 2.48+0.1(3.29)=2.81

0.9 2.81 3.71 2.81+0.1(3.71)=3.18

1.0 3.18

62
Thus the required approximate value of y = 3.18
dy y  x
Example: Given  with initial condition y = 1 at x = 0; find y for x = 0.1 by Euler’s
dx y  x
method.
Solution:
We divide the interval (0, 0.1) into five steps, i.e., we take n = 5 and h = 0.02. The various
calculations are arranged as follows:

x y dy/dx old y + 0.02 (dy/dx) = new y

0.00 1.0000 1.0000 1.0000+0.02(1.0000)=1.0200

0.02 1.0200 0.9615 1.0200+0.02(0.9615)=1.0392

0.04 1.0392 0.926 1.0392+0.02(0.926)=1.0577

0.06 1.0577 0.893 1.0577+0.02(0.893)=1.0756

0.08 1.0756 0.862 1.0756+0.02(0.862)=1.0928

0.10 1.0928

Hence the required approximate value of y = 1.0928.

3.8 RUNGE-KUTTA METHOD

The Taylor’s series method of solving differential equations numerically involved in finding the
higher order derivatives. However there is a class of methods known as RungeKutta methods
which do not require the calculations of higher order derivatives and give greater accuracy. The
Runge-Kutta formulae possess the advantage of requiring only the function values at some
selected points.

3.8.1 First order R-K method:

From Euler’s method, we have

63
 y1=y0+hf(x0, y0)= y0 +hy0 [y = f(x,y)]

Expanding by Taylor’s series

h2
y1=y(x0+h)= y0 +hy0+ y0 " +.......
2

It follows that the Euler’s method agrees with the Taylor’s series solution upto the term in h.

Hence, Euler’s method is the Runge-Kutta method of the first order.

3.8.2 Second order R-K method:

The modified Euler’s method gives

h
y1=y+ [f(x0, y0)+ f(x0+h,y1)] (1)
2

Substituting y1 = y0 + +hf(x0, y0) on the right hand side of (1), we obtain

h
y1 = y0+ [f0+f(x0+h, y0+ hf0)], (2)
2

where f0 = f(x0,y0).

Expanding L.H.S. by Taylor’s series, we get

h2 h3
y1=y(x0+h)=y0 + hy0 + y0+ y0+..... (3)
2! 3!

Expanding f(x0+h, y0+hf0) by Taylor’s series for a function of two variables, (2) gives

h   f   f   
y1= y0+  f 0   f ( x 0 , y 0 )  h    hf 0    O(h 2 ) 
2    x  0  y  0  

1  f   f   
= y0+ hf 0  hf 0  h 2      f 0   O(h 3 ) 
2   x 0  y 0  

64
 h2  df ( x, y ) f f 
= y0 + hf0 + f 0 'O(h3 )  dx  x  f y 
2  

h2
y1 = y0 + hy0 + y0 "O(h 3 ) (4)
2!

Comparing (3) and (4), it follows that this method agrees with the Taylor’s series solution upto
the term in h2 .

The Runge-Kutta method of the second order is same as modified Euler’s method.

 The second order Runge-Kutta formula is

1
y1 = y0 + ( k1  k 2 ) , where
2

ki = hf (x0,y0)

k2 = hf(x0 + h,y0 + k1).

Similarly, the third order Runge-Kutta formula is

1
y1 = y0 + ( k 1  4 k 2  k 3 ) , where
6

k1 = hf(x0, y0)

1 1
k2 = hf(x0 + h, y 0  k 1 )
2 2

and k3 = hf(x0+h, y0 +k), where

k= hf(x0+h, y0+k1).

3.8.3 Fourth order R-K method :

This method is most commonly used and working rule for finding the increment k of y
corresponding to an increment h of x by Runge-Kutta method for

65
 dy
 f ( x, y ), y ( x0 )  y0
dx

is as follows:

Calculate successively

k1 = hf(x0, y0)

1 1
k2 = hf(x0 + h, y 0  k 1 ) ,
2 2

1 1
k3 = hf(x0 + h, y 0  k 2 ) ,
2 2

k4= hf(x0 + h, y0 +k3)

1
Finally compute k= (k1 + 2k2+ 2k3+k4)
6

and y1 = y0 + k .

dy
Example: Given  1  y 2 , where y = 0 when x = 0, find y(0.2), y(0.4) and y(0.6) using
dx
Runge-Kutta method.

Solution:

we take h = 0.2 with x0 = y0 =0, we obtain

k1= 0.2,

k2 = 0.2 (1.01) = 0.202,

k3 = 0.2 (1+0.010201) = 0.20204,

k4 = 0.2 (1+0.040820) = 0.20816,

1
and y(0.2) = 0 + (k1+2k2+2k3+k4)
6

66
 = 0.2027, which is correct to four decimal places.

To compute y (0.4), we take x0 = 0.2 y0 = 0.2027 and h = 0.2. With these values, we get

k1= 0.2 [1+(0.2027)2] = 0.2082,

k2 = 0.2 [1+(0.3068)2] = 0.2188,

k3 = 0.2 [1+(0.3121)2] = 0.2195,

k4 = 0.2 [1+(0.4222)2] = 0.2356,

and y(0.4) = 0.2027+0.2201

= 0.4228, correct to four decimal places.

Finally, taking x0 = 0.4, y0 = 0.4228 and h = 0.2, and proceeding as above, we obtain y(0.6) =
0.6841.

3.9 PREDICTOR-CORRECTOR METHODS

In the methods described so far, to solve a differential equation over a single interval, say from x
= xn to x = xn+1, we require information only at the beginning of the interval, i.e., at x = xn.
Predictor-corrector methods are methods which require function values atxn, xn-1, xn-2,...... for the
computation of the function value at xn+1. A predictor formula is used to predict the value of y at
xn+1 and then a corrector formula is used to improve the value of yn+1.

Here we describe two such methods:

1. Milne’s method
2. Adams-Bashforth method
3.9.1Milne’s method:
Given dy/dx = f(x,y) and y = y0, x = x0; to find an approximate value of y for x = x0 + nh by
Milne’s method, we proceed as follows:
The value y0 = y(x0) being given, we compute

67
 y1 = y(x0+h), y2 = (x0+2h), y3 = y (x0+3h),

by Taylor’s series method.

Next we calculate

f0 = f(x0, y0), f1 = f(x0 + h, y1), f2 = f(x0 + 2h, y2), f3 = f(x0 + 3h, y3)

Then to find y4 = y(x0 +4h), we substitute Newton’s forward interpolation formula

n( n  1) 2 n(n  1)(n  2) 3
f(x, y) = f0 + nf0 +  f0+  f0+.......
2 6

in the relation

x0  4 h
y4 = y0 +  x0
f ( x, y )dx.

Therefore

x0  4 h  n(n  1) 2 
y4 = y0 +  x0
 f 0  nf 0 
 2
 f 0  ..........dx


[put x = x0+ nh, dx = hdn]

4  n(n  1) 2 
= y0 + h   f
0
0  nf 0 
2
 f 0  ..........dn


 20 8 
y4 = y0 +h  4 f 0  8f 0  2 f 0  3 f 0  ..........
 3 3 

Neglecting fourth and higher order differences and expressing f0, 2f0 and 3f0 in terms of the
function values, we get

4h
y4= y0 + ( 2 f1  f 2  2 f 3 )
3

which is called a ‘predictor’.

Having found y4, we obtain a first approximation to

68
 f4 = f(x0 + 4h, y4).

Then a better value of y4 is found by Simpson’s rule as

h
y4 = y2 + ( f 2  4 f3  f 4 )
3

which is called a ‘corrector’.

Then an improved value of f4 is computed and again the corrector is applied to find a still better
value of y4. We repeat this step until y4 remains unchanged.

Once y4 and f4 are obtained to desired degree of accuracy, y5 = y(x0 + 5h) is found from the
predictor as

4h
y5 = y1 + (2 f 2  4 f 3  2 f 4 )
3

and f5 = f(x0 + 5h, y5) is calculated. Then a better approximation to the value of y5 is obtained
from the corrector as

h
y5 = y3 + ( f3  4 f4  f5 ) .
3

We repeat this step till y5 becomes stationary and we, then proceed to calculate y6 as before.

This is Milne’s Predictor – Corrector method.

Example: Using Runge-Kutta method of order 4, find y for x = 0.1, 0.2, 0.3 given that dy/dx = xy
+ y2, y(0) = 1. Find the solution at x = 0.4 using Milne’s method.

Solution:

We have f(x,y) = xy +y2.

To find y(0.1), here x0 = 0, y0 = 1, h = 0.1.

 k1 = h f(x0,y0) = (0.1) f(0,1) = 0.1000

69
 1 1
k2 = h f(x0 + h ,y0+ k1 ) = (0.1) f(0.05,1.05) = 0.1155
2 2

1 1
k3 = h f(x0 + h ,y0+ k2 ) = (0.1) f(0.05,1.0577) = 0.1172
2 2

k4 = h f(x0+h ,y0+k3) = (0.1) f(0.1,1.1172) = 0.1360

1
k= (k1 + 2k2 +2k3 + k4)
6

1
= (0.1 + 0.231 + 0.2343 + 0.1360) = 0.1169
6

Thus y(0.1) = y1 = y0 +k = 1.1169,

To find y (0.2), here x1 = 0.1, y1 = 1.1169, h = 0.1.

k1= hf(x1, y1) = (0.1) f(0.1, 1.1169) = 0.1359

1 1
k2 = h f(x1+ h , y1+ k1 ) = (0.1) f(0.15,1.1848) = 0.1581
2 2

1 1
k3 = h f(x1+ h , y1+ k2 ) = (0.1) f(0.15,1.1959) = 0.1609
2 2

k4 = h f(x1+h, y1+k3) = (0.1) f(0.2,1.2778) = 0.1888

1
k= (k1 + 2k2 +2k3 + k4) = 0.1605
6

Thus y(0.2) = y2 = y1 +k = 1.2773.

To find y(0.3), here x2 = 0.2, y2 = 1.2773, h = 0.1.

k1= hf(x2, y2) = (0.1) f(0.2, 1.2773) = 0.1887

1 1
k2 = h f(x2+ h , y2+ k1 ) = (0.1) f(0.25,1.3716) = 0.2224
2 2

70
 1 1
k3 = h f(x2+ h , y2+ k2 ) = (0.1) f(0.25,1.3885) = 0.2275
2 2

k4 = h f(x2+h, y2+ k3) = (0.1) f(0.3,1.5048) = 0.2716

1
k= (k1 + 2k2 +2k3 + k4) =0.2267
6

Thus y(0.3) = y3 = y2 + k = 1.504.

Now the starting values for the Milne’s method are

x0= 0.0 y0= 1.0000 f0= 1.0000

x1= 0.1 y1= 1.1169 f1= 1.3591

x2= 0.2 y2= 1.2773 f2= 1.8869

x3= 0.3 y3= 1.5049 f3= 2.7132

Using the predictor,

4h
y4 = y0 + ( 2 f1  f 2  2 f 3 ) .
3

x4 = 0.4 y4 = 1.8344 f4 = 4.0988

and the corrector,

h
y4 = y2 + ( f2  4 f3  f4 )
3

0 .1
y4 = 1.2773 + [1.8869  4( 2.7132 )  4.098]
3

= 1.8386, f4 = 4.1159

Again using the corrector,

0 .1
y4 = 1.2773 + [1.8869  4( 2.7132)  4.1159])
3

71
 = 1.8391, f4 = 4.1182 (1)

Again using the corrector,

0 .1
y4 = 1.2773 + [1.8869  4( 2.7132 )  4.1182 ]
3

= 1.8392 which is same as(1)

Hence y (0.4) = 1.8392

3.9.2 Adams –Bashforth method:

dy
Given  f ( x, y ) and y0 = y(x0), we compute
dx

y-1= y(x0 – h), y-2 = y(x0 –2h), y-3 = y(x0 – 3h),

by Taylor’s series or Euler’s method or Rugne –Kutta method. Next we calculate

f-1=f(x0 – h, y-1), f-2 = f(x0 – 2h, y-2), f-3 =f(x0 – 3h, y-3).

Then to find y1, we substitute Newton’s backward interpolation formula

n ( n  1) 2 n ( n  1)( n  2 ) 3
f(x, y) =f0+ n  f o   f0   f 0  ......
2 6

x0 h
in y1 = y0 +  x0
f ( x, y )dx. (1)

x1  n(n  1) 2 
 y1 = y0 +  x0
 f 0  nf 0 
 2
 f 0  ........dx.


[put x = x0 +nh, dx = hdn]

1  n(n  1) 2 
= y0 +h  0
 f 0  nf 0 
 2
 f 0  ........dn


 1 5 3 
= y0 +h  f 0  f 0   2 f 0   2 f 0  ........
 2 12 8 

72
Neglecting fourth and higher order differences and expressing f0 ,  2 f 0 and 3 f 0 in terms of
function values, we get

h
y1= y0 + (55 f 0  59 f 1  37 f  2  9 f 3 ) (2)
24

This is called Adams-Bashforth predictor formula. Having found y1, we find f1 = f(x0+h1, y1).
Then to find a better value of y1, we derive a corrector formula by substituting Newton’s
backward formula at f1 ,i.e.,

n( n  1) 2 n( n  1)( n  2) 3
f(x,y) = f1+n f 1   f1   f 1  ...... in (1).
2 6

x1  n(n  1) 2 
 y1 = y0 +  x0
 f1  nf1 
 2
 f1  ........dx.


[put x = x1 +nh, dx = h dn]

0  n(n  1) 2 
= y0 +h  1
 f1  nf1 
 2
 f1  ........dn


 1 1 1 
or y1 = y0 +h  f1  f1   2 f1   2 f1  ........
 2 12 24 

Neglecting fourth and higher order differences andexpressing f1 ,  2 f1 and 3 f1 in terms of
function values, we obtain

h
y1= y0 + (9 f1  19 f 0  5 f 1  f  2 ) (3)
24

which is called Adams-Moulton corrector formula.

Then an improved value of f1 is calculated and again the corrector (3) is applied to find a still
better value y1. This step is repeated till y1 remains unchanged and then we proceed to calculate
y2 as above.

73
 dy
Example:Given  x 2 (1  y ) and y(1) = 1, y(1.1) = 1.233, y(1.2) = 1.548, y(1.3) = 1.979,
dx
evaluate y(1.4) by Adams-Bashforth method.

Solution:

Here f(x, y) = x2(1+y).

Starting values of the Adams-Bashforth method with h = 0.1 are

x = 1.0, y-3 = 1.000, f-3 = (1.0)2 (1+1.000) = 2.000

x = 1.1, y-2 = 1.233, f-2 = 2.702

x = 1.2, y-1 = 1.548, f-2 = 3.669

x = 1.3, y0 = 1.979, f0 = 5.035

Using the predictor,

h
y1= y0 + (55 f 0  59 f 1  37 f  2  9 f  3 )
24

x = 1.4, y1 = 2.573, f1 = 7.004

Using the corrector,

h
y1 = y0 + (9 f1  19 f 0  5 f 1  f  2 )
24

0 .1
y1 = 1.979 + (9  7.004  19  5.035  5  3.669  2.702 )
24

= 2.575

Hence y(1.4) = 2.575.

74
 3.10 CHECK YOUR PROGRESS

1. Solve the equations x + 4y – z = -5; x + y –6 z = -12;

3x - y – z = 4, by using Gauss elimination method.

2. Solve the following by Gauss – Seidel iteration method:

(i) 2x + y + 6z = 9; 8x + 3y + 22 = 13; x + 5y + z = 7

(ii) 5x + 2y + z = 12; x + 4y + 22 = 15; x + 2y + 5z = 20

3. Using Taylor’s series method, evaluate y(0.1) if y(x) satisfies

dy
 xy  1 , y(0) = 1
dx

4. Given
dy y  x
 , with y(0) = 1,
dx y  x

find y for x = 0.1

5. Apply Runge – Kutta method to find an approximate value of y for x = 0.2 in steps of
0.1, if
dy
 x  y 2 , where y(0) = 1.
dx

6. Apply Milne’s method to find a solution of the differential equation

dy
 x  y2 , in the range 0 < x < 1 for y = 0 at x = 0.
dx

75
 3.11 SUMMARY

Linear Equations: -In mathematics, a linear equation is an equation that may be put in the
form where are the variables, and are the coefficients, which are often real numbers. The
coefficients may be considered as parameters of the equation, and may be arbitrary
expressions, provided they do not contain any of the variables.

Ordinary Differential Equations:- A most general form of an ordinary differential equation

(ode) is given by

f( x, y, y', . . ., y(m) ) = 0

where x is the independent variable and y is a function of x. y', y'' . . . y(m) are
respectively, first, second and mth derivatives of y with respect to x.

Some definitions:

 the highest derivative in the ode is called the order of the ode.

 the degree of the highest order derivative in the ode is called the degree of the ode.

 the differential equation is linear if no product of the dependent variable y(x) with itself
or with any one of its derivatives occur in the equation otherwise is non-linear.

3.12 KEYWORDS

1. Gaussian Elimination Method:- The process of row reduction makes use of elementary row
operations, and can be divided into two parts. The first part (sometimes called forward
elimination) reduces a given system to row echelon form, from which one can tell whether
there are no solutions, a unique solution, or infinitely many solutions. The second part
(sometimes called back substitution) continues to use row operations until the solution is
found; in other words, it puts the matrix into reduced row echelon form.

76
2. Gauss – Seidel Iteration Method:- Gauss–Seidel method. In numerical linear algebra,
the Gauss–Seidel method, also known as the Liebmann method or the method of
successive displacement, is an iterative method used to solve a linear system of equations.
3. Euler’s Method:- yn+1 = yn+ hf(xn,yn), where xn = x0 + nh.
4. Runge – KuttaMethod:-RungeKutta methods which do not require the calculations of
higher order derivatives and give greater accuracy. The Runge-Kutta formulae possess the
advantage of requiring only the function values at some selected points.
5. Taylor’s Series Method:-
( x  x0 ) 2 ( x  x0 ) 3
y = y0+(x-x0)(y)0+ (y)0 + (y)0 + ....
2! 3!

3.13 SELF ASSESSMENT TEST

1. Solve the following system of equations using Gaussian elimination method

1x + 2y + 1z = 3
3x + 2y + 1z = 3
1x - 2y – 5z =1
2. Find the solution of the following system of equations.

3x - y + z = 1
3x + 6y + 2z = 0
3x + 3y + 7z = 4

using Gauss-Seidel method and perform the first four iterations.

3. Find y(0.5) for y’=-2x-y, y(0)=-1, with step length 0.1,Euler’s method.

4. Find y(0.2) for y’=-y, y(0)=1, with step length 0.1

5. Find y(0.3) for y’=(x.y2+y), y(0)=1, with step length 0.1, Using R-K fourth order method.

6. Find y(0.4) for y’=y-x3, y(0)=1, by Milen’s predictor corrector method.

77
 3.14 ANSWER TO CHECK YOUR PROGRESS

(1) x = 117/71, y = -81/71, z = 148/71

(2) (i) x = 1, y = 1, z = 1 (ii) x = 1, y = 2, z = 3

(3) 1.1053425

(4) 1.0928

(5) 1.2736

(6) 0.4555

3.15 REFERENCES/ SUGGESTED READINGS

1 Computer Oriented Numeical Methods, V. Rajaraman, PHI.

2 Introductory Methods of Numerical Analysis, S.S. Sastry, PHI.
3 Numerical Methods for Scientific and Engineering Computation, M.K.Jain, S.R.Lyenger,
R.K.Jain, Wiley Eastern Limited.
4 Numerical Methods with Programs in C, T. Veerarajan, T. Ramachandarn., Tata McGraw
Hill.

78
 SUBJECT: Computer Oriented Numerical Methods

COURSE CODE: BCA-PC(L)-122

AUTHOR: Prof. Dr. Aseem Miglani
LESSON NO. 4

INTERPOLION

REVISED/UPDATED SLM BY AMIT

STRUCTURE
4.0 Objective
4.1 Introduction
4.2 Errors in Polynomial Interpolation
4.3 Finite Differences
4.4 Detection of Error by use of Difference Tables
4.5 Differences of a Polynomial
4.6 Newton’s Forward Interpolation Formula
4.7 Newton’s Backward Interpolation Formula
4.8 Check your progress
4.9 Summary
4.10 Keywords
4.11 Self-Assessment Test

4.12 Answer to check your progress

4.13 References/ SuggestedReadings

78
 4.0 OBJECTIVE

In this lesson, the objective is to introduce the concepts of finite differences, thereby, to derive
the interpolation formulae using the forward and backward difference operators and tables, for
the given equi-spaced set of tabular values.

4.1 INTRODUCTION

The problem of approximating a given function f(x) by polynomials P9x) is generally used for
the construction of the function f(x), when it is not given in the form and only the values of f(x)
are given at a set of distinct points. The deviation of f(x) from P(x), i.e., f(x) – P(x), x є [a, b], is
called the error of approximation.

Let y = f(x), x0≤ x ≤ xn, be a function, then corresponding to every value of x in the range x0≤
x≤ xn, there exists one or more values of y. If the function f(x) is single-valued and continuous
and that it is known explicitly, then the values of f(x) corresponding to certain given values of x,
say x0, x1…. xn can easily be computed and tabulated. Conversely, given a set of tabulated values
(x0, y0), (x1, y1), (x2, y2),......, (xn, yn) satisfying the relation y = f(x) where the explicit nature of
f(x) is not known, then it is required to find a function, say (x), such that f(x) and (x) agree at
the set of tabulated points. Such a process is called interpolation. If (x) is a polynomial, then the
process is called polynomial interpolation and (x) is called the interpolating polynomial.
Similarly, different types of interpolation arise depending on whether (x) is a finite
trigonometric series, series of Bessel functions, etc. But, we shall be concerned with polynomial
interpolation only.

The study of interpolation is based on the calculus of finite differences, which deals with the
changes that take place in the value of the function (dependent variable) due to finite changes in
the independent variable.

79
In this lesson, we derive two important interpolation formulae by means of forward and
backward differences of a function.

4.2 ERRORS IN POLYNOMIAL INTERPOLATION

Let the function y = f(x) be defined by the (n+1) points (x i, yi), i = 0, 1, 2,….., n, and is
continuous and differentiable (n+1) times. Let y = f(x) be approximated by a polynomial n(x) of
degree not exceeding n such that

n(xi) = yi, i = 0, 1,2,……… n. ---(1)

Now, if we use n(x) to obtain approximate values of f(x) at some points other than those defined
by (1), what would be the accuracy of this approximation? Since the expression f(x) - n(x)
vanishes for x = x0, x1,….. xn, so we consider

f(x) - n(x) = Ln+1(x), ---(2)

where n+1(x) = (x – x0) (x – x1).....(x – xn), ---(3)

and L is to be determined such that equation (2) holds for any intermediate value of x, say x = x,
x0<x<xn. Clearly, we have
f ( x' )   n ( x' )
L= . ---(4)
 n 1 ( x' )
We construct a function F(x) such that

F(x) = f(x) - n(x) – L n+1(x), ---(5)

where L is given by (4) above. From the definition of F(x), it is clear that

F(x0) = F(x1) = … = F(xn) =F(x) = 0,

that is, F(x) vanishes at (n+2) points in the interval x0 ≤ x ≤ xn. Consequently, by the repeated
application of Rolle’s theorem F(x) must vanish (n+1) times, F(x) must vanish n times, etc., in

80
the interval x0 ≤x≤ xn. In particular, F(n+1)(x) must vanish once in the interval. Let this point be
given by x=, x0<< xn. Differentiating (5), (n +1) times with respect to x and put x = , we get

0 = f(n+1) ()- L. (n+1)! ---(6)

Expression (6) implies

f ( n 1) ( )
L= ---(7)
(n  1)!

Comparison of (4) and (7) gives us

f ( n 1) ( )
f(x) - n(x) = n+1 (x). ---(8)
(n  1)!

Since x = x is any intermediate value of x, so dropping the prime on x, we obtain

 n1( x)
f(x) - n(x) = f(n+1)(), x0<<xn, ---(9)
(n  1)!

which is the required expression for the error. It is extremely useful in theoretical work in
different branches of numerical analysis. In particular, we will use it to determine errors in
Newton’s interpolation formulae.

4.3 FINITE DIFFERENCES

Suppose that the function y = f(x) be tabulated for equally spaced set of values say xi = x0 + ih, i
= 0, 1,2,….n, and we have a table of values (xi, yi), i = 0, 1,2,3…..n. Finding the values of f(x)
for some intermediate values of x, or the derivative of f(x) for some x in the range x0≤ x ≤ xn, is
based on the concept of the differences of a function. The following three types of differences
are found useful.

81
Forward Differences:
If y0, y1, y2, …, yn denote a set of values of y, then y1 – y0, y2 – y1,… , yn – yn-1 are called the
differences of y. These differences of y denoted by y0, y1,…, yn-1 respectively are called the
first forward differences, and we have

y0 = y1 – y0, y1 = y2 – y1, ....., yn-1 = yn – yn-1,

or yr = yr+1 – yr, r = 0, 1,2,……… n – 1,

where  is called the forward difference operator. The differences of the first forward differences
are called second forward differences and are denoted by 2y0, 2y1,….. Similarly, one can
define third forward differences, fourth forward differences etc. Thus, we have

2y0 = y1 - y0 = y2 – y1 – (y1 – y0) = y2 – 2y1 + y0,

3y0 = 2y1 - 2y0 = y3 -2y2+ y1 – (y2 – 2y1 + y0)

= y3 - 3y2 + 3y1 – y0

4y0 = 3y1 - 3y0

= y4 - 3y3 + 3y2 – y1 – (y3 – 3y2 + 3y1 – y0)

= y4 - 4y3 + 6y2 – 4y1 + y0

And, in general

pyr = p-1 yr+1 - p-1 yr = yn – nc1 yn-1+ nc2 yn-2 +……+(-1)ny0.

It is, therefore, clear that any higher order difference can easily be expressed in terms of the
ordinates, since the coefficients occurring on the right side are the binomial coefficients.
These differences are systematically set out as follows in what is called a forward difference
table.

82
 Forward Difference Table
2 3
x y    4 5 6

x0 y0

y0

x1 y1 2y0

y1 3y0

x2 y2 2y1 4y0

y2 3y1 5y0

x3 y3 2y2 4y1 6y0

y3 3y2 5y1

x4 y4 2y3 4y2

y4 3y3

x5 y5 2y4

y5

x6 y6

Backward Differences:
The differences y1– y0, y2 – y1, .....,yn – yn-1 when denoted by y1, y2, …., yn
respectively, are called the first backward differences. So that y1 = y1 – y0, y2 = y2 –
y1,…..,yn = yn – yn-1, where  is called the backward difference operator. One can define
backward differences of higher orders as

2y2 = y2 - y1 = y2 – y 1 – (y1- y0) = y2 – 2y1 + y0,

3y3 = 2y3 - 2y2 = y3 – 3y2 + 3y1 – y0 .

These differences are exhibited in the form of a backward difference table as given below

83
 Backward Difference Table
x y  2 3 4 5 6

x0 y0

x1 y1 y1

x2 y2 y2 2y2

x3 y3 y3 2y3 3y3

x4 y4 y4 2y4 3y4 4y4

x5 y5 y5 2y5 3y5 4y5 5y5

x6 y6 y6 2y6 3y6 4y6 5y5 6y6

Central Differences:

The central difference operator  is defined by the relations

y1 – y0 = y1/2, y2 – y1 = y3/2,….., yn – yn-1 = yn – ½ .

The higher order central differences are defined as

y3/2 – y½ = 2y1, y5/2 – y3/2 = 2y2,…….
2y2 – 2y1 = 3y3/2, and so on.
These differences are shown in the following table:

Central Difference Table

x y  2 3 4 5 6

x0 y0

y1/2

x1 y1 2y1

84
 y3/2 3y3/2

x2 y2 2y2 4y2

y5/2 3y5/2 5y5/2

x3 y3 2y3 4y3 6y3

y7/2 3y7/2 5y7/2

x4 y4 2y4 4y4

y9/2 3y9/2

x5 y5 2y5

y11/2

x6 y6

We see from this table that the central differences on the same horizontal line have the
same suffix. Also, the differences of odd order are know only for half values of the suffix and
those of even order for only integral values of the suffix.
Observation: It is noted that in all the three tables, the same numbers occur in the same position
and it is only the notation which changes, e.g.,

y1 – y0 = y0 = y1 = y1/2

Other Difference Operators:

The operators ,  and  have already been defined. Besides these, there are operators E and ,
which are defined as follows :

Shift operator E is defined as the operation of increasing the argument x by h, so that

Ef(x) = f(x + h), E2f(x) = f(x + 2h), E3 f(x)= f(x+3h) etc.

The inverse operator E-1 is defined by

E-1 f(x) = f(x – h)

85
If yx is the function f(x), then

Eyx = yx+h, E-1yx = yx-h, Enyx = yx+nh,

where n may be any real number.

Averaging operator  is defined by the relation

1
yx = (y h+y h )
2 x 2 x 2

Remark: In the difference calculus E is regarded as the fundamental operator and , , ,  can
be expressed in terms of E.

Relations between the operators: Show that

(i) =E–1 (ii)  = 1 – E-1

(iii)  = E1/2- E-1/2 (iv)  = ½(E1/2 + E-1/2)

(v)  = E = E = E1/2 (vi) E = ehD

Proof:

i) yx = yx+h – yx = Eyx – yx = (E – 1)yx

This shows that the operators  and E are connected by the symbolic relation
 = E – 1 or E = 1 +

ii) yx = yx – yx-h = yx – E-1yx

= (1 – E-1) yx

  = 1 – E-1

(iii) yx = y h – y h
x x
2 2

= E1/2yx – E-1/2yx = (E1/2 – E-1/2)yx

86
  = E1/2 – E-1/2

1
(iv) yx = (y h+y h )
2 x 2 x 2

1 1/2 1
= (E yx + E-1/2yx) = (E1/2 + E-1/2)yx
2 2

1 1/2
  = (E + E-1/2)
2

(v) Eyx = E(yx – yx-h) = Eyx – Eyx – h

= yx+h – yx = yx

 E = 

 Eyx = yx+h = yx+h – yx = yx

 E = 

E1/2yx =  y h = y h h –y h h
x x  x 
2 2 2 2 2

= yx+h – yx = yx

 E1/2 = 

Hence  = E = E = E1/2

(vi) Ef(x) = f(x + h)

h2
= f(x) + hf  (x) + f (x) + ………….
2!

[by Taylor’s series]

h2 2
= f(x) + hDf(x) + D f(x) + ……
2!

87
  h 2 D 2 h3 D 3 
= 1  hD    .......... ....  f ( x )
 2! 3! 

= ehDf(x)

E = ehD

Example: Prove that

 2  x Ee x
ex =   e . 2 x ,
 E  e

the interval of differencing being h.

Solution:

 2  x
Since   e = 2. E-1 ex = 2 ex-h
 E 

= 2ex e-h = e-h2ex

 2  x Ee x Ee x
   e . 2 x = e-h2ex . 2 x = e-hEex
 E  e e

= e-h . ex+h = ex.

Example: Prove with the usual notations, that

i. (E1/2 + E-1/2) ( 1 + )1/2 = 2 + 

ii.  = ½ 2 + (1+2/4)
iii. 3y2 = 3y5.
Solution:

i. (E1/2 + E-1/2) (1+)1/2

= (E1/2 + E-1/2) E1/2

88
 =E+1=1++1

= 2 + .

ii. ½2 +  1   2 / 4

= ½(E1/2 - E-1/2)2 + (E1/2 – E-1/2) 1  ( E 1 / 2  E 1 / 2 ) 2 / 4

= ½( E+E-1 – 2) + (E1/2 – E-1/2) ( E  E 1  2) / 4

=½(E + E-1 – 2) + (E1/2 – E-1/2) (E1/2 + E-1/2)/2

=½[(E+ E-1 – 2) + (E – E-1)] = ½(2E – 2)

=E – 1 = .

iii. 3y2 = (E – 1)3y2 [ = E – 1]

= (E3 – 3E2 + 3E – 1) y2

= y5 – 3y4 + 3y3 – y2 ---(1)

3y5 = (1 – E-1)3y5 [ = 1 – E-1]

= (1 – 3E-1 + 3E-2 – E-3)y5

= y5 – 3y4 +3y3 – y2 ---(2)

From (1) and (2),

3y2 = 3y5.

Example: Using the method of separation of symbols, prove that

(i) u1x + u2x2 + u3x3 + …………

2 3
x  x   x  2
= u1     u1     u1  .......
1 x 1 x  1 x 

89
 u1 x u 2 x 2 u 3 x 3
(ii) u0 +    ..........
1! 2! 3!

 x2 2 x3 3 
= e  u 0  xu 0 
x
 u0   u 0  ........ 
 2! 3! 

Solution:

(i) L.H.S. = xu1 + x2 Eu1 +x3E2u1 +…… [ux+h = Ehux]

= x(1+xE + x2E2+…..)u1

1
=x. u1 , [taking sum of infinite G.P.]
1  xE

 1 
= x u1 [E= 1 + ]
1  x(1  ) 

1
 1  x  x 
= x  u1  1   u1
 1  x  x )  1 x  1 x 

x  x x 2 2 
= 1    ........ u1
1  x  1  x (1  x ) 2


x x2 x3
= u1   u 1  2 u1  .......... ..
1 x (1  x ) 3
2
(1  x )

= R.H.S.

x x2 2 x3 3
(ii) L.H.S. = u0 + Eu0 + E u0 + E u0+……….
1! 2! 3!

 xE x 2 E 2 x 3 E 3 
= 1     ........ u 0
 1! 2! 3! 

= exE u0 = ex(1+) u0

= ex. exu0

90
  x x 2 2 x 3 3 
= ex 1     ........ u 0
 1! 2! 3! 

 x x2 2 x3 
= ex  u 0  u 0   u 0   3 u 0  ...... 
 1! 2! 3! 

= R.H.S.

4.4 DETECTION OF ERRORS BY USE OF DIFFERENCE

TABLES

Difference tables can be used to check errors in tabular values. Suppose there is an error of +1
unit in a certain tabular value. As higher differences are formed, this error spreads out and is
considerably magnified. It affects the difference table as shown below in the table.
Error Difference Table shows the following characteristics:

i. The effect of the error increases with the order of the differences.
ii. The errors in any one column are the binomial coefficients with alternating signs.

Table

y  2 3 4 5

0 0

0 0

0 0 0

0 0 1

0 0 1

91
 0 1 -5

0 1 -4

1 -3 10

1 -2 6

-1 3 -10

0 1 -4

0 -1 5

0 0 1

0 0 -1

0 0 0

0 0

0 0

(iii) The algebraic sum of the errors in any difference column is zero, and
(iv) The maximum error occurs opposite the function value containing the error.
These facts can be used to detect errors by difference tables.

Example:

x y  2 3 4

1 3010

414

2 3424 -36

378 -39

92
3 3802 -75 +178

303 +139

4 4105 +64 -271

367 -132

5 4472 -68 +181

299 +49

6 4771 -19 -46

280 +3

7 5051 -16

264

8 5315

The term – 271 in the fourth difference column has fluctuations of 449 and 452 on either side of
it. Comparison with the error difference table suggest that there is an error of -45 in the entry for
x = 4. The correct value of y is therefore 4105 + 45 = 4150, which shows that the last two digits
have been transposed, a very common form of error. The reader is advised to form a new
difference table with this correction, and to check that the third differences are now practically
constant.
If an error is present in a given data, the differences of some order will become alternating in
sign. Hence, higher order differences should be formed till the error is revealed as in the above
example. If there are errors in several tabular values, then it is not easy to detect the errors by
differencing.

4.5 DIFFERENCES OF A POLYNOMIAL

The nth differences of a polynomial of the nth degree are constant and all higher order
differences are zero.

93
Let the polynomial of the nth degree in x, be

f(x) = axn + bx n-1 + cx n-2 + ……… + kx + l

 f(x) = f(x +h) – f(x)

= a[(x +h)n – xn] + b[(x +h)n-1- xn-1] + …… +kh

= anhxn-1 + bxn-2 + cxn-3 + ……..+kx + l, ---(1)

where b, c, …….. l are the new constant co-efficients.

Thus the first difference of a polynomial of the nth degree is a polynomial of degree (n –
1).

Similarly 2f(x) = [f(x +h) – f(x)]

= f(x+h) - f(x)

= anh[(x+h)n-1 – xn-1] + b[(x+h)n-2 – xn-2] +…..+kh

= an(n-1)h2xn-2 + bxn-3 + cxn-4 +.….+k

 The second differences represent a polynomial of degree (n – 2). Continuing this process,
for the nth differences we get a polynomial of degree zero, i.e.,

nf(x) = an(n – 1)(n – 2)……. 1hn

= an!hn , --- (2)

which is a constant. Hence the (n+1)th and higher differences of a polynomial of nth degree will
be zero.
Remark: The converse of this theorem is also true, i.e., if the nth differences of a function
tabulated at equally spaced intervals are constant, the function is a polynomial of degree n. This
fact is important in numerical analysis as it enables us to approximate a function by a polynomial
of nth degree, if its nth order differences become nearly constant.

Example: Evaluate

94
 10[(1-ax)(1-bx2)(1-cx3)(1-dx4)]

Solution:

10[(1-ax)(1-bx2)(1-cx3)(1-dx4)]

= 10[abcdx10 + ( )x9 + ( )x8 + ……+1]

= abcd 10(x10) [10(xn) = 0 for n <10]

= abcd (10!) [by (2) above ]

4.6 NEWTON’S FORWARD INTERPOLATION FORMULA

Let the function y = f(x) take the values y0, y1, y2,…. corresponding to the values x0, x0+h,
x0+2h,…. of x. Suppose it is required to find f(x) for x = x0+ph, where p is any real
number.

For any real number p, we have defined E such that

Epf(x) = f(x+ph)

 yp = f(x0+ph) = Epf(x0) = (1+)py0 [  E=1+]

p( p  1) 2 p( p  1)( p  2) 3 …….
= 1 + p +  + + y0,
2! 3!

[Using Binomial theorem]

p( p  1) 2 p( p  1)( p  2) 3
i.e., yp = y0 + py0 +  y0 +  y0 +….. ---(1)
2! 3!

If y = f(x) is a polynomial of the nth degree, then n+1y0 and higher differences will be zero.
Hence (1) will become

95
 p( p  1) 2 p( p  1)( p  2) 3
yp = y0 + py0 +  y0 +  y0 +…..
2! 3!

p( p  1).....( p  n  1) n
……+  y0 ---(2)
n!

Remark: This formula is used for interpolating the values of y near the beginning of a set of
tabulated values and extrapolating values of y a little backward (i.e. to the left) of y0.

4.7 NEWTON’S BACKWARDINTERPOLATION FORMULA

Let the function y = f(x) take the values y0, y1, y2,…. corresponding to the values x0, x0+h,
x0+2h,…. of x. Suppose it is required to evaluate f(x) for x = xn+ph, where p is any real number.
Then we have

yp = f(xn+ph) = Epf(xn) = (1-)-pyn [ E-1 = 1-]

p ( p  1) 2 p( p  1)( p  2) 3 …….
= 1 + p +  + + yn,
2! 3!

[Using Binomial theorem]

p ( p  1) 2 p( p  1)( p  2) 3
i.e., yp = yn + pyn +  yn +  yn+ ….. ---(1)
2! 3!

It is called Newton’s backward interpolation formula as it contains yn and backward

differences of yn.

Remark: This formula is used for interpolating the values of y near the end of a set of tabulated
values and also for extrapolating values of y a little ahead (to the right) of yn.

Example: The table gives the distance in nautical miles of the visible horizon for the given
heights in feet above the earth’s surface:

x = height: 100 150 200 250 300 350 400

96
y = distance: 10.63 13.03 15.04 16.81 18.42 19.90 21.27

Find the values of y when (i) x = 218 (ii) x = 410.

Solution

The difference table is as under

x y  2 3 4

100 10.63

2.40

150 13.03 -0.39

2.01 0.15

200 15.04 -0.24 -0.07

1.77 0.08

250 16.81 -0.16 -0.05

1.61 0.03

300 18.42 -0.13 -0.01

1.48 0.02

350 19.90 -0.11

1.37

400 21.27

(i) If we take x0 = 200, then y0 = 15.04, y0= 1.77, 2y0 = -0.16,

3 = 0.03 etc.

x  x 0 18
Since x = 218 and h = 50, p =  = 0.36
h 50

 Using Newton’s forward interpolation formula, we get

97
 p( p  1) 2 p( p  1)( p  2) 3
y218 = y0 + py0 +  y0 +  y0
2! 3!

p( p  1)( p  2)( p  3) 4
+  y 0 +…..
4!

0.36(0.64)
f(218) = 15.04 + 0.36(1.77) + (-0.16)
2

0.36(0.64)(1.64) 0.36(0.64)(1.64)(2.64)
+ (0.03)  (0.01)
6 24

= 15.04 + 0.637 + 0.018 + 0.002 + 0.0004

= 15.697, i.e., 15.7 nautical miles.

(ii) Since x = 410 is near the end of the table, we use Newton’s backward interpolation
formula.

x  x n 10
 Taking xn = 400, p = = = 0.2
h 50

Using the line of backward difference

yn = 21.27, yn = 1.37, 2yn = -0.11, 3yn = 0.02 etc.

Newton’s backward formula gives

p ( p  1) 2 p ( p  1)( p  2) 3
y410 = y400 + py400 +  y400 +  y400 +
2! 3!
p( p  1)( p  2)( p  3) 4
+  y 400
4!

0.2(1.2)
= 21.27 + 0.2 (1.37) + (-0.11)
2!

0.2(1.2)(2.2) 0.2(1.2)(2.2)(3.2)
+ (0.02)  (0.01)
3! 4!

= 21.27 + 0.274 – 0.0132 + 0.0018 – 0.0007

98
 = 21.53 nautical miles.

Example: From the following table, estimate the number of students who obtained marks
between 40 and 45:

Marks : 30-40 40-50 50-60 60-70 70-80

No. of student:31 42 51 35 31

Solution:

First we prepare the cumulative frequency table, as follows:

Marks less than (x): 40 50 60 70 80

No. of students (yx): 31 73 124 159 190

Now the difference table is

x yx yx 2yx 3yx 4yx

40 31

42

50 73 9

51 -25

60 124 -16 37

35 12

70 159 -4

31

80 190

We shall find y45, i.e., number of student with marks less than 45. Taking x 0 = 40, x = 45, we
have
99
 x  x0 5
p= = = 0.5 [ h = 10]
h 10

 Using Newton’s forward interpolation formula, we get

p( p  1) 2 p( p  1)( p  2) 3
y45 = y40 + p y40 +  y40 +  y40 +
2! 3!
p( p  1)( p  2)( p  3) 4
 y40.
3!

0.5(0.5) 0.5(0.5)(1.5)
= 31 + 0.5 x 42 + x9+ x (-25)
2 6
0.5(0.5)(1.5)(2.5)
+ x 37
24

= 31+21+1.125-1.5625-1.4453

= 47.87

The number of students with marks less than 45 is 47.87, i.e., 48. But the number of students
with marks less than 40 is 31.

Hence the number of students getting marks between 40 and 45=48-31=17.

Example: Find the cubic polynomial which takes the following values :

x: 0 1 2 3

f (x) : 1 2 1 10

Hence or otherwise evaluate f(4).

Solution:

The difference table is

x f(x) f(x) 2f(x) 3f(x)

0 1
100
 1

1 2 -2

-1 12

2 1 10

3 10

x0
We take x0= 0 and p = =x [ h = 1]
h

 Using Newton’s forward interpolation formula, we get

x x( x  1) 2 x( x  1)( x  2) 3
f(x) = f(0) + f(0) +  f(0)+  f(0)
1 1.2 1.2.3

x( x  1) x( x  1)( x  2)
= 1 + x(1) + (-2) + (12)
2 6

= 2x3 – 7x2 + 6x +1,

which is the required polynomial.

x  xn
To compute f(4), we take xn = 3, x = 4 so that p = =1
h

Using Newton’s backward interpolation formula, we get

p ( p  1) 2 p( p  1)( p  2) 3
f(4) = f(3) + pf(3) +  f(3) +  f(3)
1 .2 1.2.3

= 10 + 9 +10 +12 = 41

Which is the same value as that obtained by substituting x = 4 in the cubic polynomial above.

101
 4.8 CHECK YOUR PROGRESS

1. Show that
 1   f ( x )
  
 f ( x)  f ( x) f ( x  1)
2. Evaluate
 1  1
(i) 2  2  (ii) n  
 x  5x  6  x
3 With usual notations, show that
1
(i) (1  )(1  )  1 (ii)   (  )
2
(iii)  r f k  r f k r
4. Evaluate
4 (1  x)(1  2 x)(1  3 x)(1  4 x) , h=1
5. Estimate the value of f(22) and f(42)from the following data:
x: 20 25 30 35 40 45
f(x): 354 332 291 260 231 204
6. Find the number of men getting wages between Rs. 10 and 15 from the following data
Wages (Rs.): 0-10 10-20 20-30 30-40
Frequency: 9 30 35 42
7. Construct Newton’s forward interpolation polynomial for the following data:
x: 4 6 8 10
y: 1 3 8 16

102
 4.9 SUMMARY

1. We start to learn about the error in polynomial interpolation and detection of that error by
using difference table.
2. This chapter tells us how to find value near the top and near the end of table.
3. We learn Newton’s forward and backward interpolation formula for detecting the value at
any intermediate level.

4.10 KEYWORDS

Polynomial Interpolation:- In numerical analysis, polynomial interpolation is

the interpolation of a given data set by the polynomial of lowest possible degree that passes
through the points of the dataset.

Polynomial interpolation is a method of estimating values between known data points. When
graphical data contains a gap, but data is available on either side of the gap or at a few specific
points within the gap, an estimate of values within the gap can be made by interpolation.
Newton’s Forward Interpolation Formula:-
p( p  1) 2 p( p  1)( p  2) 3
yp = y0 + py0 +  y0 +  y0 +…..
2! 3!
Newton’s Backward Interpolation Formula:-
p ( p  1) 2 p( p  1)( p  2) 3
yp = yn + pyn +  yn +  yn+ …..
2! 3!

103
 4.11 SELF ASSESSMENT TEST

1. Find solution using Newton’s forward difference formula

x: 1891 1901 1911 1921 1931

f(x): 46 66 81 93 101
X= 1895

2. Find solution using Newton’s forward difference formula

x: 0 1 2 3
f(x): 1 0 1 10
X= -1

3. In the table below the value of y are consecutive terms of a series of which the number
21.6 is the 6th term. Find the 1st and the 10th term of the series.

X: 3 4 5 6 7 8 9
Y: 2.7 6.4 12.5 21.6 34.3 51.2 72.9
4. Construct difference table for following data
X: 1 2 4 5 6 7 8
Y: 23 33 45 56 65 77 89

4.12 ANSWER TO CHECK YOUR PROGRESS

(2) (i) -2/(x +2) (x +3) (x +4)

n(n  1)(n  2)......2.1
(ii) (1) n
x( x  1)( x  2)......( x  n)
(4) 576
(5) 352; 219
(6) 24
(7) 1.625

104
 4.13 REFERENCES/ SUGGESTEDREADINGS

1. Computer Oriented Numerical Methods, V. Rajaraman, PHI.

2. Introductory Methods of Numerical Analysis, S.S. Sastry, PHI.
3. Numerical Methods for Scientific and Engineering Computation, M.K.Jain, S.R.Lyenger,
R.K.Jain, Wiley Eastern Limited.
4. Numerical Methods with Programs in C, T. Veerarajan, T. Ramachandarn., Tata McGraw
Hill.

105
 SUBJECT: Computer Oriented Numerical Methods

COURSE CODE: BCA-PC(L)-122

AUTHOR: Prof. Dr. Aseem Miglani
LESSON NO. 5

INTERPOLION-II

REVISED/UPDATED SLM BY Amit

STRUCTURE
5.0 Objective
5.1 Introduction
5.2 Central Difference Interpolation Formulae
5.2.1 Gauss’s Forward Interpolation Formula
5.2.1 Gauss’s Backward Interpolation Formula
5.2.3 Stirling’s Formula
5.2.4 Bessel’s Formula
5.3 Interpolation with Unequal Intervals
5.3.1 Lagrange’s Interpolation Formula
5.4 Inverse Interpolation
5.4.1 Lagrange’s Method
5.4.2 Iterative Method
5.5 Check Your Progress
5.6 Summary
5.7 Keywords
5.8 Self Assessment Questions
5.9 Answer to Check Your Progress
5.10 References/ Suggested Readings

106
 5.0 OBJECTIVE

The objective of this lesson is to continue with the objective of the previous lesson by deriving
some more interpolation formulae such as, central difference interpolation formulae, which are
based on central difference table and operators; Interpolation formulae for unequi-spaced set of
values and also the inverse interpolation.

5.1 INTRODUCTION

The Newton’s forward and backward interpolation formulae are applicable for interpolation near
the beginning and end of the tabulated values. Now we shall develop central difference formulae
which are best suited for interpolation near the middle of the table, for equi-spaced set of values.

All the interpolation formulae discussed, have the disadvantage of requiring the values of the
independent variable to be equally spaced, i.e., equi-spaced set of values. It is therefore desirable
to have interpolation formulae with unequally spaced values of the argument and therefore, we
discuss one such formula, Lagrange’s interpolation formula. Also, the inverse interpolation is
discussed, in which we can find the value of the argument for a given value of the function, from
the set of tabulated values.

5.2 CENTRAL DIFFERENCE INTERPOLATION FORMULA

If x takes the values x0 – 2h, x0 – h, x0, x0 + h, x0 +2h and the corresponding values of y = f(x)
are y–2, y-1, y0, y1, y2, then we can write the difference table in the two notations as follows :

107
x y 1st diff. 2nd diff. 3rd diff. 4th diff.

x0 – 2h y– 2

y-2 (=y–3/2)

x0 –h y–1 2y-2(=2y-1)

y-1(=y-1/2) 3y-2(=3y – ½)

x0 y0 2y-1(=2y0) 4y-2(=4y0)

y0(=y1/2) 3y-1(=3y1/2)

x0+h y1 2y0(=2y1)

y1(=y3/2)

x0 + 2h y2

5.2.1 Gauss’s Forward Interpolation Formula:

The Newton’s forward interpolation formula is
p( p  1) 2 p( p  1)( p  2) 3
yp = y0 + p y0 +  y0 +  y0+ ---(1)
1.2 1.2.3

We have 2y0 - 2y-1 = 3 y–1,

i.e., 2y0 = 2y-1 + 3y-1. ---(2)
Similarly, 3y0 = 3y-1 + 4y-1, ---(3)
4y0 = 4y-1 + 5y-1 etc. ---(4)
Also, 3y-1 - 3y-2 = 4y-2,
i.e. 3y-1 = 3y-2 + 4y-2.
Similarly, 4y-1 = 4y-2 + 5y-2 etc. ---(5)
Substituting for 2y0, 3y0, 4y0 …… from (2), (3), (4), (5) in (1), we get
p( p  1) 2
yp = y0 + py0 + ( y-1 +3y-1)
1.2

108
 p( p  1)( p  2) 3
+ ( y-1+4y-1)
1.2.3
p( p  1)( p  2)( p  3) 4
+ ( y-1 +5y-1) + …..
1.2.3.4
Hence
p( p  1) 2 ( p  1) p( p  1) 3
yp= yo+ py0 +  y-1+  y-1
2! 3!
( p  1) p( p  1)( p  2) 4
+  y-2 +….., ---(6)
4!
which is known as Gauss’s forward interpolation formula.
In the central difference notations, this formula will be
p( p  1) 2 ( p  1) p( p  1) 3
yp= y0 + py1/2+  y0 +  y1/2
2! 3!
( p  1) p( p  1)( p  2) 4
+  y0 +….. ---(7)
4!
Remark :This formula is used to interpolate the values of y for p (0<p<1) measured forwardly
from the origin.

5.2.2 Gauss’s Backward Interpolation Formula:

The Newton’s forward interpolation formula is
p( p  1) 2 p( p  1)( p  2) 3
yp = y0 + p y0 +  y0 +  y0+ ---(1)
1.2 1.2.3

We have y0 - y-1 = 2 y–1,

i.e., y0 = y-1 + 2y-1. ---(2)
Similarly, 2y0 = 2y-1 + 3y-1, ---(3)
3y0 = 3y-1 + 4y-1 etc. ---(4)
Also, 3y-1 - 3y-2 = 4y-2,
i.e. 3y-1 = 3y-2 + 4y-2 ---(5)
Similarly, 4y-1 = 4y-2 + 5y-2 etc. ---(6)
Substituting for y0, 2y0, 3y0 …… from (2), (3), (4), (5), (6) in (1), we get

109
 p( p  1) 2
yp= y0 + p( y-1 +2y-1) + ( y-1 +3y-1)
1.2
p( p  1)( p  2) 3
+ ( y-1+4y-1)
1.2.3
p( p  1)( p  2)( p  3) 4
+ ( y-1 +5y-1) +..…
1.2.3.4
( p  1) p 2 ( p  1) p( p  1) 3
= y0+ py-1 +  y-1+  y-1
1.2 1.2.3
( p  1) p( p  1)( p  2) 4
+  y-1
1.2.3.4
( p  1) p( p  2)( p  3) 5
+  y-1 +.....
1.2.3.4
( p  1) p 2 ( p  1) p( p  1) 3
= y0+ py-1 +  y-1+ ( y-2+4y-2)
1.2 1.2.3
( p  1) p( p  1)( p  2) 4
+ ( y-2+5y-2)+.....
1.2.3.4
[Using (5) and (6)]
Hence
( p  1) p 2 ( p  1) p( p  1) 3
yp= y0+ py-1 +  y-2+  y-2
2! 3!
( p  2)( p  1) p( p  1) 4
+  y-2 +……., ---(7)
4!
which is called Gauss’s backward interpolation formula.
In the central difference notations, this formula will be
( p  1) p 2 ( p  1) p( p  1) 3
yp= y0 + py-1/2+  y0 +  y-1/2
2! 3!
( p  2)( p  1) p( p  1) 4
+  y0 +……. ---(8)
4!
Remark: It is used to interpolate the values of y for a negative value of p lying between –1 and
0.
Gauss’s forward and backward formulae are themselves not of much practical use. However,
these serve as intermediate steps for obtaining the following two important formulae:

110
5.2.3Sterling’s Formula:
Gauss’s forward interpolation formula is
p( p  1) 2 ( p  1) p( p  1) 3
yp= y0 + py0+  y-1 +  y-1
2! 3!
( p  1) p( p  1)( p  2) 4
+  y-2 +……. ---(1)
4!
Gauss’s backward interpolation formula is
( p  1) p 2 ( p  1) p( p  1) 3
yp= y0 + py-1+  y-1 +  y-2
2! 3!
( p  2)( p  1) p( p  1) 4
+  y-2 +……. ---(2)
4!
Taking the mean of (1) and (2), we obtain

 y0  y1  p 2 2 p ( p 2  1)  3 y 1  3 y  2 
yp=y0 +p  +  y-1 +  
 2  2! 3!  2 

p 2 ( p 2  1) 4
+  y -2 +…… ---(3)
4!
which is called Stirling’s formula
In the central difference notations, (3) takes the form
p2 2 2 2
yp= y0 + py0+  y0 + p ( p  1 ) 3y0
2! 3!
2 2 2
+ p ( p  1 ) 4y0 +…… , ---(4)
4!

1 1
since (y0 +y-1) = (y1/2 +y-1/2) = y0,
2 2
1 3 1
( y-1 +3y-2) = (3y1/2 +3y-1/2) = 3y0 etc.
2 2
Remark: This formula involves means of the odd differences just above and below the central
line and even differences on this line as shown below :

 y   3 y 2   5 y 3  6
……y0…..  1  …..2y-1…..  3 .....4 y 2 ......  ... y 3 ...... central line
  5 y 
 y 0    y 1   2 

111
5.2.4 Bessel’s Formula:
Guass’s forward interpolation formula is
p( p  1) 2 ( p  1) p( p  1) 3
yp= y0 + py0+  y-1 +  y-1
2! 3!
( p  1) p( p  1)( p  2) 4
+  y-2 +…… ---(1)
4!
We have 2y0 - 2y-1 = 3y– 1 ,
i.e. 2y–1 = 2y0 - 3y-1 ---(2)
Similarly, 4y–2 = 4y–1 - 5y– 2, etc. ---(3)
Now (1) can be written as
p( p  1) 1 2 1 2  p( p 2  1) 3
yp= y0 + py0+   y1   y1  +  y-1
2! 2 2  3!

p ( p 2  1)( p  2) 1 4 1 4 
+   y 2   y 2  +……
4! 2 2 
1 p ( p  1) 2 1 p ( p  1) 2
= y0 + py0 +  y-1 + ( y0 - 3y-1)
2 2! 2 2!
p ( p 2  1) 3 1 p( p 2  1)( p  2) 4
+  y-1 +  y-2
3! 2 4!
1 p( p 2  1)( p  2)
+ (4y-1 - 5y-2) + …..
2 4!
[using (2), (3) etc.]
p( p  1) 2 y1  2 y0 p( p  1)  p  1 1  3
= y0 + py0 + . +     y-1
2! 2 2!  3 2

p ( p 2  1)( p  2) 4 y 2  4 y1
+ . + ……
4! 2
Hence
1
p( p  1) 2 y1  2 y0 ( p  2 ) p( p  1) 3
yp = y0 + py0 + . +  y-1
2! 2 3!
( p  1) p( p  1)( p  2) 4 y 2  4 y1
+ . + ....., ---(4)
4! 2
which is known as the Bessel’s formula.

112
 In the central difference notations, (4) becomes
1
p( p  1) 2 ( p  ) p( p  1)
yp = y0 + py1/2+  y1/2+ 2 3y1/2
2! 3!
( p  1) p( p  1)( p  2) 4
+  y1/2+....., ---(5)
4!
1 2 1
since ( y1  2 y0 )  2y1/2 , (4 y 2  4 y1 )  4y1/2, etc.
2 2
Remark: This is a very useful formula for practical purposes. It involves odd differences below
the central line and means of even differences of and below his line.

Choice of an interpolation formula :

The right choice of an interpolation formula, depends on the position of the interpolated value in
the given data.
The following rules will be found useful :
1. To find a tabulated value near the beginning of the table, use Newton’s forward
formula.
2. To find a value near the end of the table, use Newton’s backward formula.
3. To find an interpolated value near the centre of the table, use either Stirling’s or
Bessel’s formula.]
1 1
If interpolation is required for p lying between – and , use Stirling’s formula. If
4 4
1 3
interpolation is desired for p laying between and , use Bessel’s formula
4 4

Example: Employ Stirling’s formula to compute y12.2 from the following table (yx = 1+ log10 sin
x) :
x°: 10 11 12 13 14
105ux: 23,967 28,060 31,788 35,209 38,368

Solution:
Taking the origin at x = 12°, h = 1 and p = x – 12 , we have the following central table :

113
p yx yx 2yx 3yx 4yx

-2 0.23967
0.04093
-1 0.28060 -0.00365
0.03728 .00058
0 0.31788 -0.00307 -0.00013
0.03421 -0.00045
1 0.35209 -0.00062
0.03159
2 0.38368

1 1
At x = 12.2, p = 0.2. (As p lies between – and , the use of Stirling’s formula will be quite
4 4
suitable.)
Stirling’s formula is
p y 1  y 0 p 2 2 p ( p 2  1) 3 y 2  3 y1
yp= y0 + . +  y-1 +
1 2 2! 3! 2
p 2 ( p 2  1) 4
+  y-2 + ......
4!
When p = 0.2, we have
2
 0.03728 0.03421 (0.2)
y0.2 = 0.31788 + 0.2   (-0.00307)
 2  2
2 2
(0.2)[(0.2) 2  1]  0.00058  0.00045 (0.2) [(0.2)  1]
+    (- 0.00013)
6  2  24
= 0.31788+0.00715-0.00006-0.000002+0.0000002
= 0.32497.
Example: Apply Bessel’s formula to obtain y25, given y20 = 2854, y24 = 3162, y28 = 3544, y32 =
3992.
Solution:
1
Taking the origin at x0 = 24, h = 4, we have p = (x-24).
4
 The central difference table is

114
p y y 2 y 3 y

-1 2854
308
0 3162 74
382 -8
1 3544 66
448
2 3992

At x = 25, p =(25-24)/4 = ¼ . (As p lies between ¼ and ¾, the use of Bessel’s formula will yield
accurate result.)
Bessel’s formula is
1
( p  ) p( p  1)
p( p  1) 2 y1  2 y0 2
yp = y0 + py0+ . + 3y-1 +......
2! 2 3!
---(1)
When p = 0.25, we have
0.25(0.75)  74  66 
yp = 3162 + 0.25 x 382 +  
2  2 
(0.25)(0.25)(0.75)
+ (-8)
6
= 3162 + 95.5 – 6.5625 – 0.0625 = 3250.875 approx.

5.3 INTERPOLATION WITH UNEQUAL INTERVALS

The various interpolation formulae derived so far being applicable only to equally spaced values
of the argument. It is, therefore, desirable to develop interpolation formulae for unequally spaced
values of x. Lagrange’s interpolation formula is one such formula and is as follows:

5.3.1 Lagrange’s Interpolation Formula:

If y = f(x) takes the value y0, y1,……yn corresponding to x = x0, x1, ..…, xn then

115
 ( x  x1 )( x  x2 ).....(x  xn ) ( x  x0 )( x  x2 ).....(x  xn )
f(x) = y0 + y1
( x0  x1 )( x0  x2 )....(x0  xn ) ( x1  x0 )( x1  x2 )....(x1  xn )

( x  x0 )(x  x1 ).....(x  xn1 ) _

+…...+ yn ---(1)
( xn  x0 )(xn  x1 )....(xn  xn1 )
This is known as Lagrange’s interpolation formula for unequal intervals.
Proof: Let y = f(x) be a function which takes the values (x0, y0) (x1, y1),…. (xn, yn).
Since there are n +1 pairs of values of x and y, we can represent f(x) by a polynomial in x of
degree n. Let this polynomial be of the form
y = f(x) = a0(x-x1) (x-x2) …… (x-xn) + a1(x-x0) (x-x2) …… (x-xn)
+ a2(x-x0) (x-x1) (x-x3) ……(x-xn)+……
+ an(x-x0) (x-x1) …… (x-xn-1) --- (2)
Putting x= x0, y = y0, in (2), we get
y0 = a0(x0-x1) (x0-x2) ……….. (x0-xn)
 a0= y0/[(x0-x1) (x0-x2) ……….. (x0-xn)]
Similarly putting x = x1, y = y1 in (2), we have
a1 = y1/[(x1-x0) (x1-x2) ……….. (x1-xn)]
Proceeding the same way, we find a2, a3…… an. Substituting the values of a0, a1, ….., an in (2),
we get (1).
Remark: Lagrange’s formula can be applied whether the values xiare equally spaced or not. It is
easy to remember but quite cumbersome to apply.
Example: Given the values
x: 5 7 11 13 17
f(x): 150 392 1452 2366 5202
Evaluate f(9), using Lagrange’s formula
Solution:
Here x0 = 5, x1 = 7, x2 = 11, x3 = 13, x4 = 17
y0 = 150, y1 = 392, y2 = 1452, y3 = 2366, y4 = 5202
Putting x = 9 and substituting the above values in Lagrange’s formula, we get
(9  7)(9  11)(9  13)(9  17) (9  5)(9  11)(9  13)(9  17)
f(9) = 150 +  392
(5  7)(5  11)(5  13)(5  17) (7  5)(7  11)(7  13)(7  17)

116
 (9  5)(9  7)(9  13)(9  17) (9  5)(9  7)(9  11)(9  17)
+ 1452 +  2366
(11  5)(11  7)(11  13)(11  17) (13  5)(13  7)(13  11)(13  17)
(9  5)(9  7)(9  11)(9  13)
+  5202
(17  5)(17  7)(17  11)(17  13)
50 3136 3872 2366 578
=-      810
3 15 3 3 5

5.4 INVERSE INTERPOLATION

So far, we have been finding the value of y corresponding to a certain value of x from a given set
of values of x and y. On the other hand, the process of finding the value of x for a value of y is
called the inverse interpolation. When the values of x are unequally spaced, Lagrange’s method
is used and when the values of x are equally spaced, the Iterative method should be used.

5.4.1 Lagrange’s Method:

This is similar to Lagrange’s interpolation formula, the only difference being that x is
assumed to be expressible as a polynomial in y.
Lagrange’s formula is merely a relation between two variables either of which may be
taken as the independent variable. Therefore, on inter-changing x and y in the Lagrange’s
formula, we obtain.
( y  y1 )( y  y 2 ).......( y  y n ) ( y  y 0 )( y  y 2 ).......( y  y n )
x= x0  x1
( y 0  y1 )( y 0  y 21 ).......( y 0  y n ) ( y1  y 0 )( y1  y 2 ).......( y1  y n )

( y  y0 )( y  y1 ).......( y  yn 1 )
+ ..........  xn ,
( yn  y0 )( yn  y1 ).......( yn  yn 1 )
which is used for inverse interpolation.
Example: The following table gives the values of x and y :
x: 1.2 2.1 2.8 4.1 4.9 6.2
y: 4.2 6.8 9.8 13.4 15.5 19.6
Find the value of x corresponding to y = 12, using Lagrange’s technique.
Solution:

117
Here x0 = 1.2, x1 = 2.1, x2 = 2.8, x3 = 4.1 x4 = 4.9, x5 = 6.2
and y0 = 4.2, y1 = 6.8, y2 = 9.8, y3 = 13.4, y4 = 15.5, y5 = 19.6
Taking y = 12, the above formula gives
(12  6.8)(12  9.8)(12  13.4)(12  15.5)(12  19.6)
x=  1.2
(4.2  6.8)(4.2  9.8)(4.2  13.4)(4.2  15.5)(4.2  19.6)
(12  4.2)(12  9.8)(12  13.4)(12  15.5)(12  19.6)
+  2.1
(6.8  4.2)(6.8  9.8)(6.8  13.4)(6.8  15.5)(6.8  19.6)
(12  4.2)(12  6.8)(12  13.4)(12  15.5)(12  19.6)
+  2.8
(9.8  4.2)(9.8  6.8)(9.8  13.4)(9.8  15.5)(9.8  19.6)
(12  4.2)(12  6.8)(12  9.8)(12  15.5)(12  19.6)
+  4.1
(13.4  4.2)(13.4  6.8)(13.4  9.8)(13.4  15.5)(13.4  19.6)
(12  4.2)(12  6.8)(12  9.8)(12  13.4)(12  19.6)
+  4.9
(15.5  4.2)(15.5  6.8)(15.5  9.8)(15.5  13.4)(15.5  19.6)
(12  4.2)(12  6.8)(12  9.8)(12  13.4)(12  15.5)
+  6.2
(19.6  4.2)(19.6  6.8)(19.6  9.8)(19.6  13.4)(19.6  15.5)
= 0.022-0.234+1.252+3.419-0.964+0.055
= 3.55.

5.4.2 Iterative Method or Method of Successive Approximations:

Newton’s forward interpolation formula is
p( p  1) 2 p( p  1)( p  2) 3
yp= y0 +py0+  y0+  y0+……
2! 3!
From this, we obtain
1 p( p  1) 2 p( p  1)( p  2) 3
p= [yp – y0 -  y0-  y0 - …….]
y0 2! 3!
---(1)
Neglecting the second and higher order differences, we obtain the first approximation to p as
p1 = (yp – y0)/y0 ---(2)
To find the second approximation, retaining the term with second differences in (1) and
replacing p by p1, we get

118
 1 p ( p  1) 2
p2 = [yp – y0 - 1 1  y0] ---(3)
y0 2!

To find the third approximation, retaining the term with third differences in (1) and replacing
every p by p2, we have
1 p ( p  1) 2 p ( p  1)( p2  2) 3
p3 = [yp – y0 - 2 2  y0- 2 2  y0],
y0 2! 3!
and so on. This process is continued till two successive approximations of p agree with each
other, to the desired accuracy.
Remark: This technique can equally well be applied by starting with any other interpolation
formula.
Example: The following values of y = f(x) are given
x: 10 15 20
y: 1754 2648 3564

Find the value of x for y = 3000 by iterative method.

Solution:

Taking x0 = 10 and h = 5, the difference table is

x y y 2 y

10 1754
894
15 2648 22
916

20 3564

Here yp = 3000, y0 =1754, yo = 894 and 2yo=22.

 The successive approximations to p are

1
p1 = 3000  1754  1.39
894

119
 1  1.39(1.39  1) 
p2 = 3000  1754   22  1.387 ,
894  2 

1  1.387(1.387  1) 
p3 = 3000  1754   22 ,
894  2 

= 1.3871

We, therefore, take p = 1.387 correct to three decimal places.

Hence the value of x (corresponding to y =3000) = x0 +ph

= 10 + 1.387 x 5 =16.935.

Example: Using inverse interpolation, find the real root of the equation x3+ x –3= 0, witch is
close to 1.2.

Solution:
The difference table is
x v y(=x3 + x-3) y 2 y 3 y 4 y

1 -0.2 -1
0.431
1.1 -0.1 -0.569 0.066
0.497 0.006
1.2 0.0 -0.072 0.072 -0.00013
0.569 -0.006
1.3 0.1 0.497 0.078
0.647
1.4 0.2 1.144

Clearly the root of the given equation lies between 1.2 and 1.3.
Assuming the origin at x =1.2 and using Stirling’s formula.

y = yo+ x
y o  y 1 x 2 2
  y 1 

x( x 2  1) 3 y 1  3 y  2,
2 2 6 2
we get

120
 0.569  0.467 x 2
0 = -0.072 + x  0.072
2 2


 .

x x 2  1 0.006  0.006
[y = 0]
6 2
or 0 = -0.072 + 0.532x+0.036x2 + 0.001x3
This equation can be written as
0.072 0.036 2 0.001 3
x  x  x …..(i)
0.532 0.532 0.532
0.072
 First approximation x(1) =  0.1353
0.532
Putting x = x(1) on R.H.S. of (i), we get second approximation x(2) as
x(2) = 0.1353-0.067(0.1353)2-1.8797(0.1353)3
= 0.134
Hence the desired root = 1.2 +0.1 x 0.134 = 1.2134

5.5 CHECK YOUR PROGRESS

1. Use Stirling’s formula to evaluate f(1.22) from the following

x: 1.0 1.1 1.2 1.3 1.4
f(x): 0.841 0.891 0.932 0.963 0.985
2. Use Bessel’s formula to obtain y25, given
y20 = 24 y24 = 32 y28 = 35 y32 = 40
3. Use Lagrange’s interpolation formula to find the value of y when x = 10, if the following
values of x and y are given:
x: 5 6 9 11
f(x): 12 13 11 16
4. From the following data:
x: 1.8 2.0 2.2 2.4 2.6
f(x): 2.9 3.6 4.4 5.5 6.7
find x when y = 5 using the Iterative method.

121
5. The equation x3 – 15x + 4 = 0 has a root close to 0.3, obtain this root upto 4 decimal
places using inverse interpolation.

5.6 SUMMARY

1. In this chapter we learn to find the value which lies near the center of table.
2. We learn about the center difference interpolation method. For this we use gauss’s forward
and gauss’s backward interpolation formulas.
3. We learn to solve problem having unequal interval. For this purpose we use Lagrange’s
Interpolation Formula.

5.7 KEYWORDS

Gauss’s Forward Interpolation Formula: -

p( p  1) 2 ( p  1) p( p  1) 3
yp= yo+ py0 +  y-1+  y-1
2! 3!
( p  1) p( p  1)( p  2) 4
+  y-2 +…..,
4!
Gauss’s Backward Interpolation Formula:
( p  1) p 2 ( p  1) p( p  1) 3
yp= y0+ py-1 +  y-2+  y-2
2! 3!
( p  2)( p  1) p( p  1) 4
+  y-2 +…….,
4!
Stirling’sformula:-

 y0  y1  p 2 2 p ( p 2  1)  3 y 1  3 y  2 
yp=y0 +p  +  y-1 +  
 2  2! 3!  2 

p 2 ( p 2  1) 4
+  y -2 +……
4!

122
 Bessel’s formula:-
1
p( p  1) 2 y1  2 y0 ( p  2 ) p ( p  1) 3
yp = y0 + py0 + . +  y-1
2! 2 3!
( p  1) p( p  1)( p  2) 4 y 2  4 y1
+ . + .....,
4! 2
Lagrange’s Interpolation Formula:-
( x  x1 )( x  x 2 ).....( x  x n ) ( x  x0 )( x  x 2 ).....( x  x n )
f(x) = y0 + y1
( x0  x1 )( x0  x 2 )....( x0  x n ) ( x1  x0 )( x1  x 2 )....( x1  x n )

( x  x 0 )( x  x1 ).....( x  x n 1 ) _
+…...+ yn
( x n  x 0 )( x n  x1 )....( x n  x n 1 )
This is known as Lagrange’s interpolation formula for unequal intervals.

5.8 SELF ASSESSMENT TEST

1. Use Stirling’s formula to evaluate f(22.2) from the following

x: 2.0 2.1 2.2 2.3 2.4
f(x): 0.23967 0.28062 0.31788 0.35209 0.38368
2. Using gauss forword interpolation formula, find f(337.5)

X 310 320 330 340 350 360

Y 2.4914 2.5052 2.5185 2.5315 2.5441 2.5563

3. Using gauss forword interpolation formula, find f(30)

X 21 25 29 33 37
Y 18.4708 17.8144 17.1070 16.3432 15.5154
4. Using gauss backward interpolation formula, find f(1966)

X 1931 1941 1951 1961 1971 1981

Y 12 15 20 27 39 52

123
5. Using gauss backward interpolation formula, find f(1974)

X 1939 1949 1959 1969 1979 1989

Y 12 15 20 27 39 52
6. Using Bessel’s interpolation formula, find f(25)

X 20 24 28 32
Y 2854 3162 3544 3992
7. Using Bessel’s interpolation formula, find f(25)

X 20 24 28 32
Y 24 32 35 40

5.9 ANSWER TO CHECK YOUR PROGRESS

(1) 0.934

(2) 32.945

(3) 14.63

(4) 2.3

(5) 0.2679

5.10 REFERENCES/ SUGGESTEDREADINGS

1. Computer Oriented Numerical Methods, V. Rajaraman, PHI.

2. Introductory Methods of Numerical Analysis, S.S. Sastry, PHI.
3. Numerical Methods for Scientific and Engineering Computation, M.K.Jain, S.R.Lyenger,
R.K.Jain, Wiley Eastern Limited.
4. Numerical Methods with Programs in C, T. Veerarajan, T. Ramachandarn., Tata McGraw
Hill.

124
 SUBJECT: Computer Oriented Numerical Methods

COURSE CODE: BCA-PC(L)-122

AUTHOR: Prof. Dr. Aseem Miglani
LESSON NO. 6

NUMERICAL DIFFERENTIATION AND INTEGRATION

REVISED/UPDATED SLM BY AMIT

STRUCTURE
6.0 Objective
6.1 Introduction
6.2 Numerical Differentiation
6.2.1 Derivatives using Forward Difference Formula
6.2.2 Derivatives using Backward Differences Formula
6.2.3 Derivatives using Central Differences Formulae
6.3 Maxima and Minima of a Tabulated Function
6.4 Numerical Integration
6.4.1 Trapezoidal Rule
6.4.2 Simpson’s one-third Rule
6.4.3 Simpson’s three-eighth Rule
6.5 Errors in Quadrature Formulae
6.6 Check Your Progress
6.7 Summary
6.8 Keywords
6.9 Self Assessment Questions
6.10 Answer to Check Your Progress
6.11 References/ SuggestedReading

125
 6.0 OBJECTIVE

Objective of this lesson is to develop numerical techniques for finding the differentiation and
integration of the functions, which are given in the tabulated forms, using the interpolation
formulae. In the end, the errors in numerical integration formulae are also obtained.

6.1 INTRODUCTION

If a function y = f(x) be defined at a set of n+1 distinct points x0, x1,…..,xn-1,xn lying in some
interval [a, b] such that a = x0<x1<x2<…..<xn = b. From the given tabulated data (set of values of
x and y), we require to find differentiation of different orders at tabular or non tabular points.
b
Also, we require to find the values of the definite integral  f ( x)dx , where f(x) is either given
a

explicitly or defined by a tabulated data.

6.2 NUMERICAL DIFFERENTIATION

It is the process of calculating the value of the derivative of a function at some assigned value of
x from the given set of values (xi, yi). To compute dy/dx, we first replace the exact relation y=
f(x) by the best interpolating polynomial y = (x) and then differentiate the latter as many times
as we desire. The choice of the interpolation formula to be used, will depend on the assigned
value of x at which dy/dx is desired.

If the values of x are equi-spaced and dy/dx is required near the beginning of the table, we
employ Newton’s forward formula. If it is required near the end of the table, we use Newton’s
backward formula. For values near the middle of the table, dy/dx is calculated by means of
Stirling’s or Bessel’s formula.

126
If the values of x are not equi-spaced, we use Newton’s divided difference formula to represent
the function.

Hence corresponding to each of the interpolation formulae, we can derive a formula for finding
the derivative.

Consider the function y = f(x) which is tabulated for the values xi (= x0 + ih), i= 0, 1,2, ....n.

6.2.1 Derivatives using Forward Difference Formula:

Newton’s forward interpolation formula is

p( p  1) 2 p( p  1)( p  2) 3
y  y 0  py 0   y0   y 0  ...... ---(1)
2! 3!

Differentiating both sides w.r.t. p, we have

dy 2p 1 2 3p2  6 p  2 3
 y 0   y0   y 0  ......
dp 2! 3!

( x  x0 ) dp 1
Since p = , therefore  .
h dx h

Now

dy dy dp 1  2 p 1 2 3p2  6 p  2 3
 .  yo   y0   y0
dx dp dx h  2! 3!

4 p 3  18 p 2  22 p  6 4 
  y0  ..... ---(2)
4! 

At x = x0, p = 0. Hence putting p = 0,

 dy  1 1 2 1 3 1 4 
   y o   y 0   y 0   y 0  ......... ---(3)
 dx  x0 h  2 3 4 

Again differentiating (2) w.r.t. x, we get

127
 d 2 y d  dy  dp
  
dx 2 dp  dp  dx

1 2 2 6p 6 3 12 p 2  36 p  22 4 1
   y0   y0   y 0  ....
h  2! 3! 4! h

Putting p = 0, we obtain

d2y 1  2 11 4 
 2   2  y 0   y 0  12  y 0  .....
3
---(4)
 dx  x0 h  

6.2.2 Derivatives using Backward Differences Formula:

Newton’s backward interpolation formula is

P( p  1) 2 P( P  1)( p  2) 3
y = yn + py n   yn   y n  ......
2! 3!

Differentiating both sides w.r.t. p, we get

dy 2 p 1 2 3p2  6 p  2 3
 y n   yn   y n  .......
dx 2! 3!

x  xn dp 1
Since p  , therefore 
h dx h

Now dy  dy . dp
dx dp dx

 2 p 1 2 3p3  6 p  2 2 
= 1 y n   yn   y n  ...... ---(5)
h 2! 3! 

At x = xn, p = 0. Hence putting p = 0, we get

 dy  1 1 2 1 3 1 4 
   y n   y n   y n   y n .......... ---(6)
 dx  xn h  2 3 4 

Again differentiating (5) w.r.t. x, we have

128
 d 2 y d  dy  dp
  
dx 2 dp  dx  dx

1 2 6P  6 3 6 p 2  18 p  11 4 
  y n   y n   y n  ....
h 3! 12 

Putting p = 0, we obtain

d2y
 2  = 2  2 y n   3 y n   4 y n ....
1 11
---(7)
 dx  xn h  12 

6.2.3 Derivatives using Central Difference Formulae:

Stirling’s formula is

p  y 0  y 1  p 2 2
yP= y0+    y 1
1!  2  2!

 3 y 1  3 y  2  p 2 ( p 2  12 ) 4
+ p( p  1 ) 
2 2
  y  2  .... ---(8)
 
3!  2  4!

Differentiating both sides w.r.t. p, we get

dy  y 0  y 1  2 p 2 3 p 2  1  3 y 1  3 y  2  4 p3  2p 4
 +  y 1 +    y  2  .....
dx  2  2! 3!  2 
 4!

x  x0 dp 1
Since p = ,   .
h dx h

Now

dy dy dp
 .
dx dp dx

1  y 0  y 1  3 p 2  1  3 y 1  3 y  2  2p3  p 4 
=   +p 2
y +    y  .....
-1
   2
h  2  6  2  12 

At x = x0, p = 0. Hence putting p = 0, we get

129
  dy  = 1  y 0  y 1 1  y 1   y  2 1  y  2   y 3 
3 3 5 5

      ........ --(9)
 dx  x0 h  2 6 2 30 2 

d2y 1  2 1 4 1 6 
similarly  2   2  y 1  12  y  2  90  y 3 ..... ---(10)
 dx  x0 h  

Similarly, we can use any other interpolation formula for computing the derivatives.

Example : Given that

x: 1.0 1.1 1.2 1.3 1.4 1.5 1.6

y: 7.989 8.403 8.781 9.129 9.451 9.750 10.031,

dy d2y
Find and at (a ) x  1.1 (b) x = 1.6
dx dx 2

Solution:

The difference table is

x y  23 4 5 6

1.0 7.989

0.414

1.1 8.403 -0.036

0.378 0.006

1.2 8.781 -0.030 -0.002

0.348 0.004 0.002

1.3 9.129 -0.026 0.000 -0.003

0.322 0.004 -0.001

1.4 9.451 -0.023 -0.001

130
 0.299 0.005

1.5 9.750 -0.018

0.281

1.6 10.031

a) For x = 1.1, we take

 dy  1  1 2 1 3 1 4 1 5 1 6 
   y0   y0   y0   y0   y0   y0  ........ --(i)
 dx  x0 h  2 3 4 5 6 

 d2y  1  11 5 137 6 
and  2   2 2 y0  3 y0  4 y0  5 y0   y0  ....... ---(ii)
 dx  x0 h  12 6 180 

Here h = 0.1, x0 = 1.1, y0 = 0.378, 2y0 = -0.03 etc.

Substituting these values in (i) and (ii), we get

 dy  1  1 1 1 1 1 
   0.378  (0.03)  (0.004)  (0)  (0.001)  (0.003)
 dx 1.1 0.1  2 3 4 5 6 
= 3.946

 d2y  1  11 5 137 
 2    0.03  (0.004)  12 (0)  6 (0.001)  180 (0.003)
 
2
 dx 1.1 (0.1)

= -3.545

b) For x = 1.6, we use the above difference table and the backward difference operator 
instead of , and we take

 dy  1  1 2 1 3 1 4 1 5 1 6 
   yn  2  yn  3  yn  4  yn  5  yn  6  yn  ...... ---(iii)
 dx  xn h  

and

131
  d2y  1  11 5 137 6 
 2   2  2 yn  3 yn   4 yn  5 yn   yn  ...... ---(iv)
 dx  xn h  12 6 180 

Here h = 0.1, xn = 1.6,  yn = 0.281,  2 yn = -0.018 etc.

Putting these values in (iii) and (iv), we get

 dy  1  1 1 1 1 1 
   0.281  ( 0.018)  (0.005)  ( 0.001)  ( 0.001)  (0.003)
 dx 1.6 0.1  2 3 4 5 6 

= 2.727

 d2y  1  11 5 137 
 2    0.018  0.005  12 (0.001)  6 (0.001)  180 (0.003)
 
2
 dx 1.6 (0.1)

= -1.703

dx d 2 x
Example: From the following table find & at t = 0.3
dt dt 2

T: 0 0.1 0.2 0.3 0.4 0.5 0.6

X 30.13 31.62 32.87 33.64 33.95 33.81 33.24

Solution :

The difference table is:

t x  2 3 4 5 6

0 30.13

1.49

0.1 31.62 -0.24

1.25 -0.24

0.2 32.87 -0.48 0.26

132
 0.77 0.02 -0.27

0.3 33.64 -0.46 -0.01 0.29

0.31 0.01 0.02

0.4 33.95 -0.45 0.01

-0.14 0.02

0.5 33.81 -0.43

-0.57

0.6 33.24

As the derivatives are required near the middle of the table, we use Stirling’s formulae:

 dx  = 1  x0  x1  1   x1   x2  1   x2   x3 

3 3 5 5
           ....... --(i)
 dt t0 h  2  6 2  30  2 

 d 2x  1  2 1 4 1 6 
 2   2  x1  12  x 2  90  x 3  ..... ---(ii)
 dt t 0 h  

Here h = 0.1, t0 = 0.3, x0 = 0.31, x-1 = 0.77, 2x-1= -0.46 etc.

Putting these values in (i) and (ii), we get

 dx  = 1  0.31  0.77  1  0.01  0.02  1  0.02  0.27  

         .......
 dt  0.3 0 .1  2  6 2  30  2  

= 5.33

 d 2x  1  1 1 
 2    0.46  12 (0.01)  90 (0.29)  .... = -45.6
 
2
 dt 0.3 (0.1)

133
 6.3 MAXIMA AND MINIMA OF A TABULATED FUNCTIO N

Newton’s forward interpolation formula is

p( p  1) 2 p( p  1)( p  2) 3
y  y o  py o   yo   y o  ......
2! 3!

Differentiating both sides w.r.t. p, we get

dy 2p 1 2 3p2  6 p  2 3
 y o   yo   y o  ...... ---(1)
dp 2 2

For maxima or minima, dy/dp = 0. Hence equating the right hand side of (1) to zero and retaining
only upto third differences, we obtain

2p 1 2 (3 p 2  6 p  2) 3
y o   yo   yo  0 ,
2! 6

i.e. ( ½ 3y0)p2 + (2y0 – 3 y0) p + (y0 – ½ 2y0 + 1/3 3y0) = 0

Substituting the values of y0, 2y0, 3y0 from the difference table, we solve this quadratic for p.
Then the corresponding values of x are given by x = x0 + ph at which y is maximum or
minimum.

Example: From the table given below, find the value of x for which y is minimum. Also find
this value of y.

x: 3 4 5 6 7 8

y: 0.205 0.240 0.259 0.262 0.250 0.224

Solution:

The difference table is

134
x y  2 3

3 0.205

0.035

4 0.240 -0.016

0.019 0.000

5 0.259 -0.016

0.003 0.001

6 0.262 -0.015

-0.012 0.001

7 0.250 -0.014

-0.026

8 0.224

Taking x0 = 3, we have y0 = 0.205, y0 = 0.035, 2y0 = -0.016 and 3y0 =0.

 Newton’s forward difference formula gives

p ( p  1)
y = 0.205 + p(0.035) + (-0.016) ---(i)
2!

Differentiating it w.r.t. p, we have

dy 2 p 1
= 0.035 + (-0.016)
dp 2!

For y to be minimum, dy/dp = 0

135
 0.035 – 0.008 (2p – 1) = 0

Which gives p = 2.6875

 x = x0 + ph = 3 +2.6875 x 1 = 5.6875.

Hence y is minimum when x = 5.6875.

Putting p = 2.6875 in (i), the minimum value of y

= 0.205 + 2.6875 x 0.035 + ½ (2.6875 x 1.6875) (-0.016)

= 0.2628.

6.4 NUMERICAL INTEGRATION

The process of evaluating a definite integral from a set of tabulated values of the integrand f(x) is
called numerical integration. This process when applied to a function of a single variable, is
known as quadrature.

The problem of numerical integration, like that of numerical differentiation, is solved by

representing f(x) by an interpolation formula and then integrating it between the given limits. In
this way, we can derive quadrature formulae for approximate integration of a function defined by
a set of numerical values only.

Newton-Cotes quadrature formula:

b
Let I =  f ( x)dx ,
a

where f(x) takes the values yo, y1, y2, …..,yn for x = x0, x1, x2,….,xn.

Fig – 5.1
136
Let us divide the interval (a,b) into n sub- intervals of width h so that x0 = a, x1= x0+h, x2 = x0 +
2h, ……,xn = x0+nh = b. Then

x0 nh
I =  x0
f ( x)dx [Put x = x0 +rh, dx = hdr]

n
= h  f ( x0  rh)dr
0

n r (r  1) 2 r (r  1)( r  2) 3 
= h   yo  ry0   y0   y0  ...dr
0
 2! 3! 

[Approximated by Newton’s forward interpolation formula]

Integrating term by term, we obtain

x0 nh  n n(2n  3) 2 n(n  2)2 3

x0
f ( x)dx =nh  yo  y0 
 2 12
 y0 
24
 y0

 n 4 3n3 11n 2  4 y0 
     3n   ....... ---(1)
5 2 3  24 

This is known as Newton-Cotes quadrature formula. From this general formula, we deduce the
following important quadrature rules by taking n = 1, 2, 3,…..

6.4.1 Trapezoidal Rule

Putting n = 1 in (1) and taking the curve through (x0, y0) and (x1,y1) as a straight line, i.e., a
polynomial of first order so that differences of order higher than first become zero, we get

x0 h  1  h
x0
f ( x)dx = h y0  y0   ( y0  y1 )
 2  2

Similarly

x0  2 h  1  h
x0  h
f ( x)dx = h y1  y1   ( y1  y 2 )
 2  2

………………………………………
137
 h
 yn1  yn 
x0  nh
x0  ( n 1) h
f ( x)dx =
2

Adding these n integrals, we obtain

h
 y 0  y n   2( y1  y 2  ....  y n1 )
x0 nh
x0
f ( x)dx =
2
---(2)

This is known as the trapezoidal rule.

In this rule, the area of each strip (trapezium) is found separately. Then the area under the curve
and the ordinates at x0 and x0 +nhis approximately equal to the sum of the areas of the n
trapeziums.

6.4.2 Simpson’s one-third Rule

Putting n = 2 in (1) above and taking the curve through (x0, y0), (x1, y1) and (x2, y2) as a parabola,
i.e., a polynomial of second order so that differences of order higher than second vanish, we get

x0  2 h  1  h
x0
f ( x)dx = 2h y0  y0  2 y0   ( y0  4 y1  y 2 )
 6  3

Similarly

h
 y 2  4 y3  y 4 
x0  4 h
x0  2 h
f ( x)dx =
3

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

h
 yn2  4 yn1  yn  , n being even.
x0  nh
x0  ( n  2 ) h
f ( x)dx =
3

Adding all these integrals, we have when n is even

h
 y0  y n   4( y1  y3  ....  y n1 )  2( y 2  y 4  .....  y n2 )
x0 nh
x0
f ( x)dx =
3

---(3)

138
This is known as the Simpson’s one-third rule or simply Simpson’s rule and is most commonly
used.

Remark : While applying (3), the given interval must be divided into even number of equal sub-
intervals, since we find the area of two strips at a time.

6.4.3 Simpson’s three-eighth Rule

Putting n = 3 in (1) above and taking the curve through (xi, yi): i = 0, 1, 2, 3, as a polynomial of
third order so that differences above the third order vanish, we get

x 0  3h  3 3 1 
x0
f ( x)dx = 3h y0  y0  2 y0  2 y0 
 2 2 8 

3h
= ( y0  3 y1  3 y2  y3 )
8

Similarly,

x0 5 h 3h
x0 3 h
f ( x)dx =
8
( y3  3 y 4  3 y5  y6 ) and so on.

Adding all such expressions from x0 to x0 +nh, where n is a multiple of 3, we obtain

3h
 y 0  y n   3( y1  y 2  y 4  y 5  ..  y n1 )  2( y 3  y 6  ..  y n3 ) ,
x0 nh
x0
f ( x)dx =
8

---(4)

which is known as Simpson’s three-eighth rule.

In this rule, the number of sub-intervals should be taken as multiple of three.

6 dx
Example: Evaluate  0 1  x2
by using (i) Trapezoidal rule, (ii) Simplson’s 1/3 rule, (iii)

Simpson’s 3/8 rule, and compare the results with its actual value.

Solution:

139
 1
Divide the interval (0,6) into six parts each of width h = 1. The values of f(x) = are given
1 x2
below

x 0 1 2 3 4 5 6

f(x) 1 0.5 0.2 0.1 0.0588 0.0385 0.027

y y0 y1 y2 y3 y4 y5 y6

(i) By Trapezoidal rule,

6 dx h
 0 1 x 2
= [( y0  y6 ) +2(y1+y2+y3+y4+y5)]
2

1
= [(1  0.027)  2(0.5  0.2  0.1  0.0588  0.0385)]
2

= 1.4108

(ii) By Simpson’s 1/3 rule,

6 dx h
 1 x
0 2
=
3
[( y0  y6 ) +4(y1+y3+y5)+2(y2+y4)]

1
= [(1  0.027) +4(0.5+0.1+0.0385)+2(0.2+0.0588)]
3

= 1.3662.

(iii) By Simpson’s 3/8 rule,

6 dx 3h
 1 x
0 2
=
8
[( y0  y6 ) +3(y1+y2+y4+y5)+2y3)]

3
= [(1  0.027) +3(0.5+0.2+0.0588+0.0385)+2(0.1)]
8

140
 = 1.3571

6 dx
 1 x
6
Also 2
= tan 1 x  tan 1 6  1.4056
0 0

This shows that the value of the integral found by Simpson’s 1/3 rule is the nearest to the actual
value.

Example: The velocity v(km/min) of a moped which starts from rest, is given at fixed intervals
of time t (min) as follows :

t: 2 4 6 8 10 12 14 16 18 20

v: 10 18 25 29 32 20 11 5 2 0

Estimate approximately the distance covered in 20 minutes.

Solution:

If s (km) be the distance covered in t (min), then

ds
v
dt

20 h
s t 0   vdt  [ X  4.O  2.E ].
20
 (By Simpson’s rule)
0 3

Here h = 2, v0 = 0, v1 = 10, v2 = 18, v3 = 25 etc.

X = v0 + v10 = 0 + 0 = 0

O = v1 + v3 +v5 + v7 +v9 = 10 + 25 + 32 + 11+2 = 80

E = v2 + v4 + v6 + v8 = 18 + 29 + 20+ 5 = 72

Hence the required distance

2
= s t 0  (0  4  80  2  72)  309.33 km.
20

141
Example:A solid of revolution is formed by rotating about the x-axis, the area between the x-
axis, the lines x = 0 and x = 1 and a curve through the points with the following co-ordinates.

x: 0.00 0.25 0.50 0.75 1.00

y: 1.0000 0.9896 0.9589 0.9089 0.8415

Estimate the volume of the solid formed using Simpson’s rule.

Solution:

Here h = 0.25 y0 = 1, y1= 0.9896, y2 = 0.9589 etc.

 Required volume of the solid generated

1 h
 y dx   . [( y 0  y 4 )  4( y1  y3 )  2 y 2 ]
2 2 2 2 2 2
=
0 3

0.25
= [{1  (0.8415) 2 }  4{( 0.9896) 2  (0.9089) 2 }  2(0.9589) 2 ]
3

0.25  3.1416
= [1.7081  7.2216  1.839]
3

= 0.2618(10.7687) = 2.8192.

Example: Evaluate

1 1
I = dx , correct to three decimal places.
0 1 x

Solution:

We solve this example by both the trapezoidal and Simpson’s rules with h = 0.5, 0.25 and 0.125
respectively.

142
 1
(i) h = 0.5. The values of x and y  are tabulated below
1 x

x: 0 0.5 1.0

y: 1.0000 0.6667 0.5

(a) Trapezoidal rule gives

1
I= [1.0000  2(0.6667)  0.5]
4

= 0.7084

(b) Simpson’s rule gives

1
I= [1.0000  4(0.6667)  0.5]
6

= 0.6945

1
(ii) h = 0.25. The tabulated values of x and y  are given below
1 x

x: 0 0.25 0.50 0.75 1.00

y: 1.0000 0.8000 0.6667 0.5714 0.5

(a) Trapezoidal rule gives

1
I= [1.0  2(0.8000  0.6667  0.5714)  0.5]
8

= 0.6970

(b) Simpson’s rule gives

1
I= [1.0  4(0.8000  0.5714)  2(0.6667)  0.5]
12

= 0.6932

143
(iii) Finally, we take h = 0.125. The tabulated values of x and y are

x: 0 0.125 0.250 0.375 0.5 0.625 0.750 0.875 1.0

y: 1.0 0.8889 0.8000 0.7273 0.6667 0.6154 0.5714 0.5333 0.5

(a) Trapezoidal rule gives

1
I= [1.0  2(0.8889  0.8000  0.7273  0.6667
16

+ 0.6154  0.5714  0.5333)  0.5]

= 0.6941

(b) Simpson’s rule gives

1
I= [1.0  4(0.8889  0.7273  0.6154  0.5333)
24

 2(0.8000  0.6667  0.5714)  0.5]

= 0.6932

Hence the value of I may be taken to be equal to 0.693, correct to three decimal places. The exact
value of Iis loge2, which is equal to 0.693147…. This example demonstrates that, in general,
Simpson’s rule yields more accurate results than the trapezoidal rule.

6.5 ERRORS IN QUADRATURE FORMULAE

The error in the quadrature formulae is given by

b b
E =  ydx  P( x)dx ,
a a

where P(x) is the polynomial representing the function y = f(x), in the interval [a, b].

144
1) Error in the Trapezoidal rule
Expanding y = f(x) around x = x0 by Taylor’s series, we get

( x  x0 ) 2
y  y 0  ( x  x 0 ) y 0 ' y 0 "..... ---(1)
2!

x0  h x0  h ( x  x0 ) 2
  x0
ydx  
x0
[ y0  ( x  x0 ) y0 '
2!
y0 "......]dx

h2 h3
= y0 h  y0 ' y0 "..... ---(2)
2! 3!

Also A1 = area of the first trapezium in the interval [x0, x1]

1
 h( y0  y1 ) ---(3)
2

Putting x = x0+h and y =y1 in (1), we get

h2
y1  y 0  hy 0 ' y 0 ".....
2!

Substituting this value of y1 in (3), we get

1  h2 
A1  h  y0  y0  hy0 ' y0 ".....
2  2! 

h2 h3
= hy 0  y 0 ' y 0 "..... ---(4)
2 2.2!

 Error in the interval [x0, x1]

x1
=  x0
ydx  A1

1 1  3 h3
=  h y0 "...   y0 "..... ,
 3! 2.2!  12

145
 h3
i.e., Principal part of the error in [x0, x1] = - y0 "
12

h3
Similarly principal part of the error in [x1, x2]= - y1 " and so on.
12

h3
Hence the total error E = - [y0+ y1+……+yn-1]
12

Assuming that y (X) is the largest of the n quantities y0y1,..., yn-1, we obtain

nh 3 (b  a)h 2
E y" ( X )   y" ( X ) , [ nh  b  a ] --- (5)
12 12

Hence the error in the trapezoidal rule is of the order h2.

2) Error in Simpson’s one-third rule

Expanding y = f(x) around x = x0 by Taylor’s series, we get (1).

 Over the first double strip, we get

x2 x0  2 h  ( x  x0 ) 2 
 x0
ydx  
x0
 0

( y  ( x  x 0 ) y 0 '
2!
y0 "...... dx


4h 2 8h 3 16h 4 32h 5 iv
= 2hy0  y0 ' y0 " y0 "' y0  ...... --- (6)
2! 3! 4! 5!

Also A1 = area over the first double strip by Simpson’s 1/3 rule

1
= h( y0  4 y1  y2 ) --- (7)
3

Putting x = x0 + h and y = y1 in (1), we get

h2 h3
y1= y0 + hy0 + y 0 ' ' y0 ' ' '...........
2! 3!

Again putting x = x0 + 2h and y = y2 in (1), we have

146
 4h 2 8h 3
y 2  y 0  2hy 0 ' y 0 ' ' y 0 ' ".....
2! 3!

Substituting these values of y1 and y2 in (7), we get

h   h2 
A1 =  0
y  4
 y 0  hy 0 ' y 0 " '...... 
3   2! 

 4h 2 8h 3 
+  y 0  2hy 0 ' y 0 " y 0 " '......
 2! 3! 

4h 3 2h 2 5h 5 iv
= 2hy0 + 2h y0 +
2
y0 " y0 " ' y0  ..... ---(8)
3 3 18

x2
 Error in the interval [x0, x2] =  x0
y dx  A1

4 5
=   h 5 y 0 ,
iv
[(6)-(8)],
 15 18 

i.e., Principal part of the error in [x0, x2]

4 5 h 5 iv
=   h 5 y 0 
iv
y0
 15 18  90

h 5 iv
Similarly principal part of the error in [x2, x4] =  y2 and so on.
90

h5
Hence the total error E =  [ y 0  y 2  .....  y iv 2 ( n 1) ]
iv iv

90

Assuming the yiv (X) is the largest of y0iv, y2iv, …., yiv2n-2, we get

nh 5 iv (b  a)h 4 iv
E < y0 ( X )   y (X ) , [ 2nh = b-a], ---(9)
90 180

1
i.e., the error in Simpson’s rule is of the order h4.
3

147
(3) Error in Simpson’s 3/8 rule. Proceeding as above, here the principal part of the error in the
interval [x0, x3]

3h 5 iv
=  y --- (10)
80

6.6 CHECK YOUR PROGRESS

Find the first and second derivatives of f(x) at x = 1.5 if

x: 1.5 2.0 2.5 3.0 3.5 4.0

f(x): 3.375 7.00 13.625 24.00 38.875 59.00

1. Find the first and second derivatives of the function tabulated below at x = 1.1,
x: 1.00 1.2 1.4 1.6 1.8 2.0

f(x): 0.0 0.128 0.544 1.296 2.432 4.00

2. Find the first derivative at x = 4 from the following values of x and y,

x: 1 2 4 8 10

f(x): 0 1 5 21 27


3. Calculate the value of  2 sin x dx by Simpson’s one-third rule using 10 sub-intervals.
0

4. A curve is drawn to pass through the points given in the following table:
x: 1 1.5 2 2.5 3 3.5 4

y: 2 2.4 2.7 2.8 3 2.6 2.1

Find the area bounded by the curve, x-axis and the lines x = 1, x = 4.

148
 6.7 SUMMAERY

NUMERICAL DIFFERENTIATION:- Methods based on interpolation uses the ponomial

approximation obtained by nterpolation to find the derivative of the function, which is known at
discrete points in the interval [a, b].Derivatives are find by using Forward difference Formula
,Backward difference formula, central difference formula.
NUMERICAL INTEGRATION: -Numerical Integration is nothing but finding an approximate
b
value to I =  f ( x)dx .
a

There are two different strategies to develop numerical integration formulae. One is similar to
what we have adopted to numerical differentiation. That is, we approximate a polynomial for the
given function and integrate that polynomial within the limits of the integration. This restricts us
to integrate a function known at discrete tabular points. If these points are uniformly spaced then
the corresponding integration formulae are called as Newton - Cotes formulae for numerical
integration. On the other hand if we know the function explicitly but could not integrate in the
usual means because of the nature of the function then we can use the concept called quadrature
rule to find an approximate value to the integration.
In this method the function is evalualted at some predetermined abscissa (nodal) points and then
these values are added after multiplying them with some weights which are again
predeternmined, to find an approximate value to the given integral.

149
 6.8 KEYWORDS

1. Derivatives using Forward Difference Formula:

 dy  1 1 2 1 3 1 4 
   y o   y 0   y 0   y 0  .........
  x0
dx h  2 3 4 

d2y 1  2 11 4 
 2   2  y 0   y 0  12  y 0  .....
3

 dx  x0 h  

2. Derivatives using Backward Differences Formula:-

 dy  1 1 2 1 3 1 4 
   y n   y n   y n   y n ..........
 dx  xn h  2 3 4 

d2y 1  11 
 2  = 2  2 y n   3 y n   4 y n ....
 dx  xn h  12 

3. Derivatives using Central Difference Formulae:

 dy  = 1  y 0  y 1 1  y 1   y  2 1  y  2   y 3 
3 3 5 5

      ........
 dx  x0 h  2 6 2 30 2 

d2y 1  2 1 4 1 6 
 2   2  y 1  12  y  2  90  y 3 .....
 dx  x0 h  

h
 y 0  y n   2( y1  y 2  ....  y n1 )
x0 nh
4. Trapezoidal Rule:- x0
f ( x)dx =
2

5. Simpson’s one-third Rule:-

h
 y0  y n   4( y1  y3  ....  y n1 )  2( y 2  y 4  .....  y n2 )
x0 nh
x0
f ( x)dx =
3

150
 5. Simpson’s three-eighth Rule

3h
 y 0  y n   3( y1  y 2  y 4  y 5  ..  y n1 )  2( y 3  y 6  ..  y n3 )
x0 nh
x0
f ( x)dx =
8

6.9 SELF ASSESSMENT TEST

1. From the following table of values of X and Y, obtain first and second Derivatives for
X=1.2

x 1.0 1.2 1.4 1.6 1.8 2.0 2.2

y 2.7183 3.3201 4.0552 4.9530 6.0496 7.3891 9.0250

2. From the following table of values of X and Y, obtain first and second Derivatives for
X=1.4

x 1.4 1.6 1.8 2.0 2.2

y 4.0552 4.9530 6.0496 7.3891 9.0250

3. From the following table of values of X and Y, obtain first and second Derivatives for
X=2.2

x 1.0 1.2 1.4 1.6 1.8 2.0 2.2

y 2.7183 3.3201 4.0552 4.9530 6.0496 7.3891 9.0250

4. From the following table of values of X and Y, obtain first and second Derivatives for
X=1.4

x 1.4 1.6 1.8 2.0 2.2

y 4.0552 4.9530 6.0496 7.3891 9.0250

151
5. From the following table of values of X and Y, obtain first and second Derivatives for
X=5

x 2 4 9 10

y 4 56 711 980

6. From the following table, find the area bounded by curve and x axis from X=7.47 to
X=7.52 using trapezoidal, simpson’s 1/3, simpson’s 3/8 rule.

x 7.47 7.48 7.49 7.50 7.51 7.52

f(x) 1.93 1.95 1.98 2.01 2.03 2.06

6.10 ANSWER TO CHECK YOUR PROGRESS

(1) 4.75;9

(2) 0.63;6.6

(3) 2.8326

(4) 0.9985

(5) 7.78

152
 6.11 REFERENCES/ SUGGESTEDREADINGS

Computer Oriented Numerical Methods, V. Rajaraman, PHI.

1. Introductory Methods of Numerical Analysis, S.S. Sastry, PHI.

2. Numerical Methods for Scientific and Engineering Computation, M.K.Jain,
S.R.Lyenger, R.K.Jain, Wiley Eastern Limited.
3. Numerical Methods with Programs in C, T. Veerarajan, T. Ramachandarn., Tata
McGraw Hill.

153