Numerical Analysis

References:
• Dr. Nanigopal Dutta and Dr. Rabindra Nath Jana: Introductory Numerical
Analysis, Shreedhar Prakashani.
• Amritava Gupta and Subhash Chandra Bose: Introduction to Numerical
Analysis, Academic Publisher.

ge
le
1. Numerical Error

ol
1.1. Numerical Error. The difference between the exact value and the computed value

C
of a numerical computation is called the numerical error of the computed value.

aj
Example: If we like to compute the decimal representation of 1/3 then it comes to

R
0.33333.... So in reality we consider the value upto some decimal places. If we consider
5 decimal places then the computed value of 1/3 becomes 0.33333. Hence the error in

l
da
this computation is 31 − 0.33333.

ha
Note: In any numerical computation, if error occurs, we can not measure the exact
value of it. Because if we can measure the exact value of error then adding it with the
is
computed value we can get the exact value of the computation and hence the computation
ah
involves no error. In the above example we can not determine the exact value of the
M

error- 13 − 0.33333. But we can determine an upper bound of the error of any numerical
::

computation.
ti
ai

1.1.1. Classification of error. The errors in numerical computation may occurs in three
difference ways and accordingly the error can be classified as
M

Inherent error: Error occurs due to erroneous input data in a numerical com-
ar

putation is called inherent error.

um

Example: Let A=0.333333, B=0.666666 instead of 1/3 and 2/3 respectively. Then
computed value of A+B is 0.999999, but exact value is 1. Hence an error occurs
K

in this computation which is inherent error.

Round-off error: Error occurs due to finite precision of an infinite decimal
as

representation during a numerical computation is called round-off error.

an

Example: Let us compute the value of 1/3. Its value is 0.33333..., which consists
of infinite decimal representation. So if we take its value upto 6 decimal places
M

it becomes 0.333333. The error occurs due to this computation is 31 − 0.333333,

which is the round-off error.
r.

Truncation error: Error occurs due to the replacement of an infinite or infini-

D

tesimal process by a finite one is called truncation error.

Example: Let us compute the value of e using Taylor series expansion. We know
2 3 4
that ex = 1 + x + x2! + x3! + x4! + ... and hence e = 1 + 1 + 2!1 + 3!1 + 5!1 + 6!1 + .... Since
it is an infinite series, we shall continue to sum different terms and when value
of a term becomes negligible we stop computation and final value of the sum is
considered as an approximate value of e. If we neglect the terms if its value less
than 0.000001 then
e = 1 + 1 + 0.5 + 0.16666667 + 0.04166667 + 0.00833333 + 0.00138889 +
0.00019841 + 0.00002480 + 0.0000028 + 0.00000028 = 2.71828185
1
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 2

The error, e-2.71828185, of this computation is a truncation error.

1.1.2. Rules for rounding off a number. The goal of any numerical computation is to
reduce the error. Hence a number is rounded off in such a way that the error at
the last decimal place should not exceed 1/2. Hence to round off a number N =
dm dm−1 ...d2 d1 d0 .d−1 d−2 ...d−p , upto n decimal places, discard all the digits after n-th dec-
imal place. So it becomes R = dm dm−1 ...d2 d1 d0 .d−1 d−2 ...d−n . Put a decimal point before

ge
the discarded digits and get the discarded number D = 0.d−(n+1) d−(n+2) ...d−p .

le
• If D < 0.5, then the rounded-off number is Ro = R = dm dm−1 ...d2 d1 d0 .d−1 d−2 ...d−n

ol
• If D > 0.5, then add 1 to dn to get the rounded-off number, i.e., rounded-off

C
number is Ro = dm dm−1 ...d2 d1 d0 .d−1 d−2 ...(d−n + 1)
• If D = 0.5, then add 1 to dn if it is odd to get the rounded-off number, otherwise

aj
R is the rounded-off number., i.e.,

R

dm dm−1 ...d2 d1 d0 .d−1 d−2 ...(d−n + 1) if d−n is odd

l
Ro =

da
dm dm−1 ...d2 d1 d0 .d−1 d−2 ...d−n if d−n is even

ha
Examples: The following table illustrates the rounding-off of some numbers upto 4
decimal places.
is
ah
Input Discarded Remaining Rounded off
Remark
Number(N ) Number(D) Number(R) Number(Ro )
M

Since D > 0.5, 1 is added to

::

325.0234553 0.553 325.0234 325.0235

the last digit of R to get Ro
ti
ai

Since D = 0.5, and last digit of

325.02345 0.5 325.0234 325.0234
M

R, 4 is even so Ro = R
Since D = 0.5, and last digit of
ar

325.02355 0.5 325.0235 325.0236 R, 5 is odd, 1 is added to the

um

last digit of R to get Ro

325.0235455 0.455 325.0235 325.0235 Since D < 0.5, we have Ro = R
K
as

1.2. Significant digit/figure. The digits which are used to represent a number are
called significant digits or significant figures where 0 will be considered as a significant
an

digit if it is not used to fix the position of the decimal point.

M

The following statements describe the identification of significant digits

r.

• All non-zero digits are significant.

D

• Zeros before first non-zero digit are not significant.

• If decimal point exists, then zeros after first non-zero digit are all significant.
• If decimal point does not exists, then zeros after the last non-zero digit will not
be considered as significant.
Examples:
• In the number 12.3450 all the digits are significant digit. Here 0 is a significant
digit because it indicates the number is correct upto 4 decimal places, i.e., if it is
a computed value of a quantity then maximum error in it is 0.00005.
• The number 0.003450 have four significant digits. Here the zeros between decimal
point and 3 are not significant as they are used to fix the position of decimal point.
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 3

• The number 2500 has two significant digits as it can be written as 25 × 102 and
hence zeros are used to represent the position of the decimal point.
• The number, 2500., has four significant digits as here zeros are used to indicate
that the maximum error in the value is 0.5.

1.2.1. Rules for rounding off a number upto n significant figures. To round off a number
N = dm dm−1 ...d2 d1 d0 .d−1 d−2 ...d−p , upto n significant figures at first determine the n-

ge
th significant digit. let it be d−l . Then discard all the digits after d−l . So it becomes

le
R = dm dm−1 ...d2 d1 d0 .d−1 d−2 ...d−l . Put a decimal point before the discarded digits and

ol
get the discarded number D = 0.d−(l+1) d−(n+2) ...d−p .

C
• If D < 0.5, then the rounded-off number is Ro = R = dm dm−1 ...d2 d1 d0 .d−1 d−2 ...d−l
• If D > 0.5, then add 1 to d−l to get the rounded-off number, i.e., rounded-off

aj
number is Ro = dm dm−1 ...d2 d1 d0 .d−1 d−2 ...(d−l + 1)

R
• If D = 0.5, then add 1 to dl if it is odd to get the rounded-off number, otherwise

l
R is the rounded-off number., i.e.,

da


ha
dm dm−1 ...d2 d1 d0 .d−1 d−2 ...(d−l + 1) if d−l is odd
Ro =
dm dm−1 ...d2 d1 d0 .d−1 d−2 ...d−l if d−l is even
is
ah
Examples: The following table illustrates the rounding-off of some numbers upto 5
significant figures.
M
::

Input Discarded Remaining Rounded off

Remark
ti

Number(N ) Number(D) Number(R) Number(Ro )

ai

Since D > 0.5, 1 is added to

5.0234553 0.553 5.0234 5.0235
M

the last digit of R to get Ro

Since D = 0.5, and last digit of
5.02345 0.5 5.0234 5.0234
ar

R, 4 is even so Ro = R
um

Since D = 0.5, and last digit of

5.02355 0.5 5.0235 5.0236 R, 5 is odd, 1 is added to the
K

last digit of R to get Ro

5.0235455 0.455 5.0235 5.0235 Since D < 0.5, we have Ro = R
as

Since D = 0.5, and last digit of

an

0.00354495 0.5 0.0035449 0.0035450 R, 9 is odd, 1 is added to the

last digit of R to get Ro
M

Since D = 0.5, and last digit of

0.00354465 0.5 0.0035446 0.0035446
R, 6 is even so Ro = R
r.
D

1.2.2. Absolute error. The absolute error in a computed value, Vc , is denoted by EA and
is the difference between the exact value, Ve , of the quantity and the computed value, i.e.,
EA = |Ve − Vc |

1.2.3. Relative error. The relative error in a computed value Vc , is denoted by ER and is
obtained by dividing the absolute error in Vc by the magnitude of the exact value, Ve , of
the quantity, i.e.,
Ve − Vc
ER =
Ve
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 4

1.2.4. Percentage error. The percentage error in a computed value (Vc ) is denoted by Ep
and is obtained by multiplying the relative error in Vc by 100, i.e.,
Ve − Vc
EP = × 100 %
Ve
where Ve is the exact value of the quantity.

ge
Example: If the quantity 1/64 is rounded to 4 significant figures then determine the
absolute, relative and percentage errors in the computed value.

le
ol
Ans: Here exact value of the quantity is Ve = 1/64 = 0.015625. Computed value, i.e.,

C
rounded-off value Vc = 0.01562 and hence

aj
Absolute error EA = |Ve − Vc | = 0.015625 − 0.01562 = 0.000005

R
Ve − Vc 0.015625 − 0.01562
Relative error ER = = = 0.00032
Ve 0.015625

l
da
Ve − Vc 0.015625 − 0.01562
Percentage error ER = × 100 % = × 100 % = 0.032%

ha
Ve 0.015625

is
Example: Three approximate values of the number 1/3 are given as 0.332, 0.333 and
ah
0.334. Which of these three is the best approximation.
M

Ans: Here exact value of the quantity is Ve = 1/3 = 0.33333. Let three approximate
::

values be Vc1 = 0.332, Vc2 = 0.333, Vc3 = 0.334 and respective absolute errors be EA1 ,
ti

EA2 , EA3 . Then

ai
M

EA1 = |Ve − Vc1 | = 0.33333 − 0.332 = 0.00133 = 133 × 10−5

EA2 = |Ve − Vc2 | = 0.33333 − 0.333 = 0.00033 = 33 × 10−5
ar

EA3 = |Ve − Vc1 | = 0.33333 − 0.334 = 0.00067 = 67 × 10−5

um

Since absolute error in Vc2 (= 0.333), EA2 (= 33 × 10−5 ) is least, so 0.333 is the best
K

approximation of 1/3.
as

Example: Find the relative error in the number 6.84, if all the digits are correct.
an

Ans: In the quantity 6.84, absolute error EA = 0.005 and hence relative error ER =
EA /6.84 = 0.005/6.84 = 0.000731. So 731 × 106 is the the relative error in 6.84.
M

1.3. Propagation of Round-off error. A round-off error is generated when a real

r.

number is approximated/converted to a machine number with finite size. Let x∗ and y ∗

D

be two real numbers whose approximate representation x and y involves round-off errors
1 and 2 respectively, i.e., x∗ = x + 1 and y ∗ = y + 2 . Propagation of round-off error
under binary operations on these two numbers are presented below.
Addition. We have
x∗ + y ∗ = x + 1 + y + 2
⇒ (x∗ + y ∗ ) − (x + y) = 1 + 2
So, the propagated round-off error in the sum of two approximate numbers is the
sum of the round-off errors in the individual numbers.
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 5

Multiplication. We have
x∗ .y ∗ = (x + 1 ).(y + 2 )
⇒ x∗ .y ∗ − x.y = y1 + x2 + 1 2
= y1 + x2 [Since 1 and 2 are small neglecting 1 .2 ]
x∗ .y ∗ − x.y 1 2
⇒ = +
xy x y

ge
So, the propagated round-off error in the product of two approximate numbers is

le
the sum of the relative errors in the individual numbers.

ol
Division. We have

C
x∗ x x + 1 x
− = −

aj
y∗ y y + 2 y

R
 
1 2
y1 − x2 x −
x y

l
= =

da
y(y + 2 ) (y + 2 )
x∗ x 1 2

ha
y∗
− y x
− y 1 2 2
⇒ x = 2 ' − [assuming is negligible ]
1+ y x y y
y
is
ah
So, the propagated relative error in the ratio of two approximate numbers is the
difference of the relative errors in the individual numbers.
M

Theorem: If a number is correct upto n significant figures and the first significant digit
::

1
ti

of the number is k then the relative error is less than .

k × 10n−1
ai

Proof: Let N be a number whose approximate value Na has n correct significant digits
M

but has m correct decimal figures. Now three cases may arise (i) n¿m, (ii)
ar
um
K
as
an
M
r.
D
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 6

2. Interpolation
Interpolation literally means reading in between the lines. The problem of interpolation
can be stated as follows
Let a function f (x) be not known but a table of values of it is known, i.e., it assumes values y0 ,
y1 , y2 , ..., yn at n + 1 respective values x0 , x1 , x2 , ..., xn of x. Goal is to estimate the value
of f (x) at a non-tabulated point x ∈ (x0 , xn ) under a error bound between the estimated and the

ge
true value.

le
To achieve the goal we estimate a polynomial function φ(x) through the points (x0 , y0 ),

ol
(x1 , y1 ), (x2 , y2 ), ..., (xn , yn ), which are called interpolating points, in such a manner that

C
the function φ(x) passes through the said interpolating points. Since φ(x) passes through
n + 1 known points so we can determine at most n + 1 unknown parameters of φ(x)

aj
and hence degree of φ(x) is n [since a general n-degree polynomial has n + 1 unknown

R
parameters]. Then φ(α) is considered as the estimate of f (x) at x = α. The following

l
theorem justified the existence of such a polynomial.

da
Weirstrass Theorem. If a function f (x) is continuous in an interval [a,b] then for any

ha
preassigned  > 0, there exists a polynomial φ(x) such that
is
|f (x) − φ(x)| < , ∀ x ∈ (a, b)
ah
:: M
ti
ai
M
ar
um
K

Figure 1. Pictorial representation of Weirstrass Theorem

as
an

2.1. Lagrange’s Interpolation Formula. Let f (x) be a function which has values y0 ,
y1 , y2 , ..., yn at n + 1 interpolating points x0 , x1 , x2 , ..., xn of x. We like to derive a
M

polynomial φ(x) of degree n which agrees with f (x) at the interpolating points, i.e.,
r.

φ(xi ) = f (xi ) = yi , i = 0, 1, 2, ..., n (1)

D

Let us consider
φ(x) = a0 (x − x1 )(x − x2 )...(x − xn ) + a1 (x − x0 )(x − x2 )(x − x3 )...(x − xn ) +
+a2 (x − x0 )(x − x1 )(x − x3 )(x − x4 )...(x − xn ) + ...
+ai (x − x0 )...(x − xi−1 )(x − xi+1 )...(x − xn ) + ...
+an−1 (x − x0 )(x − x1 )...(x − xn−2 )(x − xn ) + an (x − x0 )(x − x1 )...(x − xn−1 )
Xn
= ai (x − x0 )(x − x1 )...(x − xi−1 )(x − xi+1 )...(x − xn ) (2)
i=0
where ai , i = 0, 1, 2, ..., n are n + 1 unknown parameters to be determined
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 7

Putting x = xi in (2) for i = 0, 1, 2, , ..., n successively we get

φ(xi ) = ai (xi − x0 )(xi − x1 )...(xi − xi−1 )(xi − xi+1 )...(xi − xn )
⇒ yi = ai (xi − x0 )(xi − x1 )...(xi − xi−1 )(xi − xi+1 )...(xi − xn )
yi
⇒ ai = for i = 0, 1, 2, , ..., n
(xi − x0 )(xi − x1 )...(xi − xi−1 )(xi − xi+1 )...(xi − xn )
Putting these values of ai in (2) we get

ge
n
(x − x0 )(x − x1 )...(x − xi−1 )(x − xi+1 )...(x − xn )

le
X
φ(x) = yi (3)
(x − x )(x − x )...(x − x )(x − x )...(x − x )

ol
i=0 i 0 i 1 i i−1 i i+1 i n

C
which is called Lagrange’s interpolation formula of f (x) for n + 1 interpolating points x0 ,
x1 , x2 , ..., xn and for x ∈ (x0 , xn ).

aj
R
2.1.1. Lagrange’s Interpolation Formula using Lagrangian function. The coefficient of yi
in (3) is called Lagrangian function and is denoted by li (x), i.e.,

l
da
(x − x0 )(x − x1 )...(x − xi−1 )(x − xi+1 )...(x − xn )
li (x) = for i = 0, 1, 2, ..., n (4)

ha
(xi − x0 )(xi − x1 )...(xi − xi−1 )(xi − xi+1 )...(xi − xn )
For ψ(x) = (x − x0 )(x − x1 )...(x − xn ) we have
is
ah
0
ψ (x) = (x − x1 )(x − x2 )...(x − xn ) + (x − x0 )(x − x2 )...(x − xn ) + ... +
M

+(x − x0 )(x − x1 )...(x − xn−2 )(x − xn ) + (x − x0 )(x − x1 )...(x − xn−1 )

::

n
ti

X
= (x − x0 )(x − x1 )...(x − xi−1 )(x − xi+1 )...(x − xn−2 )(x − xn ) (5)
ai

i=0
M

From (5) we have

0
ar

ψ (xi ) = (xi − x0 )(xi − x1 ) ... (xi − xi−1 )(xi − xi+1 ) ... (xi − xn−2 )(xi − xn ) (6)
um

for i = 0, 1, 2, ..., n
From (4) and (6) we have
K

ψ(x)
li (x) = for i = 0, 1, 2, ..., n (7)
as

(x − xi )ψ 0 (xi )
an

From (3) and (7) we have

M

n
X ψ(x)yi
φ(x) = (8)
(x − xi )ψ 0 (xi )
r.

i=0
D

2.1.2. Lagrange’s interpolation formula for equally spaced points. Let f (x) be a function
which has values y0 , y1 , y2 , ..., yn at n + 1 equispaced interpolating points x0 , x1 , x2 , ...,
xn of x, where xi − xi−1 = h for i = 1, 2, ..., n. We like to derive a polynomial φ(x) of
degree n which agrees with f (x) at the interpolating points, i.e.,
φ(xi ) = f (xi ) = yi , i = 0, 1, 2, ..., n (1)
Let us consider the interpolating polynomial φ(x) in the following form:
φ(x) = a0 (x − x1 )(x − x2 )...(x − xn ) + a1 (x − x0 )(x − x2 )(x − x3 )...(x − xn ) +
+a2 (x − x0 )(x − x1 )(x − x3 )(x − x4 )...(x − xn ) + ...
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 8

+ai (x − x0 )...(x − xi−1 )(x − xi+1 )...(x − xn ) + ...

+an−1 (x − x0 )(x − x1 )...(x − xn−2 )(x − xn ) + an (x − x0 )(x − x1 )...(x − xn−1 )
Xn
= ai (x − x0 )(x − x1 )...(x − xi−1 )(x − xi+1 )...(x − xn )
i=0
n  
X (x − x0 )(x − x1 )...(x − xi−1 )(x − xi )(x − xi+1 )...(x − xn )
= ai

ge
i=0
(x − xi )

le
ai ψ(x)
= (2)

ol
(x − xi )

C
where ai , i = 0, 1, 2, ..., n are n + 1 unknown parameters to be determined

aj
and ψ(x) = (x − x0 )(x − x1 )...(x − xn ) (3)

R
Let x = x0 + sh. Then we have

l
xi = x0 + ih = x − sh + ih

da
⇒ x − xi = (s − i)h for i = 0, 1, 2, ..., n (4)

ha
and ψ(x) = hn s(s − 1)(s − 2)...(s − n) (5)
is
Putting x = xi in (2) for i = 0, 1, 2, , ..., n successively we get
ah
φ(xi ) = ai (xi − x0 )(xi − x1 )...(xi − xi−1 )(xi − xi+1 )...(xi − xn )
M

⇒ yi = ai (xi − x0 )(xi − x1 )...(xi − xi−1 )(xi − xi+1 )...(xi − xn )

::

yi
ti

⇒ ai =
ai

(ih)((i − 1)h)...(h)(−h)(−2h)...{−(n − i)h}

M

= (−1)n−i hn i!(n − i)!for i = 0, 1, 2, , ..., n (6)

Using (4), (5), and (6) from (2) we have
ar

n n
um

X hn s(s − 1)(s − 2)...(s − n) X s(s − 1)(s − 2)...(s − n)

φ(x) = yi = yi (7)
i=0
(s − i)hn (−1)n−i i!(n − i)! i=0
(s − i)(−1)n−i i!(n − i)!
K

which is called Lagrange’s interpolation formula of f (x) for n + 1 equispaced interpolating

as

points x0 , x1 , x2 , ..., xn and for x ∈ (x0 , xn ).

an

2.1.3. Uniqueness of interpolating formula. Let φ(x) be the interpolating polynomial of

M

degree n of a function f (x) which assumes values y0 , y1 , ..., yn at n+1 interpolating points
x0 , x1 , ..., xn respectively. If φ(x) is not unique then let φ̄(x) is another interpolating
r.

polynomial of f (x) for the same interpolating points. From these assumptions we have
D

φ(x) and φ̄(x) are distinct interpolating polynomials of degree n for f (x)
⇒ dφ(x) = φ(x) − φ̄(x) is a polynomial of degree n
also φ(xi ) = φ̄(xi ) = yi for i = 0, 1, 2, , ..., n
⇒ dφ(xi ) = 0 for i = 0, 1, 2, , ..., n
⇒ The n degree polynomial dφ(x) vanishes at n + 1 points
⇒ dφ(x) = 0, ∀ x ⇒ φ(x) = φ̄(x), ∀ x
⇒ The interpolating polynomial φ(x) of f (x) is unique.
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 9

2.1.4. Error in an interpolating formula. Let a real valued function f (x) be defined on
0 00
an interval I. Also let f (x) and its derivatives f (x), f (x), ..., f n+1 (x) are all continuous
on I and I contains the interpolating points x0 , x1 , ..., xn . For these assumptions, let
φ(x) be the interpolating polynomial of degree n of f (x) for the interpolating points
x0 , x1 , ..., xn ∈ I and E(x) be the error in approximating f (x) by φ(x) at a point x ∈ I,
i.e., E(x) = f (x) − φ(x) and let ψ(x) = (x − x0 )(x − x1 )...(x − xn ).
f (x) − φ(x) E(x)

ge
Let F (t) = f (t) − φ(t) − kψ(t), where k = = (1)
ψ(x) ψ(x)

le
⇒ F (t) = 0 for t = x, x0 , x1 , ..., xn

ol
C
⇒ F (t) vanishes at n + 2 points in I, i.e., F (t) has at least n + 2 zeros in I
0 00
⇒ F (t) has at least n + 1 zeros in I, F (t) has at least n zeros in I, ..., (2)

aj
n+1
and F (t) has at least one zero in I [by Rolle’s theorem]

R
⇒ ∃ ξ ∈ I, such that F n+1 (ξ) = 0

l
da
⇒ f n+1 (ξ) − φn+1 (ξ) − kψ n+1 (ξ) = 0

ha
E(x)
⇒ f n+1 (ξ) − 0 − (n + 1)! = 0
ψ(x)
is
ah
[Since φ(x) and ψ(x) are polynomials of degree n and n + 1 respectively.]
ψ(x)f n+1 (ξ)
M

⇒ E(x) = (3)
(n + 1)!
::
ti

Theorem: Let f (x) be real-valued, continuous and have continuous derivatives of order
ai

n + 1 for all x in an interval I containing the interpolating points x0 , x1 , ..., xn . Then at

M

any point x ∈ I the error in approximating f (x) by the interpolating polynomial φ(x) is
given by
ar

ψ(x)f n+1 (ξ)

E(x) = f (x) − φ(x) =
um

(n + 1)!
where ξ is a point in the interval I and ψ(x) = (x − x0 )(x − x1 )...(x − xn ).
K

Example: From the following table of data determine Lagrange’s interpolation formula
as

for f (x) and hence find f (1.5).

an

x: 0 1 2 3
M

f (x) : 1 2 11 34
r.

Ans: Let x0 = 0, x1 = 1, x2 = 2, x3 = 3 and y0 = 1, y1 = 2, y2 = 11, y3 = 34, then we

D

have Lagrange’s interpolation formula for f (x) as below:

(x − x1 )(x − x2 )(x − x3 ) (x − x0 )(x − x2 )(x − x3 )
f (x) = y0 + y1
(x0 − x1 )(x0 − x2 )(x0 − x3 ) (x1 − x0 )(x1 − x2 )(x1 − x3 )
(x − x0 )(x − x1 )(x − x3 ) (x − x0 )(x − x1 )(x − x2 )
+ y2 + y3
(x2 − x0 )(x2 − x1 )(x2 − x3 ) (x3 − x0 )(x3 − x1 )(x3 − x2 )
(x − 1)(x − 2)(x − 3) (x − 0)(x − 2)(x − 3)
= 1+ ×2
(0 − 1)(0 − 2)(0 − 3) (1 − 0)(1 − 2)(1 − 3)
(x − 0)(x − 1)(x − 3) (x − 0)(x − 1)(x − 2)
+ × 11 + × 34
(2 − 0)(2 − 1)(2 − 3) (3 − 0)(3 − 1)(3 − 2)
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 10

1 11 34
= − (x − 1)(x − 2)(x − 3) + x(x − 2)(x − 3) − x(x − 1)(x − 3) + x(x − 1)(x − 2)
6 2 6
1 11
= (x − 1)(x − 2){34x − (x − 3)} + x(x − 3){(x − 2) − (x − 1)}
6 2
1 1
= (x − 1)(x − 2)(11x + 1) + x(x − 3)(7 − 9x)
2 2
1

ge
= {(x − 1)(x − 2)(11x + 1) − x(x − 3)(9x − 7)}
2

le
1
= {(x2 − 3x + 2)(11x + 1) − x(9x2 − 34x + 21)}

ol
2

C
1
= (11x3 − 32x2 + 19x + 2 − 9x3 + 34x2 − 21x)
2

aj
1
(2x3 + 2x2 − 2x + 2) = x3 + x2 − x + 1

R
=
2

l
da
Hence putting x = 1.5 in the above formula we have f (1.5) = 5.125.

ha
Example: The function y = f (x) is passes through the points (7, 3), (8, 1), (9, 1),
is
(10, 9). Find the value of y at x=9.5 using Lagrange’s interpolation formula.
ah
Ans: The given data can be tabulated as below
:: M

x : 7 8 9 10
ti

y = f (x) : 3 1 1 9
ai
M

Let x0 = 7, x1 = 8, x2 = 9, x3 = 10 and so f (x0 ) = 3, f (x1 ) = 1, f (x2 ) = 1, f (x3 ) = 9,

then we have Lagrange’s interpolation formula for f (x) as below:
ar
um

(x − x1 )(x − x2 )(x − x3 ) (x − x0 )(x − x2 )(x − x3 )

f (x) = f (x0 ) + f (x1 )
(x0 − x1 )(x0 − x2 )(x0 − x3 ) (x1 − x0 )(x1 − x2 )(x1 − x3 )
K

(x − x0 )(x − x1 )(x − x3 ) (x − x0 )(x − x1 )(x − x2 )

as

+ f (x2 ) + f (x3 )
(x2 − x0 )(x2 − x1 )(x2 − x3 ) (x3 − x0 )(x3 − x1 )(x3 − x2 )
an

(9.5 − 8)(9.5 − 9)(9.5 − 10) (9.5 − 7)(9.5 − 9)(9.5 − 10)

So f (9.5) = ×3+ ×1
M

(7 − 8)(7 − 9)(7 − 10) (8 − 7)(8 − 9)(8 − 10)

(9.5 − 7)(9.5 − 8)(9.5 − 10) (9.5 − 7)(9.5 − 8)(9.5 − 9)
r.

+ ×1+ ×9
(9 − 7)(9 − 8)(9 − 10) (10 − 7)(10 − 8)(10 − 9)
D

0.375 −0.625 1.875 1.875

= + + + × 9 = 3.625
2 2 2 6

Example. A certain function f (x) defined on the interval (0,1) is such that f (0) = 0,
f (1/2) = −1, and f (1) = 0. Find the quadric polynomial p(x) which agrees with f (x)
3
for x = 0, 1/2, 1. Moreover, if | ddxf3 | ≤ 1 for 0 ≤ x ≤ 1 prove that

1
|f (x) − p(x)| ≤ for 0 ≤ x ≤ 1
12
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 11

Ans. The Lagrange interpolation formula p(x) of f (x) through three interpolating points
0, 1/2, and 1 is given by
(x − 1/2)(x − 1) (x − 0)(x − 1) (x − 0)(x − 1/2)
p(x) = f (0) + f (1/2) + f (1)
(0 − 1/2)(0 − 1) (1/2 − 0)(1/2 − 1) (1 − 0)(1 − 1/2)
(x − 0)(x − 1)
= × (−1) = 4x(x − 1)
(1/2 − 0)(1/2 − 1)

ge
Error in approximating f (x) by p(x) is given by

le
000

ol
f (ξ)
E(x) = f (x) − p(x) = (x − 0)(x − 1/2)(x − 1)

C
3!
000
|f (ξ)| 11 1

aj
⇒ |f (x) − p(x)| = |x||(x − 1/2)||(x − 1)| ≤ = , since for 0 ≤ x ≤ 1
6 26 12

R
d3 f
we have |x| ≤ 1, |(x − 1/2)| ≤ 1/2, |(x − 1)| ≤ 1 and | 3 | ≤ 1

l
dx

da
ha
Algorithm 1 Algorithm for Lagrange interpolation formula
is
ah
1: Input: Data set {(x0 , y0 ), (x1 , y1 ), ..., (xn , yn )} and unknown point xu .
2: Output: yu , the value of y at x = xu .
M

3: Sum=0
::

4: for i = 0 to n with increment 1 do

ti

5: P =1
ai

6: for j = 0 to i − 1 with increment 1 do

M

P × (xu − xj )
7: P =
(xi − xj )
ar

8: end for
um

9: for j = i + 1 to n with increment 1 do

P × (xu − xj )
10: P =
K

(xi − xj )
11: end for
as

12: Sum=Sum+P
an

13: end for

M

14: Output: Sum

r.
D

2.1.5. Advantage and limitations of Lagrangian interpolation formula. Major advantages

and limitations of Lagrange’s interpolation formula are listed here.
Advantages. The main advantage of the Lagrangian interpolation formula is
that it is applicable for both equispaced and un-equispaced points. Moreover, the
order of interpolating point is not required in this formula.
Limitations. The major limitation of the Lagrangian interpolation formula is
that it is not able to incorporate additional interpolating point, i.e., if one or more
interpolating points are added to increase the degree of the polynomial then the
the old formula could not be extended, i.e., the interpolation formula to be derived
newly.
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 12

2.2. Divided differences. Let a function f (x) be known at n + 1 points x0 , x1 , ..., xn

then the divided differences of different orders are defined as below:
• Zeroth order divided difference at a point x0 is denoted by f [x0 ] and is defined as
f [x0 ] = f (x0 ) (1)
• First order divided difference for the points x0 , x1 is denoted by f [x0 , x1 ] and is
defined as

ge
f (x0 ) − f (x1 ) f (x0 ) f (x1 )

le
f [x0 , x1 ] = = + (2)
x0 − x 1 x0 − x1 x1 − x0

ol
• second order divided difference for the points x0 , x1 , x2 is denoted by f [x0 , x1 , x2 ]

C
and is defined as

aj
 
f [x0 , x1 ] − f [x1 , x2 ] 1 f (x0 ) − f (x1 ) f (x1 ) − f (x2 )
f [x0 , x1 , x2 ] = = −

R
x 0 − x2 x0 − x2 x0 − x1 x1 − x2

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

da
= + + (3)
(x0 − x1 )(x0 − x2 ) (x1 − x0 )(x1 − x2 ) (x2 − x0 )(x2 − x1 )

ha
• Similarly, the n-th order divided difference for the points x0 , x1 , x2 , ..., xn is de-
noted by f [x0 , x1 , x2 , ..., xn ] and is defined as
is
ah
n
X f (xi )
M

f [x0 , x1 , ..., xn ] = (4)

i=0
(xi − x0 )(xi − x1 )...(xi − xi−1 )(xi − xi+1 )...(xi − xn )
::
ti

2.2.1. Divided difference with equal component. The first order divided difference at the
ai

point x is denoted by f [x, x] and is derived as

M

f (x + ) − f (x1 ) 0 0
f [x, x] = lim f [x, x + ] = lim = f (x) = f [x] (1)
ar

→0 →0 
um

The second order divided difference at the point x is denoted by f [x, x, x] and is derived
as
K

0 f [x + , x + ) − f [x, x]
f [x, x] = lim

as

→0
f [x + , x + ] − f [x + , x] f [x + , x) − f [x, x]
an

= lim + lim
→0  →0 
M

= lim f [x + , x + , x + ] + lim f [x, x, x + ] = 2f [x, x, x] (2)

→0 →0
r.

Similarly, the n + 1-th order divided difference at the point x is denoted by f [x, x, x, ...n +
D

1 times] and is derived as

0
f [x, x, x, ...n times] = nf [x, x, x, ...n + 1 times] (3)
Again second order derivative of f [x] is derived as
d2 f [x]
 
d df [x] d
2
= = f [x, x] = 2f [x, x, x] (4)
dx dx dx dx
Third order derivative of f [x] is derived as
d3 f [x] d d2 f [x]
 
d
3
= 2
= 2f [x, x, x] = 3!f [x, x, x, x] (5)
dx dx dx dx
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 13

Similarly, the r-th order derivative of f [x] is derived as

dr f [x]
= = r!f [x, x, x, ..., r + 1 times] (6)
dxr
2.3. Newton’s fundamental interpolation formula. Let y = f (x) be a function
which takes values y0 , y1 , y2 , ..., yn at n + 1 interpolating points x0 , x1 , x2 , ..., xn , of x.
So xi = x0 + ih, i = 0, 1, 2, ...n. We like to derive a polynomial φ(x) of degree n which

ge
agrees with f (x) at the interpolating points, i.e.,

le
φ(xi ) = f (xi ) = yi , i = 0, 1, 2, ..., n (1)

ol
Let us consider

C
φ(x) = a0 + a1 (x − x0 ) + a2 (x − x0 )(x − x1 ) + a3 (x − x0 )(x − x1 )(x − x2 )

aj
+... + an (x − x0 )(x − x1 )(x − x2 )...(x − xn−1 ) (2)

R
where ai , i = 0, 1, 2, ..., n are n + 1 unknown parameters to be determined

l
Putting x = x0 in (2) we have

da
φ(x0 ) = a0 ⇒ a0 = f (x0 ) = f [x0 ] [using (1)]

ha
Again putting x = x1 in (2) we have
is
ah
φ(x1 ) = a0 + a1 (x1 − x0 )
⇒ f (x1 ) = f (x0 ) + a1 (x1 − x0 ) [using (1)]
M

f (x1 ) − f (x0 )
::

⇒ a1 = = f [x0 , x1 ]
ti

x1 − x0
ai

Again putting x = x2 in (2) we have

M

φ(x2 ) = a0 + a1 (x2 − x0 ) + a2 (x2 − x0 )(x2 − x1 )

ar

⇒ f (x2 ) = f (x0 ) + f [x0 , x1 ](x2 − x0 ) + a2 (x2 − x0 )(x2 − x1 )[using (1)]

f (x2 ) − f (x0 )
um

⇒ = f [x0 , x1 ] + a2 (x2 − x1 )
x2 − x0
K

f [x2 , x0 ] − f [x0 , x1 ]
⇒ a2 = = f [x0 , x1 , x2 ]
x2 − x1
as

Similarly we have
an

ai = f [x0 , x1 , x2 , ..., xi ] for i = 0, 1, 2, ...n

M

Putting these values of ai in (2) we get

r.

φ(x) = f (x0 ) + (x − x0 )f [x0 , x1 ] + (x − x0 )(x − x1 )f [x0 , x1 , x2 ] +

D

(x − x0 )(x − x1 )(x − x2 )f [x0 , x1 , x2 , x3 ] + ...

+(x − x0 )(x − x1 )(x − x2 )...(x − xn−1 )f [x0 , x1 , x2 , ..., xn ]
which is called Newton’s fundamental interpolation formula of f (x) for n+1 interpolating
points x0 , x1 , x2 , ..., xn and for x ∈ (x0 , xn ).
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 14

2.4. Finite difference operators.

2.4.1. Forward difference. The forward difference of a function f (x) with step length h is denoted
by ∆f (x) and is defined as
∆f (x) = f (x + h) − f (x)
The r-th order forward difference of f (x) with step length h is denoted by ∆r f (x) and is defined
as

ge
∆r f (x) = ∆r−1 f (x + h) − ∆r−1 f (x) for r = 1, 2, 3, ...., where
∆1 f (x) = ∆f (x) = f (x + h) − f (x), and hence

le
ol
∆2 f (x) = ∆f (x + h) − ∆f (x) = {f (x + 2h) − f (x + h)} − {f (x + h) − f (x)}

C
= f (x + 2h) − 2f (x + h) + f (x)
∆3 f (x) = ∆2 f (x + h) − ∆2 f (x)

aj
= {f (x + 3h) − 2f (x + 2h) + f (x + h)} − {f (x + 2h) − 2f (x + h) + f (x)}

R
= f (x + 3h) − 3f (x + 2h) + 3f (x + h) − f (x)

l
da
... ... ... ... ... ... ... ...
n
f (x + nh) − n C1 f (x + (n − 1)h) + n C2 f (x + (n − 2)h) − ...

ha
∆ f (x) =
+(−1)r n Cr f (x + (n − r)h) + ... + (−1)n n Cn f (x + nh)
is
ah
Let x0 , x1 , ..., xn be a set of equispaced values of x for which values of f (x) are y0 , y1 ,...,yn
respectively. Then forward differences are calculated in a tabular form as presented below
M

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

::

x0 y0
ti

∆y0 = y1 − y0
ai

x1 y1 ∆2 y0 = ∆y1 − ∆y0
∆y1 = y2 − y1 ∆ 3 y0 = ∆ 2 y1 − ∆ 2 y0
M

x2 y2 ∆2 y1 = ∆y2 − ∆y1 ∆4 y0 = ∆3 y1 − ∆3 y0
∆y2 = y3 − y2 ∆ 3 y1 = ∆ 2 y2 − ∆ 2 y1
ar

x3 y3 ∆2 y2 = ∆y3 − ∆y2 ∆ 4 y1 = ∆ 3 y2 − ∆ 3 y1
... ... ... ... ... ...
um

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

xn−3 yn−3 ∆2 yn−4 = ∆yn−3 − ∆yn−4 ∆4 yn−3 = ∆3 yn−2 − ∆3 yn−3
∆yn−3 = yn−2 − yn−3 ∆3 yn−2 = ∆2 yn−1 − ∆2 yn−2
K

xn−2 yn−2 ∆2 yn−3 = ∆yn−2 − ∆yn−3 ∆4 yn−4 = ∆3 yn−3 − ∆3 yn−4

∆yn−2 = yn−1 − yn−2 ∆3 yn−3 = ∆2 yn−2 − ∆2 yn−3
as

xn−1 yn−1 ∆2 yn−2 = ∆yn−1 − ∆yn−2

∆yn−1 = yn − yn−1
an

xn yn
M

Algorithm 2 Algorithm for determining forward difference table

r.

1: Input: Number of function values n, and the function values {y0 , y1 , y2 , ..., yn }.
D

2: for i = 0 to n with increment 1 do

3: di0 = yi
4: end for
5: for j = 1 to n with increment 1 do
6: for i = 0 to n − j with increment 1 do
7: dij = di+1j−1 − yij−1
8: end for
9: end for
10: Output: dij , for i = 0, 1, 2, ..., n and i-th row has n − i + 1 column suffixes, j = 0, 1, ..., n − i.
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 15

2.4.2. Backward difference. The back difference of a function f (x) with step length h is denoted
by 5f (x) and is defined as
5f (x) = f (x) − f (x − h)
The r-th order backward difference of f (x) with step length h is denoted by ∆r f (x) and is
defined as
5r f (x) = 5r−1 f (x) − 5r−1 f (x − h) for r = 1, 2, 3, ..., where
51 f (x) = 5f (x) = f (x) − f (x − h), and hence

ge
52 f (x) = 5f (x) − 5f (x − h)

le
= {f (x) − f (x − h)} − {f (x − h) − f (x − 2h)}

ol
= f (x) − 2f (x − h) + f (x − 2h)

C
53 f (x) = 52 f (x) − 52 f (x − h)

aj
= {f (x) − 2f (x − h) + f (x − 2h} − {f (x − h) − 2f (x − 2h) + f (x − 3h)}

R
= f (x) − 3f (x − h) + 3f (x − 2h) − f (x − h)

l
... ... ... ... ... ... ... ...

da
n
5 f (x) = f (x) − n C1 f (x − h) + n C2 f (x − 2h) − n C3 f (x − 3h) + ...

ha
+(−1)r n Cr f (x − rh) + ... + (−1)n n Cn f (x − nh)

is
Let x0 , x1 , ..., xn be a set of equispaced values of x for which values of f (x) are y0 , y1 ,...,yn
ah
respectively. Then backward differences are calculated in a tabular form as presented
M

below
::

x y = f (x) 5y = ∆f (x) 52 y = 52 f (x) 53 y = 53 f (x) 54 y = 54 f (x)

ti

x0 y0
ai

5y1 = y1 − y0
x1 y1 52 y2 = 5y1 − 5y0
M

∆y2 = y2 − y1 53 y3 = 52 y3 − 52 y2
x2 y2 52 y3 = 5y3 − 5y2 54 y4 = 53 y4 − 53 y3
5y3 = y3 − y2 53 y4 = 52 y4 − 52 y3
ar

x3 y3 52 y4 = 5y4 − 5y3 54 y5 = 53 y5 − 53 y4
um

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

... ... ... ... ... ...
xn−3 yn−3 52 yn−2 = 5yn−2 − 5yn−3 54 yn−1 = 53 yn−1 − 53 yn−2
K

5yn−2 = yn−2 − yn−3 53 yn−1 = 52 yn−1 − 52 yn−2

xn−2 yn−2 52 y n−1 = 5yn−1 − 5yn−2 54 yn = 53 yn − 53 yn−1
as

5yn−1 = yn−1 − yn−2 53 yn = 52 yn − 52 yn−1

xn−1 yn−1 52 yn = 5yn − 5yn−1
an

5yn = yn − yn−1
xn yn
M

Algorithm 3 Algorithm for determining backward difference table

r.
D

1: Input: Number of function values n and the function values {y0 , y1 , y2 , ..., yn }.
2: for i = 0 to n with increment 1 do
3: di0 = yi
4: end for
5: for j = 1 to n with increment 1 do
6: for i = j to n with increment 1 do
7: dij = dij−1 − yi−1j−1
8: end for
9: end for
10: Output: dij , for i = 0, 1, 2, ..., n and i-th row has i column suffixes, j = 0, 1, ..., i.
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 16

2.4.3. Relation between forward and backward differences. The forward and backward dif-
ferences of a function f (x) with step length h is denoted by ∆f (x) and 5f (x) respectively
and are defined as
∆f (x) = f (x + h) − f (x) and 5 f (x) = f (x) − f (x − h)
⇒ ∆f (x) = f (x + h) − f (x) = 5f (x + h)
⇒ ∆2 f (x) = ∆f (x + h) − ∆f (x) = f (x + 2h) − f (x + h) − {f (x + h) − f (x)}

ge
⇒ ∆2 f (x) = f (x + 2h) − 2f (x + h) + f (x) = 52 f (x + h)

le
⇒ ∆3 f (x) = ∆f (x + 2h) − 2∆f (x + h) + ∆f (x)

ol
= {f (x + 3h) − f (x + 2h)} − 2{f (x + 2h) − f (x + h)} + {f (x + h) − f (x)}

C
= {f (x + 3h) − 3f (x + 2h) + 3f (x + h) − f (x)} = 53 f (x + 3h)

aj
So after r − 1-th step we have ∆r f (x) = 5r f (x + rh) for r = 1, 2, ...

R
⇒ ∆r f (x0 ) = 5r f (x0 + rh) for r = 1, 2, ...

l
∆r y0 = 5r yr for r = 1, 2, ..., where yi = f (xi ) = f (x0 + ih) for r = 1, 2, ...

da
⇒

ha
is
ah
::M
ti
ai
M
ar
um
K
as
an
M
r.
D
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 17

2.5. Gregory-Newton’s forward interpolation formula. Let y = f (x) be a function

which takes values y0 , y1 , y2 , ..., yn at n+1 equally spaced (spacing h) interpolating points
x0 , x1 , x2 , ..., xn , of x. So xi = x0 + ih, i = 0, 1, 2, ...n. We like to derive a polynomial
φ(x) of degree n which agrees with f (x) at the interpolating points, i.e.,

φ(xi ) = f (xi ) = yi , i = 0, 1, 2, ..., n (1)

ge
Let us consider

le
φ(x) = a0 + a1 (x − x0 ) + a2 (x − x0 )(x − x1 ) + a3 (x − x0 )(x − x1 )(x − x2 )

ol
+... + an (x − x0 )(x − x1 )(x − x2 )...(x − xn−1 ) (2)

C
where ai , i = 0, 1, 2, ..., n are n + 1 unknown parameters to be determined

aj
R
Putting x = x0 in (2) we have

l
da
φ(x0 ) = a0 ⇒ a0 = y0 [using (1)]

ha
Again putting x = x1 in (2) we have

φ(x1 ) = a0 + a1 (x1 − x0 ) is
ah
⇒ y 1 = y 0 + a1 h [using (1)]
M

y1 − y0 ∆y0
⇒ a1 =
::

=
h h
ti
ai

Again putting x = x2 in (2) we have

M

φ(x2 ) = a0 + a1 (x2 − x0 ) + a2 (x2 − x0 )(x2 − x1 )

ar

∆y0
⇒ y2 = y0 + 2h + a2 .2h.h [using (1)]
um

h
y2 − y0 − 2∆y0 y2 − y0 − 2(y1 − y0 )
⇒ a2 =
K

2
=
2!h 2!h2
2
y2 − 2y1 + y0
as

∆ y0
= 2
=
2!h 2!h2
an
M

Again putting x = x3 in (2) we have

φ(x3 ) = a0 + a1 (x3 − x0 ) + a2 (x3 − x0 )(x3 − x1 ) + a3 (x3 − x0 )(x3 − x1 )(x3 − x2 )

r.
D

∆y0 ∆2 y0
⇒ y3 = y0 + 3h + 3h.2h + a3 .3h.2h.h [using (1)]
h 2!h2
y3 − y0 − 3∆y0 − 3∆2 y0
⇒ a3 =
3!h3
y3 − y0 − 3(y1 − y0 ) − 3(y2 − 2y1 + y0 )
=
3!h3
y3 − 3y1 + 3y2 − y0 ∆3 y0
= =
3!h3 3!h3
i
∆ y0
Similarly ai = for i = 1, 2, 3, ...n
i!hi
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 18

Putting these values of ai in (2) we get

∆y0 ∆2 y0 ∆3 y0
φ(x) = y0 + (x − x0 ) + (x − x0 )(x − x1 ) + (x − x 0 )(x − x 1 )(x − x 2 )
h 2!h2 3!h3
n
∆ y0
+... + (x − x0 )(x − x1 )(x − x2 )...(x − xn−1 )
n!hn
which is called Newton’s forward interpolation formula of f (x) for n + 1 interpolating

ge
points x0 , x1 , x2 , ..., xn and for x ∈ (x0 , xn ).

le
Let x − x0 = sh, then x − xi = x0 + sh − (x0 + ih) = h(s − i), for i = 1, 2, ..., n, then the

ol
above formula reduces to

C
∆2 y0 ∆3 y0
φ(x) = y0 + s∆y0 + s(s − 1) + s(s − 1)(s − 2) + ...
2! 3!

aj
∆n y0

R
+s(s − 1)(s − 2)...(s − n − 1)
n!

l
Note: Newton’s forward interpolation formula is suitable for the estimation of the value

da
of f (x) for x near x0 .

ha
Algorithm 4 Algorithm for Newton’s forward interpolation formula
is
ah
1: Input: Number of interpolating point n, data set {(x0 , y0 ), (x1 , y1 ), ..., (xn , yn )} and un-
known point xu .
M

2: Output: yu , the value of y at x = xu .

::

3: for i = 0 to n with increment 1 do . Calculation of difference table starts

ti

4: di0 = yi
ai

5: end for
M

6: for j = 1 to n with increment 1 do

7: for i = 0 to n − j with increment 1 do
ar

8: dij = di+1j−1 − yij−1

9: end for
um

10: end for . Calculation of difference table ends

11: h = x1 − x0 , Sum=y0 , P = 1, s = (xu − x0 )/h
K

12: for i = 1 to nh with increment

i 1 do h i
−xi−1
P = P × xu h×i or P = P × s−(i−1)
as

13: i . Calculation of successive terms

14: Sum=Sum+P × d0i
an

15: end for

M

16: Output: Sum

r.
D

2.6. Gregory-Newton’s backward interpolation formula. Let y = f (x) be a func-

tion which takes values y0 , y1 , y2 , ..., yn at n + 1 equally spaced (spacing h) interpolating
points x0 , x1 , x2 , ..., xn , of x. So xi = x0 + ih, i = 0, 1, 2, ...n. We like to derive a
polynomial φ(x) of degree n which agrees with f (x) at the interpolating points, i.e.,
φ(xi ) = f (xi ) = yi , i = 0, 1, 2, ..., n (1)
Let us consider
φ(x) = a0 + a1 (x − xn ) + a2 (x − xn )(x − xn−1 ) + a3 (x − xn )(x − xn−1 )(x − xn−2 )
+... + an (x − xn )(x − xn−1 )(x − xn−2 )...(x − x1 ) (2)
where ai , i = 0, 1, 2, ..., n are n + 1 unknown parameters to be determined
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 19

Putting x = xn in (2) we have

φ(xn ) = a0 ⇒ a0 = yn [using (1)]
Again putting x = xn−1 in (2) we have
φ(xn−1 ) = a0 + a1 (xn−1 − xn )
⇒ yn−1 = yn + a1 (−h) [using (1)]

ge
yn − yn−1 5yn
⇒ a1 = =
h h

le
ol
Again putting x = xn−2 in (2) we have

C
φ(xn−2 ) = a0 + a1 (xn−2 − xn ) + a2 (xn−2 − xn )(xn−2 − xn−1 )
5yn

aj
⇒ yn−2 = yn + (−2h) + a2 .(−2h).(−h) [using (1)]
h

R
yn−2 − yn + 2 5 yn yn−2 − yn + 2(yn − yn−1 )
⇒ a2 = =

l
da
2!h2 2!h2
2
yn−2 − 2yn−1 + yn 5 yn

ha
= 2
=
2!h 2!h2
Again putting x = xn−3 in (2) we have is
ah
φ(xn−3 ) = a0 + a1 (xn−3 − xn ) + a2 (xn−3 − xn )(xn−3 − xn−1 )
M

+a3 (xn−3 − xn )(xn−3 − xn−1 )(xn−3 − xn−2 )

::

5yn 52 y n
ti

⇒ yn−3 = yn + (−3h) + (−3h).(−2h) + a3 .(−3h).(−2h).(−h) [using (1)]

ai

h 2!h2
yn − yn−3 − 3 5 yn + 3 52 yn
M

⇒ a3 =
3!h3
ar

yn − yn−3 − 3(yn − yn−1 ) + 3(yn − 2yn−1 + yn−2 )

=
um

3!h3
yn − 3yn−1 + 3yn−2 − yn−3 53 yn
= =
K

3!h3 3!h3
5i yn
as

Similarly ai = for i = 1, 2, 3, ...n

i!hi
an

Putting these values of ai in (2) we get

M

5yn 52 y n 53 y n
φ(x) = yn + (x − xn ) + (x − xn )(x − xn−1 ) + (x − x n )(x − x n−1 )(x − x n−2 )
h 2!h2 3!h3
r.

n
5 yn
D

+... + (x − xn )(x − xn−1 )(x − xn−2 )...(x − x1 )

n!hn
which is called Newton’s backward interpolation formula of f (x) for n + 1 interpolating
points x0 , x1 , x2 , ..., xn and for x ∈ (x0 , xn ).
Let x − xn = sh, then x − xi = xn + sh − (x0 + ih) = nh + sh − ih = h(s + n − i), for
i = 1, 2, ..., n, then the above formula reduces to
52 y n 53 y n 5n y n
φ(x) = yn + s 5 yn + s(s + 1) + s(s + 1)(s + 2) + ... + s(s + 1)(s + 2)...(s + n − 1)
2! 3! n!
Note: This formula is suitable for the estimation of value of f (x) for x near xn .
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 20

Algorithm 5 Algorithm for Newton’s backword interpolation formula

1: Input: Number of interpolating point n, data set {(x0 , y0 ), (x1 , y1 ), ..., (xn , yn )} and un-
known point xu .
2: Output: yu , the value of y at x = xu .
3: for i = 0 to n with increment 1 do . Calculation of difference table starts
4: di0 = yi
5: end for

ge
6: for j = 1 to n with increment 1 do

le
7: for i = j to n with increment 1 do
8: dij = dij−1 − yi−1j−1

ol
9: end for

C
10: end for . Calculation of difference table ends

aj
11: h = x1 − x0 , Sum=yn , P = 1, s = (xu − xn )/h.
12: for i = 1 to nh with increment 1 do

R
xu −xn−(i−1)
i h i
s+(i−1)
13: P =P × h×i or P = P × i . Calculation of successive terms

l
Sum=Sum+P × dni

da
14:
15: end for

ha
16: Output: Sum

is
ah
Example: Using the following table of tan x estimate its value at 0.71.
M

x: 0.70 0.72 0.74 0.76

tan x : 0.84299 0.87707 0.91309 0.95045
::
ti

Ans: As the point x = 0.71 is near the initial Tabular values of x, and values of x are
ai

equidistant, here we have to use Newton’s forward interpolation formula. The forward
M

difference table for the given set of data is

∆2 y ∆3 y
ar

x y = tan x ∆y
0.70 0.84299
um

0.03478
0.72 0.87707 0.00124
K

0.03602 0.0001
as

0.74 0.91309 0.00134

an

0.03736
0.76 0.95045
M

Here h = 0.72 − 0.70 = 0.02, n = 4 and hence s = x−x h

0
= 0.71−0.70
0.02
= 0.5. So using
r.

Newton’s forward interpolation formula of degree 3 we have

D

y = tan x = y0 + s∆y0 + s(s − 1)∆2 y0 + s(s − 1)(s − 2)∆3 y0

0.00124 0.0001
tan 0.71 = 0.84299 + 0.5 × 0.03478 + 0.5(0.5 − 1) + 0.5(0.5 − 1)(0.5 − 2)
2 6
0.0001
= 0.84299 + 0.5 × 0.03478 + 0.5(−0.5) × 0.00062 + 0.5(−0.5)(−1.5)
6
= 0.859535
Example: Using the following table of tan x estimate its value at 0.75.
x: 0.70 0.72 0.74 0.76
tan x : 0.84299 0.87707 0.91309 0.95045
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 21

Ans: As the point x = 0.75 is near the end of the Tabular values of x, and values of
x are equidistant, here we have to use Newton’s backward interpolation formula. The
backward difference table for the given set of data is
x y = tan x 5y 52 y 53 y
0.70 0.84299
0.03478

ge
0.72 0.87707 0.00124
0.03602 0.0001

le
0.74 0.91309 0.00134

ol
0.03736

C
0.76 0.95045

aj
Here h = 0.72 − 0.70 = 0.02, n = 4 and hence s = x−x n
= 0.75−0.76 = −0.5. So using

R
h 0.02
Newton’s backward interpolation formula of degree 3 we have

l
y = yn + s 5 yn + s(s + 1) 52 yn + s(s + 1)(s + 2) 53 yn

da
0.00134

ha
tan 0.75 = 0.95045 + (−0.5) × 0.03736 + (−0.5)(−0.5 + 1) +
2
(−0.5)(−0.5 + 1)(−0.5 + 2)
0.0001
is
ah
6
0.0001
M

= 0.95045 + (−0.5) × 0.03736 + (−0.5)(0.5) × 0.00674 + (−0.5)(0.5)(1.5)

6
::

= 0.930079
ti
ai

Example: Using the following table of tan x estimate its value at 0.73.
M

x: 0.70 0.72 0.74 0.76

tan x : 0.84299 0.87707 0.91309 0.95045
ar

Ans: As the point x = 0.73 is near the middle of the Tabular values of x, though
um

the values of x are equidistant, here we have to use Lagrange’s interpolation formula.
Let x0 = 0.70, x1 = 0.72, x2 = 0.74, x3 = 0.76, and y0 = tan x0 = 0.84299, y1 =
K

tan x1 = 0.87707, y2 = tan x2 = 0.91309, y3 = tan x3 = 0.95045. Now using Lagrange’s

as

interpolation formula through these four points we have at x = 0.73

an

(0.73 − 0.72)(0.73 − 0.74)(0.73 − 0.76)

tan 0.73 = × 0.84299 +
(0.70 − 0.72)(0.70 − 0.74)(0.70 − 0.76)
M

(0.73 − 0.70)(0.73 − 0.74)(0.73 − 0.76)

× 0.87707 +
r.

(0.72 − 0.70)(0.72 − 0.74)(0.72 − 0.76)

D

(0.73 − 0.70)(0.73 − 0.72)(0.73 − 0.76)

× 0.91309 +
(0.74 − 0.70)(0.74 − 0.72)(0.74 − 0.76)
(0.73 − 0.70)(0.73 − 0.72)(0.73 − 0.74)
× 0.95045
(0.76 − 0.70)(0.76 − 0.72)(0.76 − 0.74)
(0.01)(−0.01)(−0.03) (0.03)(−0.01)(−0.03)
= × 0.84299 + × 0.87707 +
(−0.02)(−0.04)(−0.06) (0.02)(−0.02)(−0.04)
(0.03)(0.01)(−0.03) (0.03)(0.01)(−0.01)
+ × 0.91309 + × 0.95045
(0.04)(0.02)(−0.02) (0.06)(0.04)(0.02)
= −0.052686875 + 0.493351875 + 0.513613125 − 0.059403125 = 0.894875
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 22

Example: The area of a circle of diameter d id given for the following values:

d: 80 85 90 95
A : 5026 5674 6362 7088

Find approximate area of a circle having diameter 82.

Ans: As the point d = 82 is near the initial Tabular values of d, and the values of d are

ge
equidistant, here we have to use Newton’s forward interpolation formula. The forward

le
difference table for the given set of data is

ol
C
d A = f (d) ∆y ∆2 y ∆3 y
80 5026

aj
648

R
85 5674 40

l
688 −2

da
90 6362 38

ha
726
95 7088
is
ah
Here h = 85 − 80 = 5, n = 4 and hence s = d−d h
0
= 82−80
5
= 0.4. So using Newton’s
M

forward interpolation formula of degree 3 we have

::
ti

A = f (d) = y0 + s∆y0 + s(s − 1)∆2 y0 + s(s − 1)(s − 2)∆3 y0

ai

40 −2
f (82) = 5026 + 0.4 × 648 + 0.4(0.4 − 1) + 0.4(0.4 − 1)(0.4 − 2)
M

2 6
= 5026 + 0.4 × 648 + 0.4(−0.6) × 20 + 0.4(0.2)(−1.6)
ar

= 5280.272
um

Example: The area of a circle of diameter d id given for the following values:
K

d: 80 85 90 95
as

A : 5026 5674 6362 7088

an

Find approximate area of a circle having diameter 91.

M

Ans: As the point d = 91 is near the end of the Tabular values of d, and the values
r.

of d are equidistant, here we have to use Newton’s backward interpolation formula. The
D

backward difference table for the given set of data is

d A = f (d) 5y 52 y 53 y
80 5026
648
85 5674 40
688 −2
90 6362 38
726
95 7088
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 23

Here h = 85 − 80 = 5, n = 4 and hence s = d−d h

3
= 91−95
5
= −0.8. So using Newton’s
backward interpolation formula of degree 3 we have
A = f (d) = y3 + s 5 y3 + s(s + 1) 52 y3 + s(s + 1)(s + 2) 53 y3
38 −2
f (82) = 7088 − 0.8 × 726 − 0.8(−0.8 + 1) − 0.8(−0.8 + 1)(−0.8 + 2)
2 6
= 7088 − 0.8 × 726 − 0.8(0.2) × 19 + 0.8(0.2)(1.2)/3

ge
= 6504.224

le
Example: The area of a circle of diameter d id given for the following values:

ol
d: 80 85 90 95

C
A: 5026 5674 6362 7088
Find approximate area of a circle having diameter 86.

aj
R
Ans: As the point d = 86 is near the middle of the Tabular values of d, though values of
d are equidistant, here we have to use Lagrange’s interpolation formula. Using Lagrange

l
da
interpolation formula we have area for d = 86 is
(86 − 85)(86 − 90)(86 − 95) (86 − 80)(86 − 90)(86 − 95)

ha
A = × 5026 + × 5674
(80 − 85)(80 − 90)(80 − 95) (85 − 80)(85 − 90)(85 − 95)
(86 − 80)(86 − 85)(86 − 95) is
(86 − 80)(86 − 85)(86 − 90)
ah
+ × 6362 + × 7088
(90 − 80)(90 − 85)(90 − 95) (95 − 80)(95 − 85)(95 − 90)
M

(1)(−4)(−9) (6)(−4)(−9)
::

= × 5026 + × 5674
(−5)(−10)(−15) (5)(−5)(−10)
ti
ai

(6)(1)(−9) (6)(1)(−4)
+ × 6362 + × 7088
M

(10)(5)(−5) (15)(10)(5)
= −241.248 + 4902.336 + 1374.192 − 226.816 = 5808.464
ar
um
K
as
an
M
r.
D
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 24

3. Numerical Integration
Numerical integration is used when the function is unknown but a table of values taken by
the function is known. In that case we determine an appropriate interpolating polynomial
using the tabulated values and then integrating the interpolating polynomial over the
specified interval an approximate value of the required integration is obtained. Numerical
integration is also used when function is known but is not integrable analytically. In that

ge
case we divide the interval of integration by some equidistant points and then determine a
table of values of the function at those points. Then the function is replaced by a suitable

le
interpolating polynomial using the table of values and do the same process as discussed

ol
for the earlier case.

C
3.1. Newton-Cotes’ quadrature formula. Let us consider the definite integral

aj
R
Zb
I= f (x)dx (1)

l
da
a

ha
Divide the interval [a,b] into n equal subintervals by the equidistant points x0 , x1 , ..., xn ,
where xi = a + ih, h = (b − a)/n. Now Lagrangian interpolation polynomial φ(x) of f (x)
through n + 1 interpolating points x0 , x1 , ..., xn is given by is
ah
n
(x − x0 )(x − x1 )...(x − xi−1 )(x − xi+1 )...(x − xn )
M

X
φ(x) = f (xi ) (2)
(xi − x0 )(xi − x1 )...(xi − xi−1 )(xi − xi+1 )...(xi − xn )
::

i=0
ti

and error in approximating f (x) by φ(x) is

ai
M

(x − x0 )(x − x1 )..(x − xn ) n+1

E(x) = f (ξ) where x0 < ξ < xn (3)
(n + 1)!
ar

Let x − x0 = sh then x − xi = x0 + sh − xi = sh − ih = h(s − i) for i = 0, 1, 2, ..., n then

um

(2) reduces to
n
K

X hn s(s − 1)...(s − (i − 1))(s − (i + 1))...(s − n)

φ(x) = f (xi )
hn i.(i − 1)...1.(−1)(−2)...(−(n − i))
as

i=0
n
an

X s(s − 1)(s − 2)...(s − n)

= f (xi ) (4)
(−1)n−i i!(n − i)!(s − i)
M

i=0

Also (5) reduces to

r.
D

hn+1 s(s − 1)...(s − n) n+1

E(x) = f (ξ) where x0 < ξ < xn (5)
(n + 1)!
Now replacing f (x) by φ(x) in (1) we have
Zb Zn Zn
b−a
I ' φ(x)dx = φ(s)hds = φ(s)hds
n
a 0 0
Zn X
n
b−a s(s − 1)(s − 2)...(s − n)
= f (xi )ds
n i=0
(−1)n−i i!(n − i)!(s − i)
0
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 25
 n 
n
 s(s − 1)(s − 2)...(s − n) ds f (xi )
X Z
= (b − a)
i=0
(−1)n−i ni!(n − i)!(s − i)
0
n
X
= (b − a) Hi f (xi ) (6)
i=0
Rn

ge
1 (−1)n−i s(s−1)(s−2)...(s−n)
where Hi = n i!(n−i)! (s−i)
ds, i = 0, 1, 2, ..., n, are called Cotes’ coefficients
0

le
and formula (6) is called Newton-Cotes’ quadrature formula. Error in (6) is given by

ol
Zxn Zxn n+1

C
h s(s − 1)...(s − n) n+1
E(x)dx = f (ξ)dx
(n + 1)!

aj
x0 x0

R
Zn
hn+2 s(s − 1)...(s − n) n+1
= f (ξ)ds (6)

l
da
(n + 1)!
0

ha
3.2. Trapezoidal rule of integration. Let us consider the definite integral
Zb is
ah
I = f (x)dx (1)
M

a
::

Let the Lagrange’s interpolation polynomial of f (x) through the points (a, f (a)), (b, f (b))
ti

be φ(x). Then
ai

x−b x−a
M

φ(x) = f (a) + f (b) (2)

a−b b−a
ar

Now replacing f (x) by φ(x) in (1) we get

um

Zb Zb
I = f (x)dx = φ(x)dx
K

a a
as

Zb  
x−b x−a
an

= f (a) + f (b) dx
a−b b−a
M

a
b
(x − b)2 (x − a)2

= f (a) + f (b)
r.

2(a − b) 2(b − a) a
D

(a − b)2 (b − a)2
 
= − f (a) + f (b)
2(a − b) 2(b − a)
b−a
= [f (a) + f (b)] (-1)
2
This formula is called trapezoidal rule of integration of the integral (1).
The integral (1) physically represents the area bounded by the curve y = f (x) and the
lines x = a, x = b, y = 0. So above formula replaces the average function value of f (x)
in the interval [a, b] by average of f (a) and f (b). Hence it involves significant amount
of error. To reduce the error, the interval [a, b] is normally divided into n subintervals
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 26

by n + 1 equidistant points such that length of each subinterval h is sufficiently small.

Then applying Trapezoidal rule for integrating f (x) in each of the subintervals and then
summing all the integrals we determine the value of (1).
Let us divide [a, b] into n equal subintervals by the equidistant points x0 , x1 , ..., xn , where
xi = a + ih, h = (b − a)/n. Then we have
Zb Zxn

ge
I = f (x)dx = f (x)dx

le
a x0

ol
Zx1 Zx2 Zx3 Zxn

C
= f (x)dx + f (x)dx + f (x)dx + ... + f (x)dx

aj
x0 x1 x2 xn−1
x1 − x0 x2 − x1

R
= [f (x0 ) + f (x1 )] + [f (x1 ) + f (x2 )]
2 2

l
da
x3 − x 2 xn − xn−1
+ [f (x2 ) + f (x3 )] + ... + [f (xn−1 ) + f (xn )]
2 2

ha
h h h h
= [f (x0 ) + f (x1 )] + [f (x1 ) + f (x2 )] + [f (x2 ) + f (x3 )] + ... + [f (xn−1 ) + f (xn )]
2 2 2
is 2
ah
h
= [f (x0 ) + 2{f (x1 ) + f (x2 ) + f (x3 ) + ... + f (xn−1 )} + f (xn )] (-5)
M

2
This formula is called composite Trapezoidal rule of integration of the integral (1).
::
ti

3.3. Geometrical interpretation of Trapezoidal formula. The Trapezoidal formula

ai

Rb
M

of deriving f (x)dx is
a
ar

Zb
(b − a)
um

f (x)dx = [f (a) + f (b)]

2
a
K

(b−a)
The expression [f (a)
+ f (b)] represents the area of the trapezoid bounded by x-axis
as

2
(i.e., y = 0) the lines x = a, x = b, and the line joining the points (a, f (a)) and (b, f (b)).
an

Rb
On the other hand f (x)dx represents the area bounded by x-axis (i.e., y = 0) the lines
M

a
x = a, x = b, and the curve y = f (x). So the trapezoidal formula replaces the curve
r.

y = f (x) by the line joining the points (a, f (a)) and (b, f (b)).
D

Figure 2. Geometrical representation of Trapezoidal formula

 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 27

3.4. Error in Trapezoidal rule of integration. The Trapezoidal rule of deriving

Rb
f (x)dx is
a

Zb
(b − a) h
f (x)dx = [f (a) + f (b)] = [f (a) + f (a + h)], where h = b − a
2 2
a

ge
dF (x)
Let F (x) be the anti-derivative of f (x), i.e., dx
= f (x), then we have

le
Zb a+h

ol
Z
dF (x)
f (x)dx = dx = F (a + h) − F (a)

C
dx
a a

aj
h2 00 0 h3 000 h4 0000

R
= F (a) + hF (a) + F (a) + F (a) + F (a) + ..... − F (a)
2! 3! 4!

l
h2 0 h3 00 h4 000

da
= hf (a) + f (a) + f (a) + f (a) + .....
2! 3! 4!

ha
So error E is calculated as
Zb
h is
ah
E = f (x)dx − [f (a) + f (a + h)]
2
M

a
h2 0 h3 00 h4 000
::

= hf (a) + f (a) + f (a) + f (a) + ..... −

ti

2! 3! 4!
ai

h 0 h2 00 h3 000
[f (a) + {f (a) + hf (a) + f (a) + f (a) + .....}]
M

2 2! 3!
3
h 00
ar

= − f (a) + .....
12
um

h3 00
' − f (ξ), where a < ξ < b [Since h is small neglecting its higher order.]
12
K
as

Algorithm 6 Algorithm for numerical integration using Trapezoidal rule

an

1: Input: Interval of integration [a, b], and number of subintervals n.

M

2: Define: Integral function f (x)

3: h = (b − a)/n
r.

4: Sum=0
D

5: for i = 0 to n − 1 with increment 1 do

6: Sum=Sum+(h/2){f (a + ih) + f (a + (i + 1)h)}
7: end for
8: Output: Sum
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 28

3.5. Simpson’s 1/3-rd rule of integration. Let us consider the definite integral
Zb
I= f (x)dx (1)
a

Let us divide the interval of integration [a,b] into two subintervals by three points x0 , x1 ,
x2 , where x0 = a, x1 = a + h, x2 = a + 2h = b, h = (b − a)/2. The Newtons’s forward

ge
interpolation polynomial of f (x) through the points (x0 , f (x0 )), (x1 , f (x1 )), (x2 , f (x2 )) be

le
φ(x). Then

ol
∆2 f (x0 )

C
φ(x) = f (x0 ) + s∆f (x0 ) + s(s − 1) (2)
2

aj
where s = (x − x0 )/h

R
Now replacing f (x) by φ(x) in (1) we get

l
da
Zb Zb

ha
I = f (x)dx = φ(x)dx

is
a a
Z2 
ah
∆2 f (x0 )

= h f (x0 ) + s∆f (x0 ) + s(s − 1) ds
M

2
0
::

2
∆2 f (x0 ) s2 ∆2 f (x0 ) s3
  
ti

= h f (x0 )s + ∆f (x0 ) − +
ai

2 2 2 3
 0
M

2 2
  
∆ f (x0 ) ∆ f (x0 ) 8
= h 2f (x0 ) + 2 ∆f (x0 ) − +
2 2 3
ar

 
4
um

= h 2f (x0 ) + 2∆f (x0 ) − ∆2 f (x0 ) + ∆2 f (x0 )


3
K

 
1
= h 2f (x0 ) + 2∆f (x0 ) + ∆2 f (x0 )
3
as

 
1
an

= h 2f (x0 ) + 2{f (x1 ) − f (x0 )} + {f (x2 ) − 2f (x1 ) + f (x0 )}

3
M

h
= [f (x0 ) + 4f (x1 ) + f (x2 )]
3
r.
D

This formula is called Simpson’s 1/3-rd rule of integration of the integral (1).
Let us divide [a, b] into n equal subintervals by the equidistant points x0 , x1 , ..., xn , where
xi = a + ih, h = (b − a)/n, where n is even. Then we have
Zb Zxn
I= f (x)dx = f (x)dx
a x0
Zx2 Zx4 Zxn
= f (x)dx + f (x)dx + ... + f (x)dx
x0 x2 xn−2
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 29

h h
= [f (x0 ) + 4f (x1 ) + f (x2 )] + [f (x2 ) + 4f (x3 ) + f (x4 )] + ... +
3 3
h
[f (xn−2 ) + 4f (xn−1 ) + f (xn )]
3
h
= [f (x0 ) + 4{f (x1 ) + f (x3 ) + ... + f (xn−1 )} + 2{f (x2 ) + f (x4 ) + ... + f (xn−2 )}]
3
+f (xn )]

ge
This formula is called composite Simpson’s 1/3-rd rule of integration of the integral (1).

le
ol
3.6. Geometrical interpretation of Simpson’s 1/3-rd formula. The Simpson’s1/3-

C
Rb
rd formula of deriving f (x)dx is

aj
a

R
Zb
(b − a)
f (x)dx = [f (a) + f ((a + b)/2) + f (b)]

l
3

da
a

ha
(b−a)
The expression 3
[f (a) + f ((a + b)/2) + f (b)]
represents the area of the region bounded
by x-axis (i.e., y = 0) the lines x = a, x = b, and the parabola through the points
is Rb
ah
(a, f (a)), ((a + b)/2, f ((a + b)/2)), and (b, f (b)). On the other hand f (x)dx represents
M

a
the area bounded by x-axis (i.e., y = 0) the lines x = a, x = b, and the curve y = f (x).
::

So the Simpson’s1/3-rd formula replaces the curve y = f (x) by the parabola through the
ti

points (a, f (a)), ((a + b)/2, f ((a + b)/2)), and (b, f (b)).
ai
M
ar
um
K
as
an
M
r.

Figure 3. Geometrical representation of Simpson’s1/3-rd formula

D

3.7. Error in Simpson’s 1/3-rd rule of integration. The Simpson’s 1/3-rd rule of
Rb
deriving f (x)dx is
a

Zb
(b − a)
f (x)dx = [f (a) + 4f ((a + b)/2) + f (b)]
3
a
h
= [f (a) + 4f (a + h) + f (a + 2h)], where h = (b − a)/2
3
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 30

dF (x)
Let F (x) be the anti-derivative of f (x), i.e., dx
= f (x), then we have
Zb a+2h
Z
dF (x)
f (x)dx = dx = F (a + 2h) − F (a)
dx
a a
0 (2h)2 00 (2h)3 000 (2h)4 iv
= F (a) + 2hF (a) + F (a) + F (a) + F (a) +
2! 3! 4!

ge
(2h)5 v
F (a) + ..... − F (a)

le
5!

ol
0 4h3 00 2h4 000 4h5 iv
= 2hf (a) + 2h2 f (a) +

C
f (a) + f (a) + f (a) + .....
3 3 15

aj
So error E is calculated as

R
Zb
h
E = f (x)dx − [f (a) + 4f (a + h) + f (a + 2h)]

l
da
3
a

ha
2 0 4h3 00 2h4 000 4h5 iv
= 2hf (a) + 2h f (a) + f (a) + f (a) + f (a) + ..... −
3 3
is
15
ah
h 0 h2 00 h3 000 h4
[f (a) + 4{f (a) + hf (a) + f (a) + f (a) + f iv (a) + .....}
3 2! 3! 4!
M

2 3
0 (2h) 00 (2h) 000 (2h)4 iv
::

+{f (a) + 2hf (a) + f (a) + f (a) + f (a) + .....}]

2! 3! 4!
ti

h5
ai

= − f iv (a) + .....
M

90
h5
' − f iv (ξ), where a < ξ < b [Since h is small neglecting its higher order.]
ar

90
um

Algorithm 7 Algorithm for numerical integration using Simpson 1/3 rd rule

K

1: Input: Interval of integration [a, b], and number of subintervals n(multiple of two).
as

2: Define: Integral function f (x)

h = (b − a)/n
an

3:
4: Sum=0
M

5: for i = 0 to n with increment 2 do

6: Sum=Sum+(h/3){f (a + ih) + 4f (a + (i + 1)h) + f (a + (i + 2)h)}
r.

7: end for
D

8: Output: Sum

3.8. Simpson’s 3/8-th rule of integration. Let us consider the definite integral
Zb
I= f (x)dx (1)
a

Let us divide the interval of integration [a,b] into three subintervals by four points x0 , x1 ,
x2 , x3 , where x0 = a, x1 = a + h, x2 = a + 2h = b, x3 = a + 3h = b, h = (b − a)/3.
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 31

The Newtons’s forward interpolation polynomial of f (x) through the points (x0 , f (x0 )),
(x1 , f (x1 )), (x2 , f (x2 )), and (x3 , f (x3 )) be φ(x). Then
∆2 f (x0 ) ∆3 f (x0 )
φ(x) = f (x0 ) + s∆f (x0 ) + s(s − 1) + s(s − 1)(s − 2) (2)
2 6
where s = (x − x0 )/h
Now replacing f (x) by φ(x) in (1) we get

ge
Zb Zb

le
ol
I = f (x)dx = φ(x)dx

C
a a
Z3 

aj
∆2 f (x0 ) ∆3 f (x0 )

= h f (x0 ) + s∆f (x0 ) + s(s − 1) + s(s − 1)(s − 2) ds

R
2 6
0

l
da
3
∆2 f (x0 ) ∆3 f (x0 ) s2 ∆2 f (x0 ) − ∆3 f (x0 ) s3 ∆3 f (x0 ) s4
  
= h f (x0 )s + ∆f (x0 ) − + + +
2 3 2 2 3 6 4 0

ha
2 3 2 3 3
   
∆ f (x0 ) ∆ f (x0 ) 9 ∆ f (x0 ) − ∆ f (x0 ) 27 ∆ f (x0 ) 81
= h 3f (x0 ) + ∆f (x0 ) −
2
+
3 2
+
is 2 3
+
6 4
ah
   
3 9 2 27 3
M

2 3 3
= h 3f (x0 ) + 6∆f (x0 ) − 3∆ f (x0 ) + 2∆ f (x0 ) + {∆ f (x0 ) − ∆ f (x0 )} + ∆ f (x0 )
4 2 8
::

 
3h
ti

2 3 2 3 3
= 8f (x0 ) + 2{6∆f (x0 ) − 3∆ f (x0 ) + 2∆ f (x0 )} + 12{∆ f (x0 ) − ∆ f (x0 )} + 9∆ f (x0 )
ai

8
M

 
3h 2 3
= 8f (x0 ) + 12∆f (x0 ) + 6∆ f (x0 ) + ∆ f (x0 )
8
ar


3h
um

= 8f (x0 ) + 12{f (x1 ) − f (x0 )} + 6{f (x2 ) − 2f (x1 ) + f (x0 )} +

8

K

{f (x3 ) − 3f (x2 ) + 3f (x1 ) − f (x0 )}

as

 
3h
an

= f (x0 ) + 3f (x1 ) + 3f (x2 ) + f (x3 )

8
M

This formula is called Simpson’s 3/8-th rule of integration of the integral (1).
r.

Let us divide [a, b] into n equal subintervals by the equidistant points x0 , x1 , ..., xn , where
D

xi = a + ih, h = (b − a)/n, where n is multiple of 3. Then we have

Zb Zxn
I = f (x)dx = f (x)dx
a x0
Zx3 Zx6 Zxn
= f (x)dx + f (x)dx + ... + f (x)dx
x0 x3 xn−3
3h 3h
= [f (x0 ) + 3f (x1 ) + 3f (x2 ) + f (x3 )] + [f (x3 ) + 3f (x4 ) + 3f (x5 ) + f (x6 )]
8 8
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 32

3h
+... + [f (xn−3 ) + 3f (xn−2 ) + 3f (xn−1 ) + f (xn )]
8
3h
= [f (x0 ) + 3{f (x1 ) + f (x2 ) + f (x4 ) + f (x5 ) + ... + f (xn−2 ) + f (xn−1 )}
8
+2{f (x3 ) + f (x6 ) + ... + f (xn−3 ) + f (xn )} + f (xn )]
This formula is called composite Simpson’s 3/8-th rule of integration of the integral (1).

ge
le
3.9. Error in Simpson’s 3/8-th rule of integration. The Simpson’s 3/8-th rule of

ol
Rb
deriving f (x)dx is

C
a

aj
Zb
3h

R
f (x)dx = [f (a) + 3f (a + h) + 3f (a + 2h) + f (a + 3h)], where h = (b − a)/3
8

l
a

da
dF (x)
Let F (x) be the anti-derivative of f (x), i.e., = f (x), then we have

ha
dx

Zb
is
a+3h
Z
dF (x)
ah
f (x)dx = dx = F (a + 3h) − F (a)
dx
M

a a
(3h)2 00 (3h)3 000 (3h)4 iv
::

0
= F (a) + 3hF (a) + F (a) + F (a) + F (a) +
ti

2! 3! 4!
ai

(3h)5 v
F (a) + ..... − F (a)
M

5!
9h2 0 9h3 00 27h4 000 81h5 iv
ar

= 3hf (a) + f (a) + f (a) + f (a) + f (a) + .....

2 2 8 40
um

So error E is calculated as
K

Zb
3h
E = f (x)dx − [f (a) + 3f (a + h) + 3f (a + 2h) + f (a + 3h)]
as

8
an

a
9h2 0 9h3 00 27h4 000 81h5 iv
M

= 3hf (a) + f (a) + f (a) + f (a) + f (a) + ..... −

2 2 8 40
3h h2 00 h3 000 h4
r.

0
[f (a) + 3{f (a) + hf (a) + f (a) + f (a) + f iv (a) + .....}
8 2! 3! 4!
D

0 (2h)2 00 (2h)3 000 (2h)4 iv

+3{f (a) + 2hf (a) + f (a) + f (a) + f (a) + .....}
2! 3! 4!
0 (3h)2 00 (3h)3 000 (3h)4 iv
+{f (a) + 3hf (a) + f (a) + f (a) + f (a) + .....]
2! 3! 4!
3h5 iv
= − f (a) + .....
80
3h5 iv
' − f (ξ), where a < ξ < b [Since h is small neglecting its higher order.]
80
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 33

Algorithm 8 Algorithm for numerical integration using Simpson’s 3/8-th rule

1: Input: Interval of integration [a, b], and number of subintervals n(multiple of three).
2: Define: Integral function f (x)
3: h = (b − a)/n
4: Sum=0
5: for i = 0 to n with increment 3 do
Sum=Sum+(h/3){f (a + ih) + 3f (a + (i + 1)h) + 3f (a + (i + 2)h) + f (a + (i + 3)h)}

ge
6:
7: end for

le
8: Output: Sum

ol
C
3.10. Weddle’s rule of integration. Let us consider the definite integral

aj
Zb

R
I= f (x)dx (1)

l
da
a

ha
Let us divide the interval of integration [a,b] into six subintervals by four points x0 , x1 ,
x2 , x3 , x4 , x5 , x6 , where x0 = a, x1 = a + h, x2 = a + 2h = b, x3 = a + 3h, x4 = a + 4h,
is
x5 = a + 5h, x6 = a + 6h = b, h = (b − a)/6. The Newtons’s forward interpolation
ah
polynomial of f (x) through the points (x0 , f (x0 )), (x1 , f (x1 )), (x2 , f (x2 )), (x3 , f (x3 )),
M

(x4 , f (x4 )), (x5 , f (x5 )), and (x6 , f (x6 )) be φ(x). Then
::

∆2 f (x0 ) ∆3 f (x0 )
ti

φ(x) = f (x0 ) + s∆f (x0 ) + s(s − 1) + s(s − 1)(s − 2)

ai

2 6
4
∆ f (x0 ) ∆5 f (x0 )
M

+s(s − 1)(s − 2)(s − 3) + s(s − 1)(s − 2)(s − 3)(s − 4)

24 120
ar

∆6 f (x0 )
+s(s − 1)(s − 2)(s − 3)(s − 4)(s − 5) , where s = (x − x0 )/h
um

720
Now replacing f (x) by φ(x) in (1) we get
K

Zb Zb
as

I = f (x)dx = φ(x)dx
an

a a
M

Z6 
∆2 f (x0 ) ∆3 f (x0 )
= h f (x0 ) + s∆f (x0 ) + s(s − 1) + s(s − 1)(s − 2)
2 6
r.

0
D

∆4 f (x0 ) ∆5 f (x0 )
+s(s − 1)(s − 2)(s − 3) + s(s − 1)(s − 2)(s − 3)(s − 4)
24 120
∆6 f (x0 )

+s(s − 1)(s − 2)(s − 3)(s − 4)(s − 5) ds
720
3h
= [f (x0 ) + 5f (x1 ) + f (x2 ) + 6f (x3 ) + f (x4 ) + 5f (x5 ) + f (x6 )]
10
This formula is called Weddle’s rule of integration of the integral (1).
h 7
Error of this formula is − 140 f vi (ξ), where a < ξ < b.
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 34

Algorithm 9 Algorithm for numerical integration using Weddle’s rule

1: Input: Interval of integration [a, b], and number of subintervals n (multiple of six).
2: Define: Integral function f (x)
3: h = (b − a)/n
4: Sum=0
5: for i = 0 to n with increment 6 do
Sum=Sum+ 3h [f (a + ih) + 5f (a + (i + 1)h) + f (a + (i + 2)h) + 6f (a + (i + 3)h) +

ge
6: 10
f (a + (i + 4)h) + 5f (a + (i + 5)h) + f (a + (i + 6)h)]

le
7: end for

ol
8: Output: Sum

C
aj
3.11. Gauss quadrature formula. Gauss quadrature formula was introduced by Friedrich

R
Gauss in 1815. It is a quadrature rule constructed to provide an exact result for poly-
nomials of degree 2n − 1 or less by a suitable choice of the nodes xi and weights wi for

l
da
i = 1, ..., n. The general form of this formula is

ha
Z1 n
X
I= F (u)du = wi F (ui ) (1)
is
ah
−1 i=1

where the weights w1 , w2 , ..., wn and interpolating points u1 , u2 , ..., un are selected in such
M

a manner that (1) provides an exact result for a polynomial of degree ≤ 2n − 1. Let F (u)
::

be polynomial of degree 2n − 1 as given in (2).

ti
ai

2n−1
X
ai ui
M

F (u) = (2)
i=0
ar

Then from (1) we get

um

Z1 "2n−1
X
# n
X 2n−1
X
ai ui du = wi aj uji
K

−1 i=0 i=1 j=0

as

n 2n−1
2 2 2a2n−2 X X
an

⇒ 2a0 + a2 + a4 + ... + = wi aj uji (2)

3 5 2n1
M

i=1 j=0

Equating coefficients of ai from both sides of (2) we have the following system of equations
r.


w1 + w2 + ... + wn = 2 
D


u1 w1 + u2 w2 + ... + un wn = 0 


u21 w1 + u22 w2 + ... + u2n wn = 32 (3)
... ... ... ... ... ... 


u2n−1 2n−1
+ un2n−1 wn = 0 

1 w1 + u 2 w 2 + ...
For n = 2 we have

w1 + w2 = 2 

u1 w 1 + u2 w 2 = 0

(4)
u21 w1 + u22 w2 = 32 

u31 w1 + u32 w2 = 0 
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 35

Now solving (4) we have

−1 1
w1 = 1, w2 = 1, x1 = √ , x2 = √
3 3
So we have two point Gaussian interpolation formula as
Z1    
−1 1

ge
I= F (u)du = F √ +F √ (5)
3 3

le
−1

ol
For n = 3 from (3) we have

C
w1 + w2 + w3 = 2


aj

u1 w1 + u2 w2 + u3 w3 = 0



R


u21 w1 + u22 w2 + u23 w3 = 32

(6)
u31 w1 + u32 w2 + u33 w3 = 0

l
da

u41 w1 + u42 w2 + u43 w3 = 52




u51 w1 + u52 w2 + u53 w3 = 0


ha
Now solving (6) we have
is
ah
r r
8 5 3 3
M

w1 = , w2 = w3 = , x1 = 0, x2 = − , x3 =
9 9 5 5
::

So we have three point Gaussian interpolation formula as

ti
ai

Z1 r ! r !
8 5 3 5 3
M

I = F (u)du = F (0) + F − + F (7)

9 9 5 9 5
ar

−1
um

Let us consider the definite integral

Zb
K

I= f (x)dx (8)
as

a
an

To solve it using Gauss quadrature transform the interval of integration [a,b] to [-1,1] by
M

the transformation
b−a a+b
x= u+
r.

2 2
D

Which transform (8) to (9) as presented below.

Z1  
b−a b−a a+b
I= f u+ du (9)
2 2 2
−1

Now using n point Gaussian quadrature formula we have

Z1   n  
b−a b−a a+b b−aX b−a a+b
I= f u+ du = wi f ui + (10)
2 2 2 2 i=1 2 2
−1
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 36

3.11.1. Composite Gauss quadrature formula. Let us consider the definite integral
Zb
I= f (x)dx (1)
a

To derive composite two-point Gauss quadrature formula for (1) first divide the interval
[a,b] into n intervals by taking n + 1 equidistant points x0 , x1 ,...,xn , where xi = a + ih,

ge
h = (b − a)/n. Then we have

le
ol
Zb Zx1 Zx2 Zxn

C
I = f (x)dx = f (x)dx + f (x)dx + ... + f (x)dx

aj
a x0 x1 xn−1
xi+1

R
n−1 n−1 Z1  
xi+1 − xi xi+1 − xi
Z
X X xi + xi+1
= f (x)dx = f u+ dx

l
2 2 2

da
i=0 x i=0 −1
i

ha
n−1 Z1   
X h hu a + ih + a + (i + 1)h
= f + du
2 2 2
is
ah
i=0 −1
Z1  
M

n−1 
X h hu h
= f + a + ih + du
::

2 2 2
ti

i=0 −1
ai

n−1     
X h hu1 h hu2 h
M

= f + a + ih + +f + a + ih +
i=0
2 2 2 2 2
ar

−1 √1 .
where u1 = √ and u2 =
um

3 3

Similarly, to derive composite three-point Gauss quadrature formula for (1) first divide
K

the interval [a,b] into even number of intervals, say n, by taking n + 1 equidistant points
x0 , x1 ,...,xn , where xi = a + ih, h = (b − a)/n. Then we have
as
an

Zb Zx2 Zx4 Zxn

M

I = f (x)dx = f (x)dx + f (x)dx + ... + f (x)dx

a x0 x2 xn−2
r.

n
−1 x2(i+1) n
−1 1
D

2 2
X x2(i+1) − x2i Z x2(i+1) − x2i
 
X Z x2i + x2(i+1)
= f (x)dx = f u+ dx
i=0 i=0
2 2 2
x2i −1
n
2
−1
X  
x2i + x2(i+1)
 
x2i + x2(i+1)

= h w1 f hu1 + + w2 f hu2 +
i=0
2 2
 
x2i + x2(i+1)
+w3 f hu3 +
2
q q
where w1 = 9 , w2 = w3 = 9 , u1 = 0, u2 = − 5 , u3 = 35
8 5 3
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 37

Algorithm 10 Numerical integration using two-point composite Gauss quadrature

1: Input: Interval of integration [a, b], and number of subintervals n.
2: Define: Integral function f (x)
3: h = (b − a)/n
4: Sum=0
−1
5: u1 = √ 3
, u2 = √13
for i = 0 to n − 1with increment 1 do

ge
6:
Sum=Sum+ h2 f hu2 1 + a + ih + h2 + f hu2 2 + a + ih + h2



7:

le
8: end for

ol
9: Output: Sum

C
aj
Algorithm 11 Numerical integration using three-point composite Gauss quadrature

R
1: Input: Interval of integration [a, b], and number of subintervals n (even number).

l
2: Define: Integral function f (x)

da
3: h = (b − a)/n

ha
4: Sum=0 q q
5: 8 5
w1 = 9 , w2 = w3 = 9 , u1 = 0, u2 = − 5 , u3 = 35
3

for i = 0 to n2 − 1 with increment 1 do is

ah
6:
7: x2i = a + 2ih,h x2(i+1)
 = a + 2(i + 1)h
M

    i
x +x x +x x +x
8: Sum=Sum+ w1 f hu1 + 2i 22(i+1) + w2 f hu2 + 2i 22(i+1) + w3 f hu3 + 2i 22(i+1)
::

end for
ti

9:
Output: Sum
ai

10:
M

3.12. Midpoint rule. Let us consider the definite integral

ar
um

Zb
I= f (x)dx (1)
K

a
as

Let us divide the interval of integration [a, b] into n subintervals by n + 1 equidistant

an

points x0 , x1 , x2 , ..., xn , where x0 = a, xi = a + ih, i = 0, 1, 2, ..., n. The midpoint rule of

integration is given by
M

      
x0 + x1 x1 + x 2 xn−1 + xn
I = h f +f + ... + f (2)
r.

2 2 2
D

3.13. Boole’s rule. Let us consider the definite integral

Zb
I= f (x)dx (1)
a

Let us divide the interval of integration [a,b] into four subintervals by five equidistant
points x0 , x1 , x2 , x3 , x4 , where x0 = a, x1 = a + h, x2 = a + 2h = b, x3 = a + 3h = b,
x4 = a + 4h = b, h = (b − a)/4. The Newtons’s forward interpolation polynomial of
f (x) through the points (x0 , f (x0 )), (x1 , f (x1 )), (x2 , f (x2 )), (x3 , f (x3 )), and (x3 , f (x3 ))
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 38

be φ(x). Then
∆2 f (x0 ) ∆3 f (x0 )
φ(x) = f (x0 ) + s∆f (x0 ) + s(s − 1) + s(s − 1)(s − 2)
2 6
∆4 f (x0 )
+s(s − 1)(s − 2)(s − 3) , where s = (x − x0 )/h
24
Now replacing f (x) by φ(x) in (1) we get

ge
Zb Zb

le
I = f (x)dx = φ(x)dx

ol
C
a a
Z4

aj
∆2 f (x0 ) ∆3 f (x0 )

= h f (x0 ) + s∆f (x0 ) + s(s − 1) + s(s − 1)(s − 2)

R
2 6
0

l
∆4 f (x0 )

da

+s(s − 1)(s − 2)(s − 3) ds
24

ha
2h
=
45 is
[7f (x0 ) + 32f (x1 ) + 12f (x2 ) + 32f (x3 ) + 7f (x4 )]
ah
R2
M

Example: Evaluate e−x/2 dx using four intervals and compare the result with actual
::

1
ti

value(or determine the error).

ai

Ans: Here the interval of integration is [1, 2] and is to be divided by four subintervals.
M

So length of each subinterval is h = (2 − 1)/4 = 0.25 and hence points of division are
x0 = 1, x1 = 1.25, x2 = 1.5, x3 = 1.75 and x4 = 2. Then values of the integrating function
ar

f (x) = e−x/2 at the points of subdivisions are tabulated below:

um

x: 1 1.25 1.50 1.75 2

−x/2
f (x) = e : 0.60653 0.53526 0.47237 0.41686 0.36788
K

Then using Trapezoidal rule value of the integral is

as
an

Z2
h
e−x/2 dx = [f (x0 ) + 2{f (x1 ) + f (x2 ) + f (x3 )} + f (x4 )]
M

2
1
r.

0.25
= [0.60653 + 2{0.53526 + 0.47237 + 0.41686} + 0.36788]
D

2
= 0.47792375
Then using Simpson’s 1/3-rd rule value of the integral is
Z2
h
e−x/2 dx = [f (x0 ) + 4{f (x1 ) + f (x3 )} + 2f (x2 ) + f (x4 )]
3
1
0.25
= [0.60653 + 4{0.53526 + 0.41686} + 2 × 0.47237 + 0.36788]
3
= 0.4773025
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 39

Actual value of integral is

Z2
2
e−x/2 dx = −2 e−x/2 1 = −2 e−1 − e−0.5 = 0.4773024 [correct upto 7 significant figures]
  

Then error in the evaluated value using Trapezoidal rule is |0.47792375 − 0.4773024| =
0.00062135 and error in the evaluated value using Simpson’s 1/3-rd rule is |0.4773025 −

ge
0.4773024| = 0.0000001.

le
0.7
R
Example: Evaluate (ex + 2x)dx using h=0.1

ol
0.1

C
Ans: Here the interval of integration is [0.1, 0.7] and length of each subinterval is 0.1

aj
and hence number of subintervals is n = (0.7 − 0.1)/0.1 = 6. So the points of subdivision

R
are x0 = 0.1, x1 = 0.2, x2 = 0.3, x3 = 0.4, x4 = 0.5, x5 = 0.6, x6 = 0.7. Then values of
the integrating function f (x) = ex + 2x at the points of subdivisions are tabulated below:

l
da
x: 0.1 0.2 0.3 0.4 0.5 0.6 0.7
f (x) = ex + 2x : 1.30517 1.62140 1.94986 2.29182 2.64872 3.02212 3.41375

ha
Then using Trapezoidal rule value of the integral is
is
ah
Z0.7
h
(ex + 2x)dx = [f (x0 ) + 2{f (x1 ) + f (x2 ) + f (x3 ) + f (x4 ) + f (x5 )} + f (x6 )]
M

2
0.1
::

0.1
ti

= [1.30517 + 2{1.62140 + 1.94986 + 2.29182 + 2.64872 + 3.02212} + 3.41375]

ai

2
= 1.389338
M

Then using Simpson’s 1/3-rd rule value of the integral is

ar

Z2
um

h
e−x/2 dx = [f (x0 ) + 4{f (x1 ) + f (x3 ) + f (x5 )} + 2{f (x2 ) + f (x4 )} + +f (x6 )]
3
K

1
0.1
= [1.30517 + 4{1.62140 + 2.29182 + 3.02212} + 2{1.94986 + 2.64872} + 3.41375]
as

3
an

= 1.38858
M
r.
D
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 40

4. Solution of algebraic and transcendental equations

4.1. Bisection method. Let us consider the following algebraic or transcendental equa-
tion
f (x) = 0 (1)
To solve (1) using bisection method we first search an interval [a, b] such that f (a)f (b) < 0,
which implies, there is a real root of (1) in [a, b]. We rename the interval [a, b] as [a1 , b1 ]

ge
and divide it into two equal subintervals [a1 , c1 ] and [c1 , b1 ] by the point c1 = (a1 + b1 )/2.

le
If f (c1 ) = 0 or |f (c1 )| is negligible then we consider c1 as the root of (1). Otherwise two

ol
cases will arise- either f (a1 )f (c1 ) < 0 or f (c1 )f (b1 ) < 0. In the first case root lies inthe

C
interval [a1 , c1 ] and we take it as search interval of next iteration and rename it as [a2 , b2 ].
In the second case root lies in the interval [c1 , b1 ] and we take it as search interval of next

aj
iteration and rename it as [a2 , b2 ]. Again we divide the interval [a2 , b2 ] by two equal subin-

R
tervals by the point c2 = (a2 + b2 )/2. If f (c2 ) = 0 or |f (c2 )| is negligible then we consider

l
c2 as the root of (1). Otherwise proceeding same as the previous iteration get the next

da
search interval [a3 , b3 ]. We continue this process and get a sequence of points c1 , c2 , c3 , ...

ha
and we continue the process until we get a point cm such that |f (cm )| is negligible. Then
cm is considered as the real root of (1). Bisection method always converges.
is
ah
To get a result upto k significant figures the iteration stops when k significant figures
M

of two successive approximations becomes equal and final approximated value is taken as
::

the root of the equation with desired accuracy.

ti

4.1.1. Order of convergence of an iterative method. Let xn and xn+1 be n-th and n + 1-th
ai

approximations of a real root ξ of the equation f (x) = 0. Let error in xn and xn+1 are n
M

and n+1 respectively, i.e., n = ξ − xn and n+1 = ξ − xn+1 . Then the iterative process is
said to have a order of convergence k if ∃ a real number M > 0 such that the following
ar

relation holds:
um

n+1
< M, ∀ n ∈ N
kn
K

4.2. Fixed-point iteration method. Let us consider the following algebraic or tran-
as

scendental equation
an

f (x) = 0 (1)
M

To solve (1) using fixed-point iteration method we first search a small interval [a, b] such
that f (a)f (b) < 0, which implies, there is a real root of (1) in [a, b]. Now, x0 = (a + b)/2
r.

can be considered as an initial guess of the real root of (1) in the said interval. Then (1)
D

is rewritten in the form

x = φ(x) (2)
Then successive approximations x1 , x2 , ..., xn of the real root near x0 are obtained using
the formula (3).
xn = φ(xn−1 ), i = 1, 2, ..., n (3)
Now, if the sequence of successive approximations {xn } converges, then after some itera-
tion k, we should have
xk ' xk+1 ⇒ xk ' φ(xk ) ⇒ f (xk ) ' 0 ⇒ xk is an approximate real root of (1)
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 41

In reality we continue to evaluate successive approximations of the roots of (1) using (3)
and when difference between two successive approximations become negligible, then final
approximate value is considered as the root of (1).

To get a result upto k significant figures the iteration stops when k significant figures
of two successive approximations becomes equal and final approximated value is taken as
the root of the equation with desired accuracy.

ge
le
4.2.1. Convergence criterion of fixed-point iteration method. Let us consider the following

ol
algebraic or transcendental equation

C
f (x) = 0 (1)

aj
The fixed-point iterative formula for finding a real root ξ of (1) in an interval [a, b] for an

R
initial guess x0 ∈ [a, b] is

l
da
xn = φ(xn−1 ), n = 1, 2, ... (2)

ha
where x = φ(x) is an alternative representation of the equation (1). Since ξ is a real root
of (1), i.e., of x = φ(x), we have
is
ah
ξ = φ(ξ) (3)
M

Now using (2) and (3) we have

::
ti

ξ − xn+1 = φ(ξ) − φ(xn ), n = 1, 2, ... (4)

ai

Let φ(x) is continuous and derivable in [a, b], then it is continuous and derivable in [xn , ξ].
M

So using mean value theorem ∃ ξn ∈ (xn , ξ) such that the relation (5) holds.
ar

0
ξ − xn+1 = (ξ − xn )φ (ξn ), n = 1, 2, ... (5)
um

Using successive relations of (5) we have

K

0
ξ − xn+1 = (ξ − xn )φ (ξn )
as

0 0 0
= (φ(ξ) − φ(xn−1 ))φ (ξn ) = (ξ − xn−1 )φ (ξn−1 )φ (ξn )
an

0 0 0 0 0
= (φ(ξ) − φ(xn−2 ))φ (ξn−1 )φ (ξn ) = (ξ − xn−2 )φ (ξn−2 )φ (ξn−1 )φ (ξn )
M

0 0 0 0 0
= ........ = (ξ − x0 )φ (ξ0 )φ (ξ1 )...φ (ξn−2 )φ (ξn−1 )φ (ξn )
r.

0
Now, if |φ (x)| < l < 1, ∀ x ∈ [a, b] then it follows from the above equation that
D

|ξ − xn+1 | < |ξ − x0 |ln+1 → 0 as n → ∞

⇒ xn+1 → ξ as n → ∞, i.e., {xn } → ξ as n → ∞
Theorem. If φ(x) is a continuous function in the interval [a,b] that contains a real root
0
of f (x) = 0, which has been written as x = φ(x), and if |φ (x)| < l < 1, ∀ x ∈ [a, b] then
for any choice of x0 ∈ [a, b], the sequence {xn } determined from
xn+1 = φ(xn ), n = 0, 1, 2, ...
converges to the root of x = φ(x).
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 42

4.2.2. Order of convergence of fixed-point iteration method. Let us consider a fixed-point

iterative formula xn = φ(xn−1 ) in finding successive approximations a real root ξ of an
equation f (x) = 0 in an interval [a, b]. Then x = φ(x) is an alternative representation of
0
the equation f (x) = 0. A sufficient condition for convergence of the process is φ (x)| <
1, ∀ x ∈ [a, b]. Assuming φ(x) is continuous and derivable on [a, b], using mean value
theorem we have
0 en+1 0

ge
|ξ − xn+1 | = |φ(ξ) − φ(xn )| = |ξ − xn ||φ (ξn )| ⇒ = |φ (ξn )| < 1
en

le
where en = |ξ − xn | is the error if the exact root ξ is approximated by xn , i.e., error in

ol
any step is proportional to the error in the previous step. So order of convergence of the

C
iterative process is linear.

aj
Algorithm 12 Solution of a transcendental equation by fixed-point iteration method

R
1: Input: Search interval [a, b] and accuracy of significant figures k.

l
da
2: Define: Function f (x) and φ(x)
3: if f (a)f (b) > 0 then

ha
4: Output: A real root may not exists in [a, b]
5: Stop execution.
is
ah
6: end if
7:  ← 10−k
M

8: x0 ← (a + b)/2
::

9: x1 ← φ(x0 )
ti

10: while |x1 − x0 | ≥  do

ai

11: x0 ← x1
M

12: x1 ← φ(x0 )
13: end while
ar

14: Output: x1 , f (x1 )

um
K

4.3. Newton-Raphson method. Let us consider the following algebraic or transcen-

dental equation
as

f (x) = 0 (1)
an

To solve (1) using Newton-Rapshon method we first search a small interval [a, b] such
M

that f (a)f (b) < 0, which implies, there is a real root of (1) in [a, b]. Now we consider
x0 = (a + b)/2 as initial guess of the real root of (1) in the said interval. Let h be the
r.

error in x0 . Then x1 = x0 + h is actual root of (1). Hence we have

D

f (x0 + h) = 0 and so using Taylor series expansion we have

0 h2 00
f (x0 ) + hf (x0 ) + f (x0 ) + ... = 0. Now assuming h is small and neglecting its 2nd and
2!
0 −f (x0 ) f (x0 )
higher order we have f (x0 ) + hf (x0 ) = 0 ⇒ h = 0 ⇒ x1 = x0 − 0 .
f (x0 ) f (x0 )
x1 is the first approximation of the real root of (1). Continuing this process we shall get
the successive approximations using the following formula
f (xn−1 )
xn = xn−1 − 0 for n = 1, 2, 3, ... (-1)
f (xn−1 )
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 43

Now, if the sequence of successive approximations {xn } converges, then after some itera-
tion k, we should have
f (xk )
xk = xk+1 ⇒ xk = xk − 0 ⇒ f (xk ) = 0 ⇒ xk is a real root of (1)
f (xk )
In reality we continue to evaluate successive approximations of the roots of (1) using (2)
and when difference between two successive approximations become negligible, then final

ge
approximate value is considered as the root of (1). To get a result upto k significant
figures the iteration stops when k significant figures of two successive approximations

le
becomes equal and final approximated value is taken as the root of the equation with

ol
desired accuracy.

C
4.3.1. Convergence criterion of Newton-Raphson method. Let us consider the following

aj
algebraic or transcendental equation

R
f (x) = 0 (1)

l
da
The Newton-Raphson iterative formula for finding a real root ξ of (1) in an interval [a, b]
for an initial guess x0 ∈ [a, b] is

ha
f (xn )
xn+1 = xn −
f 0 (xn )
for n = 1, 2, ...
is (2)
ah
Comparing this equation with the fixed-point iterative formula
M

xn = φ(xn−1 ) for n = 1, 2, 3, ...

::
ti

f (x)
we have φ(x) = x − f 0 (x)
. Since the condition of convergence of a fixed-point iterative
ai

0
method is |φ (x)| < 1,i∀ x ∈ [a, b], the condition of convergence of Newton-Raphson
M

h 00
00 0
method is dx x − ff0(x)
d
(x)
< 1, i.e., f (x)f (x)
0
[f (x)]2
< 1, i.e., |f (x)f (x)| < [f (x)]2 , ∀ x ∈ [a, b].
ar

4.3.2. Order of convergence of Newton-Raphson method. The iterative formula for finding
um

a real root ξ of f (x) = 0 in an interval [a, b] in Newton-Raphson method is xn+1 =

xn − ff0(x n)
for an initial guess x0 ∈ [a, b]. Since ξ is a real root of f (x) = 0 and xn is n-th
K

(xn )
order approximation for initial guess x0 ∈ [a, b], so we have
as

f (ξ) = 0 (2)
an

⇒ f (ξ) = f (ξ − xn + xn ) = f (en + xn ) = 0
M

0 e2 00
⇒ f (xn ) + en f (xn ) + n f (xn ) + ... = 0 by Taylor series expansion of f about xn
2
r.

2 00
f (xn ) e f (xn )
D

⇒ 0 + en + n 0 + ... = 0
f (xn ) 2 f (xn )
00
f (xn ) e2n f (xn )
⇒ + (ξ − x n ) + '0
f 0 (xn ) 2 f 0 (xn )
00
f (xn ) e2 f (xn )
⇒ xn − 0 −ξ ' n 0
f (xn ) 2 f (xn )
00 00
e2 f (xn ) en+1 1 f (xn )
⇒ xn+1 − ξ ' n 0 ⇒ '
2 f (xn ) 2
en 2 f 0 (xn )
So, the order of convergence of the Newton-Raphson method is two.
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 44

4.3.3. Geometrical representation of Newton-Raphson method. Let us consider the follow-

ing algebraic or transcendental equation
f (x) = 0 (1)
The iterative formula for finding a real root ξ of (1) in an interval [a, b] in Newton-Raphson
method is xn+1 = xn − ff0(x n)
(xn )
for an initial guess x0 ∈ [a, b]. Now equation of the tangent
to (1) at (x0 , f (x0 )) is

ge
0
y − f (x0 ) = f (x0 )(x − x0 ) (2)

le
If this tangent cuts the x-axis at x = x1 then we have

ol
C
0 f (x0 )
0 − f (x0 ) = f (x0 )(x1 − x0 ) ⇒ x1 = x0 −
f 0 (x0 )

aj
Here x1 represents the first approximation of the real root of (1) for the initial guess

R
x0 ∈ [a, b]. So an approximation of the root is obtained from the point of intersection of

l
the tangent to (1) at the previous approximation and the x-axis, i.e., x2 is the point of

da
intersection of the tangent to (1) at x = x1 and x-axis, x3 is the point of intersection of

ha
the tangent to (1) at x = x2 and x-axis, and so on.

is
ah
:: M
ti
ai
M
ar
um

Figure 4. Successive approximations of a real root of f (x) = 0 in Newton-

Raphson method
K

4.4. Regula-Falsi method. Let us consider the following algebraic or transcendental

as

equation
an

f (x) = 0 (1)
M

To solve (1) using Regula-Falsi method we first search a small interval [x0 , x1 ] such that
r.

f (x0 )f (x1 ) < 0, which implies, there is a real root of (1) in [x0 , x1 ]. Now draw a straight
D

line (2) through the points (x0 , f (x0 )) and (x1 , f (x1 )).
y − f (x0 ) x − x0
= (2)
f (x1 ) − f (x0 ) x1 − x0
Now (2) cuts the x-axis at x = x2 , whose expression is presented in (3).
f (x0 )(x1 − x0 )
x2 = x0 − (3)
f (x1 ) − f (x0 )
x2 is the first approximation of the real root of (1) in [x0 , x1 ]. Now, three cases may arises:
1: The value of |f (x2 )| is near to zero. In this case the process stops and x2 is taken
as the approximate real root of (1).
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 45

Algorithm 13 Solution of a transcendental equation by Newton Raphson method

1: Input: Search interval [a, b] and accuracy of significant figures k.
0
2: Define: Function f (x) and f (x)
3:  ← 110−k
4: if f (a)f (b) > 0 then
5: Output: A real root may not exists in [a, b]
Stop execution.

ge
6:
7: end if

le
8: x0 ← (a + b)/2

ol
9: x1 ← x0 − ff0(x 0)

C
(x0 )
10: while |x1 − x0 | ≥  do

aj
11: x0 ← x1
x1 ← x0 − ff0(x 0)

R
12:
(x0 )
13: end while

l
da
14: Output: x1 , f (x1 )

ha
is
2: The value of f (x0 )f (x2 ) < 0. In this case [x0 , x2 ] is taken as new search interval
ah
and is renamed as [x0 , x1 ], i.e., the value of x1 is replaced by the value of x2 . The
search process continues with the new search interval by finding new value of x2
M

using (3).
::

3: The value of f (x1 )f (x2 ) < 0. In this case [x2 , x1 ] is taken as new search interval
ti

and is renamed as [x0 , x1 ], i.e., the value of x0 is replaced by the value of x2 . The
ai

search process continues with the new search interval by finding new value of x2
M

using (3).
This iterative process always converges to the desired root if initial search interval is
ar

sufficiently small so that supremum is less than twice of infimum in that interval. A
um

pictorial representation of finding successive approximations of a real root of f (x) = 0 is

presented in the following figure. To get a result upto k significant figures the iteration
K
as
an
M
r.
D

Figure 5. Successive approximations of a real root of f (x) = 0 in Regula Falsi method

stops when k significant figures of two successive approximations becomes equal and final
approximated value is taken as the root of the equation with desired accuracy.

4.4.1. Convergence analysis of Regula-Falsi method. Let us consider the following alge-
braic or transcendental equation
f (x) = 0 (1)
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 46

The Regula-Falsi iterative formula for finding a real root ξ of (1) in a small interval [x0 , x1 ]
such that f (x0 )f (x1 ) < 0, is
f (xn )(xn − xn−1 )
xn+1 = xn − for n = 1, 2, ...
f (xn ) − f (xn−1 )
(f (xn ) − f (xn−1 ))
⇒ (xn+1 − xn ) = −f (xn−1 )
(xn − xn−1 )

ge
0
⇒ (xn+1 − xn )f (ξn ) = −f (xn ), ξn ∈ [xn−1 , xn ] [Applying mean value theorem]

le
0
⇒ (xn − xn+1 )f (ξn ) = f (xn )

ol
0
⇒ (xn − ξ + ξ − xn+1 )f (ξn ) = f (xn )

C
0 0 0 0
⇒ (−n + n+1 )f (ξn ) = f (xn ) − f (ξ) = (xn − ξ)f (ξn ), ξn ∈ [xn−1 , xn ]

aj
0 0 0 0 0
f (ξ ) f (ξn ) − f (ξn ) M −m

R
⇒ n+1 = n − n 0 n = n 0 ≤ n
f (ξn ) f (ξn ) m

l
da
where M and m are the supremum and infimum of f (x) in [x0 , x1 ]
2 n

ha
 
M −m M −m
⇒ n+1 ≤ n−1 ≤ ... ≤ 0
m m
is
ah
 −m 
If the initial interval [x0 , x1 ] is sufficiently small so that Mm < 1, i.e., M < 2m then
lim n = 0, i.e., lim xn = ξ.
M

n→∞ n→∞
::

4.4.2. Order of convergence of Regula-Falsi method. Let us consider the following alge-
ti

braic or transcendental equation

ai
M

f (x) = 0 (1)
The Regula-Falsi iterative formula for finding a real root ξ of (1) in a small interval [x0 , x1 ]
ar

such that f (x0 )f (x1 ) < 0, is

um

f (xn )(xn − xn−1 )

xn+1 = xn − for n = 1, 2, ...
f (xn ) − f (xn−1 )
K

(f (xn ) − f (xn−1 ))
as

⇒ (xn+1 − xn ) = −f (xn−1 )
(xn − xn−1 )
an

0
⇒ (xn+1 − xn )f (ξn ) = −f (xn ), ξn ∈ [xn−1 , xn ] [Applying mean value theorem]
M

0
⇒ (xn − xn+1 )f (ξn ) = f (xn )
0
⇒ (xn − ξ + ξ − xn+1 )f (ξn ) = f (xn )
r.
D

0 0 0 0
⇒ (−n + n+1 )f (ξn ) = f (xn ) − f (ξ) = (xn − ξ)f (ξn ), ξn ∈ [xn−1 , xn ]
0 0 0 0 0
f (ξ ) f (ξn ) − f (ξn ) M −m
⇒ n+1 = n − n 0 n = n 0 ≤ n
f (ξn ) f (ξn ) m
So, order of cnvergence of Regula falsi method is 1.
4.5. Secant method. Let us consider the following algebraic or transcendental equation
f (x) = 0 (1)
To solve (1) using Secant method we first search a small interval [x0 , x1 ] such that
f (x0 )f (x1 ) < 0, which implies, there is a real root of (1) in [x0 , x1 ]. Now draw a straight
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 47

Algorithm 14 Solution of a transcendental equation by Regula Falsi method

1: Input: Search interval [x0 , x1 ] and accuracy of significant figures k.
2: Define: Function f (x)
3: if f (x0 )f (x1 ) > 0 then
4: Output: A real root may not exists in [x0 , x1 ]
5: Stop execution.
end if

ge
6:
7:  ← 10−k

le
0 )(x0 −x1 )
8: x2 ← x0 − ff(x

ol
(x0 )−f (x1 )
9: while |x1 − x2 | ≥  and |x0 − x2 | ≥  do

C
10: if f (x0 )f (x2 ) < 0 then

aj
11: x1 ← x2
12: else

R
13: x0 ← x2

l
end if

da
14:
0 )(x0 −x1 )
15: x2 ← x0 − ff(x (x0 )−f (x1 )

ha
16: end while
17: Output: x2 , f (x2 )
is
ah
M

line (2) through the points (x0 , f (x0 )) and (x1 , f (x1 )).
::

y − f (x1 ) x − x1
= (2)
ti

f (x1 ) − f (x0 ) x1 − x0
ai

Now (2) cuts the x-axis at x = x2 , whose expression is presented in (3).

M

f (x1 )(x1 − x0 )
x2 = x1 − (3)
ar

f (x1 ) − f (x0 )
um

x2 is the first approximation of the real root of (1) in [x0 , x1 ]. Second approximation of
the real root is obtained using (4) by taking the new search interval [x1 , x2 ] or [x2 , x1 ]
K

according as x2 > x1 or x1 > x2 .

as

f (x2 )(x2 − x1 )
x3 = x2 − (4)
an

f (x2 ) − f (x1 )
Similarly, successive approximations of the real root is obtained using (5).
M

f (xn )(xn − xn−1 )

xn+1 = xn − , n = 1, 2, 3, ... (5)
r.

f (xn ) − f (xn−1 )
D

To get a result upto k significant figures the iteration stops when k significant figures of
two successive approximations becomes equal and final approximated value is taken as
the root of the equation with desired accuracy.
Example: Using Newton-Raphson method find the real root of x sin x + cos x = 0, near
x = π, correct upto 4 decimal places.
0
Ans: Let f (x) = x sin x + cos x, then f (x) = sin x + x cos x − sin x = x cos x and hence
iterative formula is
xn−1 sin xn−1 + cos xn−1
xn = xn−1 − , where x0 = π as we have to find the real root near π.
xn−1 cos xn−1
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 48

ge
le
Figure 6. Successive approximations of a real root of f (x) = 0 in Secant method

ol
Now successive approximations are tabulated below:

C
0
i xi f (xi ) f (xi ) xi+1

aj
0 3.142857 −1.003973 −3.142854 2.823412

R
1 2.823412 −0.066533 −2.681694 2.798602

l
3 2.798602 −0.000373 −2.635044 2.798460

da
4 2.798460 −0.000195 −2.635325 2.798386

ha
5 2.798386 −0.00000012 −2.635185 2.798386

is
So the required root is 2.7984 correct to four decimal places.
ah
Example: Using Bisection method find the real root of xex = 1, correct upto 4 significant
M

figures.
Ans: Let f (x) = xex − 1. Now f (0) = −1 < 0 and f (1) = e − 1 = 1.718 > 0, then there
::
ti

is a real root of f (x) = 0 in [0, 1]. Let a1 = 0, b1 = 1, then successive approximations of

ai

the root of (1) are tabulated below:

M

i ai bi f (ai ) f (bi ) ci f (ci )

ar

1 0 1 −1 1.718 0.5 −0.1756

2 0.5 1 −0.1756 1.718 0.75 2.8
um

3 0.5 0.75 −0.1756 2.8 .625 0.167

4 0.5 0.625 −0.1756 0.167 .5625 −0.0128
K

5 0.5625 0.625 −0.0128 0.167 0.59375 0.075

as

6 0.5625 0.59375 −0.0128 0.075 0.578125 0.031

an

7 0.5625 0.578125 −0.0128 0.031 0.5703125 0.009

8 0.5625 0.5703125 −0.0128 0.009 0.56640625 0.002
M

9 0.5625 0.56640625 −0.0128 0.002 0.564453125 −0.007

10 0.564453125 0.56640625 −0.007 0.002 0.5654578125 −0.00465
r.

11 0.5654578125 0.56640625 −0.00465 0.002 0.56593203125 −0.0033

D

12 0.56593203125 0.56640625 −0.0033 0.002 0.56616328125 −0.0027

13 0.56616328125 0.56640625 −0.0027 0.002 0.56628476563 −0.0024
14 0.56628476563 0.56640625 −0.0024 0.002 0.56637587891 −0.0021
15 0.56637587891 0.56640625 −0.0021 0.002 0.56639106446 −0.002
16 0.56639106446 0.56640625 −0.002 0.002 0.56639865723 −0.002
So the required root is 0.5664 correct to four decimal places.
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 49

5. Solution of ordinary differential equation

5.1. Euler’s method. Let us consider the following differential equation
dy
= f (x, y) (1)
dx
where y(a) = y0 , and we have to determine y(b)
To solve it using Euler’s method we first divide the interval [a, b] into n equal subintervals

ge
by equidistant points, x0 = a, x1 , x2 ,...,xn = b, where xi = a + ih, i = 0, 1, 2, ..., n,

le
h = (b − a)/n. We choice n sufficiently large so that second and higher power of h is

ol
negligible. Let values of y at x = x1 , x = x2 ,...,x = xn , be y1 , y2 , ...,yn respectively. Now
using Taylor’s infinite series we have

C
0 h2 00

aj
y(x1 ) = y(x0 + h) = y(x0 ) + hy (x0 ) + y (x0 ) + ...
2!

R
0
y1 = y(x0 ) + hy (x0 ) [Neglecting second and higher power of] h

l
y1 = y(x0 ) + hf (x0 , y0 )

da
Hence yi+1 = yi + hf (xi , yi ) for i = 1, 2, ..., n (-2)

ha
So using values of x0 and y0 we can determine y1 , the value of y at x = x1 . Then using
is
values of x1 , y1 we can determine y2 , the value of y at x = x2 and so on. Finally after n
ah
steps we shall get yn the value of y at x = xn = b. This process is called Euler’s method.
M

dy
Algorithm 15 Algorithm to find y(b) for dx
= f (x, y), y(a) = y0 by Euler method
::

1: Input: a, b, n, y0 .
ti

Define: Function f (x, y)

ai

2:
3: h ← (b − a)/n
M

4: x0 ← a
ar

5: for i = 1 to n do
6: y1 ← y0 + hf (x0 , y0 )
um

7: x0 ← a + ih
8: y0 ← y1
K

9: Output: x0 , y0
as

10: end for

an

5.2. Modified-Euler method. Let us consider the following differential equation

M

dy
= f (x, y) (1)
r.

dx
D

where y(a) = y0 , and we have to determine y(b)

To solve it using Modified-Euler method we first divide the interval [a, b] into n equal
subintervals by equidistant points, x0 = a, x1 , x2 ,...,xn = b, where xi = a + ih, i =
0, 1, 2, ..., n, h = (b − a)/n. We choice n sufficiently large so that second and higher power
of h is negligible. Let values of y at x = x1 , x = x2 ,...,x = xn , be y1 , y2 , ...,yn respectively.
Now using Taylor’s infinite series we have
0 h2 00
y(x1 ) = y(x0 + h) = y(x0 ) + hy (x0 ) + y (x0 ) + ...
2!
0
y1 = y(x0 ) + hy (x0 ) [Neglecting second and higher power of] h
y1 = y(x0 ) + hf (x0 , y0 )
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 50

This value of y1 is taken as the initial approximation of y1 and is denoted by y10 , i.e.,
y10 = y(x0 ) + hf (x0 , y0 )
This formula is called predictor formula of y1 . Then successive corrections are obtained
using the following formula
h
y1i = y(x0 ) + [f (x0 , y0 ) + f (x1 , y1i−1 )]
2

ge
The process continued until difference between two successive corrections becomes negli-

le
gible and final approximated value is taken as the value of y1 . Then using values of x1 , y1

ol
we can determine y2 , the value of y at x = x2 and so on. Finally after n steps we shall

C
get yn the value of y at x = xn = b. This process is called Modified-Euler method.

aj
dy
Algorithm 16 Find y(b) for = f (x, y), y(a) = y0 using modified Euler method

R
dx
1: Input: a, b, n, y0 and degree of precession k.

l
da
2: Define: Function f (x, y)
3:  ← 10−k

ha
4: h ← (b − a)/n
5: x0 ← a
is
ah
6: for i = 1 to n do
7: y1p ← y0 + hf (x0 , y0 )
M

8: y1c ← y0 + hf (x0 , y0 ) + f (x0 + h, y1p )

::

9: while |y1p − y1c | >  do

ti

10: y1p ← y1c

ai

11: y1c ← y0 + (h/2)[f (x0 , y0 ) + f (x0 + h, y1p )]

M

12: end while

13: x0 ← a + ih
ar

14: y0 ← y1c
um

15: Output: x0 , y0
16: end for
K
as

5.3. 2nd order Runge-Kutta method. Let us consider the following differential equa-
an

tion
M

dy
= f (x, y) (1)
dx
r.

where y(a) = y0 , and we have to determine y(b)

D

To solve it using 2nd order Runge-Kutta method we first divide the interval [a, b] into n
equal subintervals by equidistant points, x0 = a, x1 , x2 ,...,xn = b, where xi = a + ih, i =
0, 1, 2, ..., n, h = (b − a)/n. We choice n sufficiently large so that second and higher power
of h is negligible. Let values of y at x = x1 , x = x2 ,...,x = xn , be y1 , y2 , ...,yn respectively.
Now using 2nd order Runge-Kutta method successive values of y1 is determined using the
following rule:
k1 = hf (x0 , y0 )
k2 = hf (x0 + h, y0 + k1 )
y1 = y0 + (k1 + k2 )/2
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 51

So using the values of x0 and y0 we can determine y1 , the value of y at x = x1 . Similarly

successive values of y are determined as
k1 = hf (xi , yi )
k2 = hf (xi + h, yi + k1 )
yi+1 = yi + (k1 + k2 )/2 for i = 1, 2, ..., n (-4)
So after n steps we shall get yn the value of y at x = xn = b. This process is called 2nd

ge
order Runge-Kutta method.

le
ol
dy
Algorithm 17 Find y(b) for dx
= f (x, y), y(a) = y0 using 2nd order R-K method

C
1: Input: a, b, n, y0 .

aj
2: Define: Function f (x, y)
h ← (b − a)/n

R
3:
4: x0 ← a

l
da
5: for i = 1 to n do
6: k1 ← hf (x0 , y0 )

ha
7: k2 ← hf (x0 + h, y0 + k1 )
8: y1 ← y0 + [k1 + k2 ]/2
is
ah
9: x0 ← a + ih
10: y0 ← y1
M

11: Output: x0 , y0
::

12: end for

ti
ai
M

5.4. 4-th order Runge-Kutta method. Let us consider the following differential equa-
tion
ar

dy
um

= f (x, y) (1)
dx
where y(a) = y0 , and we have to determine y(b)
K

To solve it using 4-th order Runge-Kutta method we first divide the interval [a, b] into n
as

equal subintervals by equidistant points, x0 = a, x1 , x2 ,...,xn = b, where xi = a + ih, i =

an

0, 1, 2, ..., n, h = (b − a)/n. We choice n sufficiently large so that second and higher power
M

of h is negligible. Let values of y at x = x1 , x = x2 ,...,x = xn , be y1 , y2 , ...,yn respectively.

Now using 2nd order Runge-Kutta method successive values of y1 is determined using the
r.

following rule:
D

k1 = hf (x0 , y0 ), k2 = hf (x0 + h/2, y0 + k1 /2)

k3 = hf (x0 + h/2, y0 + k2 /2), k4 = hf (x0 + h, y0 + k3 )
y1 = y0 + (k1 + 2k2 + 2k3 + k4 )/6
So using the values of x0 and y0 we can determine y1 , the value of y at x = x1 . Similarly
successive values of y are determined as
k1 = hf (xi , yi ), k2 = hf (xi + h/2, yi + k1 /2)
k3 = hf (xi + h/2, yi + k2 /2), k4 = hf (xi + h, yi + k3 )
yi+11 = y0 + (k1 + 2k2 + 2k3 + k4 )/6 for i = 1, 2, ..., n
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 52

So after n steps we shall get yn the value of y at x = xn = b. This process is called 4-th
order Runge-Kutta method.

dy
Algorithm 18 Find y(b) for dx
= f (x, y), y(a) = y0 using 4-th order R-K method
1: Input: a, b, n, y0 .
2: Define: Function f (x, y)

ge
3: h ← (b − a)/n

le
4: x0 ← a

ol
5: for i = 1 to n do
k1 ← hf (x0 , y0 )

C
6:
7: k2 ← hf (x0 + h/2, y0 + k1 /2)

aj
8: k3 ← hf (x0 + h/2, y0 + k2 /2)

R
9: k4 ← hf (x0 + h, y0 + k3 )
10: y1 ← y0 + [k1 + 2k2 + 2k3 + k4 ]/6

l
da
11: x0 ← a + ih
y0 ← y1

ha
12:
13: Output: x0 , y0
14: end for
is
ah
M

Example: Solve the following differential equation by second order Runge-Kutta method:
::
ti
ai

dy
= x + y2 (1)
M

dx
where y(0) = 1, and we have to determine y(0.2)
ar
um

Here f (x, y) = x + y 2 , x0 = 0, y0 = 1 and we take h=0.1. Then using 2-th order

Runge-Kutta method we have
K

k1 = hf (x0 , y0 ) = 0.1 × f (0, 1) = 0.1{0 + 12 } = 0.1

as

k2 = hf (x0 + h, y0 + k1 ) = 0.1 × f (0.1, 1.1) = 0.1{0.1 + (1.1)2 } = 0.131

an

y1 = y0 + (k1 + k2 )/2 = 1 + (0.1 + 0.131)/2 = 1.1155

M

So value of y at x = x1 = x0 + h = 0 + 0.1 = 0.1 is y1 = 1.1155. Now we shall determine

y2 using values of x1 and y1 as below.
r.
D

k1 = hf (x1 , y1 ) = 0.1 × f (0, 1) = 0.1{0.1 + (1.1155)2 } = 0.1344

k2 = hf (x1 + h, y1 + k1 ) = 0.1 × f (0.2, 1.1155 + 0.1344) = 0.1{0.2 + (1.2499)2 } = 0.1762
y2 = y1 + (k1 + k2 )/2 = 1.1155 + (0.1344 + 0.1762)/2 = 1.2708

Therefore y(0.2)=1.2708.
Example: Solve the following differential equation by 4-th order Runge-Kutta method:
dy
= x+y (1)
dx
where y(0) = 1, and we have to determine y(0.2)
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 53

Here f (x, y) = x + y, x0 = 0, y0 = 1 and we take h=0.1. Then using 4-th order Runge-
Kutta method we have
k1 = hf (x0 , y0 ) = 0.1 × f (0, 1) = 0.1{0 + 1} = 0.1
k2 = hf (x0 + h/2, y0 + k1 /2) = 0.1 × f (0.1/2, 1 + 0.1/2) = 0.1{0.05 + 1.05} = 0.110
k3 = hf (x0 + h/2, y0 + k2 /2) = 0.1 × f (0.1/2, 1 + 0.110/2) = 0.1{0.05 + 1.055} = 0.1105
k4 = hf (x0 + h, y0 + k3 ) = 0.1 × f (0.1, 1 + 0.1105) = 0.1{0.05 + 1.1105} = 0.1205

ge
y1 = y0 + (k1 + 2k2 + 2k3 + k4 )/6 = 1 + (0.1 + 2 × 0.11 + 2 × 0.1105 + 0.1205)/6 = 1.11025

le
So value of y at x = x1 = x0 + h = 0 + 0.1 = 0.1 is y1 = 1.11025. Now we shall determine

ol
y2 using values of x1 and y1 as below.

C
k1 = hf (x1 , y1 ) = 0.1 × f (0.1, 1.11025) = 0.1{0.1 + 1.11025} = 0.121025

aj
k2 = hf (x1 + h/2, y1 + k1 /2) = 0.1 × f (0.1 + 0.1/2, 1.11025 + 0.121025/2)

R
= 0.1{0.15 + 1.1707625} = 0.13207625

l
da
k3 = hf (x1 + h/2, y1 + k2 /2) = 0.1 × f (0.1 + 0.1/2, 1.11025 + 0.13207625/2)

ha
= 0.1{0.15 + 1.176288125} = 0.1326288125
k4 = hf (x1 + h, y1 + k3 ) = 0.1 × f (0.1 + 0.1, 1.11025 + 0.1326288125)
is
ah
= 0.1{0.2 + 1.12428788125} = 0.14428788125
y2 = y1 + (k1 + 2k2 + 2k3 + k4 )/6
M

= 1.11025 + (0.121025 + 2 × 0.13207625 + 2 × 0.1326288125 + 0.14428788125)/6

::
ti

= 1.24270
ai

Therefore y(0.2)=1.24270.
M

Example: Solve the following differential equation by Euler method:

ar

dy
= xy (1)
um

dx
where y(0) = 1, and we have to determine y(1) taking h = 0.2
K

Here f (x, y) = xy, x0 = 0, y0 = 1 and we take h=0.2. Then using Euler method method
as

we have
an

y(0.2) = y1 = y0 + hf (x0 , y0 ) = 1 + 0.2 × f (0, 1) = 1 + 0.2{0 × 1} = 1

y(0.4) = y2 = y1 + hf (x1 , y1 ) = 1 + 0.2 × f (0.2, 1) = 1 + 0.2{0.2 × 1} = 1.04
M

y(0.6) = y3 = y2 + hf (x2 , y2 ) = 1.04 + 0.2 × f (0.4, 1.04)

r.

= 1 + 0.2{0.4 × 1.04} = 1.1232

D

y(0.8) = y4 = y3 + hf (x3 , y3 ) = 1.1232 + 0.2 × f (0.6, 1.1232)

= 1 + 0.2{0.6 × 1.1232} = 1.25798
y(1.0) = y5 = y4 + hf (x4 , y4 ) = 1.25798 + 0.2 × f (0.8, 1.25798)
= 1 + 0.2{0.8 × 1.25798} = 1.45926
Hence y(1) = 1.45926
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 54

6. Solution of system of linear equations

6.1. Gauss Elimination Method. Consider the following system of linear equations

a11 x1 + a12 x2 + a13 x3 + ... + a1n−1 xn−1 + a1n xn = b1

a21 x1 + a22 x2 + a23 x3 + ... + a2n−1 xn−1 + a2n xn = b2
a31 x1 + a32 x2 + a33 x3 + ... + a3n−1 xn−1 + a3n xn = b3

ge
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...

le
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...

ol
an−11 x1 + an−12 x2 + an−13 x3 + ... + an−1n−1 xn−1 + an−1n xn = bn−1

C
an1 x1 + an2 x2 + an3 x3 + ... + ann−1 xn−1 + ann xn = bn

aj
The process consists of the following steps:

R
• If the system of linear equations is not diagonally dominant rearrange the equations

l
da
in a diagonally dominant form. Let the system is diagonally dominant.
• Determine the multipliers M1i = aa11 i1
for equations i = 2, 3, ..., n.

ha
• Multiply the first equation by M1i and subtract the same from the i-th equation
for i = 2, 3, ..., n. is
ah
• After the completion of this operation the system of equation reduces to
M

a11 x1 + a12 x2 + a13 x3 + ... + a1n−1 xn−1 + a1n xn = b1

::

a122 x2 + a123 x3 + ... + a12n−1 xn−1 + a12n xn = b12

ti
ai

a132 x2 + a133 x3 + ... + a13n−1 xn−1 + a13n xn = b13

M

... ... ... ... ... ... ... ... ... ... ... ... ... ...
... ... ... ... ... ... ... ... ... ... ... ... ... ...
ar

a1n−12 x2 + a1n−13 x3 + ... + a1n−1n−1 xn−1 + a1n−1n xn = b1n−1

um

a1n2 x2 + a1n3 x3 + ... + a1nn−1 xn−1 + a1nn xn = b1n

K

where a1ij = aij −a1j M1i for i, j = 2, 3, ..., n and b1i = bi −b1 M1i for i = 3, 4, ..., n.
as

a1
• Determine the multipliers M2i = a1i2 for equations i = 3, 4, ..., n.
an

22
• Multiply the second equation by M2i and subtract the same from the i-th equation
M

for i = 3, 4, ..., n.
• After the completion of this operation the system of equation reduces to
r.
D

a11 x1 + a12 x2 + a13 x3 + ... + a1n−1 xn−1 + a1n xn = b1

a122 x2 + a123 x3 + ... + a12n−1 xn−1 + a12n xn = b12
a233 x3 + ... + a23n−1 xn−1 + a23n xn = b23
... ... ... ... ... ... ... ... ... ... ...
... ... ... ... ... ... ... ... ... ... ...
2 2 2
an−13 x3 + ... + an−1n−1 xn−1 + an−1n xn = b2n−1
a2n3 x3 + ... + a2nn−1 xn−1 + a2nn xn = b2n

where a2ij = a1ij −a12j M2i for i, j = 3, 4, ..., n and b2i = b1i −b12 M2i for i = 3, 4, ..., n.
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 55

• Proceeding this way after (n − 1)-th iteration we have

a11 x1 + a12 x2 + a13 x3 + a14 x4 + ... + a1n−1 xn−1 + a1n xn = b1
a122 x2 + a123 x3 + a124 x4 + ... + a12n−1 xn−1 + a12n xn = b12
a233 x3 + a234 x4 + ... + a23n−1 xn−1 + a23n xn = b23
a344 x4 + ... + a33n−1 xn−1 + a33n xn = b34

ge
... ... ... ... ... ... ... ... ...
... ... ... ... ... ... ... ... ...

le
n−2 n−2 n−2

ol
an−1n−1 xn−1 + an−1n xn = bn−1

C
an−1 n−1
nn xn = bn

aj
n−1
• Then from the last equation we have xn = abnn−1 . Putting this value in the second
nn

R
equation from last we shall get value of xn−1 . Then putting the values of xn and
xn−1 in the third equation from last we shall get the value of xn−2 . In this way we

l
da
shall get the values of successive variables.

ha
Algorithm 19 Solution of a system of linear equations by Gauss Elimination method
1: Input: n and augmented matrix A = (aij )n×n+1 . is
ah
2: for k = 1 to n − 1 do
M

3: for i = k + 1 to n do
M = aakk
::

ik
4:
ti

5: for j = k to n + 1 do
ai

6: aij = aij − M akj

M

7: end for
8: end for
ar

9: end for
um

10: for i = n to 1 do
11: xi = ain+1
K

12: for j = i + 1 to n do
13: xi = xi − aij xj
as

14: end for

an

15: xi = axiii
M

16: end for

17: Output: x1 , x2 , ..., xn
r.
D

6.2. Gauss Jacobi iteration method. Consider the system of linear equations (1)
a11 x1 + a12 x2 + a13 x3 + ... + a1n−1 xn−1 + a1n xn = b1 

a21 x1 + a22 x2 + a23 x3 + ... + a2n−1 xn−1 + a2n xn = b2 



a31 x1 + a32 x2 + a33 x3 + ... + a3n−1 xn−1 + a3n xn = b3

(1)
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... 

... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... 



an1 x1 + an2 x2 + an3 x3 + ... + ann−1 xn−1 + ann xn = bn
The process to solve (1) using Gauss Jacobi iteration method consists of the following
steps:
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 56

• If the system of linear equations is not diagonally dominant rearrange the equations
in a diagonally dominant form. Let the system is diagonally dominant.
• Now rewrite the equation in the form
x1 = [b1 − (a12 x2 + a13 x3 + ... + a1n−1 xn−1 + a1n xn )]/a11


x2 = [b2 − a21 x1 − (a23 x3 + ... + a2n−1 xn−1 + a2n xn )]/a22




x3 = [b3 − (a31 x1 + a32 x2 ) − (a34 x4 + ... + a3n−1 xn−1 + a3n xn )]/a33



... ... ... ... ... ... ... ... ... ... ... ... ... ... (2)

ge
xi = bi − (ai1 x1 + ai2 x2 + ... + aii−1 xi−1 ) − (ai+11 xi+1 + ... + ain xn )]/aii 



le

... ... ... ... ... ... ... ... ... ... ... ... ... ... 


ol


xn = [bn − (an1 x1 + an2 x2 + an3 x3 + ... + ann−1 xn−1 )]/ann

C
• Let x1 = x01 , x2 = x02 , ..., xn = x0n be the initial guess of the solution of the system

aj
of equation (1).

R
• Then putting these values in the right hand side of the i-th equation of (2) we get
the first approximation x1i of xi for i = 1, 2, ..., n as presented below.

l
da
x11 = [b1 − (a12 x02 + a13 x03 + ... + a1n−1 x0n−1 + a1n x0n )]/a11


ha

x12 = [b2 − a21 x01 − (a23 x03 + ... + a2n−1 x0n−1 + a2n x0n )]/a22




1 0 0 0 0 0
x3 = [b3 − (a31 x1 + a32 x2 ) − (a34 x4 + ... + a3n−1 xn−1 + a3n xn )]/a33

is


... ... ... ... ... ... ... ... ... ... ... ... ... ... (3)
ah
x1i = bi − (ai1 x01 + ai2 x02 + ... + aii−1 x0i−1 ) − (ai+11 x0i+1 + ... + ain x0n )]/aii 


M


... ... ... ... ... ... ... ... ... ... ... ... ... ... 


1 0 0 0 0

xn = [bn − (an1 x1 + an2 x2 + an3 x3 + ... + ann−1 xn−1 )]/ann
::
ti

• Putting the values of first approximations in the right hand side of the successive
ai

equations of (2) we shall get the second approximations and so on. In this way
M

k-th approximations of the solutions of (1) for k = 1, 2, 3, ... are given by

ar


xk1 = [b1 − (a12 xk−1 2 + a13 xk−1
3 + ... + a1n−1 xk−1 k−1
n−1 + a1n xn )]/a11 

k−1 k−1
xk2 = [b2 − a21 x1 − (a23 x3 + ... + a2n−1 xn−1k−1 + a2n xn )]/a22 k−1
um





k−1 k−1 k−1 k−1
xk3 = [b3 − (a31 x1 + a32 x2 ) − (a34 x4 + ... + a3n−1 xn−1 + a3n xn )]/a33 k−1 


K

... ... ... ... ... ... ... ... ... ... ... ... ... ... (4)
k−1 k−1 k−1 k−1
xki k−1
= bi − (ai1 x1 + ai2 x2 + ... + aii−1 xi−1 ) − (aii+1 xi+1 + ... + ain xn )]/aii 


as


... ... ... ... ... ... ... ... ... ... ... ... ... ...



an

k−1 k−1 k−1 k−1


xkn = [bn − (an1 x1 + an2 x2 + an3 x3 + ... + ann−1 xn−1 )]/ann

M

• The process continues until the sum of the absolute differences of the solutions in
two successive differences becomes negligible, i.e., when
r.

" n #
D

X
k k−1
|xi − xi | < 
i=1

for a predefined small value . Then x1 = xk1 , x2 = xk2 , ..., xn = xkn is taken as
solution of (1).
• To get a result correct upto d decimal places, the iteration continues until the
absolute differences of the each variables in two successive iterations becomes less
than  = 10−d , i.e., when
|xki − xik−1 | <  for i = 1, 2, ..., n
Then x1 = xk1 , x2 = xk2 , ..., xn = xkn is taken as the solution of (1).
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 57

Algorithm 20 Solution of a system of linear equations by Gauss-Jacobi iteration method

1: Input: n, k and augmented matrix A = (aij )n×n+1 .
2:  = 10−k
3: for i = 1 to n do
4: x1i = 0 . Initial guess
5: end for
Continue = 1

ge
6:
7: while Continue do

le
8: for i = 1 to n do

ol
9: x2i = ain+1

C
10: for j = 1 to i − 1 do
x2i = x2i − aij x1j

aj
11:
12: end for

R
13: for j = i + 1 to n do

l
14: x2i = x2i − aij x1j

da
15: end for

ha
16: i = |x2i − x1i |
17: end for
Continue = 0 is
ah
18:
19: for i = 1 to n do
M

20: if i >  then

::

21: Continue = 1
ti

22: Break
ai

23: end if
M

24: end for

25: end while
ar

26: Output: x1 , x2 , ..., xn

um
K

Theorem. A sufficient condition for the convergence of the Gauss-Jacobi iteration

method for the system of linear equations (1)
as

a11 x1 + a12 x2 + a13 x3 + ... + a1n−1 xn−1 + a1n xn = b1 


an

a21 x1 + a22 x2 + a23 x3 + ... + a2n−1 xn−1 + a2n xn = b2 



M


a31 x1 + a32 x2 + a33 x3 + ... + a3n−1 xn−1 + a3n xn = b3

(1)
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... 
r.


... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... 

D



an1 x1 + an2 x2 + an3 x3 + ... + ann−1 xn−1 + ann xn = bn
n
X
is |aij | < |aii | for i = 1, 2, ..., n.
i=1,i6=j

Proof. For an initial guess x1 = x01 , x2 = x02 , ..., xn = x0n the k-th and k + 1-th approxi-
mation of the solution of (1) are given by (2) and (3) respectively for i = 1, 2, ..., n.
xki = bi − (ai1 x1k−1 + ai2 xk−1
2 + ... + aii−1 xk−1 k−1 k−1
i−1 ) − (ai+11 xi+1 + ... + ain xn )]/aii (2)
xk+1
i = bi − (ai1 xk1 + ai2 xk2 + ... + aii−1 xki−1 ) − (ai+11 xki+1 + ... + ain xkn )]/aii (3)
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 58

Let k+1
i be the difference of the k-th approximation and k + 1-th approximation of xi ,
k+1
i.e., i = |xk+1i − xki | for i = 1, 2, ..., n and k = 1, 2, .... Then we have from (2) and (3)
" i−1 n
#
1 X X
xki − xk+1
i = aij (xjk−1 − xkj ) + aij (xjk−1 − xkj )
aii j=1 j=i+1
" n #
1 X
⇒ |xki − xk+1 |≤ |aij ||(xk−1 − xkj )|

ge
i j
|aii | j=1, j6=i

le
" n #
1

ol
X
⇒ k+1
i ≤ |aij |kj
|aii | j=1, j6=i

C
aj
" n #
1 X
⇒ k+1 ≤ |aij |k where k = sup{ki , i = 1, 2, ..., n} for k = 1, 2, ...

R
|aii | j=1, j6=i

l
n

da
k+1 k 1 X
⇒  ≤  K where K = |aij |
|aii | j=1, i6=j

ha
k+1 ≤ K k 1 → 0 as k → ∞ [as K < 1]
⇒
is
ah
⇒ {xki } is a convergent sequence for i = 1, 2, ..., n
M

⇒ The iterative process is convergent

::

6.3. Gauss-Seidel iteration method. Consider the system of linear equations (1)
ti
ai

a11 x1 + a12 x2 + a13 x3 + ... + a1n−1 xn−1 + a1n xn = b1 


M

a21 x1 + a22 x2 + a23 x3 + ... + a2n−1 xn−1 + a2n xn = b2 




a31 x1 + a32 x2 + a33 x3 + ... + a3n−1 xn−1 + a3n xn = b3

ar

(1)
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... 
um


... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... 



an1 x1 + an2 x2 + an3 x3 + ... + ann−1 xn−1 + ann xn = bn
K

The process to solve (1) using Gauss-Seidel iteration method consists of the following
as

steps:
an

• If the system of linear equations is not diagonally dominant rearrange the equations
in a diagonally dominant form. Let the system is diagonally dominant.
M

• Now rewrite the equation in the form

r.

x1 = [b1 − (a12 x2 + a13 x3 + ... + a1n−1 xn−1 + a1n xn )]/a11



D

x2 = [b2 − a21 x1 − (a23 x3 + ... + a2n−1 xn−1 + a2n xn )]/a22





x3 = [b3 − (a31 x1 + a32 x2 ) − (a34 x4 + ... + a3n−1 xn−1 + a3n xn )]/a33



... ... ... ... ... ... ... ... ... ... ... ... ... ... (2)
xi = bi − (ai1 x1 + ai2 x2 + ... + aii−1 xi−1 ) − (aii+1 xi+1 + ... + ain xn )]/aii 



... ... ... ... ... ... ... ... ... ... ... ... ... ... 



xn = [bn − (an1 x1 + an2 x2 + an3 x3 + ... + ann−1 xn−1 )]/ann
• Let x1 = x01 , x2 = x02 , ..., xn = x0n be the initial guess of the solution of the system
of equation (1).
• Then putting these values in the right hand side of the first equation of (2) we
get the first approximation x11 of x1 . Putting x1 = x11 , x3 = x03 , ..., xn = x0n in the
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 59

second equation of (2) we get the first approximation x12 of x2 . Similarly putting
x1 = x11 , x2 = x12 , ..., xi−1 = x1i−1 , xi+1 = x0i+1 ..., xn = x0n in the second equation
of (2) we get the first approximation x1i of xi for i = 1, 2, ..., n as presented below.
x11 = [b1 − (a12 x02 + a13 x03 + ... + a1n−1 x0n−1 + a1n x0n )]/a11


x12 = [b2 − a21 x11 − (a23 x03 + ... + a2n−1 x0n−1 + a2n x0n )]/a22




1 1 1 0 0 0
x3 = [b3 − (a31 x1 + a32 x2 ) − (a34 x4 + ... + a3n−1 xn−1 + a3n xn )]/a33




ge
... ... ... ... ... ... ... ... ... ... ... ... ... ... (3)
1 1 1 1 0 0
xi = bi − (ai1 x1 + ai2 x2 + ... + aii−1 xi−1 ) − (aii+1 xi+1 + ... + ain xn )]/aii  

le


... ... ... ... ... ... ... ... ... ... ... ... ... ... 

ol


1 1 1 1 1

xn = [bn − (an1 x1 + an2 x2 + an3 x3 + ... + ann−1 xn−1 )]/ann

C
• Then putting these values in the right hand side of the first equation of (3) we

aj
get the second approximation x21 of x1 . Putting x1 = x21 , x3 = x13 , ..., xn = x1n in

R
the second equation of (3) we get the second approximation x22 of x2 . Similarly

l
putting x1 = x21 , x2 = x22 , ..., xi−1 = x2i−1 , xi+1 = x1i+1 ..., xn = x1n in the second

da
equation of (2) we get the second approximation x2i of xi for i = 1, 2, ..., n, and

ha
so on. In this way k-th approximations of the solutions of (1) for k = 1, 2, 3, ...
are given by
is
ah

xk1 = [b1 − (a12 xk−1 2 + a13 xk−1
3 + ... + a1n−1 xn−1k−1
+ a1n xk−1
n )]/a11 

k k k−1 k−1
M

x2 = [b2 − a21 x1 − (a23 x3 + ... + a2n−1 xn−1k−1 + a2n xn )]/a22





k k k k−1 k−1 k−1
x3 = [b3 − (a31 x1 + a32 x2 ) − (a34 x4 + ... + a3n−1 xn−1 + a3n xn )]/a33

::



(4)
ti

... ... ... ... ... ... ... ... ... ... ... ... ... ...
k−1
ai

k k k k k−1
xi = bi − (ai1 x1 + ai2 x2 + ... + aii−1 xi−1 ) − (aii+1 xi+1 + ... + ain xn )]/aii 



M

... ... ... ... ... ... ... ... ... ... ... ... ... ...




k k k k k
xn = [bn − (an1 x1 + an2 x2 + an3 x3 + ... + ann−1 xn−1 )]/ann

ar

• The process continues until the sum of the absolute differences of the solutions in
um

two successive differences becomes negligible, i.e., when

K

" n #
X
|xki − xik−1 | < 
as

i=1
an

for a predefined small value . Then x1 = xk1 , x2 = xk2 , ..., xn = xkn is taken as the
M

solution of (1).
• To get a result correct upto d decimal places, the iteration continues until the
r.

absolute differences of the each variables in two successive iterations becomes less
D

than  = 10−d , i.e., when |xki − xk−1

i | <  for i = 1, 2, ..., n. Then x1 = xk1 , x2 =
k k
x2 , ..., xn = xn is taken as the solution of (1).
Theorem. A sufficient condition for the convergence of the Gauss-Seidel iteration method
for the system of linear equations (1)

a11 x1 + a12 x2 + a13 x3 + ... + a1n−1 xn−1 + a1n xn = b1 
a21 x1 + a22 x2 + a23 x3 + ... + a2n−1 xn−1 + a2n xn = b2 



a31 x1 + a32 x2 + a33 x3 + ... + a3n−1 xn−1 + a3n xn = b3

(1)
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... 

... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... 



an1 x1 + an2 x2 + an3 x3 + ... + ann−1 xn−1 + ann xn = bn

 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 60

n
X
is |aij | < |aii | for i = 1, 2, ..., n.
i=1, i6=j

Proof. For an initial guess x1 = x01 , x2 = x02 , ..., xn = x0n the k-th and k + 1-th approxi-
mation of the solution of (1) are given by (2) and (3) respectively for i = 1, 2, ..., n.
xki = bi − (ai1 xk1 + ai2 xk2 + ... + aii−1 xki−1 ) − (ai+11 xk−1 k−1
i+1 + ... + ain xn )]/aii (2)

ge
xk+1 = bi − (ai1 xk+1 + ai2 xk+1 + ... + aii−1 xk+1
i−1 ) − (ai+11 xki+1 + ... + ain xkn )]/aii (3)

le
i 1 2

Let k+1

ol
i be the difference of the k-th approximation and k + 1-th approximation of xi ,
k+1
i.e., i = |xk+1 − xki | for i = 1, 2, ..., n and k = 1, 2, .... Then we have from (2) and (3)

C
i
 

aj
i−1 n
1 X X
xki − xk+1 = aij (xk+1 − xkj ) + aij (xkj − xk−1 )

R

i j j
aii
j=1 j=i+1

l
 

da
i−1 n
1 X X
⇒ |xki − xk+1
i |≤  |aij ||(xk+1
j − xkj )| + |aij ||(xkj − xk−1
j )|
|aii |

ha
j=1 j=i+1

is
 
i−1 n
k+1 1 X k+1
X
|aij |kj 
ah
⇒ i ≤ |aij |j +
|aii |
j=1 j=i+1
M

 
i−1 n
1 
::

X X
⇒ k+1 ≤ |aij |k+1 + |aij |k 
ti

|aii |
j=1 j=i+1
ai

k
where  = sup{ki , i = 1, 2, ..., n} for k = 1, 2, ...
M

 
i−1 i−1 n
1  X X X
ar

⇒ k+1 ≤ |aij |k+1 − |aij |k + |aij |k 

|aii |
um

j=1 j=1 j=1, j6=i

i−1 n
1 X 1 X
⇒ k+1 ≤ k+1 Ki − k Ki + k K, where Ki = |aij |, K = |aij |
K

aii |aii |
j=1 j=1, i6=j
as

K − Ki k
⇒ k+1 ≤  = αk , where α < 1 as K < 1
an

1 − Ki
⇒ k+1 ≤ αk 1 → 0 as k → ∞
M

⇒ {xki } is a convergent sequence for i = 1, 2, ..., n

r.

⇒ The iterative process is convergent

D

6.4. LU decomposition method or Matrix factorization method. Consider the

system of linear equations (1)
a11 x1 + a12 x2 + a13 x3 + ... + a1n−1 xn−1 + a1n xn = b1 

a21 x1 + a22 x2 + a23 x3 + ... + a2n−1 xn−1 + a2n xn = b2 



a31 x1 + a32 x2 + a33 x3 + ... + a3n−1 xn−1 + a3n xn = b3

(1)
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... 

... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... 



an1 x1 + an2 x2 + an3 x3 + ... + ann−1 xn−1 + ann xn = bn
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 61

Algorithm 21 Solution of a system of linear equations by Gauss-Seidal iteration method

1: Input: n, k and augmented matrix A = (aij )n×n+1 .
2:  = 10−k
3: for i = 1 to n do
4: x1i = 0 . Initial guess
5: end for
Continue = 1

ge
6:
7: while Continue do

le
8: for i = 1 to n do

ol
9: x2i = ain+1

C
10: for j = 1 to i − 1 do
x2i = x2i − aij x2j

aj
11:
12: end for

R
13: for j = i + 1 to n do

l
14: x2i = x2i − aij x1j

da
15: end for

ha
16: i = |x2i − x1i |
17: end for
Continue = 0 is
ah
18:
19: for i = 1 to n do
M

20: if i >  then

::

21: Continue = 1
ti

22: Break
ai

23: end if
M

24: end for

25: end while
ar

26: Output: x1 , x2 , ..., xn

um

The process to solve (1) using decomposition method consists of the following steps:
K

• Write the system of equation (1) in the form

as

AX = B (2)
an

a11 a12 a13 ... a1n−1 a1n x1 b1

M

     
 a21 a22 a23 ... a2n−1 a2n   x2   b2 
r.

a31 a32 a33 ... a3n−1 a3n x3  b3 

     
where A =  , X = , B =
   
D

... ... ... ... ... ... ...  ... 

    
   
 an−11 an−12 an−13 ... an−1n−1 an−1n   ...   ... 
an1 an2 an3 ... ann−1 ann xn bn
• Write the coefficient matrix A in the form A = LU where
l11 0 0 ... 0 0 1 u12 u13 ... u1n−1 u1n
   
 l21 l22 0 ... 0 0   0 1 u23 ... u2n−1 u2n 
 l31 l32 l33 ... 0 0   0 0 1 ... u3n−1 u3n
   
L= U = 

 ... ... ... ... ... ...   ... ... ... ... ... ...


 l l l ... ln−1n−1 0   ... ... ... ... 1 un−1n 
n−11 n−12 n−13
ln1 ln2 ln3 ... lnn−1 lnn 0 0 0 ... 0 1
Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 62

• Then equation (2) reduces to

LU X = B (3)

• From LU = A, equating the elements of the first column of both matrices we have

li1 = ai1 , for i = 1, 2, ..., n

ge
• From LU = A, equating the elements of the first row of both matrices we have

le
ol
a1j
l11 = a11 , l11 u1j = a1j ⇒ u1j = , for j = 2, 3, ..., n

C
l11

aj
• From LU = A, equating the elements of the second column of both matrices we

R
have
li1 u12 + li2 = ai2 ⇒ li2 = ai2 − li1 u12 , for i = 2, 3, ..., n

l
da
• From LU = A, equating the elements of the second row of both matrices we have

ha
l21 u1j + l22 u2j = a2j ⇒ u2j = is
a2j − l21 u1j
, for j = 3, 4, ..., n
ah
l22
M

• From LU = A, equating the elements of the third column of both matrices we

::

have
ti
ai

li1 u13 + li2 u23 + li3 = ai3 ⇒ li3 = ai3 − li1 u13 − li2 u23 , for i = 3, 4, ..., n
M

• From LU = A, equating the elements of the third row of both matrices we have
ar

a3j − l31 u1j − l32 u2j

um

l31 u1j + l32 u2j + l33 u3j = a3j ⇒ u3j = , for j = 4, 5, ..., n
l33
K

• Take U X = Z = [z1 z2 ... zn ]T and then (3) reduces to

as

l11 0 0 ... 0 0 z1 b1
    
an

 l21 l22 0 ... 0 0  z2   b2 

M

l31 l32 l33 ... 0 0 z3 b3

    
= (4)
    
... ... ... ... ... ... ... ...
  
    
r.

 l
n−11 ln−12 ln−13 ... ln−1n−1 0 ... ...
   
D

ln1 ln2 ln3 ... lnn−1 lnn zn bn

• System of equation (4) can be solved by forward substitution as presented pre-

sented here. From the i-th equation of (4) we have

li1 z1 + li2 z2 + ... + lii−1 zi−1 + lii zi = bi (5)

⇒ lii zi = bi − [li1 z1 + li2 z2 + ... + lii−1 zi−1 ] (6)
i−1
P
bi − lij zj
j=1
⇒ zi = f ori = 1, 2, ..., m (7)
lii
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 63

After getting values of z1 , z2 , ..., zn from the relation U X = Z we have the system
of equations (8).
1 u12 u13 ... u1n−1 u1n x1 z1
    
 0 1 u23 ... u2n−1 u2n  x2   z2 
 0 0 1 ... 0 0 x3 z3
    
= (8)
   
... ... ... ... ... ... ... ...
  
    

ge
 ... ... ... ... 1 un−1n  ...   ... 
0 0 0 ... 0 1 xn zn

le
ol
Now, from the i-th equation of (8) we have

C
xi + uii+1 xi+1 + ... + uim−1 xm−1 + uii+1 xm = zi (9)

aj
⇒ xi = zi − [uii+1 xi+1 + ... + uim−1 xm−1 + uim xm ]f ori = m, m − 1, ..., 1 (10)

l R
da
Algorithm 22 Solution of a system of linear equations by LU decomposition method

ha
1: Input: n and augmented matrix A = (aij )n×n+1 .
2: for k = 1 to n do
for i = k to n do lik = aik is
ah
3:
4: for j = 1 to k − 1 do
M

5: lik = lik − lij ujk

::

6: end for
ti

7: end for
ai

8: for j = k + 1 to n do
M

9: ukj = akj
10: for i = 1 to k − 1 do
ar

11: ukj = ukj − lki uij

um

12: end for

ukj
13: ukj = lkk
K

14: end for

15: end for
as

16: for i = 1 to m do
an

17: zi = bi
M

18: for j = 1 to i − 1 do
19: zi = zi − lij zj
r.

20: end for

D

21: zi = zi /lii
22: end for
23: for i = m to 1 do
24: xi = zi
25: for j = i + 1 to m do
26: xi = xi − uij xj
27: end for
28: end for
29: Output: x1 , x2 , ..., xn
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 64

7. Determination of dominant Eigen value & vector by Power method

The eigen values λ1 , λ2 , ..., λn of a square matrix
 
a11 a12 a13 ... a1n−1 a1n
 a21 a 22 a 23 ... a 2n−1 a 2n

 
A=  ... ... ... ... ... ... 

 an−11 an−12 an−13 ... an−1n−1 an−1n

ge

an1 an2 an3 ... ann−1 ann

le
are obtained by solving the characteristic equation

ol
C
an λn + an−1 λn−1 + ... + a1 λ + a0 = 0

aj
Dominant eigen value and vector. Let λ1 , λ2 , ..., λn be the eigenvalues of an n × n

R
matrix A. λ1 is called the dominant eigen value of A if

l
|λ1 | > |λi |, i = 2, 3, ..., n

da
and corresponding eigen vector is called dominant eigen vector.

ha
Theorem. is
If x is an eigenvector of a matrix A, then its corresponding eigenvalue
ah
is given by
M

Ax.x
λ=
x.x
::
ti

This quotient is called the Rayleigh quotient.

ai

Proof. Let x be an eigen vector of A for the eigen value λ then we have
M

Ax.x
Ax = λx ⇒ Ax.x = λx.x ⇒ λ =
x.x
ar

The power method for approximating eigenvalues is a iterative process. First we assume
um

that the matrix A has a dominant eigenvalue with corresponding dominant eigenvectors.
Then we choose an initial approximation x0 of one of the dominant eigenvectors of A.
K

This initial approximation must be a non-zero vector in <n . Finally we form the sequence
as

of approximations given by
an

x1 = Ax0 , x2 = Ax1 = A2 x0 , x3 = Ax2 = A3 x0 , ..., xk = Ak x0

M

For large powers of k, and by properly scaling this sequence, we will see that we obtain a
good approximation of the dominant eigenvector of A.
r.
D

Example. Consider the matrix

 
1 2 0
A =  −2 1 2 
1 3 1
with initial guess of the dominant eigen vector
 
1
x0 = 1 

1
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 65

Then successive approximations of the eigen vectorand eigen value are obtained as
      
1 2 0 1.00 3.00 0.60
Ax0 =  −2 1 2   1.00  =  1.00  = 5.00  0.20  = λ11 x1
1 3 1 1.00 5.00 1.00
      
1 2 0 0.60 1.00 0.45
Ax1 = −2 1 2
   0.20 = 1.00
   = 2.20  0.45  = λ21 x2

ge
1 3 1 1.00 2.20 1.00

le
      
1 2 0 0.45 1.35 0.48

ol
Ax2 = −2 1 2
   0.45 = 1.55
   = 2.80  0.55  = λ31 x3

C
1 3 1 1.00 2.80 1.00

aj
      
1 2 0 0.48 1.58 0.50
 = 3.13  0.51  = λ41 x4

R
Ax3 =  −2 1 2   0.55  = 1.59

1 3 1 1.00 3.13 1.00

l
da
      
1 2 0 0.50 1.52 0.50
 = 3.03  0.50  = λ51 x5

ha
Ax4 =  −2 1 2   0.51  =  1.51
1 3 1 1.00 3.03 1.00
    
is
 
ah
1 2 0 0.50 1.50 0.50
Ax5 =  −2 1 2   0.50  =  1.50  = 3.00  0.50  = λ61 x6
M

1 3 1 1.00 3.00 1.00

::

 
0.50
ti

So dominant eigen value λ = 3.00 and eigen vector= x =  0.50 

ai

1.00
M
ar
um
K
as
an
M
r.
D
 Prepared by: Dr. Manas Kumar Maiti, Mahishadal Raj College. For any error/ambiguity M:7908210880. 66

Algorithm 23 Dominant eigen value by power method

1: Input: n and the matrix A = (aij )n×n .
2: for i = 1 to n do
3: xi = 1
4: end for
5: λ0 = 1
D = 100

ge
6:
7:  = 0.000001

le
8: while D >  do

ol
9: for i = 1 to n do . Set initial guess of the eigen vector

C
10: yi = 0
for j = 1 to n do

aj
11:
12: yi = yi + aij xj

R
13: end for

l
14: end for

da
15: λ = y1

ha
16: for i = 2 to n do . Determine the element with maximum magnitude
17: if |λ| < |yi | then
λ = yi is
ah
18:
19: end if
M

20: end for

::

21: for i = 1 to n do . Set new approximation of eigen vector

ti

22: yi = yi /λ
ai

23: xi = y i
M

24: end for

25: D = |λ − λ0 |
ar

26: λ0 = λ
um

27: end while

28: Output: x1 , x2 , ..., xn
K
as
an
M
r.
D