MANONMANIAM SUNDARANAR UNIVERSITY

DIRECTORATE OF DISTANCE & CONTINUING EDUCATION

TIRUNELVELI 627012, TAMIL NADU

B.Sc. MATHEMATICS - II YEAR

DJM2C - NUMERICAL METHODS

(From the academic year 2016-17)

Most Student friendly University - Strive to Study and Learn to Excel

For more information visit: http://www.msuniv.ac.in

 NUMERICAL METHODS
Unit I : Finite differences – difference table – operators E,Δ and - Relations between these
operations – Factorial notation – Expressing a given polynomial in factorial notation –
Difference equation – Linear difference equations – Homogeneous linear difference equation
with constant coefficients.
Unit II: Interpolation using finite differences – Newton – Gregory formula for forward
interpolation – Divided differences – Properties – Newton’s formula for unequal intervals -
Lagrange’s formula – Relation between ordinary differences and divided differences –
inverse interpolation.
Unit III: Numerical differentiation and integration – General Quadrature formula for
equidistant ordinates – Trapezoidal Rule – Simpson’s one third rule – Simpson’s three eight
rule – Waddle’s rule – Cote’s method.
Unit IV: Numerical solution of ordinary differential equations of first and second orders –
Piccards method. Eulers method and modified Euleis method – Taylor’s series method –
Milne’s method – Runge – Kutta method of order 2 and 4 – Solution of algebraic and
transcendent equations. Finding the initial approximate value of the root – Iteration method –
Newton Raphson’s method.
Unit V: Simultaneous linear algebraic equations – Different methods of obtaining the
solution – The elimination method by Gauss – Jordan method – Grouts’ method – Method of
factorization .
Books:
Calculus of finite differences and Numerical Analysis, P.P. Gupta & G.S. Malik, Krishna
PrakashamMardin, Mecrutt.

1
 UNIT I
FINITE DIFFERENCES
1.1 Introduction
1.2 Difference Operations
1.3 Factorial Function
1.4 Difference Equations
1.5 Linear Difference Equations

1.1 Introduction
We introduce the idea of finite differences and associated concepts, which have
important applications in numerical analysis.
For example, interpolation formulae are based on finite differences. Through finite
differences, we study the relations that exist between the values that are assumed by the
functions whenever the independent variables change by finite jumps.
1.2 Difference Operations
There are three difference operators namely forward, backward and central difference
operators.
Forward Difference Operator
Consider the function y = f(x). Suppose we are given a table of values of the function
at the points
x0 , x1  x0  h, x2  x0  2h, , xn  x0  nh .

Let f x0   y0 , f x1   y1 , f xn   y n .

We define

 f x   f x  h  f x 

Thus y0  f x0  h  f x   f x1   f x0   y1  y0 .

 y0  y1  y0 .

Further, x0 , x1 , , xn are called arguments. The corresponding values of f(x) are

called entries and h is called the interval of differencing.

Similarly,  y1  y 2  y1

  

 y n1  y n  y n1

2
 ∆ is called the forward difference operator and  y0 ,  y1 ,,  y n1 are called the first
forward differences of the function y = f(x).
The second order differences of the function are defined by

2 y0   y1   y0

2 y1   y 2   y1

  

2 y n1   y n   y n1 .

Similarly, higher order differences can be defined. In general the nth order differences
are defined by the equations

n yi  n1 yi 1  n1 yi .

These differences of the function y = f(x) can be systematically represented in the

form of a table called forward difference table. We can construct the difference table for any
number arguments and a sample difference table is given for six consecutive arguments.

x y = f(x) ∆y ∆2 y ∆3 y ∆4 y ∆5 y

x0 y0
[
 y0 [

x1  x0  h y1 2 y 0
[
 y1 3 y 0
x2  x0  2h y2
[

 y1
2

4 y 0
 y2
[[

3 y1
2 y 2 5 y 0
x3  x0  3h y3
4 y1
 y3
[

 y2
3

x4  x0  4h y4  y3
2

 y4
x5  x0  5h y5
[[

In the above forward difference table y 0 is known as the first entry and
 y0 ,  y0 ,,  y0 are called leading differences.
2 5

Note:
Since each higher order difference is defined in terms of the previous lower
differences by continuous substitution of each higher order difference can be expressed in
terms of the values of the function.
Thus
3
2 y0   y1   y0   y 2  y1    y1  y0 

 y 2  2 y1  y0

3 y0  2 y1  2 y0

  y3  2 y 2  y1    y 2  2 y1  y0 

 y3  3 y 2  3 y1  y0

4 y0  3 y1  3 y0

  y 4  3 y3  3 y 2  y1    y3  3 y2  3 y1  y0 

 y 4  4 y3  6 y 2  4 y1  y0 .

We observed that the coefficients occurring in the right hand side are simply the
binomial coefficients in (1-x)n. Hence in general, we have

n y0  y n  nc1 y n1  nc2 y n2     1 y0 .

n

Properties of the Operator ∆:

1. ∆ is linear, that is  a f x   b g x   a   f x   b  g x  where a and b are

constants.

Proof:

 a f x  b g x  a f x  h  b g x  h  a f x   b g x

 a  f x  h  f x  b g x  h  g x
  a f x   b gx   a  f x   b  gx .

2. m n  f x   m n  f x .

Proof:

m n  f x    m times  n times  f x  g

  m  n times  f x 
   f x   m n  f x .
m n

3.  f x g x  f x  h g x  g x    f x.

Proof:

 f x g x  f x  h g x  h  f x g x

  f x  h g x  h  f x  h g x   f x  h g x  f x g x
4
  f x  hg x  h  g x  g x f x  h  f x
 f x g x  f x  h g x  g x    f x.

 f x  g x    f x   f x   g x 
4.   .
 g x   g x  h  g x 

Proof:

 f x  f x  h  f x 
  
 g x   g x  h  g x 
f x  h  g x   f x  g x  h 

g x  g x  h 
f x  h  g x   f x  g x   f x  g x   f x  g x  h 

g x  g x  h 
g x  f x  h   f x   f x g x  h   g x 

g x  g x  h 
 f x  g x    f x   f x   g x 
  .
 g x   g x  h  g x 

Backward Differences
Consider the function y = f(x). Suppose we are given a table of values of the function
at the points
x0 , x1  x0  h, x2  x0  2h, , xn  x0  nh .

Let f x0   y0 , f x1   y1 , f x2   y 2 , , f xn   y n .

We define

  f x  f x  f x  h

Thus  y1  y1  y0

 y 2  y 2  y1

  

 y n  y n  y n1 .

 is called the backward difference operator and  y1 ,  y 2 ,,  y n are called the
first order backward differences of the function y = f(x).

5
 The second order differences of the function y = f(x) are defined by

 2 y 2   y 2   y1

 2 y 2   y3   y 2

  

 2 y n   y n   y n1 .

Similarly, higher order differences can be defined. In general, the nth order differences
are given by,

 n yi   n1 yi   n1 yi 1 .

These differences of the function y= f(x) can be systematically represented in the

form of a table called backward differences table.
x y  f x  y 2 y 3 y 4 y 5 y
[

x0 y0
 y0
[

x1  x0  h
[

 2 y0
y1
[[
 y1  3 y0
x2  x0  2h y2  2 y1  4 y0
 y2
[[

[  y1
3
 5 y0
x3  x0  3h y3  y2
2
[[

 4 y1
 y3 [

 3 y2
x4  x0  4h y4  2 y3
 y4

x5  x0  5h y5

As in the case of ∆ we can prove that  is also linear.

Remark:

The relation between the two difference operators is given by   f x  h   f x  .

For,   f x  h  f x  h  f x    f x 

Similarly,

 2  f x  2h    f x  2h  f x  h

  f x  2h  f x  h sin ce  is linear 

  f x  h    f x 

   f x  h  f x 

6
 2  f x  2h  2 f x .

In general,

  f x  nh  2 f x .

Hence from the forward difference table of the function f(x) we can obtain backward
differences of all orders.
Central Difference Operator
Sometimes it is convenient to employ another system of differences known as central
differences. We define central difference operator δ as

 h  h
 f x   f  x    f  x   .
 2  2

Thus if f xi   yi then we have

 y 1  y1  y 0
2

 y 3  y 2  y1
2

 

 y n 1  y n  y n1
2

Here the subscript of δy is the average of the subscripts of the two members of the
difference. The higher order differences can define similar to forward and backward
differences.

 2 y1   y 3   y 1
2 2

 2 y2   y 5   y 3
2 2

 3 y 3   2 y 2   2 y1 etc.
2

These differences of the function y = f(x) can be systematically represented in the

form of a table called central difference table.
x y  f x   2 y 3 y 4 y
[[

x0 y0
 y1
[

x1  x0  h y1  2 y1 [

 y3  3 y3
x2  x0  2h y2 2 2

[[ [
 2 y2 [

7
 x3  x0  3h y3  y5  4 y2
2
 2 y3  3 y5
[

x4  x0  4h y4 2
[
 y7
2

Fundamental theorem for Finite Differences

Let f(x) be a polynomial of degree n. Then the nth difference of f(x) is a constant and
all higher order differences are zero.

1.3 Factorial Function (or) Factorial Notation

Definition
A product of the form x(x-h) (x-2h) … [x-(n-1) h] is called a factorial function and is
denoted by x(n).

 x n   xx  hx  2hx  n  1h.

Thus x(1) = x, x(2) = x(x-1) and x(3) = x(x-1) (x-2).

We observed that x(n) is a polynomial of degree n with leading coefficient 1.
The following theorem shows that the formula for the first difference of x (n) is
obtained by the simple rule of differentiation.
Theorem:

 x n   nh x n1 . In particular when h = 1,  x n   n x n1 .

Proof:

 x n   x  h 
n 
 x n 

 x  h x x  hx  n  2h  xx  hx  2hx  n  1h

 x x  hx  2hx  n  2h x  h  x  n  1h

 x n   x n1 nh.

When h = 1,

 x n   n x n1 . (1)

Note: 1
From equation (1) we get the formula for first order difference, which is obtained by
the simple differentiation rule.
8
  
For example, 2 x n    nhx n1  n h  x n1  n n  1 h 2 x n2  proceeding like this
we get

n x n    nn  1n  21h n x 0  n!h n .

Note: 2
Any polynomial f(x) of degree n can be expressed in the form
f x   c0 x n   c1 x n1    cn1 x 1  cn . If f(x) is divided successively by x-0, x-1, x-2, …,
x-(n-1), then the remainders give the coefficients cn , cn1 ,, c1 , c0 .

Definition
The reciprocal factorial function x(n) for positive integer n is defined as
1
x  n  
x  hx  2h... x  nh
As in the case of factorial function the formula for first order difference of x (-n) is
similar to differentiation rule when h = 1.
Theorem:

 x  n    nh x n1 . In particular  x  n   n x n1 .

Proof:

 x  n   x  h 
 n 
 x  n 

1 1
 
x  2hx  3hx  n  1h x  hx  2hx  nh


1
x  h  x  n  1h
x  hx  2hx  n  1h
 nh

x  hx  2hx  n  1h
 x  n    nh x n1 .

Remark:


2 x  n     nhx n1 
 nh  n  1hx n2    1 h 2 n n  1 x n2 
2

In general r x  n    1 h r nn  1n  r  1 x n1 .

r

Example: 1
Form the forward difference table for the following data:

9
 x 0 1 2 3 4

y 8 11 9 15 6

Solution
The difference table is given below.
x y ∆y ∆2y ∆3y ∆4y
0 8
3
1 11
-5
-2
2 9 13
8 -36
6
3 15 [

-23
-15
9
4 6

Example: 2
Find ∆ (2x).
Solution
∆(2x) = 2x+h-2x
∆(2x) = 2x (2h – 1).

Example: 3
Find the nth difference of ex.
Solution

e x   e xh  e x  e x e h  1

   
2 e x   e x

  e e  1
x h

 e  1 e 
h x

 e  1e e  1
h x h


 e x eh 1 
2

  
Similarly 3 e x  e x e h  1 . 
3

 
Proceeding like this we get n e x  e x e h  1 .  
n

Other Difference Operators

10
 The shift operator E and averaging operator .
Definition:
The shift operator E is defined by Ef(x) = f(x+h). Hence E2 f(x) = Ef(x+h) = f(x+2h).
In general, for any positive integer n
En f(x)=f(x+nh)
In particular we have
Ey0 = y1
E2 y0 = y2
.............
En y0 = yn.
The inverse operator E-1 is defined as
E-1f(x) =f(x-h)
For any real number n we have
En f(x) =f(x+nh).
Note:
Em En f(x) = Em+n f(x)
Definition:
The averaging operator is defined as

 h  h
fx    fx  
f ( x )  
2  2
2
There are several relation connecting the operators , , , E,  and the differentiation
operator D.
Example 1:
E = E = .
Solution:
E = E (1 - E-1) E = E – 1 = .
Also, E = (1 – E-1) E = E – 1 = .
Example 2:

 1  
1 1
 E 2  E 2  1   2  2   .
 
 
Solution:

11
 1  
1 1  1  
1 1
 E 2  E 2  1   2   E 2  E 2  E 2
   
   
=E+1
= (1 + ) + 1

 1  
1 1
 E 2  E 2  1   2  2   .
 
 
1.4 Difference Equations
Any situation in which there exists a sequential relation at discrete values of the
independent variable, leads to difference equations. Difference equations may be thought of a
discrete counterpart of differential equations and there is a striking similarity between the
methods of solving difference equations and differential equations.
Definition
An equation involving the differences of an unknown function y = y(x) at one or more
general values of the argument n is called a difference equation.
The following are some examples of difference equations.

 yn  2 yn  n (1)

2 y n  5 y n  3 y n  0 (2)

2 u x  4u x  4u x  3 x (3)

We assume that the consecutive values of independent variable differ by unity. With
this assumption a difference equation can be written in an alternative form, as illustrated
below.
Consider the difference equation,  y n  2 y n  n .

Since  y n  E  1y n  y n1  y n the above difference equation can be written as

y n1  yn  2 y n  n (i.e) yn1  yn  n .

Consider the difference equation (3), 2 u x  4u x  4u x  3 x .

Since u x  E  1u x  u 1  u x and 2 u x  E  1 u x  E 2  2E  1u x  u x2  2u x1  u x

2

the difference equation (3) can be written as

u x2  2u x1  u x   4u x1  u x   4u x  3 x .

(i.e) u x2  6u x1  9u x  3 x .

Definition

12
 The order of a difference equation is the difference between the largest and smallest
subscripts occurring in it, when the equation is expressed in the form free of the difference
operator ∆.
The degree of a difference equation, expressed in the form free of ∆, is the higher
power of y.
Examples:
1. The difference equation y n  2 y n  n , when expressed in a form free of ∆ is
y n1  y n  n . The largest subscript in the equation is n+1 and the smallest subscript
in the equation is n. Hence the order of the difference equation is 1. The highest
power of y is 1. Hence the degree of the difference equation is also 1.
2. The difference equation 2 u x  4u x  4u x  3 x , when expressed in a form free of ∆
is u x2  6u x1  9u x  3 x . The largest subscript in the equation is x+2 and the smallest
subscript in the equation is x. Hence the order of the difference equation is 2. The
highest power of y is 1. Hence the degree of the difference equation is 1.
3. Consider the difference equation 4 y n23  2 y n y n1  y n2 y n41  0 . It is free form ∆. The
largest subscript here is n+3 and the smallest subscript is n. Hence the order of the
difference equation is 3. The highest power of y is 4 hence the degree of the
difference equation is 4.

Note:
The order of the difference equation may not be the highest power of ∆ involved in it.
For example, consider the equation 2 y n  2y n  y n  2 n . This can be written as

E  12 yn  2E  1yn  yn  2 n .

i.e. E 2 y n  2 n

i.e. y n 2  2 n

which is not even a difference equation.

Definitions
Solution
A solution of a difference equation is an expression for yn which satisfies the given
difference equation.
General Solution
A solution in which the number of arbitrary constants is equal to the order of the
difference equation is called the general solution.
13
Particular Solution
Any solution which is obtained from the general solution by giving particular values
to the arbitrary constants is called a particular solution.
Example: 1

Write the difference equation 3 y x  2 y x  y x  y x  0 in the subscript notation.

Solution
The given difference equation can be written as

E  13 y x  E  12 y x  E  1y x  y x  0

(i.e) E 3  3E 2  3E  1y x  E 2  2E  1y x  E 2  2E  1y x  E  1y x  y x  0

(i.e) E 3  2E 2  4E y x  0

(i.e) y x3  2 y x2  4 y x1  0 .

Example: 2

Find the order of the difference equation 3 y n  32 y n  2y n  y n  cos x .

Solution
The given difference equation can be written as

E  13 yn  3E  12 yn  2E  1yn  yn  cos x

i.e. E 3  3E 2  3E  1y n  3E 2  2E  1y n  2E  1y n  y n  cos x

i.e. E 3  6E 2  11E  5y n  cos x

 yn3  6 yn2  11y n1  5 yn  cos x

This difference equation is free from ∆.

 The order of the given difference equation is (n+3) - n = 3.
Example: 3
2
Show that yn  1  is a solution of the difference equation
n
n  1yn1  n yn  2n  3 .
Solution

n  1 y n1  n y n  n  11  2   2
  n1  
 n 1  n 

=n+1–2+n–2

14
n  1 yn1  n yn  2n  3 .
2
 yn  1  is a solution of the given difference equation.
n
Formation of Difference Equations
In the case of differential equations difference equation can be formed by eliminating
the constants from the given equation. We can see some examples.
Example: 1

Form the difference equation by eliminating the constant a from y n  a 3n .

Solution

y n  a 3n .

 
 y n1  a3n1  3 a3n  3y n .

 y n 1  3y n  0 is the required difference equation.

Example: 2
Form the difference equation by eliminating the arbitrary constants A and B from the
equation y n  A a n  B b n where a  b.

Solution

yn  A a n  B bn (1)

 y n 1  A a n 1  B b n1

  
(i.e) y n1  a A a n  b B b n  (2)

Similarly

  
y n 2  a 2 A a n  b 2 B b n  (3)

Eliminating Aan and Bbn from equation (1) and (3) we get

yn 1 1
y n 1 a b  0.
y n 2 a2 b2

   
(i.e) y n ab 2  a 2 b  y n1 b 2  a 2  y n2 b  a   0

 y n2 b  a   y n1 b  a b  a   y n abb  a   0

 y n2  a  by n1  aby n  0 (Since a  b).

1.5 Linear Difference Equations

15
 The difference equation of the form
a 0 y nr  a1y nr 1    a r y n  f n  (1)

where a 0 , a1 ,, a r and f(n) are functions of n is called a linear difference equation.

If a 0 , a1 ,, a r are constants then equation (1) is called a linear difference equation
with constant coefficients.
Equation (1) can also be written in the form

a E 0
r

 a1E r 1    a r 1E  a r y n  f n  .

In this section, deal with linear difference equation with constant coefficients and
discuss the methods of solving them. The methods are analogous to the methods of linear
differential equations with constant coefficients.
Definition
Consider the difference equation

a E
0
r

 a1E r 1    a r 1E  a r y n  f n  (1)

Let 1 n , 2 n , , r n  be r independent solutions of

a E
0
r
 a1E r1    a r1E  a r y n  0 (2)

Then U n  c11 n   c 22 n     c r r n  is the general solution of equation (2) and is called
the complementary function (C.F) of equation (1).
If Vn is a particular solution of equation (1) then y n  U n  Vn is the complete
solution of equation (1). Vn is called the particular integral (P.I) of equation (1).
Thus the complete solution of equation (1) is given by
y n  C.F  P.I .

Rules for finding Complementary Function

I. We first consider a linear difference equation of order one given by
y n 1  ay n  0 where a is a constant.

y n 1 y n
Dividing by an+1 we get,  0
a n 1 a n

y 
i.e.  nn 0.
a 
yn
  c where c is a constant.
an

 y n  ca n  0 is the solution of the given difference equation.

16
 II. Consider a linear difference equation of second order given by
 
E 2  aE  b y n  0 (1)
where a and b are constants.

Then the equation E2 + aE + b = 0 is called the auxiliary equation of equation (1). Let
1 , 2 be the roots of the auxiliary equation.

Case (i):
1 and 2 are real and distinct.

Equation (1) can be written as (E - 1) (E - 2) yn = 0.

We can derive two independent solutions of equation (1) by solving the equations
(E - 1)yn = 0 and (E - 2)yn = 0.

By, I the solution of these equations are y n  c11n and y n  c 2  n2 where c1, c2 are
arbitrary constants. Hence
y n  c11n  c 2  n2
is the general solution of equation (1).

Case (ii):
The roots are real and equal.

Equation (1) takes the form (E - 1)2 yn = 0.

Let y n  1n z n .
 E  1  1n z n  0 .
2

 
(i.e) E 2  2E1  1n 1n z n  0
1n2 z n2  21n 2 z n1  1n2 z n  0
z n2  2 z n1  z n  0
E 2
 2E  1  0
(i.e) E  12 z n  0
 2 z n  0 .
 z n  c1  c 2 n where c1, c2 are arbitrary constants. Hence
y n  c1  c 2 n 1n is the general solution.

Case (iii):
The roots are imaginary.

Let the roots be  + i and  - i.

 y n  c1   i  c 2   i
n n

17
  c1 rcos   i sin   c 2 rcos   i sin  [putting  = r cos  and  = r sin ]
n n

 r n c1 cos n  i sin n  c 2 cos n  i sin n

 r n A1 cos n  A 2 sin n

where A1, A2 are arbitrary constants and r   2   2 and   tan 1   .

 y n  r n A1 cos n  A 2 sin n is the general solution where r and  are given above.

III. Working Rule to find the C. F of the Difference Equation

a E
0
r

 a1E r 1    a r 1E  a r y n  f n 

From the auxiliary equation

a 0 E r  a1E r1    a r1E  a r  0

Let 1, 2, ..., r be its roots.

If 1, 2, ..., r are all distinct then the C. F is

c11n  c 2  n2    c r  nr .

If 1 = 2, the corresponding part of the C.F is c1  c 2 n 1n .


If 1 = 2 = 3, the corresponding part of the C.F is c1  c 2 n  c3n 2 1n . 
If there is a part of complex roots  + i,  - i the corresponding part of the C.F is

r n c1 cos n  c 2 sin n where r   2   2 and   tan 1   .

Rules for finding P.I:


Consider the difference equation a 0 E r  a1E r 1    a r 1E  a r y n  f n  . 
(i.e) Ey n  f n  where E  a 0 E r  a1E r1    a r1E  a r .

Type I:
F(n) = an where a is a constant.


E a n  a 0 E r  a1E r1    a r1E  a r a n 
 a 0 a nr  a1 a nr 1    a r a n


 a 0 a r  a1 a r1    a r a n 
 a a n

18
 an
if a  0 .
1 n
 P.I  a 
E  a 

Suppose (a) = 0.
Then (a) = (E – a)  (E).
1 n
Let a  bn
Ea

E  a b n  a n

b n 1  a b n  a n

b n 1 b n 1
  
a n 1 a n a

b  1
   nn  
a  a
bn n
 
an a

 b n  n a n 1

1
 a n  n a n 1 .
Ea
By similar argument, we have
1 n n  1 n 2
an  a
E  a  2
2!

1 n n  1n  2 n 3
an  a
E  a  3
3!

In general,

1 n r  a n r
an  .
E  a r r!

Type II:
F(n) is a polynomial in n.

P.I 
1
f n 
E 


1
f n 
1   

P.I  1    f n  .
1

19
We expand  1   1 in ascending powers of  and operate on f(n).

Type III:
f(n) = cos kn or sin kn.
1 ikn
Since coskn = real part of eitn and sinkn = imaginary part of eikn we can compute e by
E 
using the formula given in Type I.
Equating real and imaginary part we get the required P.I.
Type IV:

f n   a n gn .

P.I 
1
E 
f n  
1
E 

a n gn  
Now,

   
 E  a n gn   a 0 E r  a1 E r 1    a r1 E  a r a n gn  
 a 0 a nr gn  r   a1 a nr1 gn  r  1    a r a n gn 

 
 a 0 a r E r gn   a1 a r 1 E r 1 gn     a r gn  a n

 a n  a E gn 


1
a E 

a n gn  
1
 E 
 
a n gn  
1

 E 
f n 

gn  .
1
 P.I  a n
 a E 

Example: 1
Solve yn+1 – 2yn-1 + 2yn-1 = 0.
Solution
The given equation can be written as
E2yn+1 – 2Eyn-1 + 2yn-1 = 0

 
 E 2  2E  2 y n1  0

The auxiliary equation is E2 – 2E + 2 = 0.

2 4
The roots are  1  i    i say 
2

The complete integral is y n1  r n1 A cos n  1  Bsin n  1 where

20
  
r   2   2  2 and   tan 1    tan 1 1 
 4

 n  1   n  1 
 y n 1   2
n 1 
A cos 4   B sin 4  .
    

Example: 2

 
Solve 4E 2  4E  1 y n  2 n  2n .

Solution
The auxiliary equation is 4E2 – 4E + 1 = 0.
1 1
(i.e) (2E – 1)2 = 0. Hence the roots are , .
2 2
n
1
The C.F  A  B n  
2


Particular integral P.I   2
1  n
 2 2
n
 
 4E  4E  1 


Now  2
1  n
2   
2n

2n
.
 4E  4E  1  4  22  4  2  1 9

 
n
 1  n  1  1 
Also  2  2  2  
 4E  4 E  1   4E  4E  1   2 

 1   1 n

 
2  
 2 E  1  2 
n 2
n (n  1) 1
   .
2! 2

Hence the complete solution is yn = C.F + P.I.

n n  1  1 
n n 2
1 2
n
y n  A  B n       .
2 9 2 2

Exercises:
1. Construct the difference table for the following data.

x 45 50 55 60 65 70 75

y 80 95 120 100 85 70 60

2. Form the difference table for the following data:

21
 x 0 1 2 3 4 5

y 1 5 19 55 125 241

and hence find 5 y 0 .

3. Form the difference table for the following data:

x 0 1 2 3 4 5 6

y 2 5 8 20 30 10 3

and hence find 6 y 0 .

4. Form the difference table of the function y = x3 + x2 – 2x + 1 for x = -1, 0, 1, 2, 3, 4.
5. Form the difference table for the function y = x3 + x + 1 at x = 0, 1, 2, 3, 4 and hence
find (i) 2 y 2 (ii) 3 y1 (iii) 4 y 0 .
6. Find the order and degree of the following difference equations.
(i) y n1  3 y n  3n (ii) y n2  y n1  y n  0 (iii) y n  y n1  6 y n2  0 .
7. Find the order and degree of the following difference equations.
(i) E 2 y n  3Ey n 1  y n  n 2 (ii) y n3  3y n 2  y n  n 2 2 n .
8. Show that y n  A  Bn 2 n is a solution of the equation y n2  4y n1  4y n  0 .
9. Form the difference equation by eliminating a from y n  a 5n .
10. Form the difference equation by eliminating a and b from each of the relations given
below.
(i) y n  a 2 n  b 3n (ii) y x  a 2 x  b 5x (iii) y x  a 2 x  b 2x
(iv) y n  a n  b3n .
11. Solve the difference equation y n2  3y n1  2y n  5n  2n .
12. Solve the difference equation u n 2  u n1  u n  n 2  n  1 .
13. Solve the following difference equations

(i) y x2  8y x1  15y x  0 (ii) E 2  5E  6 y n  0  
(iii) 2  3  1 y n  0 . 

******

UNIT-II
22
 INTERPOLATION

2.1 Introduction
2.2 Lagrange’s interpolation formula for unequal intervals
2.3 Inverse interpolation by Lagrange’s method
2.4 Newton’s Divided Difference
2.5 Newton’s divided difference interpolation formula for unequal intervals

2.1 Introduction

Interpolation
Definition
Interpolation is the process of estimating the value of a function at an intermediate
point or the process of finding the value of the function inside the given range is called
interpolation.
Interpolation is the process of finding the most appropriate estimate for missing data.
It is the “art of reading between the lines of a table”. For making the most probable estimate
it requires the following assumptions:

(i) The frequency distribution is normal and marked by sudden ups and downs.
(ii) The changes in the series are uniform within a period.

Interpolation technique is used in various disciplines like statistics, economics,

business, population studies, price determination etc. It is used to fill in the gaps in the
statistical data for the sake of continuity of information. For example, if we know the records
for the past five years except the third year which is not available due to unforeseen
conditions, the interpolation technique helps to estimate the record for that year too under the
assumption that the changes in the records over these five years have been uniform.

Extrapolation
Definition
Extrapolation is the process of finding the values outside the given interval.

It is also possible that we may require information for future in which case the process
of estimating the most appropriate value is known as extrapolation.

23
 Given a set of tabular values of a function y=f(x) where the explicit nature of the
function is not known, then f(x) is replaced by a simpler function (x) such that f(x) and (x)
agree with the set of tabulated points. Any other value may be calculated from (x). This
function (x) is known as interpolation function. In particular if (x) is a polynomial then the
process is called polynomial interpolation and (x) is called an interpolating polynomial. The
existence of an interpolating polynomial is supported by Weierstras’s approximation theorem
which asserts that any continuous function on a closed interval can be approximated by a
polynomial.

2.2 Lagrange’s interpolation formula for unequal intervals

Let y=f(x) be the function such that f(x) is taking the values y0, y1, …, yn
corresponding to x= x0, x1, …, xn.
In the case of the values of independent variable are not equally spaced and when the
differences of dependent variable are not small, we will use Lagrange’s interpolation formula.
Let f(x) be a polynomial in x of degree n. Lagrange’s interpolation formula for
unequal intervals is

y  f x  
x  x1 x  x 2 x  x n  f x  
x 0  x1 x 0  x1 x 0  x n  0
x  x 0 x  x 2 x  x 3 x  x n  f x    
x1  x 0 x1  x 2 x1  x n  1

x  x 0 x  x1 x  x 2 x  x n 1  f x  .

x n  x 0 x n  x1 x n  x 2 x n  x n 1  n
Example: 1
Using Lagrange’s interpolation formula, find the value corresponding to x=10 from
the following table:
x 5 6 9 11

y 12 13 14 16

Solution
Given x0 = 5, x1 = 6, x2 = 9, x3 = 11, x = 10
Y0 = f(x0) = 12
Y1 = f(x1) = 13
Y2 = f(x2) = 14
Y3 = f(x3) = 16

24
 y  f x  
x  x1 x  x 2 x  x 3  f x   x  x 0 x  x 2 x  x 3  f x  
x 0  x1 x 0  x 2 x 0  x 3  0 x1  x 0 x1  x 2 x1  x 3  1
x  x 0 x  x1 x  x 3  f x   x  x 0 x  x1 x  x 2  f x 
x 2  x 0 x 2  x1 x 2  x 3  2 x 3  x 0 x 3  x1 x 3  x 2  3


10  610  910  11 12  10  510  910  11 13 
5  65  95  11 6  56  96  11
10  510  610  11 14  10  510  610  9 16
9  59  69  11 11  511  611  9

4.1.1
12  5.1. 1 13  5.4.1 14  5.4.1 14  5.4.1 16
1.4.6 1.3.5 4.3.2 4.3.2 6.5.2
= 14.63.
Y=f(x) = f(10) = 14.63.

Example: 2
Find the polynomial f(x) by using Lagrange’s formula and hence find f(3) for
x 0 1 2 5

f(x) 2 3 12 147

Solution
Given x0 = 0; x1 = 1; x2 = 2; x3 = 5
y0 = f(x0) = 2
y1 = f(x1) = 3
y2 = f(x2) = 12
y3 = f(x3) = 147
The Lagrange’s formula is

y  f (x) 
x  x1 x  x 2 x  x 3  f x   x  x 0 x  x 2 x  x 3  f x  
x 0  x1 x 0  x1 x 0  x 3  0 x1  x 0 x1  x 2 x1  x 3  1
x  x 0 x  x1 x  x 3  f x   x  x 0 x  x1 x  x 2  f x 
x 2  x 0 x 2  x1 x 2  x 3  2 x 3  x 0 x 3  x1 x 3  x 2  3

y  f (x) 
x  1x  1x  5 2  x  0x  2x  5 3 
0  10  20  3 1  01  21  5
x  0x  1x  5 12  x  0x  1x  2 147
2  02  12  5 5  05  15  2

25
 
x  1x  1x  5 2  xx  2x  5 3 
6 4
xx  1x  5
12  xx  1x  2 147
6 60
Which is the polynomial of y=f(x)

y  f (3) 
3  13  13  5 2  33  23  5 3 
6 4
33  13  5
12  33  13  2 147
6 60
y = f(3) = 44.5.
2.3 Inverse interpolation by Lagrange’s method
The process of finding a value of x for the corresponding value of y is called inverse
interpolation. In such a case, we will take y as independent variable and x as dependent
variable.
Therefore the Lagrange’s inverse interpolation formula is

x  f ( y) 
y  y1 y  y 2 y  y n  .f y  
y 0  y1 y 0  y 2 y 0  y n  0
y  y0 y  y 2 y  y n  .f y    
y1  y0 y1  y 2  y1  y n  1
y  y0 y  y1 y  y 2 y  y n1  .f y  .
y n  y0 y n  y1 y n  y 2 y n  y n1  n
Example: 1
Find the value of x, corresponding to y = 100 from the following table:
x 3 5 7 9 11

y 6 24 58 108 174

Solution
Given y0 = 6, y1 = 24, y2 = 58, y3 = 108, y4 = 174 also, y = 100
x = f(y);
x0 = f(y0) = 3, x1 = f(y1) = 5, x2 = f(y2) = 7, x3 = f(y3) = 9
The Lagrange’s formula for inverse interpolation is

x  f ( y) 
y  y1 y  y 2 y  y n  .f y  
y 0  y1 y 0  y 2 y 0  y n  0
y  y0 y  y 2 y  y n  .f y    
y1  y0 y1  y 2  y1  y n  1

26
 y  y0 y  y1 y  y 2 y  y n1  .f y 
y n  y0 y n  y1 y n  y 2 y n  y n1  n
= 0.3534 – 1.5155 + 2.8870 + 7.0676 – 0.1369
x = f(y) = 8.6556.

Example: 2
Find the value of x when y = 85 using Lagrange’s formula for the table
x 2 5 8 14

y 94.8 87.9 81.3 68.7

Solution
Given y0 = 94.8, y1 = 87.9, y2 = 81.3, y3 = 68.7, also, y = 85
x = f(y);
x0 = f(y0) = 2, x1 = f(y1) = 5, x2 = f(y2) = 8, x3 = f(y3) = 14
The Lagrange’s formula for inverse interpolation is

x  f ( y) 
y  y1 y  y 2 y  y n  .f y  
y 0  y1 y 0  y 2 y 0  y n  0
y  y0 y  y 2 y  y n  .f y    
y1  y0 y1  y 2  y1  y n  1
y  y0 y  y1 y  y 2 y  y n1  .f y 
y n  y0 y n  y1 y n  y 2 y n  y n1  n
= -0.1438778 + 3.3798011 + 3.3010599 – 0.2331532
= 6. 3038.
Therefore the value of x when y=85 is 6.3038.

2.4 Newton’s Divided Difference

Let the function y = f(x) takes the values f(x0), f(x1), …, f(xn) corresponding to the
values x1-x0, x2 – x1, x3 – x2, …, xn-xn-1 need not to be equal.
f ( x1 )  f ( x 0 )
Be first divided difference of f(x) for the argument x0, x1 is defined as
x1  x 0

it is denoted by f(x0, x1) or [x0, x1] or f(x0).

f ( x1 )  f ( x 0 )
(i.e) f ( x 0 , x1 )  .
x1  x 0

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

x0 f(x0)

f ( x 0 , x1 )
f ( x1 )  f ( x 0 )
x1 f(x1) 
x1  x 0
f ( x 0 , x1 , x 2 )
f ( x 2 x1 )  f ( x1x 0 )
 f ( x 0 , x1 , x 2 , x 3 )
f ( x1 , x 2 ) x2  x0
f ( x1x 2 x 3 )  f ( x 0 x1x 2 )
f ( x 2 )  f ( x1 ) 
x2 f(x2)  f ( x1 , x 2 , x 3 ) x3  x0
x 2  x1
f ( x 3 x 2 )  f ( x 2 x1 )
 f ( x1 , x 2 , x 3 , x 4 )
x 3  x1
f ( x 2 x 3 x 4 )  f ( x1x 2 x 0 )
f (x 2 , x 3 ) 
f (x 2 , x 3 , x 4 ) x 4  x1
f (x 2 )  f (x 3 ) f (x 4 x 3 )  f (x 3 x 2 )
x3 f(x3)  
x3  x 2 x 4  x1

f ( x 0 , x1 )
f ( x1 )  f ( x 0 )
x4 f(x4) 
x 4  x3

f ( x1 x 2 x 3 x 4 )  f ( x 0 x1 x 2 x 3 )
Fourth divided difference is 4f(x)  f (x 0 , x1 , x 2 , x 3 , x 4 )  .
x4  x0

Properties of Divided Differences

Property: 1
The value of any divided difference is independent of the order of the arguments. That
is, the divided differences are symmetrical in all their arguments.
f x1   f x 0  f x 0   f x1 
f x 0 , x1     f x1 , x 0  (1)
x1  x 0 x 0  x1

28
 f x 0  f x1  f x 0  f x1 
Again, f x 0 , x1      (2)
x 0  x1 x 0  x1 x 0  x1 x1  x 0

f x1  f x 0 
In the same way, f x1 , x 0    (3)
x1  x 0 x 0  x1

From equation (2) and equation (3), we have f x 0 , x1   f x1 , x 0  .

Similarly,
f x1 , x 2   f x 0 , x1 
f x 0 , x1 , x 2  
x2  x0

1  f x1  f x 2    f x 0  f x1  
       
x2  x0  1x  x 2 x 2  x 1  0
x  x 1 x 1  x 0 

 1 1  f x 2  f x 0  
 f x1  
1
    
x2  x0  x1  x 2 x1  x 0  x 2  x1 x 0  x1 

 x2  x0 f x 2  f x 0  
f x1  
1
   
x 2  x 0  x1  x 2  x1  x 0  x 2  x1 x 0  x1 

f x 0  f x1  f x 2 
f x 0 , x1 , x 2     (4)
x 0  x1 x 0  x 2  x1  x 0 x1  x 2  x 2  x 0 x 2  x1 
From equation (4), we find
f x 0 , x1 , x 2   f x1 , x 0 , x 2   f x1 , x 2 , x 0   

This shows that f x 0 , x1 , x 2  is independent of the order of the arguments.

By mathematical induction, we can prove that
f x 0  f x1 
f x 0 , x1 , x 2 ,, x n    
x 0  x1 x 0  x 2 x 0  x n  x1  x 0 x1  x 2 x1  x n 
f x n 
 .
x n  x 0 x n  x1 x n  x n 1 
This is symmetrical we know that any two arguments. Therefore, the divided differences are
symmetrical we know that any two arguments.

Property: 2
The operator  is linear.
Proof:
If f(x) and g(x) are two functions  and  are constants, then

  f x    gx  
 f x1    gx1  f x 0    gx 0 
x1  x 0

29
 f x1   f x 0  gx1   gx 0 
 
x1  x 0 x1  x 0

  f x    gx     f x     gx  .
Corollary: 1
Setting  =  = 1,  f x   gx    f x    gx  .
Corollary: 2
Setting  = 0,   f x     f x  .

Property: 3
The nth divided differences of a polynomial of degree n are constants.
Proof:
Taking f(x) = xn where n is a positive integer,

f x1   f x 0  x1n  x 0n
f x 0 , x1   
x1  x 0 x1  x 0

 x1n1  x 0 x1n2  x 02 x1n3    x 0n1

= a polynomial function of degree (n-1) and symmetrical in x0, x1 with leading
coefficient 1.
Again,
f x1 , x 2   f x 0 , x1 
f x 0 , x1x 2  
x2  x0


x n 1
2  
 x1 x n2 2    x1n 1  x 0n 1  x1 x 0n 2    x1n 1 
x2  x0



x n2 1  x 0n 1 x1 x n2 2  x 0n 2


x n 2 x 2  x 0 
 1
x2  x0 x2  x0 x2  x0

   
 x n22  x 0 x n23    x 0n2  x1 x n23  x 0 x n24  x 0n3    x1n2
= a polynomial of degree (n - 2) and symmetrical i x0, x1, x2 with leading
coefficient 1.
Proceeding in this way, the rth divided differences of xn will be a polynomial of degree
(n-r) and symmetrical in x0, x1, x2, ..., xr with leading coefficient 1.
Hence nth order divided differences of xn will be a polynomial of degree n – n =0 ,
with leading coefficient 1. That is, its value is 1.
That is n x n  1.
n i x n  0, for i  1,2,...

30
 
Hence, n a 0 x n  a1x n 1  a n 
 a 0 n x x  a1 n x n1    n a n

 a 0 1 0  0   0  a 0 .

Note:
Conversely, if the nth divided difference of a polynomial is constant, then the
polynomial is degree of n.

Relation between Divided Differences and Forward Differences

If the arguments x 0, x 1, x2 , ... are equally spaced, then we have,
x1  x 0  x 2  x1  x 3  x 2  x n  x n1  h .

f x1   f x 0   f x 0 
 f x 0   f x1 , x 0   
x1  x 0 h

 f x1    f x 0 
1 1
2 f x 0  
 f x 1    f x 0  h h
x2  x0 2 h

 f x 0 
1 2

2h 2
Similarly,

3 f x 0 
3 f x 0  
3!h 3

n f x 0 
n f x 0   .
n!h n
2.5 Newton’s divided difference interpolation formula for unequal intervals
y  f (x)  f (x 0 )  (x  x 0 )f (x 0 , x1 )  (x  x 0 )(x  x1 ).f (x 0 , x1 , x 2 )

 (x  x 0 )(x  x1 )(x  x 2 ).f (x 0 , x1 , x 2 x 3 )  

 (x  x 0 )(x  x1 )(x  x 2 )(x  x n 1 )f (x 0 , x1, x 2 x 3 ,x n ) .

Example: 1
Using Newton’s divided difference formula, find the values of f(2), f(8) and f(15)
given the following table:
x 4 5 7 10 11 13

f(x) 48 100 294 900 1210 2028

Solution
We form the divided difference table since the intervals are unequal.

31
 x f(x) f(x) 2f(x) 3f(x) 4f(x)

4 48

100  48
5 100  52
54 97  52
 15
294  100 74 21  15
 97 1
7 294 75 202  97 10  4 0
 21
900  294 10  5 27  21
 202 1
10 900 10  7 310  202 11  5
 27 0
1210  900 11  7 33  27
 310 1
11  10 409  310 13  7
11 1210  33
2028  1210 13  10
 409
13  11
13 2028

By Newton’s Divided Difference interpolation formula is

y  f (x)  f (x 0 )  (x  x 0 )f (x 0 , x1 )  (x  x 0 )(x  x1 ).f (x 0 , x1 , x 2 )

 (x  x 0 )(x  x1 )(x  x 2 ).f (x 0 , x1 , x 2 x 3 )

Here x0 = 4, x1 = 5, x2 = 7, x3 = 10, x4 =11, x5 = 13

Also, f(x0) = 48, f (x 0 , x1 )  52 , f (x 0 , x1 , x 2 )  15 , f (x 0 , x1 , x 2 x 3 )  1
y = f(x) = 48 +(x-4)(52)+(x-4)(x-5)(15)+(x-4)(x-5)(x-7).1
f(2) = 48+(2-4)(52)+(2-4)(2-5)(15)+(2-4)(2-5)(2-7)
f(2) = 4
f(8) = 48+(4)(52)+(4)(3)(15)+(4)(3)(1)
f(8) = 448
f(15) = 48 + 11 (52) + (11)(10)(15) + (11)(10)(8)
f(15) = 3150.
Example: 2
Using Newton’s divided difference formula, find u(3) given that u(1) = -26, u(2) = 12,
u(4) = 256, u(6) = 844.
Solution
We form the divided difference table since the intervals are unequal.

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

32
 1 -26

12  26
 38
2 12 2 1 122  38
 28
256  12 4 1 43  28
 122 3
42 6 1
4 256 294  122
 43
844  256 62
 294
6 844 64

Newton divided difference interpolation formula is

y  f (x)  f (x 0 )  (x  x 0 )f (x 0 , x1 )  (x  x 0 )(x  x1 ).f (x 0 , x1 , x 2 )

 (x  x 0 )(x  x1 )(x  x 2 ).f (x 0 , x1 , x 2 x 3 )

Here,
y = u(x) = u(x0) + (x –x0) u(x0,x1) + (x –x0) (x –x1) u(x0,x1,x2)
+(x –x0) (x –x1) (x –x2) u(x0,x1,x2,x3)

Here,
x0 = 1, x1 = 2, x2 = 4, x3 = 6
u(x0) = -26, u(x0,x1) = 38, u(x0,x1,x2) = 28, u(x0,x1,x2,x3) = 3
u(x) = -26 + (x-1)(38) + (x-1)(x-2)(28) + (x-1) (x-2)(x-3) (3)
for x = 3,
y= u(x) = -26 + (3-1) (38) + (3-1) (3-2) (28) + (3-1) (3-2) (3-4) (3)
= -26 + 76 + 56 – 6
u(3) = 100.
Newton’s Forward and Backward Interpolation Formula for Equal Intervals
Newton’s Forward Interpolation Formula
p(p  1) 2 p(p  1)(p  2) 3
y  f ( x )  y 0  py 0   y0   y0  
2! 3!
x  x0
where p  , h is the width of interval
h
x = x0 + ph.
Newton’s Forward Interpolation Formula
p(p  1) 2 p(p  1)(p  2) 3
y  f ( x )  y 0  py n   yn   yn  
2! 3!

33
 x  xn
where p  , h is the width of interval
h
x = xn + ph.
Example: 1
Using Newton’s forward interpolation formula, find the polynomial f(x) satisfying the
following data. Hence evaluate y at x = 5.
x 4 6 8 10

y 1 3 8 10

Solution
We form the difference table
x y y 2 y 3 y

4 1
3-1=2
5-2=3
6 3
8-3=5 -3-3=-6

2-5=-3
8 8
10-8=2

10 10

There are only 4 data given. Hence the polynomial will be degree 3.
Therefore Newton’s –Gregory Forward interpolation Formula is
p p(p  1) 2 p(p  1)(p  2) 3
y  f ( x )  y 0  y 0   y0   y0
1! 2! 3!
x  x0 x  4
Here y0 = 1; p  
h 2
 x 4  x  4  x  4   x  4  x  4  x  4 
     1    1  2
y  f (x)  1  
2 
(2)  
2  2  (3)   2  2  2  (6)
1! 2! 3!

 y  f (x) 
1
8

 x 3  21x 2  126x  240 
When x= 5,

 y  f (5) 
1
8
 
 (5) 3  21(5) 2  126(5)  240 =1.25

Y = 1.25 when x = 5.
Example: 2

34
 A third degree polynomial passes through the points (0,-1), (1,1), (2,1) and (3,-2)
using Newton’s forward interpolation formula find the polynomial. Hence find the value at
1.5.
Solution
We form the difference table
x y y 2y 3y

0 -1
1+1 = 2
0-2=-2
1 1
1-1=0 -3+2=-1

-3-0=-3
2 1
-2-1=-3

3 -2

There are only 4 data given. Hence the polynomial will be degree 3.
Therefore Newton’s –Gregory Forward interpolation Formula is
p p(p  1) 2 p(p  1)(p  2) 3
y  f ( x )  y 0  y 0   y0   y0
1! 2! 3!
x  x0 x  0
Here y0 = -1; p   px
h 1
x x (x  1) x (x  1)(x  2)
y  f (x )  1  2 (2)  (1)
1! 2! 3!

 y  f (x)  
6

1 3
x  3x 2  16x  6 
When x= 1.5,

 y  f (1.5)  
1
6
 
(1.5) 3  3(1.5) 2  16(1.5)  6 =1.3125

Y = 1.3125 when x = 1.5.

Example: 3
Use Newton’s backward difference formula to construct an interpolating polynomial
of degree 3 for the data: f(-0.75) = - 0.07181250; f(-0.5) = -0.024750; f(-0.25) = 0.33493750;
f(0) = 1.10100. Hence find f(-1/3).
Solution
We form the difference table

35
 x y y 2y 3y

-0.75 -0.07181250
0.0470625
0.312625
-0.50 -0.024750
0.3596875 0.09375

0.400375
-0.25 0.33493750
0.7660625

0 1.10100

Newton’s backward difference formula is

p p(p  1) 2 p(p  1)(p  3) 3
y  f ( x )  y 3  y 3   y3   y3
1! 2! 3!
x  x3
where p  ; p = 4x.
h
4x 4x (4x  1) 4x (4x  1)(4x  2)
 y  f ( x )  1.10100  (0.7660625)  (0.406375) 
1! 2! 3!
y  f (x)  x 3  4.001x 2  4.002x 1.101
3 2
 1  1  1  1
 f         4.001    4.002    1.101 = 0.174518518.
 3  3  3  3

Exercises:
1. Using Lagrange’s interpolation formula, find f(4) given that f(0) = 2, f(1) = 3, f(2) 12,
f(15) = 3587.
2. Find the third degree polynomial f(x) satisfying the following data. Also, find f(4),
f(6).
x 1 3 5 7

y 24 120 336 720

3. Using Lagrange’s interpolation find y(2) from the following data

x 0 1 3 4 5

y 0 1 81 256 625

4. Apply Lagrange’s inverse formula to obtain the root of equation f(x) = 0. Given that,
f(0) = -4; f(1) = 1; f(3) = 29; f(4) 52.

36
5. Find f(x) as a polynomial in x for the following data by Newton’s divided difference
formula:
x -4 -1 0 2 5

f(x) 1245 33 5 9 1335

6. Using Newton’s Divided difference formula, fit a polynomial to the data and hence
find y chosen x = 1.
x -1 0 2 3

y -8 3 1 12

3 1  3 1  1
1     
7. If f(x) = , find f(a, b, c, d) or  a  (or)  a   
x abcd

8. Using Newton’s forward interpolation formula, find the polynomial satisfying the
following data. Hence find f(x).
x 0 5 10 15

y 14 379 1444 3584

9. Use Newton’s forward interpolation formula find the cubic polynomial which takes
places the following values:
x 0 1 2 3

y 1 2 1 10

10. State Lagrange’s interpolation formula.

11. What is the Lagrange’s formula to find y if there sets of values (x0, y0), (x1, x2) and
(y1, y2) are given.
12. What is the assumption we make when Lagrange’s formula is used?
13. Give the inverse of Lagrange’s interpolation formula.
14. Using the Newton’s divided difference formula, find the missing value from the table:
x 1 2 4 5 6

y 14 15 5 - 9

*****

37
 UNIT – III
NUMERICAL DIFFERENTIATION AND INTEGRATION

3.1 Introduction
3.2 Numerical Differentiation
3.3 Numerical Integration
3.4 Trapezoidal Rule
3.5 Simpson’s One Third Rule
3.6 Simpson’s Three Eight Rule

38
 3.7 Waddle’s Rule
3.8 Cote’s Method

3.1 Introduction
We assume that a function f(x) is given in a tabular form at a set of n+1 distinct points
x0, x1, …, xn. From the given tabular data, we require approximations to the derivatives
f ( r ) (x), r  1 , where x  may be a tabular or a non-tabular point. We consider the cases
r = 1, 2.
In many applications of science and engineering, we require to compute the value of
b
the definite integral  f ( x ) dx , where f(x) may be given explicitly or as a tabulated data. Even
a

when f(x) is given explicitly, it may be a complicated function such that integration is not
easily carried out.
Here, we shall derive numerical methods to compute the derivatives or evaluate an
integral numerically.

3.2 Numerical Differentiation

Approximation to the derivatives can be obtained numerically using the following two
approaches
(i) Methods based on finite differences for equispaced data.
(ii) Methods based on divided differences or Lagrange interpolation for non-uniform data.
Numerical differentiation is the process of calculating the derivation of a given
function by means of a table of given values of that function. That is, if (xi, yi) are the given
dy d 2 y
set of values, then the process of computing the values of , ,, etc. at a given point is
dx dx 2
called numerical representation.
The interpolation formula depends on the particular value of x at which the
derivatives are required.
(i) If the values of x are not equally spaced, we represent the function by Newton’s
divided difference formula and the derivatives are obtained.
(ii) If the values of x are equally spaced the derivatives are calculated by using Newton’s
Forward or backward interpolation formula.

3.3 Numerical Integration

39
 b
The process of evaluating a definite integral  f x  dx from a set of tabulated values
a

(xi, yi); i=0, 1, …, n of the integrand y = f(x) is called numerical integration.

Newton cote’s formula (or) General Quadrature formula for equidistant coordinates
b
Let I   y dx where y = f(x) takes the values y0, y1, …, yn for x0, x1, …, xn. Let us
a

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. After simplification, we get

 n2n  3 2 nn  2 3 
xn 2
n
I   ydx  nh  y 0  y 0   y0   y 0   (1)
x0  2 12 24 

which is the general quadrature formula.

By putting n = 1 in equation (1), Trapezoidal rule is obtained.

1
By putting n = 2 in equation (1), Simpson’s rule is obtained.
3

3
By putting n = 3 in equation (1), Simpson’s rule is obtained.
8

3.4 Trapezoidal Rule

Putting n = 1 in equation (1), we get

x1 x0  h
 1 
 y dx 
x0

x0
y dx  h  y 0  y 0 
 2 


h
2 y0  y0 
2


h
y0   y0  y0 
2


h
y0  y1 
2

Similarly,
x2 x0  2 h
 1 
 y dx  
x1 x0  h
y dx  h  y1  y1 
 2 


h
2 y1  y1 
2
40
 
h
y1   y1  y1 
2


h
y1  y 2 
2

and so on.
xn x0  nh

 y dx   y dx 
h
y n1  y n 
xn 1 x0  n 1h
 2

 Adding all the above, we get

xn

 y dx  2  y  y n   2 y1  y 2  y3    y n 1  which is called Trapezoidal Rule.

h
0
x0

1
3.5 Simpson’s Rule (or) Simpson’s Rule
3

Putting n = 2 in equation (1) and neglecting the differences of higher order than
second order. We get,
xn

 y dx  3  y  y n   4 y1  y3    y n 1   2 y 2  y 4    y n 2  .
h
0
x0

Note:

It should be noted that for applying this rule, the interval must be divided into even
number of sub intervals of width h.

3
3.6 Simpson’s th Rule
8

Putting n = 3 in equation (1) and neglecting the higher order differences above the
third, we get
xn

 y dx  8  y  y n   3 y1  y 2  y 4  y5   y n 1   2 y3  y 6    y n 3 
3h
0
x0

Note:

This is not as accurate as Simpson’s rule. This rule is used when the number of
subintervals is a multiple of 3.

Example: 1

1
dx
Using Trapezoidal rule, evaluate 1 x
1
2
taking 8 intervals.

41
Solution

Given y  f x  
1
1 x2

Given x = -1 to 1.

 Length of interval = 2.

2 ba
h   0.25 . h 
8 number of int ervals

1   1
  0.25
8

 We form a table

x -1 -0.75 -0.5 -0.25 0 0.25 0.50 0.75 1

y 0.5 0.64 0.8 0.9412 1 0.9412 0.8 0.64 0.5

xn

 y dx  2  y  y n   2 y1  y 2  y3    y n 1 
h
 Trapezoidal rule is 0
x0

1 x
1
2
dx 
0.25
0.5  0.5  20.64  0.8  0.9412  1  0.9412  0.8  0.64
1
2


0.25
1  25.7624
2

=1.5656.

Example: 2
1
dx 1
Evaluate 1 x
0
2
with h 
6
by Trapezoidal rule.

Solution

Given y  f x  
1
1 x2

Given x = 0 to 1

1
Also h  .
6
42
The table is

x 0 1/6 2/6 3/6 4/6 5/6 1

y 1 36/37 9/10 4/5 9/13 36/61 1/2

By Trapezoidal rule,
xn

 y dx  2  y  y n   2 y1  y 2  y3    y n1 
h
0
x0

1
1
1  
6  1   36 9 4 9 36 
0 1  x 2 dx  2 1  2   2 37  10  5  13  31 

1 3 
   23.9554
12  2 

= 0.7842.

Example: 3

5.2

 log e dx by using (i) Trapezoidal rule (ii) Simpson’s rule (iii) Simpson’s
x
Evaluate
4

3
rule, given that h = 0.2.
8

Solution

Given y = f(x) = log ex

x = 4 to 5.2, h = 0.2.

The table is

x 4 4.2 4.4 4.6 4.8 5.0 5.2

y 1.737 1.824 1.910 1.997 2.084 2.171 2.258

(i) By Trapezoidal rule,

xn

 y dx  2  y  y n   2 y1  y 2  y3    y n 1 
h
0
x0


0.2
1.737  2.258  21.824  1.910  1.997  2.084  2.171
2

43
 = 0.3967.

1
(ii) By Simpson’s rule (or) Simpson’s rule,
3
xn

 y dx  3  y  y n   4 y1  y3    y n 1   2 y 2  y 4    y n 2 
h
0
x0

5.2

 log e
x
dx 
0.2
1.737  2.258  41.824  1.997  2.171  21.910  2.084
4
3

= 0.0666[3.995+4(5.992)+2(3.994)]

= 2.394.

3
(iii) By Simpson’s rule,
8
xn

 y dx  8  y  y n   3 y1  y 2  y 4  y5   y n 1   2 y3  y 6    y n 3 
3h
0
x0

30.2
 1.737  2.258  31.824  1.910  2.084  2.171  21.997
8

= 2.3967.

Romberg’s Method

b
Romberg’s method is used to evaluate the integral of the form I   y dx .
a

For Romberg’s method, let us apply Trapezoidal rule several times find the value of
I’s as follows:

h
I1: Dividing ‘h’ into 2 parts (i.e)
2

h h/2
I2: Dividing ‘h’ into 4 parts (i.e)  
4 2 

h h/4
I3: Dividing ‘h’ into 8 parts (i.e)  
8  2 

h  h/8
I4: Dividing ‘h’ into 16 parts (i.e)  
16  2 

and so on.

44
Applying Romberg’s formula,

I I 
I  I2   2 1 
 3 

For I1, I2; I2, I3; I3, I4; …

We get the values of I. This method continues till we get two successive values of I’s
are equal. The systematic refinement of the values of I’s is called Romberg’s method.

Example: 1

1
dx
Use Romberg’s method to compare  . Correct to 4 decimal places and hence
0 1  x 2

find an approximate value of π.

Solution

1
dx
Let I  
0 1 x
2

To find I1:

h
Dividing h into 2 parts (i.e) .
2

1 0
h   0.5
2

X 0 0.5 1

1 1 0.8 0.5
y=
1 x2

By Trapezoidal rule,

I1 
h
 y0  y 2   2 y1 
2


0.5
1  0.5  20.8
2

= 0.775.

To find I2:

h
Dividing h into 4 parts (i.e) .
4
45
 1 0
h   0.25
4

x 0 0.02 0.50 0.75 1

y 1 0.941 0.8 0.64 0.5

By Trapezoidal rule,

I2 
0.25
1  0.5  20.941  0.8  0.64
2

= 0.7828.

To find I3:

h/4
Dividing h into 8 parts (i.e)
2

1 0
h   0.125
8

x 0 0.125 0.25 0.375 0.5 0.625 0.75 0.875 1

y 1 0.9846 0.9412 0.8767 0.80 0.7191 0.64 0.5664 0.5

By Trapezoidal rule,

I3 
0.125
1  0.5  20.984  0.94  0.876  0.8  0.719  0.64  0.56
2

= 0.78475.

By Romberg’s Formula:

Iteration 1:

I I 
I  I2   2 1 
 3 

 0.7828  0.775 
 0.7828   
 3 

I = 0.7854.

Iteration 2:

46
 I I 
I  I3   3 2 
 3 

 0.78475  0.7828 
 0.78475   
 3 

I = 0.7854.

From the last two iterations, I = 0.7854.

To find π:

1
 1  x 
1
dx
0 1  x 2  tan  1 
0

 tan 1 1  tan 1 0


 .
4

1
dx
But 1 x
0
2
 0.7854


0.7854 
4

π = 4(0.7854)

  3.1416 .

Two points and three points Gaussian Quadrature formulas

Gaussian Quadrature
1
Gauss derived a formula which is used to evaluate the integral of the form  F u  du.
1

One point Gaussian Formula

1

The one point Gaussian formula is given by  f x  dx  2. f 0 which is exact for
1
polynomials of degree up to 1.

Two point Gaussian Formula

47
 1
 1  1
The two point Gaussian formula is 1 f  x  dx  f 
  
  f  
 3  and this is exact
 3   
for polynomial of degree up to 3.

Three point Gaussian Formula

5  3  3 
1

 f x  dx  9 f 0  9  f  
8
The three point Gaussian formula is  f  
1   5   5 

which is exact for polynomial of degree up to 5.
Note:
b 1
(i) The integral  F t  dt , can be transformed into
a
 f x  dx by the linear
1

transformation.
 b  a x  b  a 
t  .
 2
b 1
(ii) The integral  F x  dx , can be transformed into
a
 f t  dt by the linear
1

transformation.
 b  a t  b  a 
x  .
 2

Example: 1
1
dx
Evaluate 1 x
1
2
by two point and three point Gaussian formula and compare with

exact value.
Solution
1
 1  1
By two point Gaussian formula,  f x  dx  f   
3 
f 
3


1   

Given f x  
1
.
1 x2

 1 1 3
 f   
 
 3 1 4
1
3

 1 1 3
 f  
 
 3  1 1 4
3

48
 1
  f x  dx 
3 3 6
   1.5
1
4 4 4
1
 f x  
1
1 x
1
2
dx  1.5 .

5  3  3 
1

 f x  dx  9 f 0  9  f  
8
By three point Gaussian formula,  f  
1   5  
 5 

f 0 
1
1
1 0

 3 1 5
 f   
 
 5 3 8
1
5

 3 1 5
 f  
 
 5  1 3 8
5

5 5 5 
1
  f x  dx 
8
1    
1
9 9 8 8 

8 5 10
  
9 9 8
1
  f x  dx  1.5833 .
1

But exact value is

1 1
dx dx
11  x 2  20 1  x 2

 
 2 tan 1 x  0
1

 2tan 1  tan 0

1 1

 
 2  0 
4 


2
1
dx
1 x
1
2
 1.5708 .

Thus, exact value is 1.5708.

By Gaussian two point formula value is 1.5.
49
By Gaussian three point formula value is 1.5833.
Example: 2
1
x2
Using three point Gaussian Quadrature formula, evaluate  dx .
1 1  x
4

Solution

5  3  3 
1
f x  dx  f 0   f  
8
Three point Gaussian Quadrature formula is  9  
 f
5 

 5 
1
9  
1
x2
Given  dx
1 1  x
4

x2
Here f x  
1 x4

 f 0 
0
0
1 0

3
 3 15
f    5 2 
 5  3 34
1 2
5
Similarly,

 3  15
f  
 34
 5 

5  3  3 
1
f x  dx  f 0   f  
8
 9  
 f
5 

 5 
1
9  

8 5  15 15 
 0   
9 9  34 34 

1
x2
 dx  0.4902 .
11  x
4

Double integrals using Trapezoidal and Simpson’s Rules

b d
We shall evaluate double integral I    f x, y  dx dy using Trapezoidal and
a c

Simpson’s rule.

The formula for the evaluation of a double integral can be obtained by repeatedly
applying the Trapezoidal and Simpson’s rules.

50
Trapezoidal Rule for Double Integral

The formula for the evaluation of double integral using Trapezoidal rule is,

sum of the values at four corner of the box   

hk 
I  2sum of the values at the boundary of the box except the corner   
4
4sum of the remaing values  

where h – length of x values and k – length of y values.

Simpson’s Rule for Double Integral

The formula for the evaluation of double integral using Simpson’s rule is,

sum of the values at four corner of the box   

hk 
I  4sum of the values at the boundary of the box except the corner   
9
16sum of the remaing values  

where h – length of x values and k – length of y values.

Example: 1

2 2
Evaluate   f x, y  dx dy by Trapezoidal rule for the following data:
0 0

x 0 0.5 1 1.5 2

0 2 3 4 5 5

1 3 4 6 9 11

2 4 6 8 11 14

Solution

Here, h = 0.5, k = 1.

2 2
 I    f x, y  dx dy
0 0

By Trapezoidal rule,
51
 sum of the values at four corner of the box   
hk 
I  2sum of the values at the boundary of the box except the corner   
4
4sum of the remaing values  


0.5  1 2  5  14  4  23  4  5  11  11  8  6  3  44  6  9
4

I = 25.375.

Example: 2

1 1
1 1
Using Simpson’s
3
rule evaluate   1  x  y dx dy taking h = k = 0.5.
0 0

Solution

x = 0, 0.5, 1; y = 0, 0.5, 1. The table values are

x 0 0.5 1

0 1 0.6667 0.5

0.5 0.6667 0.5 0.4

1 0.5 0.4 0.3333

Let

1
gy  
1
dx dy
0
1  x  y

By Simpson’s rule,

sum of the values at four corner of the box   

hk 
I 4sum of the values at the boundary of the box except the corner   
9 
16sum of the remaing values  

g 0 
0.5
1  0.5  40.6667  1.5  2.66668
3

g 0  0.69441889 .

g 0.5 
0.5
0.6667  0.4  3.0667  40.5
3
52
g 0.5  0.51111 .

g 1 
0.5
0.5  0.3333  40.4
3

g 1  0.405538 .

Hence

1
I   g  y  dy
0


0.5
0.9441889  0.405538  40.51111
3

I  0.5241.
3.7 Weddle’s Rule
Put n = 6 in Newton-Cot’s Quadratic formula and neglecting all differences of orders
higher than sixth, we get
x 0  nh

 f (x)dx  10 y0  5y1  y 2  6y3  y 4  5y5  y6   y6  5y7  y8  6y9  y10  5y11  y12 
3h
x0

   y n6  5y n5  y n4  6y n3  y n2  5y n1  y n 

This equation is called Weddle’s rule.

x5  x 0  5h

Example: 1
1
dx
Evaluate  1  x using Weddle’s rule with 6 equal intervals.
0

Solution
1
Here n = 6,  h  ,
6
1 1 2 3 6
Let y  f ( x )   x  0, , , ,, .
1 x 6 6 6 6

x 0 1 2 3 4 5 1
6 6 6 6 6
1 1 0.8571 0.75 0.667 0.6 0.5455 0.5
y
1 x

53
By Weddle’s rule,
1
 1  x dx  10 y 0  5y1  y 2  6y3  y 4  5y5 
1 3h
0

1
3 
   1  5(0.8571)  0.75  6(0.6667)  0.6  5(0.5455  0.5
6
10
= 0.69320.
3.8 Newton’s – Cote’s Formula (or) Cote’s Formula
b
Newton’s-Cote’s formula gives a way for computing the integral  f ( x )dx
a

numerically, when y = f(x) is known at equidistance values of x, but its derivation is based
on the integration of Lagrange’s interpolation formula.
b n n
 f (x)dx  h  y k  C k (u)du
a k 0 0

Where C k (u )  (1) n k u k . n  u nk

k!(n  k )!
Here Ck (u ) is a polynomial of degree n equation is called the cote’s polynomial.
Example: 1
Find the Cote’s polynomials for n=1.
Solution
The Cote’s polynomials are C0(u) and C1(u)

C k (u )  (1) n k
u k . n  u nk
k!(n  k )!
Put n = 1 and k = 0 in Ck (u )

C 0 (u ) 
u 0 . 1  u 1  u  1.
0!(1  0)!

C1 (u ) 
u 1. 1  u 0  u.
1!(1  1)!
Example: 2
Find the Cote’s polynomials for n=2.
Solution

C k (u )  (1) n k
u k . n  u nk
k!(n  k )!

54
 C 0 (u ) 
u 0 . 2  u 2
0!(2  0)!
(2  u )(1  u )

2
2  3u  u 2

2

C (u ) 
u  . 2  u 
1 1

1!(2  1)!
1

u (2  u )

1
= 2u-u2

C 2 (u ) 
u 2 . 2  u 0
2!(2  2)!
u (u  1)

2
 u  u2
C2 (u)  .
2

Exercises:

2
dx
1. Evaluate the integral 1 x
1
2
using Trapezoidal rule with two subintervals.


2
2. Dividing the range into 10 equal parts, find the value of  sin x dx by (i) Trapezoidal rule
0

(ii) Simpson’s rule.

1
3. Using Simpson’s one third rule evaluate  xe
x
dx taking 4 intervals. Compare your result
0

with actual value.

0.7

e
x
4. Calculate x dx taking 5 ordinates by Simpson’s rule.
0.5
2
dx
5. Evaluate x
0
2
4
using Romberg’s method. Hence, obtain an approximate value of π.


6. Using Romberg’s method, Evaluate  sin x dx correct to four decimal places.
0
1
1
7. Using three point Gaussian Quadrature formula, Evaluate  dt .
0 1 t
2

55
 x 2  2x  1
2
8. Evaluate  1  x  1
0
4
dx by Gaussian three point formula.

1.5

e
 x2
9. Evaluate dx using the three point Gaussian Quadrature.
0.2
1
dx
10. Evaluate  1  x by Gaussian formula with two points.
0
5
1
11. Use Gaussian three point formula and evaluate  dx .
1
x

 3x   3x 
1 1
12. Using Gaussian three point formula, evaluate (i) 2
 5 x 4 dx (ii) 2
 5 x 4 dx .
1 0

Also compare with exact values.

2 x

13. Evaluate e
2
2
dx by Gaussian two point formula.

1 2
2 xy
14. Evaluate   1  x 1  y  dx dy by Trapezoidal rule with h = k = 0.25.
0 1
2 2

2 2
dx dy
15. Evaluate  x
1 1
2
 y2
numerically with h = 0.2 along x direction and k = 0.25 along y

direction.
3 2
16. The function f(x, y) is defined by the following table. Compute   f x, y  dx dy using
1 0

Simpson’s rule in both direction.

x 0 0.5 1 1.5 2

1 2 1.5 1.3 1.4 1.6

2 3.1 2.5 2 2.3 2.9

3 4.2 4 3.8 4.1 4.4

6
dx
17. Apply Weddle’s rule to evaluate the approximate value of the integral  1 x2 by
0

dividing the range into 6 equal parts.

18. Compute Cote’s polynomials for n = 3, 4, 5 and 6.
19. Compute Cote Numbers for n=1, 2, 3, 4, 5 and 6.
20. Verify that sum of Cote’s Numbers is 1 for n=1, 2, 3, 4, 5 and 6.

56
 *****

UNIT-IV
NUMERICAL SOLUTIONS OF DIFFERENTIAL EQUATIONS

4.1 Introduction
4.2 Picard’s Method
4.3 Initial Value Problem for Ordinary Difference Equations
4.4 Multistep Method
4.5 Runge - Kutta Method
4.6 Solution of Algebraic and Transcendent Equations

4.1 Introduction
Many problems in science and engineering can be reduced to the problem of solving
differential equations satisfying given conditions. By applying analytical methods we can
solve several standard types of differential equations. However the differential equations
appearing in physical problems are quite complex and may not posses closed form solutions.
In such cases they can be solved numerically.
We know that the general solution of a differential equation of the nth order has n
arbitrary constants. In order to compute the numerical solution of such an equation we need
n conditions. If all the n conditions are specified at the initial point only then it is called an
initial value problem. If the conditions are specified at two or more points, then it is called a
boundary value problem.

 f x, y  with the initial condition y(x0) = y0.

dy
Consider the initial value problem
dx
This problem can be solved any of the methods give the solution in one of the two
forms given below:

57
(i) A series for y in terms of powers of x, from which the value of y can be obtained by
direct substitution. The methods of Taylor and Picard belong to this type. In these
methods y is approximated by a truncated series and each term of the series is a
function of x. The information about the curve at one point is used and the solution is
not iterated. Hence these methods are called single step methods or point wise
methods. A solution of this type is called a point wise solution.
(ii) Given a set of tabulated values of x and y, we obtain y by iterative process. The
methods of Euler, Runge – Kutta, Milne, Adams – Bashforth etc. belong to this type.
In these methods, the values of y are computed by short steps ahead for equal
intervals h of the independent variable. These values are iterated till we get the
desired accuracy. Hence these methods are called step by step methods.

4.2 Picard’s Method

Consider the first order differential equation

 f x, y 
dy
(1)
dx
with initial condition y = y0 when x = x0.
We now replace equation (1) by an equivalent integral equation.
Integrating equation (1) we get
y x
 d y   f x, y dx
y0 x0

x
(i.e) y  y 0   f x, y dx (2)
x0

This is an integral equation which contains the unknown y under the integral sign.
Equation (2) is equivalent to equation (1) since any solution of equation (2) is a solution of
equation (1) and vice versa.
The first approximation y1 to the solution is obtained by putting y = y0 in f(x,y) and
from equation (2) we have
x
y1  y 0   f x, y 0 dx .
x0

Similarly for the second approximation y2, put y = y1 in f(x,y) and from equation (2)
we have
x
y 2  y 0   f x, y1  dx .
x0

Continuing this process the nth approximation is given by

58
 x
y n  y 0   f x, y n 1 dx .
x0

This is known as Picard’s iteration formula.

Note:
Picard’s method gives a sequence of approximations y1, y2, ... each giving a better
result than the preceding one. But this can be applied only to equations in which the
successive integration can be obtained easily.
Example: 1
dy
Using Picard’s method solve  1  xy with y(0) = 2. Find y(0.1), y(0.2) and y(0.3).
dx

Solution
The Picard’s iteration formula for the differential equation
x
 f x, y  is y n  y 0   f x, y n 1  dx where n = 1, 2, ...
dy
dx x0

Given f(x, y) = 1+ xy, x0 = 0 and y0 = 2.

The first approximation is
x
y1  y 0   f x, y dx
x0

x
 2   f x,2dx
0

x
 2   1  2x dx
0

y1  2  x  x 2 .
The second approximation is
x
y 2  y 0   f x, y1  dx
x0

  
x
 2   1  x 2  x  x 2 dx
0

x3 x 4
y2  2  x  x 2   .
3 4
The third approximation is

59
 x
y 3  y 0   f x, y 2  dx
x0

x   x 3 x 4 
 2  1  x  2  x  x 2
   dx


x0   3 4 

x3 x 4 x5 x6
y3  2  x  x 2     .
3 4 15 24
Putting x = 0.1, 0.2 and 0.3 in equation (1) we get
y1  y0.1  2.1104
y 2  y0.2  2.2431
y3  y0.3  2.4012 .
Example: 2
dy y  x
Find the value of y(0.1) by Picard’s method given  , y0  1 .
dx y  x

Solution
The Picard’s iterative formula for the differential equation
x
 f x, y is y n  y 0   f x, y n 1 dx where n  1,2,...
dy
dx x0

yx
Here f x, y   , x 0  0 and y 0  1 .
xy
The first approximation is
x
y1  y 0   f x, y 0 dx
0

x
1 x 
 1    dx
0 1 x 

x
 2 
 1    1   dx (By partial fraction)
0 1 x 

 1   x  2 log e 1  x 0
x

y1  1  x  2 log e 1  x  .

Putting x = 0.1 we get,

y1  y0.1  1  0.1  2 log e 1  0.1

 0.9  2  0.0953
y1  1.0906 .

60
4.3 Initial Value Problem for Ordinary Difference Equations
The differential equation together with initial conditions is called an initial value
problem. In this unit, we are going to solve numerically, the first order initial value problem
defined by,
dy
 f ( x, y); y( x 0 )  y 0 .
dx
The solution of such initial value problem can be obtained by two different methods:
1. Single step method
2. Multi step method.
The following are the single step method:
1. Euler method
2. Modified Euler method
3. Taylor series method
4. Runge - Kutta method.

All the above methods, require the information at a single point at x = x0.

The following are the multi step methods:

1. Milne’s method
2. Adam’s – Bashforth method.

Euler and Modified Euler methods

Taylor’s series method and Picard’s method are used to yield the solution of a
differential equation in the form of power series. But Euler methods are used to find the
solution in the form of table values of equally spaced points.
Euler Method
The formula is y n1  y n  hf (x n , y n ); n  0,1,2,3,
Modified Euler method Formula
 1 1 
y n 1  y n  h f x n  h, y n  h f ( x n , y n ) .
 2 2 
Example: 1
dy
Using Euler’s method, find y(0.2), y(0.4), y(0.6) from  x  y; y(0)  1 .
dx
Solution

61
 x 0 0.2 0.4 0.6
y 1 y1 y2 y3
Given y  x  y ; h=x1 – x0
Here h = 0.2
dy
y   f ( x , y)  x  y
dx
By Euler method
To find y1 = y(0.2)
y1  y 0  hf (x 0 , y 0 )

 1  0.2(x 0  y 0 )
= 1 + 0.2 (0+1) = 1.2
y1 = 1.2
To find y1 = y(0.4)
y 2  y1  hf (x1 , y1 )

 1.2  0.2(x1  y1 )
= 1.2 + 0.2 (0.2+1.2) = 1.48
y2 = 1.48

To find y3 = y(0.6)
y3  y 2  hf (x 2 , y 2 )

 1.48  0.2(x 2  y 2 )
= 1.48 + 0.2 (0.4+1.48) = 1.48
y3 = 1.856.
Example: 2
Using Euler’s method, solve y  x  y  xy , y(0) =1. Compute y at x = 0.1, by taking
h = 0.05.
Solution
x 0 0.05 0.10
y 1 y1 y2

Given y  x  y  xy  f (x, y) ; h=0.05

dy
 f ( x, y)  x  y  xy
dx
By Euler method
62
To find y1 = y(0.05)
y1  y 0  hf (x 0 , y 0 )

 1  0.05(x 0  y0  x 0 y0 )
= 1 + 0.05 (0+1+0) = 1.05
y1 = 1.05
To find y2 = y(0.10)
y2 = y1 +hf (x0, y0)
 1.05  0.05(x1  y1  x1y1 )
= 1.05 + 0.05 (0.05+1.05+0.051.05) = 1.107625
y2 = 1.107625.
Example: 3
Compute y at x = 0.25 by modified Euler method given y  2xy , y(0) = 1
Solution
dy
Given y  2xy  f (x, y)  2xy
dx
x 0 0.25

y 1 y1

To find y1 = y(0.25)
By modified Euler method,
 1 1 
y n 1  y n  h f x n  h, y n  h f ( x n , y n )
 2 2 
 h h 
y1  y 0  h f x 0  , y 0  f ( x 0 , y 0 )
 2 2 
 0.25 0.25 
y1  1  0.25 f 0  ,1  .2x 0 y 0 
 2 2 
y1 = y(0.25) = 1.625.
Example: 4
Using modified Euler’s method, compute y(0.1) with h=0.1 from
2x
y  y  , y(0)  1.
y
Solution
x 0 0.1
y 1 y1

63
 dy 2x
Given y   f ( x , y)  y 
dx y
To find y1 = y(0.1)
By modified Euler method,
 1 1 
y n 1  y n  h f x n  h, y n  h f ( x n , y n )
 2 2 
 h h 
y1  y 0  h f x 0  , y 0  f ( x 0 , y 0 )
 2 2 
 0.1 0.1  2 x 
y1  1  0.1f 0  ,1  . y 0  0 
 2 2  y 0 

y1 = y(0.1) = 1.09548.

Taylor Series Method

Consider the differential equation
dy
y   f ( x, y); y( x 0 )  y 0
dx
The solution of the above equation obtained by Taylor series as follows:

y( x )  y 0 
x  x 0  y  x  x 0 2 y  x  x 0 3 y  
0 0 0
1! 2! 3!
It is called power series solution.
In general, the Taylor’s algorithm is given as follows:
h h2 h3
y n 1  y n  y n  y n  yn   where n = 0, 1, 2, 3, …
 

1! 2! 3!
where h is the step size; h = x1 – x0.

Example: 1
Find by Taylor’s series method, the values of y at x = 0.1 and x = 0.2, correct to four
dy
decimal places from  x 2 y  1, y(0)  1 .
dx
Solution
Here x0 = 0; y0 = 1
 y  x 2 y  1 y0  x 02 y 0  1

= 0 (1) -1 = -1

y  x 2 y  1 y0  1

64
 y  x 2 y  y2x y0  x 02 y0  2x 0 y 0
y  x 2 y  2xy y0  0

y  2xy  y2  x 2 y  y2x y0  2y0  4x 0 y0  x0 y0

y  2y  4xy  x 2 y y0  2

y(iv )  2y  4(xy  y)  (x 2 y  y2x) y (iv )  2y0  4(x 0 y0  y0 )

y (iv )  6y  6xy  x 2 y y (iv )  6

Therefore Taylor series of y(x) about x0 = 0 is given by

y( x )  y 0 
x  x0 

y 
x  x0 
2


y 
x  x0 
3
y  
0 0 0
1! 2! 3!

y( x )  1 
x  0
(1) 
x  0
2
0
x  0
3
(2) 
x  0
4
(6)  
1! 2! 3! 4!
x x3 x 4
y( x )  1     ...
1! 3! 4!
Hence y(0.1) = 0.9003
Also, y(0.2) = 0.8023.
Example: 2
Solve y  y 2  x ; y(0)=1 using Taylor series method and computer y(0.1) and y(0.2).
Solution
Given y(0) = 1
Here x0 = 0; y0 = 1.

y  x 2  x y0  1

y  2yy  1 y0  3

y  2( yy  yy)  4xy  x 2 y y0  8

y (iv )  2( yy  yy  yy  yy) y (iv )  34

To find y(0.1):
Let y1 = y(x1)
h h2 h3
Taylor algorithm is y1  y 0  y0  y0  y0  
1! 2! 3!
0.1 0.01 0.001 0.0001
y(0.1)  1  (1)  (3)  (8)  (34)  
1! 2! 3! 4!
y(0.1)  1.116411
To find y(0.2):

65
 Let y2 = y(x2)
Then by Taylor’s algorithm
h h2 h3
y 2  y1  y1  y1  y1   where h = 0.1
1! 2! 3!
Given y  y 2  x

y1  y1  x1 = 1.3463
2

y1  2y1 2y1  4.006

y1  2yy  2y2  12.5696

0.01 0.001
y(0.2)  1.1164  (0.1)(1.3463)  (4.006)  (12.5696)  
2 3!
y(0.2) = 1.2732.
4.4 Multistep Method (Predictor-corrector method)
Predictor-corrector methods are methods which require function values at xn, xn-1,
xn-2, xn-3, for the computation of the function value at xn+1.
A predictor is used to find the value of y at xn+1 and then a corrector formula is used
to improve the value of yn+1.
The following two methods are multi-step methods:
1. Milne’s predictor-corrector method
2. Adam’s Predictor-corrector method.
Milne’s predictor-corrector method
(i) To use Milne’s predictor-corrector method, we need at least 4 values prior to the
required values.
(ii) Knowing four consecutive values of y namely yn-3, yn-2, yn-1, yn, we calculate yn+1
using predictor formula. (Use yn+1 on right hand formula to get, better yn+1 after
correction)
(iii)Predictor formula is used to predict the values of yn+1 at xn+1 and then a corrector
formula is used to improve the values of yn+1.
(iv) In Milne’s predictor corrector method may have been computed by Taylor’s series or
Euler’s method of modified Euler’s method or Runge - Kutta method or Picard’s
method.
Milne’s Predictor Formula:

y 4,p  y 0 
4h
2y1  y2  2y3 .
3
Milne’s Corrector Formula

66
y 4, c  y 2 
h
y2  4y3  y4  where y4  f (x 4 , y4 , p) .
3
Example: 1
dy
Using Milne’s method, compute y(0.8) given that  1  y 2 with y(0) = 1, y(0.2) =
dx
0.2027, y(0.4) = 0.4228 and y(0.6) = 0.6841.
Solution
Given

x y y  1  y 2

x0 0 (y0) y0  1  y 02 =1+12 = 2

x1 0.2 0.2027 y1  1  y12 =1+(0.2027)2 = 1.0410

x2 0.4 0.4228 y2  1  y 22 =1+(0.4228)2 =1.1787

x3 0.6 0.684 y3  1  y32  1  (0.6841) = 1.4681

x4 0.8 ? ?

To find y(0.8):
x4 = 0.8; h=0.2
By Milne’s predictor formula,

y 4,p  y 0 
4h
2y1  y2  2y3 
3
= 1.0239
y4  f (x 4 , y 4 ) = 1 + (1.02398)2 = 2.0480

y4 =2.0480.
By Milne’s Corrector formula,

y 4,c  y 2 
h
y2  4y3  y4 
3

= 0.4228 +
0.2
1.178  4(1.4681)  2.0480
3
y(0.8) = 1.0294.
Example: 2

67
 dy
Given  x 3  y , y(0) = 2, the values of y(0.2) = 2.073, y(0.4) = 2.452 and
dx
y(0.6) = 3.023 are got by Runge - Kutta method of fourth order, Find y(0.8) by Milne’s
predictor -corrector method taking h = 0.2.
Solution
Given
x y y  x 3  y

x0 0 2y0 y0  x 30  y 0 =0+2 = 2

x1 0.2 2.073 y1  x13  y1 =(0.2)3 +2.073 = 2.081

x2 0.4 2.452 y2  x 32  y 2 =(0.4)3 +2.452 = 2.516

x3 0.6 3.023 y3  x 33  y 3 = (0.6)3 + 3.023 = 3.239

x4 0.8 ? ?

To find y(0.8):
x4 = 0.8; h=0.2
By Milne’s predictor formula,

y 4,p  y 0 
4h
2y1  y2  2y3 
3

 2
4(0.2)
2(2.081)  2.516  2(3.239)
3
= 4.1664
 y4  f (x 4 , y 4,p )

= f(0.8, 4.1664)
= (0.8)3 + 4.1664 = 4.6784
y4  4.6784
By Milne’s Corrector Formula,

y 4, c  y 2 
h
y2  4y3  y4 
3

 2.452 
0.2
2.516  4(3.239)  4.6784
3
= 3.79536
The corrected value of y(0.8) = 3.79536.
4.5 Runge - Kutta method
68
 The Taylor’s series method of solving differential equations numerically in restricted
because of the evaluation of the higher order derivatives Runge - Kutta methods of solving
initial value problems do not require the calculations of higher order derivatives and give
greater accuracy.
Second – Order Runge - Kutta method
dy
Consider  f ( x, y), y( x 0 )  y 0 , then the value of y1 is obtained as follows:
dx
 y1  y 0  y where y= k2

 h k 
where k 2  hf x  , y  1 
 2 2
k1  hf (x, y) .
Third – Order Runge - Kutta method
dy
Consider  f ( x, y), y( x 0 )  y 0 , then the value of y is obtained as follows:
dx
y1  y 0  y

where y 
1
k1  4k 2  k 3 
6
where k1  hf (x, y)

 h k 
k 2  hf  x  , y  1 
 2 2
k 3  hf x  h, y  2k 2  k1  .

Fourth – Order Runge - Kutta method

This method is commonly used for solving the initial value problem
dy
 f ( x, y), y( x 0 )  y 0 .
dx
The value of y1  y(x1 ) is obtained as follows:
To find y1
y1  y 0  y

Where y 
1
k1  2k 2  2k 3  k 4 
6
where k1  hf (x 0 , y 0 )

 h k 
k 2  hf  x 0  , y 0  1 
 2 2

 h k 
k 3  hf  x 0  , y 0  2 
 2 2 
69
 k 4  hf x 0  h, y0  k 3  .

Example: 1
dy
Given  x 3  y , y(0) = 2. Compute y(0.2), y(0.4) and y(0.6) by Runge - Kutta
dx
method of fourth order.
Solution
Given y  f (x, y)  x 3  y
Also
x 0 0.2 0.4 0.6

y 2 y1 y2 y3

To find y1:
Fourth order Runge - Kutta formula is
y1  y 0  y

Where y 
1
k1  2k 2  2k 3  k 4 
6
where k1  hf (x 0 , y 0 )

 h k 
k 2  hf  x 0  , y 0  1 
 2 2

 h k 
k 3  hf  x 0  , y 0  2 
 2 2 
k 4  hf x 0  h, y0  k 3 

k1  0.2(0,2)  0.2(03  2)  0.4

 0.2 0.4 
k 2  0.2 0  ,2    0.4402
 2 2 

 0.2 0.4402 
k 3  0.2 0  ,2    0.44422
 2 2 
k 4  0.20  0.2,2  0.44422  0.490444

y 
1
k1  2k 2  2k 3  k 4 
6


1
0.4  2(0.4402)  2(0.44422)  0.490444  0.44321
6
y1 = y0 + y = 2 + 0.44321 = 2.44321
y (0.2)  2.44321 .

70
Similarly, y2 = y(0.4) = 2.99 (k1 = 6.4902, k2 = 0.5430, k3 =0.5483, k4 =0.6111, y = 0.5473)
y3 = y(0.6) = 3.68 (k1 = 0.6108, k2 = 0.6841, k3 =0.6914, k4 =0.7795, y = 0.6902).
Example: 2
dy y 2  x 2
Using Runge - Kutta method of 4th order, solve  2 with y(0) = 1 at
dx y x
x = 0.2.
Solution
y2  x 2
Given y  ; h = 0.2
y2  x

x 0 0.2

y 1 y1

To find y1:
Fourth order Runge - Kutta formula is
y1  y 0  y

where y 
1
k1  2k 2  2k 3  k 4 
6
where k1  hf (x 0 , y 0 )

 h k 
k 2  hf  x 0  , y 0  1 
 2 2

 h k 
k 3  hf  x 0  , y 0  2 
 2 2 
k 4  hf x 0  h, y0  k 3 

 y2  x 2 
k1  0.2 2  = 0.2
 y x 
 0.2 0.2 
k 2  0.2f  0  ,1    0.19672
 2 2 
k 3  0.1967

k 4  0.1891
y1 = y(0.2) = y0 + y = 1 + 0.19598 = 1.19598.
4.6 Solution of Algebraic and Transcendent Equations
In mathematics, we often come across problems of obtaining solutions of equations of
the form f(x) = 0. If f(x) is a polynomial then the equation f(x) = 0 is called an algebraic
equation.

71
 Equations which involve transcendental functions like sin x, cos x, tan x, log x and ex
etc., are called transcendental equations.
x 2  5x  6  0, 2x 2  x  4  0, x 5  x 3  2x  3  0 and some examples of algebraic
equations.
3x  sin x  2  0, log10 x  2x  10, a e x  b sin x  c cos x  d log x  x  0 are some
examples of transcendental equations.
If f(x) = 0 is a quadratic equation ax2 + bx + x = 0, we have a simple formula namely

 b  b 2  4ac
x to find its roots.
2a
However, if f(x) is a polynomial of higher degree or an expression involving
transcendental functions we have no simple formula to find roots.
Due to limitations of analytical methods, formula giving exact numerical values of the
solutions exist any in very simple cases.
Hence, we have to use approximate methods to get solutions with good degree of
accuracy.
We have different methods for obtaining approximate solutions for algebraic and
transcendental equations.
(i) Iterative Method
(ii) Aitken’s 2 Method
(iii) Bisection Method
(iv) Regula – Falsi Method
(v) Newton – Raphson Method.

Iterative Method (or) Method of Successive Approximation (or) fixed Point Method
To solve the equation f(x) = 0 by the iteration method, we start with the
approximation value of the root. The equation f(x) = 0 is expressed as x = (x) is called fixed
point equation.
If x0 is the starting approximate value to the actual root ‘’ of x = (x), be first
approximation is x1 = (x0), second approximation is x2= (x1) and so on.
In general we have xn = (xn-1), n = 1, 2, 3, ... Here xn is the nth iteration and the values
of xn gives the root of the given equation at the nth iteration.
Sufficient Condition for Convergence of Iteration (statement of fixed point Theorem)
Let x =  be a root of the equation f(x) = 0 which is equivalent to x = (x). Let I
be any interval containing the root . If x   1 for all x in I, then the sequence of
72
approximations x 0 , x1 ,, x n to the root , provided the initial approximation x0 is chosen

in I.
Example: 1
Solve the equation x3 + x2 – 1 = 0 for the root by iteration method correct to 4 decimal
places.
Solution
Let f(x) = x3 + x2 -1
f(0) = -1 (negative)
f(1) = 1 (positive)
The root lies between 0 and 1.
Let x0 = 0.5.
Express f(x) = 0 as x = (x).
x3 + x2 -1 = 0
 x3 + x2 = 1
 x2 (x+1) = 1
1
x
x 1
1
 1  x  2
1 
x
x 1
1
 x  x   1  x 

2

3
x    1  x  2
1
2

x  
1
3
21  x  2
Here I = [0, 1].

 0 
1
 0.5  1
2

 1 
1
3
 0.1768  1 .
2 22
The condition of convergence is satisfied.
Iteration formula:
xn =  (xn-1); n = 1, 2, 3, ...
Iteration 1: n = 1; x0 = initial value is 0.5

73
  x1  x 0  
1 1
  0.8165 .
1 x0 1  0.5

Iteration 2: n = 2,

 x 2  x1  
1 1
  0.742 .
1  x1 1  0.8165

Iteration 3: n = 3,

 x 3  x 2  
1 1
  0.7577 .
1 x 2 1  0.742

Iteration 4: n = 4,

 x 4  x 3  
1 1
  0.7543 .
1 x3 1  0.7577

Iteration 5: n = 5,

 x 5  x 4  
1 1
  0.7550 .
1 x 4 1  0.7543

Iteration 6: n = 6,

 x 6  x 6  
1 1
  0.7549 .
1 x5 1  0.7550

Iteration 7: n = 7,

 x 7  x 6  
1 1
  0.7549 .
1 x6 1  0.7549

From 6th and 7th iterations, x = 0.7549.

Example 2:
Find the cube root of 15 correct to 4 decimal places by the iteration method.
Solution
Let x be the cube root of 15.
X3 = 15.
f(x) = x3 – 15
The equation is x3 – 15 = 0
x3 = 15
1
 x  153
 x  (x), which is constant, convergence is not satisfied.

 x 2  x  15

 x  x  .
15 15
x2  x
x x
74
 1

 x  
15
 15 x 2 .
x
3
 1 
 x   15    x 2 .
 2

 x  
15
3
.
2 x2
f(x) = x3 -15
f(2) = 8 - 15 = - 7 (negative)
f(3) = 27 – 15 = 12 (positive)
The root lies between 2 and 3.

Let x0 = 2.5; x  

15
.
x
Iteration formula is xn = (xn-1); n = 1, 2, 3, ...

x1  x 0   2.5 
15
 2.4495 .
2.5

x 2  x1   2.4495 
15
 2.4746 .
2.4495

x 3  x 2   2.4746 
15
 2.4620 .
2.4746

x 4  x 3   2.4620 
15
 2.4683 .
2.4620

x 5  x 4   2.4683 
15
 2.4652 .
2.4683

x 6  x 5   2.4652 
15
 2.4667 .
2.4652

x 7  x 6   2.4667  
15
 2.4659 .
2.4667

x 8  x 7   2.4659 
15
 2.4663 .
2.4659

x 9  x 8   2.4663 
15
 2.4662 .
2.4663

x10  x 9   2.4662 

15
 2.4662 .
2.4662
 x  2.4662 .

75
Newton – Raphson Method (or) Newton’s Method (or) Method of Tangents
This method starts with an initial approximation to the root of the equation. A better
and closer approximation to the root can be found by using an iterative process.
Newton – Raphson formula:
 f x i  
x i1  x i    ; i  0,1,2,
 f x i 
Example: 1
Using Newton – Raphson method, find the root of x3 – 6x + 4 = 0 and correct its 4
decimal points.
Solution
Newton – Raphson formula is
 f x i  
x i1  x i    ; i  0,1,2,...
 f x i 

Let f x   x 3  6x  4 ; f x   3x 2  6
 f 0  4 (positive)
f(1) = -1 (negative)
The root lies between 0 and 1.
Let the initial approximation be x0 = 0.5.
Iteration: 1, i = 0

 f x 0    0.53  6 0.5  4 
x1  x 0     0.5     0.7143 .
 f x 0   3 0.5  6 
2

Iteration: 2, i = 1

 f x1    0.71433  6 0.7143  4 

x 2  x1     0.7143     0.7319 .
 f x1  3 0.7143  6
2
 
Iteration: 3, i = 2

 f x 2    0.73193  6 0.7319  4 
x3  x 2     0.7319     0.7320 .
 f x 2  3 0.7319  6
2
 
Iteration: 4, i = 3

 f x 3    0.73203  6 0.7320  4 
x 4  x3     0.7320     0.7320 .
 f x 3  3 0.7320  6
2
 
The value of x is 0.7320.
Example: 2

76
 Solve by Newton - Raphson method x4 – x – 10 = 0.

Solution
Let f x   x 4  x  10; f x   4x 3  1
f(0) = -10 (negative)
f(1) = -10 (negative)
f(2) = 5 (positive)
The root lies between 0 and 2.
Let the initial approximation be x0 = 1.5.
Newton – Raphson formula is,
 f x i  
x i1  x i    ; i  0,1,2,...
 f x 
i 

Iteration: 1, i = 0

 f x 0    1.54  1.5  10 
x1  x 0     1.5     2.015 .
 f x 0   4 1.5  1 
3

Iteration: 2, i = 1

 f x1    2.0154  2.015  10 

x 2  x1     2.015     1.8741 .
 f x1  4 2.015  1
3
 
Iteration: 3, i = 2

 f x 2    1.87414  1.8741  10 
x3  x 2     1.8741     1.8559 .
 f x 2  4 1.8741  1
3
 
Iteration: 4, i = 3

 f x 3    1.85594  1.8559  10 
x 4  x3     1.8559     1.8556 .
 f x 3  4 1.8559  1
3
 
Iteration: 5, i = 4

 f x 4    1.85564  1.8556  10 
x5  x 4     1.8556     1.8556 .
 f x 4  4 1.8556  1
3
 
Comparing the 4th and 5th iteration, we conclude that x = 1.8556.
Exercises:
1. Find the approximate solution for x = 0.1, x = 0.2 by Picard’s method for the equation
y  x  y, y0  1 . Check the result with exact value.

77
  x  y 2  1, y0  0 by Picard’s method.
dy
2. Find the second approximation for
dx
dy
3. Solve  y 2  x 2 with y(0) = 1. Use Taylor series at x = 0.2 and 0.4. Find x = 0.1.
dx
4. Using Taylor series method find y at x = 0.1 correct to four decimal places from
dy
 x 2  y , y(0) = 1 with h = 0.1. Compute terms up to x4.
dx
5. Using Taylor’s series method, find y(1.1) given y  x  y, y(1)  0
6. Using Taylor’s series method in the first five terms in the expansion , find y(0.1)
dy
correct tot three decimal places, given that  e x  y 2 , y(0) =1
dx
7. By Taylor’s series method find y(0.1) given that y  y  xy , y(0) = 1, y(0)  0 .
8. Using Euler’s method find y(0.3) of y(x) satisfies the initial value problem.
dy 1 2
 ( x  y) y 2 , y(0.2)  1.1114 .
dx 2
dy
9. Using Euler’s method find the solution of the initial value problem  log( x  y) ,
dx
y(0) = 2 at x = 0.2 by assuming h = 0.2.
10. Evaluate y(1.2) correct to three decimal places, by the modified Euler method, given

that
dy
dx
 
3
 y  x 2 , y(1)  0 taking h = 0.2.

11. Solve y  1  y , y(0) = 0 by modified Euler method.

dy
12. Using modified Euler’s method find f(0.1) if  x 2  y2
dx
dy
13. Consider the initial value problem,  y  x 2  1 , y(0) = 0.5, using modified Euler
dx
method, find y(0.2).
14. Using Runge - Kutta method of fourth order find y(0.1) and y(0.2) for the initial value
dy
problem  x  y 2 , y(0) = 1.
dx
xy
15. Use the fourth order Runge - Kutta method to compute y for x = 0.1 given y 
1 x2
, y(0) = 1, take h = 0.1 .

16. Find y(0.8) given that y  y  x 2 , y(0.6) = 1.737 by using Runge - Kutta method of
fourth order. Take h = 0.1

78
17. By applying the fourth order Runge - Kutta method find y(0.2) from y  y  x , given
that y(0) = 2 and h = 0.1.
18. Using Milne’s method find y(4.4) given that 5xy  y 2  2  0 given y(4) = 1,
y(4.1) = 1.0049, y(4.2) = 1.009 and y(4.3) = 1.0143.

19. Solve y  x  y 2 , 0  x  1, y(0) = 0, y(0.2) = 0.02, y(0.4) = 0.0795, y(0.6) = 0.1762

by Milne’s method to find y(0.8) and y(1).
dy 1
20. Using Milne’s method find y(2) if y(x) is the solving of  ( x  y) given y(0) = 2,
dx 2
y(0.5) = 2.636, y(1) = 3.3595 and y(1.5) = 4.968.
21. Find real root of the equation x3 + x2 – 100 = 0 correct to 5 decimal places.
22. Solve ex – 3x = 0 by iteration method.
23. Find the negative root of the equation x3 – 2x + 5 = 0.
24. Use the method of fixed point iteration to solve the equation 3x – log10 x = 6.
25. Solve x = cos x by Newton - Raphson method.
26. Solve ex = 3x by Newton - Raphson method.
27. Solve e-x = sin x by Newton - Raphson method.
28. Find the real root of x3 – 3x – 5 = 0, that lies between 2 and 3, correct to 3 decimal
places by Newton - Raphson method.

*****

79
 UNIT V
SIMULTANEOUS LINEAR ALGEBRAIC EQUATIONS
5.1 Introduction
5.2 Difference Methods of Obtaining the Solution
5.3 Gauss Elimination Method
5.4 Gauss Jordan Method
5.5 Method of Factorization
5.6 Crout’s Method

5.1 Introduction

Algebraic Equation

An expression of the form f x   a 0 x n  a1x n1    a n1x  a n where a 0 , a1 ,, a n are

constants with a0  0 and n is a positive integer is called a polynomial in x of degree n.

The polynomial f(x) = 0 is called an algebraic equation of degree n.

Transcendental Equation

If f(x) contains some other function. Such as trigonometric, logarithmic and

exponential etc. Then f(x) is called transcendental equation.

The value of x which satisfies f(x) = 0 is called its root.

Solution

The process of finding the roots of an equation is known as the solution of that
equation.

We shall discuss some numerical methods for the solution of algebraic and
transcendental equation.

5.2 Difference Methods of Obtaining the Solution

80
Simultaneous Linear Algebraic Equations

Simultaneous linear algebraic equations occur in various fields of Science and

Engineering. We know that a given system of linear equations can be solved by applying
Cramer’s rule. But this method is found to be impractical for large system of linear
equations, since the calculations are tedious. To solve such equations, there are other
numerical methods, which are particularly suited for computer operations.

To find the solution for the simultaneous linear equations, we have two types of
numerical methods.

(i) Direct Method

(ii) Indirect Method

Direct Method

The following are the direct methods.

(i) Gauss – Elimination method

(ii) Gauss – Jordan Method

Indirect Method (Iterative Method)

The following are the indirect methods.

(i) Gauss – Seidel method

(ii) Gauss – Jacobi method

5.3 Gauss – Elimination Method

This method is the most effective direct solution technique. In this method, consider
the given system of equations to be AX = B.

In Gauss elimination method, we start with the augmented matrix A|B (A with B) of
the given system and transform it to U|K (upper triangular matrix with k – rows) i.e a matrix
in which all elements below the leading diagonal elements are zero by eliminating row
operations. Finally, the solution is obtained by back substitution process.

Principles of Gauss Elimination method

[A, B] → upper triangular matrix (U|K), then find x, y, z by back substitution process.

Example: 1

Solve the equations

2x + y + 4z = 12

8x – 3y + 2z = 20

4x + 11y – z = 33
81
by Gauss elimination method.

Solution

The given system of equations is equivalent to AX = B where

2 1 4  x  12 
     
A   8  3 2  ; X   y  and B   20 
 4 11  1 z  33 
     

2 1 4 12  R1
 
  A | B    8  3 2 20  R2
 4 11  1 33  R
  3

2 1 4 12 
  R  R2  4 R1
  0  7  14  28  2
0 9 R  R3  2 R1
 9 9  3

2 1 4 12 
  9
  0  7  14  28  R3  R3  R2
 0 0  27  27  7
 

= (U|K)

Using Back substitution method,

-27z = -27

z=1

-7y -14z = -28

y=2

 2x + y + 4z = 12

x=3

 x = 3; y = 2 and z = 1.

Example: 2

Solve the system of equations

28x + 4y – z = 32

x + 3y + 10z = 24

2x + 17y + 4z = 35

82
using Gauss – elimination method.

Solution

The given system of equation is equivalent to AX = B where

 28 4  1  x  32 
     
A   1 3 10  ; X   y  and B   24 
 2 17 4  z  35 
     

 28 4  1 32  R1
 
  A | B    1 3 10 24  R2
 2 17 4 35  R
  3

 1 3 10 24 
 
  28 4  1 32  R1  R2
 2 17 4 35 
 

1 3 10 24 
  R  R2  28R1
  0  80  281  640  2
 0 11 R  R3  2 R1
  16  13  3

1 3 10 24 
  11
  0  80  281  640  R3  R3  R2
0 80
 0  54.64  101 

= (U|K).

using back substation method,

-54.64z = -101

z = 1.85

Also,

-80y – 281z = -640

y = 1.5

Also,

x + 3y + 10z = 24

x = 1.

 x = 1; y = 1.5 and z = 1.85.

83
5.4 Gauss – Jordan Method

This method is a modification of Gauss – elimination method. Here, we consider the

given system of equations to be AX = B. In Gauss – Jordan method, we start with the
augmented matrix (A|B) of the given system of equations and transform it to diagonal matrix
of unit matrix by elementary row operations. Finally, the solution is obtained directly
without back substitution process.

i.e.  A | B GJM
( D | K ) or I | K  .

Example: 1

Solve the following equations by the Gauss – Jordan method.

x + 2y + z = 8; 2x + 3y + 4z = 20; 4x + y +2z = 12.

Solution

The given equations is equivalent to AX = B where

1 2 1  x 8
     
A   2 3 4  ; X   y  and B   20 
 4 1 2 z  12 
     

1 2 1 8 
 
  A | B    2 3 4 20 
 4 1 2 12 
 

1 2 1 8 
  R  R2  2 R1
 0 1 2 4  2
 0  7  2  20  R3  R3  4 R1
 

1 0 5 16 
  R  R1  2 R2
 0 1 2 4  1
 0 0  16  48  R3  R3  7 R2
 

 1 0 5 16 
  R
  0  1 2 4  R3  3
0 0 1 3  16
 

1 0 0 1 
  R  R1  5R3
  0 1 0  2 1
 0 0 1 3  R2  R2  2 R3
 

84
 = (D|K).

 x = 1; y = 2 and z = 3.

Example: 2

Solve the following equations by Gauss Jordan method.

10x + y + z = 12

2x + 10y + z = 13 and

x + y +5z = 7.

Solution

The given equations is equivalent to be AX = B where

10 1 1   x 12 
     
A   2 10 1  ; X   y  and B  13  .
 1 1 5 z 7
     

10 1 1 12 
 
  A | B    2 10 1 13 
1 1 5 7
 

1 1 5 7
 
  2 10 1 13  R1  R3
10 1 1 12 
 

1 1 5 7 
  R  R2  2 R1
 0 8 9 1  2
 0  9  49  58  R3  R3  10 R1
 

8 0 49 57 
  R  8R1  R2
 0 8  9 1  1
 0 0  473  473  R3  8R3  9 R2
 

 8 0 49 57 
  R3
  0 8  9  1 R3 
0 0 1  473
 1 

 8 0 0 8
  R  R1  49 R3
  0 8 0 8 1
 0 0 1 1  R2  R2  9 R3
 

85
 = (D|K).

 x = 1; y = 1 and z = 1.

5.5 Method of Triangularization (or) Method of Factorization

This method is also called as decomposition method. In this method, the coefficient
matrix A of the system AX = B is decomposed or factorized into the product of a lower
triangular matrix L and an upper triangular matrix U. We will explain this method in the case
of three equations in three unknowns.

Consider the system of equations

a11x1  a12 x 2  a13 x 3  b1

a 21x1  a 22 x 2  a 23 x 3  b 2 (1)

a 31x1  a 32 x 2  a 33 x 3  b3

This system is equivalent to AX = B (2)

 a11 a12 a13   x1   b1 

     
where A   a 21 a 22 a 23  ; X   x 2  and B   b 2 
a a 33  x  b 
 31 a 32  3  3

Now we will factorize A as the product of lower triangular matrix

1 0 0
 
L   l 21 1 0
l 1 
 31 l 32

and an upper triangular matrix

 u11 u12 u13 

 
U 0 u 22 u 23 
 0 u 33 
 0

So that

LUX = B (3)

Let UX = Y (4)

and hence LY = B (5)

1 0 0   y1   b1 
    
That is,  l 21 1 0 y2    b2  (6)
l 1   y 3   b 3 
 31 l 32
86
  y1  b1 , l 21 y1  y 2  b 2 , l31 y1  l32 y 2  y3  b3

By forward substitution, y1, y2, y3 can be found out if L is known. From equation (4),

 u11 u12 u13   x1   y1 

    
 0 u 22 u 23   x 2    y 2 
 0 u 33   x 3   y 3 
 0

i.e u11x1  u12 x 2  u13 x 3  y1

u 22 x 2  u 23 x 3  y 2

u 33 x 3  y3

From these x1, x2, x3 can be solved by back substitution, since y1, y2, y3 are known if
U is known.

Now L and U can be found from LU = A.

1 0 0   u11 u12 u13   a11 a12 a 31 

    
i.e  l 21 1 0 0 u 22 u 23    a 21 a 22 a 32 
l 1   0 u 33   a 31 a 32 a 33 
 31 l 32 0

 u11 u12 u13   a11 a12 a13 

   
i.e  l 21u11 l 21u12  u 22 l 21u13  u 23    a 21 a 22 a 23 
l u l31u13  l32 u 23  u 33   a 31 a 32 a 33 
 31 11 l31u12  l32 u 22

Equating corresponding coefficients we get nine equations in nine unknowns. From

these 9 equations, we can solve for 3l’s and 6 u’s.

That is, L and U are known. Hence X is found out. Going into details, we get
u11  a11 , u12  a12 , u13  a13 . That is the elements in the first row of U are same as the
elements in the first of A.

Also, l 21 u11  a 21 , l 21 u12  u 22  a 22 , l 21 u13  u 23  a 23

a 21 a a
 l 21  , u 22  a 22  21 a12 and u 23  a 23  21 a13
a11 a11 a11

Again, l31 u11  a 31 , u12 u12  l32 u 22  a 32 and l31 u13  l32 u 23  u 33  a 33

a 31
a 32  a 12
a 31 a 11
Solving, l 31  , l 32 
a 11 a
a 22  21 a 12
a 11

87
  a 31 
 a 32  a12 
u 33  a 33 
a 31
a13  
a11   a  a 21 a 
 a 21   23 a 11 13 
a12  
a 11
 a 22 
 a11 

Therefore L and U are known.

Note:

In selecting L and U we can also take as

 l11 0 0  1 u12 u13 

   
L   l 21 l 22 0  and U   0 1 u 23 
l l33  0 0 1 
 31 l 32 

So that A = LU.

Example: 1

Solve by triangularization method, the following system x + 5y + z = 14,

2x + 5y + z = 13 and 3x + y + 4z = 17.

Solution

This is equivalent to

 1 5 1   x  14 
    
 2 1 3   y   13 
 3 1 4   z  17 
    

i.e AX = B

1 0 0   u11 u12 u13   1 5 1 

    
let LU   l 21 1 0  0 u 22 u 23    2 1 3 
l 1   0 u 33   3 1 4 
 31 l32 0

By seeing, we can write u11  1, u12  5, u13  1 .

88
 1 0 0 1 5 1  1 5 1
    
 l 21 1 0   0 u 22 u 23    2 1 3 
l 1   0 0 u 33   3 1 4 
 31 l 32

Hence, l 21  2, 5l 21  u 22  1, l 21  u 23  3 .

l 21  2, u 22   9, u 23  1

Again, l31  3, 5l31  l32 u 22  1, l31  l32 u 23  u 33  4

1  15 14 14 5
l32   ; u 33  4  3   
9 9 9 9

LUX = B implies LY = B where UX = Y

LY = B gives,

 
 1 0 0   y1  14 
 2 1 0   y   13 
  2  
1   y 3  17 
14
3
 9 

14
i.e y1  14, 2 y1  y 2  13, 3y1  y 2  y 3  17
9

5
 y1  14, y 2  15, y 3  
3

UX = Y implies,

   
1 5 1   x   14 
 0  9 1   y     15 
  
5    5 
 0 0   z    
 9  3

(i.e) x + 5y + z =14

-9y + z = -15

5 5
 z
9 3

 z= 3, y = 2, x = 1.

Example: 2

Solve the following system by triangularization method: x + y + z = 1,

4x + 3y – z = 6, 3x + 5y + 3z = 4.
89
Solution

1 1 1  x 1
     
Here A   4 3  1 , X   y  , B   6 
3 5 3  z  4
     

1 0 0   u11 u12 u13   1 1 1 

    
LU   l 21 1 0  0 u 22 u 23    4 3  1
l 1   0 u 33   3 5 3 
 31 l32 0

 u11  u12  u13  1 .

l 21u11  4, l 21u12  u 22  3, l 21u13  u 23  1

 l 21  4, u 22  1, u 23  5

l31  3, l31  l32 u 22  5, l31  l32 u 23  u 33  3

l32  2, u 33  10

Now, LUX = B implies LY = B where UX = Y

 1 0 0   y1   1 
    
 4 1 0  y2    6
 3  2 1  y   4
  3  

y1  1, 4y1  y 2  6, 3y1  2y 2  y3  4

 y1  1, y 2  2, y3  5

UX = Y gives,

1 1 1   x  1
    
 0 1  5   y    2 
 0 0  10   z   5 
    

x+y+z=1

-y -5z = 2

-10z = 5

1 1
Hence, z   , y  , x  1 .
2 2

5.6 Crout’s Method (Direct Method)

90
 This is also a direct method. Here also, we decompose the coefficient matrix A as LU
and proceed. But we will follow a different technique as suggested by Crout.

As in the previous method, we want it solve the system

AX = B (1)

 a11 a12 a13   x1   b1 

     
where A   a 21 a 22 a 23  , X   x 2  and B   b1 
a a 33  x  b 
 31 a 32  3  3

Suppose we decompose A = LU (2)

 l11 0 0  1 u12 u13 

   
where L   l 21 l 22 0  and U   0 1 u 23 
l l33  0 0 1 
 31 l 32 

Since AX = B, LUX = B

LY  B where UX  Y (3)

LU = A reduces to

 l11 0 0  1 u12 u13   a11 a12 a 31 

l 0  0 1 u 23   a 21 a 22 a 32 
 21 l 22
l31 l32 l33  0 0 1  a 31 a 32 a 33 

 l11 l11 u12 l11 u13   a11 a12 a13 


i.e l 21 l 21u12  l 22 l 21u13  l 22 u 23   a 21 a 22 a 23 
l 31 l31u12  l 32 l31u13  l32 u 23  l 33  a 31 a 32 a 33 

Equating coefficients and simplifying as in the previous method, we have

l11  a11 , l 21  a 21 , l31  a 31

a12 a
u12  , u13  13
a11 a11

l 22  a 22  l 21u12 , l32  a 32  l31u12

a 23  l 21u13
u 23  , l33  a 33  l31u13  l32 u 23
l 22

Now L and U are known.

Since LY = B, we get

91
 l11 0 0   y1   b1 
   
 l 21 l 22 0   y2   b2 
l l 33   y 3   b 3 
 31 l 32

Multiplying and equating coefficients,

l11y1  b1

l 21y1  l 22 y 2  b 2

l31y1  l32 y 2  l33 y3  b3

Therefore,

b1
y1 
a 11

b 2  l 21 y1
y2 
l 22

b 3  l31 y1  l 32 y 2
y3 
l33

 y1 
 
Knowing Y   y 2  , L and U.
y 
 3

X can be found out from UX = Y.

Note:

Computation scheme by Crout’s method:

We write down the 12 unknowns l11 , l 21 , l 22 , l31 , l32 , l33 , u12 , u13 , u 23 , y1 , y 2 , y3 as a
matrix below called auxiliary matrix or derived matrix.

 l11 u12 u13 y1 

 
Derived matrix =  l 21 l 22 u 23 y2  .
l y 3 
 31 l 32 l 33

If we know the derived matrix, we can write L, U and Y. The derived matrix is got as
explained below, using the augmented matrix (A, B).

(i) The first column of D. M (derived matrix) is the same as the first column of A.
(ii) The remaining elements of first row of D. M. Each elements of the first of D. M.
(except the first elements l11) is got by dividing the corresponding element in
(A, B) by the leading diagonal element of that row.

92
 (iii) Remaining elements of second column of D. M.
Since l 22  a 22  l 21u12 , l32  a 32  l31u12
Each element of second column except the top element = corresponding elements
in (A, B) minus the product of the first element in that row and in that column.
(iv) Remaining elements of second row.
Each element = corresponding elements in (A, B) minus sum of the inner products
of the previously column divided by diagonal element in that row.
(v) Remaining element of third column.
l33  a 33  l31u13  l32 u 23
The element = corresponding element of (A, B) – (sum of the inner products of
the previously calculated elements in the same row and column).
(vi) Remaining element of the third row.
b  l31 y1  l32 y 2 
y3  3
l33
The element = corresponding element of (A, B) – (sum of the inner products of
the previously calculated elements in the same row and column divided by the
diagonal element in that row.

Example: 1
By Crout’s method, solve the system: 2x + 3y + z = -1, 5x + y + z = 9 and
3x + 2y + 4z = 11.

Solution
2 3 1  1
Augmented matrix = (A, B) = 5 1 1 9 
3 2 4 11 
 l11 u12 u13 y1 
 
Let the derived matrix be D. M =  l 21 l 22 u 23 y2 
l y 3 
 31 l 32 l 33

Step: 1
2 . . .
Elements of the first column of D. M are 5 . . .
3 . . .
Step: 2
Elements of first row:
a 3
u12  12 
l11 2
a 1
u13  13 
l11 2
b1 1
y1  
l1 2

93
  3 1 1
2  
2 2 2
D. M. = 5 . . . 
 
3 . . . 
 
Step: 3
Elements of second column:
l 22  a 22  u12 l 21
3 13
 1 5  
2 2
l32  a 32  l31u12
3 5
 2  3 
2 2
 3 1 1
2  
2 2 2
 13 
D. M. = 5  . . 
 2 
3  5 . . 
 2 
Step: 4
Elements of second row:
a  u13l31
u 23  23
l 22
1
1 5
 2 3
13 13

2
 1
9  5  
y2   2    23
13 13

2
 3 1 1
2  
2 2 2
 13 3 23 
D. M. = 5   
 2 13 13 
3  5 . . 
 2 
Step: 5
 1   5  3  3 15 40
l33  4  3         4   
 2   2   13  2 26 13
 3 1 1
2 
2 2 2
 13 3 23 
D. M. = 5   
 2 13 13 
3  5 40
. 
 2 13 

94
Step: 6
 1   5   23 
11  3          
y3   2   2   13   21
40 8
13
 3 1 1
2 
2 2 2
 13 3 23 
D. M. = 5   
 2 13 13 
3  5 40 21 
 2 13 8 
The solution is got from UX = Y
1 u12 u13   x   y1 
i.e 0 1 u 23   y    y 2 
0 0 1   z   y 3 
 3 1  1
1 
2 2
 x   2 
 3
 y    23 
0 1 
   13 
 13 
0 0  z   21 
1
   8 
21 3 23 3y 1 1
z  , y  z   , x   
8 13 13 2 z 2
23 3  21  19
y       
13 13  8  8
3  19  1  21  1 7
x         
2 8  2 8  2 4
7 19 21
x  , y   , z  .
4 8 8
Example: 2
Solve by Crout’s method, the following: x + y + z = 3, 2x – y + 3z = 16 and
3x + y – z = -3.

Solution
1 1 1 3
Here, (A, B) = 2  1 3 16 

3 1  1  3
 l11 u12 u13 y1 
Let the derived matrix be D. M. = l 21 l 22 u 23 y 2 
l31 l32 l33 y 3 
Step: 1
1 . . .
Elements of the first column of D. M. are = 2 . . .
3 . . .
95
Step: 2
Elements of first row of D. M:
1 1 3
u12   1, u13   1, y1   3
1 1 1
1 1 1 3
 D.M.  2 . . . 
3 . . . 
Step: 3
Elements of second column:
l 22  a 22  u12 l 21  1  2  3
l32  a 32  u12 l31  1  1 3  2

1 1 1 3
 D.M.  2  3 . . 
3  2 . . 
Step: 4
Elements of second row:
3  1 2 1
u 23  
3 3
16  3  2 10
y2  
3 3
1 1 1 3 
 1 10 
 D.M.  2  3   
3  2 .3 3
. 

Step: 5
Elements of third column:
 1
l33  1  13      2  
14
 3 3
Step: 6
Elements of third row:
 10 
 3  3 33   2  
y3   3  4
14

3
 
1 1 1 3 
 1 10 
 D.M.  2  3   
 3 3
3  2  14
4 
 3 
The solution is got from UX = Y, i.e.

96
1 1 1   x   3 
 1     10 
0 1  3   y     3 
0 0 1   z   4 
    
x+y+z=3
1 10
y z  
3 3
z = 4.
By back substitution, z= 4, y = -2, x =1.

Exercises:

1. Solve the following system of equations by Gauss elimination method.

3x + 4y + 5z = 18
2x – y + 8z = 13
5x – 2y +7z = 20.
2. Solve the following system of equations by Gauss elimination method:
5x1  x2  x3  x4  4
x1  7 x2  x3  x4  12
x1  x2  6 x3  x4  5
x1  x2  x3  4 x4  6 .
3. Solve the following system of equations by Gauss elimination method.
x1  2 x2  12 x3  8x4  27
5x1  4 x2  7 x3  2 x4  4
 3x1  7 x2  9 x3  5x4  11
6 x1  12 x2  8x3  3x4  49 .
4. Solve the following equations by Gauss Jordan method.
2x + y + 4z = 12; 8x – 3y + 2z = 20; 4x + 11y – z = 33.
5. Solve the following equations by Gauss Jordan method.
2 x1  7 x2  4 x3  9
x1  9 x2  6 x3  1
 3x1  8x2  5x3  6 .
6. Solve the following system of equations by triangularization method:
(i) x + y + 5z = 16, 2x + 3y + z = 4 and 4x + y – z = 4
(ii) x – y + z = 6, 2x + 4y + z = 3 and 3x + 2y -2z = -2
(iii) 2x + y + z = 12, 8x – 3y + 2z = 20 and 4x + 11y – z = 33.
7. Using Crout’s method, solve the following muster of equation:
(i) x + y + 2z = 7, 3x + 2y + 4z =13 and 4x + 8y + 2z = 8
(ii) 2x + 4y + z = 5, 4x + 4y + 3z = 8 and 4x + 8y + z = 9
(iii) 2x – 6y + 8z = 24, 5x + 4y – 3z = 2 and 3x + y + 2z = 16.

*******

97