SHORT NOTES ON NUMERICAL METHODS

(rough draft)
Chapter 1

Finite differences

1.1 Differences of y = f (x)

Let y = f (x) be a function and let 4x = h denote the increment in the independent
variable x (4x is also known as interval or spacing). Suppose that h is a constant.
• The first finite difference of the function y = f (x) is defined by

4f (x) = f (x + 4x) − f (x)

• The second finite difference of the function y = f (x) is defined by

42 f (x) = 4(4f (x))

• The third finite difference of the function y = f (x) is defined by

43 f (x) = 4(42 f (x))

• In general, the nth finite difference of the function y = f (x) is defined by

4n f (x) = 4(4n−1 f (x)) n = 2, 3, 4 · · ·

Example 1.1.1. If h = 1 then find 4n (ex )

Solution. Let f (x) = ex and let h = 4x = 1. We know that

4f (x) = f (x + 4x) − f (x)

Therefore,

4ex = ex+4x − ex
= ex+1 − ex
= ex e − ex

2
 3

Thus,
4ex = (e − 1)ex (1.1.1)

42 ex = 4(4ex )
= 4[(e − 1)ex ]
= (e − 1)4ex
= (e − 1)(e − 1)ex
Thus,
42 ex = (e − 1)2 ex (1.1.2)

43 ex = 4(42 ex )
= 4[(e − 1)2 ex ]
= (e − 1)2 4ex
= (e − 1)2 (e − 1)ex
Thus,
43 ex = (e − 1)3 ex (1.1.3)
Continuing this we get
4n ex = (e − 1)n ex

Theorem 1.1.2. The nth difference of a polynomial of degree n is constant and all higher
order differences are zero.

Proof. Let f (x) = an xn + an−1 xn−1 + · · · + a1 x + a0 be a polynomial of degree n. Consider

4f (x) = f (x + h) − f (x))
= [an (x + h)n + an−1 (x + h)n−1 + · · · + a1 (x + h) + a0 ] − [an xn + an−1 xn−1 + · · · + a1 x + a0 ]
= an [(x + h)n − xn ] + an−1 [(x + h)n−1 − xn−1 ] + · · · + a1 [(x + h) − x]
= an nhxn−1 + a0n−2 xn−2 + · · · + a00

42 f (x) = an n(n − 1)h2 xn−2 + a00n−3 xn−3 + · · · + a000

43 f (x) = an n(n − 1)(n − 2)h3 xn−3 + a000

n−4 x
n−4
+ · · · + a000
0

Continuing this process we get

4n f (x) = an n(n − 1)(n − 2) · · · 3 · 2 · 1hn

= an n!hn

which is a constant and hence the (n + 1)st and higher differences of a polynomial of
degree n are zero.
4

1.2 Other difference operators

Besides the operator 4, we have other difference operators:

• ∇f (x) = f (x) − f (x − h)

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

• E −1 f (x) = f (x − h)

Relation between the operators:

1. 4 = E − 1

2. ∇ = 1 − E −1

3. 4 = E∇ = ∇E

Proof. 1. Consider

4f (x) = f (x + h) − f (x)
= Ef (x) − f (x)
= (E − 1)f (x)

Thus 4 = E − 1

2. Consider

∇f (x) = f (x) − f (x − h)
= f (x) − E −1 f (x)
= (1 − E −1 )f (x)

Thus ∇ = 1 − E −1

3. Consider

E∇f (x) = E(∇f (x))

= E(f (x) − f (x − h))
= Ef (x) − Ef (x − h)
= f (x + h) − f (x)
= 4f (x)

Thus
E∇ = 4 (1.2.1)
 5

Consider

∇Ef (x) = ∇(Ef (x))

= ∇f (x + h)
= f (x + h) − f (x)
= 4f (x)

Thus
∇E = 4 (1.2.2)
From (1.2.1) and (1.2.2) we get 4 = E∇ = ∇E.

1.3 Differences of functions in tabular form

Suppose a function y = f (x) is tabulated, for the equally spaced values of x, as under.

x x0 x1 x2 ··· xn
y y0 y1 y2 ··· yn

The differences y1 − y0 , y2 − y1 , y3 − y2 , · · · , yn − yn−1 , when denoted by 4y0 , 4y1 , 4y2 ,

· · · , 4yn−1 respectively, are called first forward differences where 4(read as “delta”)
is the forward difference operator. Thus, the first forward differences are

4yr = yr+1 − yr (0 ≤ r ≤ n − 1)

The second forward differences are defined by

42 yr = 4yr+1 − 4yr (0 ≤ r ≤ n − 2)

In general the pth (p ≤ n) forward differences are defined by

4p yr = 4p−1 yr+1 − 4p−1 yr (0 ≤ r ≤ n − p)

Example 1.3.1. Write the forward difference table for the following data:

x 0 1 2 3 4
y 1 2 4 8 16

Solution. The required table is as under.

6

x y 4y 42 y 43 y

x0 y0
4y0
x1 y1 42 y0
4y1 43 y0
x2 y2 42 y1
4y2
x3 y3

Table 1.1: Forward difference table

x y 4y 42 y 43 y 44 y

0 1
1
1 2 1
2 1
2 4 2 1
4 2
3 8 4
8
4 16

Example 1.3.2. Construct the forward difference table for the following data:

x 1 2 3 4 5 6 7 8
y 0 5 6 4 10 8 20 25

Solution. The required table is as under.

 7

x y 4y 42 y 43 y 44 y 45 y 46 y 47 y

1 0
5
2 5 -4
1 1
3 6 -3 10
-2 11 -37
4 4 8 -27 102
6 -16 65 -248
5 10 -8 38 -146
-2 22 -81
6 8 14 -43
12 -21
7 20 -7
5
8 25
Example 1.3.3. Find the missing term in the following data:
x 0 1 2 3 4
y 1 2 - 82 257
Solution. Let the missing term be a. The forward difference table is as shown below.
x y 4y 42 y 43 y 44 y

0 1
1
1 2 a−3
a−2 87 − 3a
2 a 84 − 2a 6a − 78
82 − a 9 + 3a
3 82 93 + a
175
4 257
As only four y values are given, the function y can be represented by a polynomial of
degree less than or equal to 3. Therefore, 44 y = 0 i.e., 6a − 78 = 0 which gives a = 13.
Example 1.3.4. Find the missing term in the following data:
x 1 2 3 4 5 6
y 2 6 12 - 30 42
Solution. Let the missing term be a. The forward difference table is as shown below.
8

x y 4y 42 y 43 y 44 y 45 y

1 2
4
2 6 2
6 a − 20
3 12 a − 18 80 − 4a
a − 12 60 − 3a 10a − 200
4 a 42 − 2a 6a − 120
30 − a 3a − 60
5 30 a − 18
12
6 42

As only five y values are given, the function y can be represented by a polynomial of
degree less than or equal to 4. Therefore, 45 y = 0 i.e., 10a − 200 = 0 which gives a = 20.
Suppose a function y = f (x) is tabulated, for the equally spaced values of x, as under.

x x0 x1 x2 ··· xn
y y0 y1 y2 ··· yn

The differences y1 − y0 , y2 − y1 , y3 − y2 , · · · , yn − yn−1 , when denoted by ∇y1 , ∇y2 , ∇y3 ,

· · · , ∇yn respectively, are called first backward differences where ∇(read as “del”) is
the backward difference operator. The second backward differences are defined by

∇2 yr = ∇yr − ∇yr−1

In general the pth backward differences are defined by

∇p yr = ∇p−1 yr − ∇p−1 yr−1

 9

x y ∇y ∇2 y ∇3 y

x0 y0
∇y1
x1 y1 ∇2 y2
∇y2 ∇3 y3
x2 y2 ∇2 y3
∇y3
x3 y3

Table 1.2: Backward difference table

Chapter 2

Interpolation

Given

x 0 2 4 6
y 0 16 128 432

then to find f (1), we first fit a curve to the given data i.e., we find a function y = f (x)
such that f (0) = 0, f (2) = 16,f (4) = 128 and f (6) = 432. The graph of such a function
pass through the points (0, 0), (2, 16), (4, 128) and (6, 432). Then we use y = f (x) to find
f (1)(this process is called “interpolation”). Note that y = 2x3 fits the given data and we
take f (1) = 2.

Note. There could be several curves that fit the given data.

x 0 1
Example 2.0.5. Curves that fit the data are
y 0 1

πx
• y = sin
2
• y = xex−1

• y=x

x 0 1
There are, in fact, infinitely many functions which fit the data . Look at
y 0 1
the function y = xn (n > 1). But there is a unique polynomial of degree less than two
which fits the above data. In this example, it is y = x(polynomial of degree 1).

Theorem 2.0.6. There exists a unique polynomial of degree less than n + 1 which fits the
x x0 x1 x2 · · · xn
data
y y0 y1 y2 · · · yn

10
 11

Figure 2.1: Three different curves passing through two points

12

The unique polynomial in the above theorem, is known as the interpolating polyno-
mial and the process of finding this polynomial is known as polynomial interpolation.
Question: How to find the interpolating polynomial for the given data?
If in the given data, x values are equally spaced then Newton’s forward difference interpo-
lation formula (NFDIF) or Newton’s backward difference interpolation formula (NBDIF)
can be used to find the interpolating polynomial.
 13

2.1 Newton’s interpolation formulae

Suppose the values of y = f (x) be given for a set of equally spaced values of x:

x x0 x1 x2 ··· xn
y y0 y1 y2 ··· yn

Let h be the spacing between x values. Then

xi = x0 + ih (0 ≤ i ≤ n) (2.1.1)

Let

F (x) = a0 + a1 (x − x0 ) + a2 (x − x0 )(x − x1 )
+ a3 (x − x0 )(x − x1 )(x − x2 ) + · · · (2.1.2)
+ an (x − x0 )(x − x1 ) · · · (x − xn−1 )

be the interpolating polynomial that fits given data. Note that yi = F (xi ) (0 ≤ i ≤ n).
Putting x = x0 in the equation (2.1.2), we get

F (x0 ) = a0

which implies
a0 = y0 (as y0 = F (x0 )) (2.1.3)
Putting x = x1 in the equation (2.1.2), we get

F (x1 ) = a0 + a1 (x1 − x0 )
= a0 + a1 h (as x1 − x0 = h)

which implies

F (x1 ) − a0
a1 =
h
y1 − y0
= (as F (x1 ) = y1 and a0 = y0 )
h
Thus,
1
a1 = 4y0 (as 4y0 = y1 − y0 ) (2.1.4)
h
Putting x = x2 in the equation (2.1.2), we get

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

= a0 + a1 (2h) + a2 (2h)h (as x2 − x0 = 2h and x2 − x1 = h)
14

which implies

F (x2 ) − a0 − 2ha1
a2 =
2h2
y2 − y0 − 2h h1 4y0 1
= 2
(as F (x2 ) = y2 , a0 = y0 and a1 = 4y0 )
2h h
y2 − y0 − 2(y1 − y0 )
= (as 4y0 = y1 − y0 )
2h2
y2 − y0 − 2y1 + 2y0
=
2h2
y2 − 2y1 + y0
=
2h2
(y2 − y1 ) − (y1 − y0 )
=
2h2
4y1 − 4y0
= (as 4yr = yr+1 − yr )
2h2

Thus,
1 2
a2 = 4 y0 (as 42 y0 = 4y1 − 4y0 ) (2.1.5)
2h2
Putting x = x3 in the equation (2.1.2), we get

F (x3 ) = a0 + a1 (x3 − x0 ) + a2 (x3 − x0 )(x3 − x1 )

+ a3 (x3 − x0 )(x3 − x1 )(x3 − x2 )
= a0 + a1 (3h) + a2 (3h)(2h) + a3 (3h)(2h)h
(as x3 − x0 = 3h, x3 − x1 = 2h and x3 − x2 = h)

which implies

F (x3 ) − a0 − 3ha1 − 6h2 a2

a3 =
6h3
y3 − y0 − 3h h1 4y0 − 6h2 2h1 2 42 y0
=
 6h3 
1 1 2
as F (x3 ) = y3 , a0 = y0 , a1 = 4y0 and a2 = 2 4 y0
h 2h
y3 − y0 − 3(y1 − y0 ) − 3(y2 − 2y1 + y0 )
=
6h3
(as 4y0 = y1 − y0 and 42 y0 = y2 − 2y1 + y0 )
y3 − y0 − 3y1 + 3y0 − 3y2 + 6y1 − 3y0
=
6h3
 15

Therefore,
y3 − 3y2 + 6y1 − y0
a3 =
6h3
(y3 − y2 ) − (y2 − y1 ) − [(y2 − y1 ) − (y1 − y0 )]
=
6h3
4y2 − 4y1 − [4y1 − 4y0 ]
= (as 4yr = yr+1 − yr )
6h3
42 y1 − 42 y0
= 3
(as 42 yr = 4yr+1 − 4yr )
6h
Thus,
1
a3 = 43 y0 (as 43 y0 = 42 y1 − 42 y0 ) (2.1.6)
3!h3
Continuing this, we get
1
a4 = 44 y0 (2.1.7)
4!h4
1
a5 = 5
45 y0 (2.1.8)
5!h
1
an = 4n y0 (2.1.9)
n!hn
1 1
F (x) = y0 + 4y0 (x − x0 ) + 42 y0 (x − x0 )(x − x1 )
h 2!h2
1
+ 43 y0 (x − x0 )(x − x1 )(x − x2 )
3!h3
1
+ ··· + 4n y0 (x − x0 )(x − x1 )(x − x2 ) · · · (x − xn−1 )
n!hn
x − x0 42 y0 x − x0 x − x1
= y0 + 4y0 +
h 2! h h
3
4 y0 x − x0 x − x1 x − x2
+
3! h h h
4n y0 x − x0 x − x1 x − x2 x − xn−1
+ ··· + ···
n! h h h h
x − x0 x − xr x − [x0 + rh] [x − x0 ] − rh x − x0
Put s = . Note that = = = −r = s − r.
h h h h h
Therefore,
42 y0
F (x) = y0 + 4y0 s + s(s − 1)
2!
43 y0
+ s(s − 1)(s − 2)
3!
4n y0
+ ··· + s(s − 1)(s − 2) · · · (s − n + 1)
n!
Xn
s
= Cr 4r y0
r=0
16

s(s − 1)(s − 2) · · · (s − r + 1)
where s Cr = and s C0 = 1.
r!
Newton’s interpolation formulae - rewritten

• (NFDIF):
n
X
s
F (x) = Cr 4r y0
r=0

x − x0 s s(s − 1)(s − 2) · · · (s − r + 1)
where s = , Cr = and s C0 = 1.
h r!
• (NBDIF):

s(s + 1) 2
F (x) = yn + s∇yn + ∇ yn
2!
s(s + 1)(s + 2) · · · (s + n − 1) n
+ ··· + ∇ yn
n!
x − xn
where s = .
h
Example 2.1.1. Find the cubic polynomial which takes the following values:

x 0 1 2 3
y 1 4 13 34

and hence find y(0.5).

Solution. The forward difference table is as under.

x y 4y 42 y 43 y

0 1
3
1 4 6
9 6
2 13 12
21
3 34

The interpolating polynomial is given by

n
X
s
F (x) = Cr 4r y0
r=0
 17

x − x0
Here n = 3, x0 = 0, h = 1. Thus, s = = x. Further y0 = 1, 4y0 = 3, 42 y0 = 6
h
and 43 y0 = 6. Thus, we have

3
X
s
F (x) = Cr 4r y0
r=0
= y0 +s C1 4y0 +s C2 42 y0 +s C3 43 y0
s(s − 1) 2 s(s − 1)(s − 2) 3
= y0 + s4y0 + 4 y0 + 4 y0
2! 3!
x(x − 1) 2 x(x − 1)(x − 2) 3
= y0 + x4y0 + 4 y0 + 4 y0
2! 3!
6 6
= 1 + 3x + (x2 − x) + (x3 − 3x2 + 2x)
2 6
= 1 + 3x + 3(x2 − x) + x3 − 3x2 + 2x
= x3 + 2x + 1

Hence the value of y at x = 0.5 is F (0.5) = (0.5)3 + 2(0.5) + 1 = 2.125.

Example 2.1.2. Find the interpolating polynomial that fits the following data:

x 0 1 2 3
y 1 2 9 28

and hence find y(0.1).

Solution. The forward difference table is as under.

x y 4y 42 y 43 y

0 1
1
1 2 6
7 6
2 9 12
19
3 28

The interpolating polynomial is given by

n
X
s
F (x) = Cr 4r y0
r=0
18

x − x0
Here n = 3, x0 = 0, h = 1. Thus, s = = x. Further y0 = 1, 4y0 = 1, 42 y0 = 6
h
and 43 y0 = 6. Thus, we have

3
X
s
F (x) = Cr 4r y0
r=0
= y0 +s C1 4y0 +s C2 42 y0 +s C3 43 y0
s(s − 1) 2 s(s − 1)(s − 2) 3
= y0 + s4y0 + 4 y0 + 4 y0
2! 3!
x(x − 1) 2 x(x − 1)(x − 2) 3
= y0 + x4y0 + 4 y0 + 4 y0
2! 3!
6 6
= 1 + x + (x2 − x) + (x3 − 3x2 + 2x)
2 6
= 1 + x + 3(x2 − x) + x3 − 3x2 + 2x
= x3 + 1

Hence the value of y at x = 0.1 is F (0.1) = (0.1)3 + 1 = 1.001

Example 2.1.3. Estimate the population in 1895 and 1925 from the following data:

Years(x) 1891 1901 1911 1921 1931

Population(y) 46 66 81 93 101

Solution. The difference table is as under.

x y 4y 42 y 43 y 44 y

1891 46
20
1901 66 -5
15 2
1911 81 -3 -3
12 -1
1921 93 -4
8
1931 101

First we find y(1895). According to NFDIF

n
X
s
F (x) = Cr 4r y0
r=0
 19

x − x0 1895 − 1891
Here n = 4, x0 = 1891, h = 10 and x = 1895. Thus, s = = = 0.4.
h 10
Further y0 = 46, 4y0 = 20, 42 y0 = −5 43 y0 = 2 and 44 y0 = −3. Thus, we have

4
X
s
F (1895) = Cr 4r y0
r=0
= y0 +s C1 4y0 +s C2 42 y0 +s C3 43 y0 +s C4 44 y0
s(s − 1) 2 s(s − 1)(s − 2) 3
= y0 + s4y0 + 4 y0 + 4 y0
2! 3!
s(s − 1)(s − 2)(s − 3) 4
+ 4 y0
4!
0.4(0.4 − 1) 2
= y0 + (0.4)4y0 + 4 y0
2!
0.4(0.4 − 1)(0.4 − 2) 3
+ 4 y0
3!
0.4(0.4 − 1)(0.4 − 2)(0.4 − 3) 4
+ 4 y0
4!

Therefore,

0.4(0.4 − 1) 0.4(0.4 − 1)(0.4 − 2)

F (1895) = 46 + (0.4)20 + (−5) + 2
2! 3!
0.4(0.4 − 1)(0.4 − 2)(0.4 − 3)
+ (−3)
4!
= 46 + 8 + 0.6 + 0.128 + 0.1248
= 54.8528

Now we find y(1925). According to NBDIF

s(s + 1) 2
F (x) = yn + s∇yn + ∇ yn
2!
s(s + 1)(s + 2) · · · (s + n − 1) n
+ ··· + ∇ yn
n!

x − xn 1925 − 1931
Here n = 4, xn = x4 = 1931, h = 10 and x = 1925. Thus, s = = =
h 10
2 2 3 3
−0.6. Further yn = y4 = 101, ∇yn = ∇y4 = 8, ∇ yn = ∇ y4 = −4 ∇ yn = ∇ y4 = −1
20

and ∇4 yn = ∇4 y4 = −3. Thus, we have

s(s + 1) 2
F (1925) = y4 + s∇y4 + ∇ y4
2!
s(s + 1)(s + 2) 3 s(s + 1)(s + 2)(s + 3) 4
+ ∇ y4 + ∇ y4
3! 4!
(−0.6)(−0.6 + 1) 2
= y4 + (−0.6)∇y4 + ∇ y4
2!
(−0.6)(−0.6 + 1)(−0.6 + 2) 3
+ ∇ y4
3!
(−0.6)(−0.6 + 1)(−0.6 + 2)(−0.6 + 3) 4
+ ∇ y4
4!
(−0.6)(0.4) 2
= y4 + (−0.6)∇y4 + ∇ y4
2!
(−0.6)(0.4)(1.4) 3
+ ∇ y4
3!
(−0.6)(0.4)(1.4)(2.4) 4
+ ∇ y4
4!

Therefore,

(−0.6)(0.4)
F (1925) = 101 + (−0.6)8 + (−4)
2!
(−0.6)(0.4)(1.4)
+ (−1)
3!
(−0.6)(0.4)(1.4)(2.4)
+ (−3)
4!
= 101 − 4.8 + 0.48 + 0.056 + 0.1008
= 96.8368

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

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

No. of students 31 42 51 35 31

Solution. We first prepare the following (cumulative frequency) table:

Marks less than (x) 40 50 60 70 80

No. of students(y) 31 73 124 159 190

The difference table is as under:

 21

x y 4y 42 y 43 y 44 y
40 31
42
50 73 9
51 -25
60 124 -16 37
35 12
70 159 -4
31
80 190
Now we find the number of students with marks less than 45. According to NFDIF
Xn
s
F (x) = Cr 4r y0
r=0
x − x0 45 − 40
Here n = 4, x0 = 40, h = 10 and x = 45. Thus, s = = = 0.5. Further
h 10
2 3 4
y0 = 31, 4y0 = 42, 4 y0 = 9 4 y0 = −25 and 4 y0 = 37. Thus, we have
4
X
s
F (45) = Cr 4r y0
r=0
= y0 +s C1 4y0 +s C2 42 y0 +s C3 43 y0 +s C4 44 y0
s(s − 1) 2 s(s − 1)(s − 2) 3
= y0 + s4y0 + 4 y0 + 4 y0
2! 3!
s(s − 1)(s − 2)(s − 3) 4
+ 4 y0
4!
0.5(0.5 − 1) 2
= y0 + (0.5)4y0 + 4 y0
2!
0.5(0.5 − 1)(0.5 − 2) 3
+ 4 y0
3!
0.5(0.5 − 1)(0.5 − 2)(0.5 − 3) 4
+ 4 y0
4!
Therefore,
0.5(0.5 − 1) 0.5(0.5 − 1)(0.5 − 2)
F (45) = 31 + (0.5)42 + (9) + (−25)
2! 3!
0.5(0.5 − 1)(0.5 − 2)(0.5 − 3)
+ (37)
4!
= 31 + 21 − 1.125 − 1.5625 − 1.4453125
= 47.8671875
Therefore, the number of students with marks less than 45 is 48. Note that the number
of students with marks less than 40 is 31. Thus the number of students getting marks
between 40 and 45 is 48 − 31 = 17.
22

2.2 Langrange’s interpolation formula

Suppose the values of y = f (x) be given for a set of (not necessarily equally spaced) values
of x:

x x0 x1 x2 ··· xn
y y0 y1 y2 ··· yn

The interpolating polynomial that fits the above data is given by

(x − x1 )(x − x2 ) · · · (x − xn )
F (x) = y0
(x0 − x1 )(x0 − x2 ) · · · (x0 − xn )
(x − x0 )(x − x2 ) · · · (x − xn )
+ y1 + · · ·
(x1 − x0 )(x1 − x2 ) · · · (x1 − xn )
(x − x0 )(x − x1 ) · · · (x − xi−1 )(x − xi+1 ) · · · (x − xn )
+ yi + · · ·
(xi − x0 )(xi − x1 ) · · · (xi − xi−1 )(xi − xi+1 ) · · · (xi − xn )
(x − x0 )(x − x1 ) · · · (x − xn−1 )
+ yn
(xn − x0 )(xn − x1 ) · · · (xn − xn−1 )

This is known as Langrange’s interpolation formula.

Example 2.2.1. Using Langrange’s interpolation formula, find f (4) given that f (0) = 6,
f (1) = 8,f (3) = 16 and f (6) = 34

Solution. Given

x0 = 0 x1 = 1 x2 = 3 x3 = 6
y0 = 6 y1 = 8 y2 = 16 y3 = 34

If F (x) is the interpolating polynomial that fits the given data, then, by Langrange’s
interpolation formula, we have

(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 − 3)(x − 6) (x − 0)(x − 3)(x − 6)
= 6+ 8
(0 − 1)(0 − 3)(0 − 6) (1 − 0)(1 − 3)(1 − 6)
(x − 0)(x − 1)(x − 6) (x − 0)(x − 1)(x − 3)
+ 16 + 34
(3 − 0)(3 − 1)(3 − 6) (6 − 0)(6 − 1)(6 − 3)
 23

Therefore,

(4 − 1)(4 − 3)(4 − 6) (4 − 0)(4 − 3)(4 − 6)

F (4) = 6+ 8
(0 − 1)(0 − 3)(0 − 6) (1 − 0)(1 − 3)(1 − 6)
(4 − 0)(4 − 1)(4 − 6) (4 − 0)(4 − 1)(4 − 3)
+ 16 + 34
(3 − 0)(3 − 1)(3 − 6) (6 − 0)(6 − 1)(6 − 3)
3(1)(−2) 4(1)(−2)
= 6+ 8
(−1)(−3)(−6) (−2)(−5)
4(3)(−2) 4(3)(−1)
+ 16 + 34
3(2)(−3) 6(5)(3)
= 2 − 6.4 + 21.33 + 4.533
= 21.463

2.3 Inverse interpolation

In (direct)interpolation we find the value of y for a given value of x. In the inverse
interpolation, we find the value of x for a given value of y. Finding a root of of y = f (x)
is an inverse interpolation problem. Interchanging the roles of x and y in the Lagrange’s
interpolation formula, we get a formula for inverse interpolation and is given by

(y − y1 )(y − y2 ) · · · (y − yn )
F (y) = x0
(y0 − y1 )(y0 − y2 ) · · · (y0 − yn )
(y − y0 )(y − y2 ) · · · (y − yn )
+ x1 + · · ·
(y1 − y0 )(y1 − y2 ) · · · (y1 − yn )
(y − y0 )(y − y1 ) · · · (y − yi−1 )(y − yi+1 ) · · · (y − yn )
+ xi + · · ·
(yi − y0 )(yi − y1 ) · · · (yi − yi−1 )(yi − yi+1 ) · · · (yi − yn )
(y − y0 )(y − y1 ) · · · (y − yn−1 )
+ xn
(yn − y0 )(yn − y1 ) · · · (yn − yn−1 )

Example 2.3.1. Apply Langrange’s method to find the value of x when f (x) = 15 from
the following data

x 5 6 9 11
y 12 13 14 16

Solution. Given

x0 = 5 x1 = 6 x2 = 9 x3 = 11
y0 = 12 y1 = 13 y2 = 14 y3 = 16
24

Now, by Langrange’s formula for inverse interpolation, we have

(y − y1 )(y − y2 )(y − y3 ) (y − y0 )(y − y2 )(y − y3 )

F (y) = x0 + x1
(y0 − y1 )(y0 − y2 )(y0 − y3 ) (y1 − y0 )(y1 − y2 )(y1 − y3 )
(y − y0 )(y − y1 )(y − y3 ) (y − y0 )(y − y1 )(y − y2 )
+ x2 + x3
(y2 − y0 )(y2 − y1 )(y2 − y3 ) (y3 − y0 )(y3 − y1 )(y3 − y2 )
(y − 13)(y − 14)(y − 16) (y − 12)(y − 14)(y − 16)
= 5+ 6
(12 − 13)(12 − 14)(12 − 16) (13 − 12)(13 − 14)(13 − 16)
(y − 12)(y − 13)(y − 16)
+ 9
(14 − 12)(14 − 13)(14 − 16)
(y − 12)(y − 13)(y − 14)
+ 11
(16 − 12)(16 − 13)(16 − 14)

(15 − 13)(15 − 14)(15 − 16) (15 − 12)(15 − 14)(15 − 16)

F (15) = 5+ 6
(12 − 13)(12 − 14)(12 − 16) (13 − 12)(13 − 14)(13 − 16)
(15 − 12)(15 − 13)(15 − 16)
+ 9
(14 − 12)(14 − 13)(14 − 16)
(15 − 12)(15 − 13)(15 − 14)
+ 11
(16 − 12)(16 − 13)(16 − 14)
(2)(1)(−1) (3)(1)(−1)
= 5+ 6
(−1)(−2)(−4) (1)(−1)(−3)
(3)(2)(−1) (3)(2)(1)
+ 9+ 11
(2)(1)(−2) (4)(3)(2)
5 27 11
= −6+ +
4 2 4
= 11.5

Example 2.3.2. Apply Langrange’s formula inversely to obtain a root of the equation
f (x) = 0, given that f (30) = −30, f (34) = −13, f (38) = 3 and f (42) = 18

Solution. Given

x0 = 30 x1 = 34 x2 = 38 x3 = 42
y0 = −30 y1 = −13 y2 = 3 y3 = 18
 25

Now, by Langrange’s formula for inverse interpolation, we have

(y − y1 )(y − y2 )(y − y3 ) (y − y0 )(y − y2 )(y − y3 )

F (y) = x0 + x1
(y0 − y1 )(y0 − y2 )(y0 − y3 ) (y1 − y0 )(y1 − y2 )(y1 − y3 )
(y − y0 )(y − y1 )(y − y3 ) (y − y0 )(y − y1 )(y − y2 )
+ x2 + x3
(y2 − y0 )(y2 − y1 )(y2 − y3 ) (y3 − y0 )(y3 − y1 )(y3 − y2 )
(y − (−13))(y − 3)(y − 18)
= 30
(−30 − (−13))(−30 − 3)(−30 − 18)
(y − (−30))(y − 3)(y − 18)
+ 34
(−13 − (−30))(−13 − 3)(−13 − 18)
(y − (−30))(y − (−13))(y − 18)
+ 38
(3 − (−30))(3 − (−13))(3 − 18)
(y − (−30))(y − (−13))(y − 3)
+ 42
(18 − (−30))(18 − (−13))(18 − 3)

(0 − (−13))(0 − 3)(0 − 18)

F (0) = 30
(−30 − (−13))(−30 − 3)(−30 − 18)
(0 − (−30))(0 − 3)(0 − 18)
+ 34
(−13 − (−30))(−13 − 3)(−13 − 18)
(0 − (−30))(0 − (−13))(0 − 18)
+ 38
(3 − (−30))(3 − (−13))(3 − 18)
(0 − (−30))(0 − (−13))(0 − 3)
+ 42
(18 − (−30))(18 − (−13))(18 − 3)
(13)(−3)(−18) (30)(−3)(−18)
= 30 + 34
(−17)(−33)(−48) (17)(−16)(−31)
(30)(13)(−18) (30)(13)(−3)
+ 38 + 42
(33)(16)(−15) (48)(31)(15)
= −0.78209 + 6.5323 + 33.6818 − 2.2016
= 37.23
Chapter 3

Numerical differentiation and

integration

Not all functions (of a single variable) can be differentiated/integrated with the rules
(or methods or techniques) of differentiation/integration that we have studied so far. To
find an approximate value of the derivatives/integral of such functions, we use numerical
methods.

3.1 Numerical differentiation

If the function y = f (x) is highly complex or given in tabular form then we can’t find its
derivatives with the rules of differentiation studied so far. In such cases we use numerical
differentiation. The idea here is very simple; we just replace the given function y = f (x)
by the interpolating polynomial F (x) and take f 0 (x) = F 0 (x), f 00 (x) = F 00 (x) etc.

3.1.1 Numerical differentiation using NFDIF

Suppose we are given y = f (x) in the tabular form

x x0 x1 x2 ··· xn
y y0 y1 y2 ··· yn

where x values are equally spaced with spacing h. Then xi = x0 + ih (0 ≤ i ≤ n). Now
replacing the given function by the interpolating polynomial (got by using NFDIF) we

26
 27

have
n
X
s x − x0
y= Cr 4r y0 (where s = )
r=0
h
= y0 + s C1 4y0 + s C2 42 y0 + s C3 43 y0 + s C4 44 y0 + · · ·
42 y0 43 y0 44 y0
= y0 + 4y0 s + s(s − 1) + s(s − 1)(s − 2) + s(s − 1)(s − 2)(s − 3) + · · ·
2! 3! 4!
42 y0 2 43 y0 3 44 y0 4
= y0 + 4y0 s + (s − s) + (s − 3s2 + 2s) + (s − 6s3 + 11s2 − 6s) + · · ·
2 6 24

Differentiating y w.r.t. s we get,

42 y0 2 43 y0 3 44 y0 4
 
dy d 2 3 2
= y0 + 4y0 s + (s − s) + (s − 3s + 2s) + (s − 6s + 11s − 6s) + · · ·
ds ds 2 6 24
d 42 y0 2 d 43 y0 3
   
d d 2
= (y0 ) + (4y0 s) + (s − s) + (s − 3s + 2s)
ds ds ds 2 ds 6
d 44 y0 4
 
3 2
+ (s − 6s + 11s − 6s) + · · ·
ds 24
ds 42 y0 d 2 43 y0 d 3
= 0 + 4y0 + (s − s) + (s − 3s2 + 2s)
ds 2 ds 6 ds
44 y0 d 4
+ (s − 6s3 + 11s2 − 6s) + · · ·
24 ds

Therefore,

42 y0 43 y0 2 44 y0 3
 
dy 2
= 4y0 + (2s − 1) + (3s − 6s + 2) + (4s − 18s + 22s − 6) + · · ·
ds 2 6 24
(3.1.1)
Now differentiating y w.r.t. x we get,

dy dy ds
= (by chain rule)
dx ds dx  
1 dy ds 1
= because =
h ds dx h

Therefore, by Equation (3.1.1), we have

42 y0 43 y0 2 44 y0 3
 
dy 1 2
= 4y0 + (2s − 1) + (3s − 6s + 2) + (4s − 18s + 22s − 6) + · · ·
dx h 2 6 24
(3.1.2)
28

Now differentiating Equation (3.1.2) w.r.t. s, we get

42 y0 43 y0 2 44 y0 3
    
d dy d 1 2
= 4y0 + (2s − 1) + (3s − 6s + 2) + (4s − 18s + 22s − 6) + · · ·
ds dx ds h 2 6 24
d 42 y0 d 43 y0 2
    
1 d
= (4y0 ) + (2s − 1) + (3s − 6s + 2)
h ds ds 2 ds 6
d 44 y0 3
  
2
+ (4s − 18s + 22s − 6) + · · ·
ds 24
42 y0 d 43 y0 d 44 y0 d
 
1 2 3 2
= 0+ (2s − 1) + (3s − 6s + 2) + (4s − 18s + 22s − 6) + · · ·
h 2 ds 6 ds 24 ds
1 42 y0 43 y0 44 y0
 
2
= (2) + (6s − 6) + (12s − 36s + 22) + · · ·
h 2 6 24
43 y0 44 y0
 
1 2 2
= 4 y0 + 6(s − 1) + 2(6s − 18s + 11) + · · ·
h 6 24
Therefore,
44 y0 2
   
d dy 1 2 3
= 4 y0 + 4 y0 (s − 1) + (6s − 18s + 11) + · · · (3.1.3)
ds dx h 12
Differentiating Equation (3.1.2) w.r.t. x, we get
d2 y
 
d dy
=
dx2 dx dx
 
d dy ds
= (by chain rule)
ds dx dx
   
1 d dy ds 1
= because =
h ds dx dx h
Therefore, by Equation (3.1.3), we have
d2 y 44 y0 2
 
11 2 3
= 4 y0 + 4 y0 (s − 1) + (6s − 18s + 11) + · · ·
dx2 hh 12
(3.1.4)
44 y0 2
 
1 2 3
= 2 4 y0 + 4 y0 (s − 1) + (6s − 18s + 11) + · · ·
h 12
By putting s = 0 in equations (3.1.2) and (3.1.4), we get the first and second derivatives
of y = f (x) at x = x0 and are given by
 
dy 1 1 2 1 3 1 4
= 4y0 − 4 y0 + 4 y0 − 4 y0 + · · ·
dx x=x0 h 2 3 4

d2 y
 
1 2 3 11 4
= 2 4 y0 − 4 y0 + 4 y0 + · · ·
dx2 x=x0 h 12
 29

Example 3.1.1. Compute y 0 (1.5) and y 00 (1.5) given that

x 1 2 3 4 5
y 3 21 91 273 651

Solution. The difference table is as under.

x y 4y 42 y 43 y 44 y

1 3
18
2 21 52
70 60
3 91 112 24
182 84
4 273 196
378
5 651
x − x0 1.5 − 1
Here x0 = 1, x = 1.5, h = 1 and s = = = 0.5.
h 1
We know that
42 y0 43 y0 2 44 y0 3
 
dy 1 2
= 4y0 + (2s − 1) + (3s − 6s + 2) + (4s − 18s + 22s − 6) + · · ·
dx h 2 6 24

Therefore,
 
0 1 52 60 2 24 3 2
y (1.5) = 18 + (2(0.5) − 1) + (3(0.5) − 6(0.5) + 2) + (4(0.5) − 18(0.5) + 22(0.5) − 6)
1 2 6 24
= [18 + 26(0) + 10(−0.25) + 1]
= 16.5

we know that
d2 y 44 y0 2
 
1 2 3
= 2 4 y0 + 4 y0 (s − 1) + (6s − 18s + 11) + · · ·
dx2 h 12

Therefore,
 
00 1 24 2
y (1.5) = 2 52 + 60(0.5 − 1) + (6(0.5) − 18(0.5) + 11)
1 12
= [52 − 30 + 2(3.5)]
= 29
30

Example 3.1.2. Compute y 0 (1) and y 00 (1) given that

x 1 1.2 1.4 1.6 1.8
y 100 158.56 265.76 441.76 710.56
Solution. The difference table is as under.
x y 4y 42 y 43 y 44 y

1 100
58.56
1.2 158.56 48.64
107.2 20.16
1.4 265.76 68.8 3.84
176 24
1.6 441.76 92.8
268.8
1.8 710.56
We know that
 
dy 1 1 2 1 3 1 4
= 4y0 − 4 y0 + 4 y0 − 4 y0 + · · ·
dx x=x0 h 2 3 4
Therefore,
 
0 1 1 1 1
y (1) = 58.56 − (48.64) + (20.16) − (3.84)
0.2 2 3 4
1
= [58.56 − 24.32 + 6.72 − 0.96]
0.2
1
= (40)
0.2
= 200
We know that
d2 y
 
1 2 3 11 4
= 2 4 y0 − 4 y0 + 4 y0 + · · ·
dx2 x=x0 h 12
Therefore,
 
00 1 11
y (1) = 48.64 − 20.16 + (3.84)
(0.2)2 12
1
= [48.64 − 20.16 + 3.52]
0.04
1
= (32)
0.04
= 800
 31

3.1.2 Numerical differentiation using NBDIF

Using Newton’s backward difference interpolation formula (NBDIF) to find the interpo-
lating polynomial and following the steps described in subsection 3.1.1 we get

∇2 yn ∇3 yn 2 ∇4 yn 3
 
dy 1 2
= ∇yn + (2s + 1) + (3s + 6s + 2) + (4s + 18s + 22s + 6) + · · ·
dx h 2 6 24
(3.1.5)
2
 
dy 1 1
2
= 2 ∇2 yn + (s + 1)∇3 yn + (6s2 + 18s + 11)∇4 yn + · · · (3.1.6)
dx h 12
By putting s = 0 in Equations (3.1.5) and (3.1.6) we get the first and second derivatives
of y = f (x) at the end point x = xn and are given by
 
dy 1 1 2 1 3 1 4
= ∇yn + ∇ yn + ∇ yn + ∇ yn + · · ·
dx x=xn h 2 3 4

d2 y
 
1 2 3 11 4
= 2 ∇ yn + ∇ yn + ∇ yn + · · ·
dx2 x=xn h 12

Example 3.1.3. Compute y 0 (4.5) and y 00 (4.5) given that

x 1 2 3 4 5
y 3 21 91 273 651

Solution. The difference table is as under.

x y 4y 42 y 43 y 44 y

1 3
18
2 21 52
70 60
3 91 112 24
182 84
4 273 196
378
5 651

x − xn 4.5 − 5
Here xn = 5, x = 4.5, h = 1 and s = = = −0.5.
h 1
We know that
 
dy 1 1 2 1 2 3 1 3 2 4
= ∇yn + (2s + 1)∇ yn + (3s + 6s + 2)∇ yn + (4s + 18s + 22s + 6)∇ yn + · · ·
dx h 2 6 24
32

Therefore,

0 1 1 1
y (4.5) = 378 + (2(−0.5) + 1)196 + (3(−0.5)2 + 6(−0.5) + 2)84
1 2 6

1 3 2
+ (4(−0.5) + 18(−0.5) + 22(−0.5) + 6)24 + · · ·
24
1
= 378 + (−0.25)84 − 1
6
= 373.5

We know that
d2 y
 
1 2 3 1 2 4
= 2 ∇ yn + (s + 1)∇ yn + (6s + 18s + 11)∇ yn + · · ·
dx2 h 12

Therefore,
 
00 1 1 2
y (4.5) = 2 196 + (−0.5 + 1)84 + (6(−0.5) + 18(−0.5) + 11)24
1 12
= 196 + 42 + 7
= 245

Example 3.1.4. Compute y 0 (5) and y 00 (5) given that

x 1 2 3 4 5
y 0.121 1.936 9.801 30.976 75.625

Solution. The difference table is as under.

x y ∇y ∇2 y ∇3 y ∇4 y

1 0.121
1.815
2 1.936 6.05
7.865 7.26
3 9.801 13.31 2.904
21.175 10.164
4 30.976 23.474
44.649
5 75.625

We know that
 
dy 1 1 2 1 3 1 4
= ∇yn + ∇ yn + ∇ yn + ∇ yn + · · ·
dx x=xn h 2 3 4
 33

Therefore,
 
0 1 1 1 1
y (5) = 44.649 + (23.474) + (10.164) + (2.904)
1 2 3 4
= 44.649 + 11.737 + 3.388 + 0.726
= 60.5

We know that
d2 y
 
1 2 3 11 4
= 2 ∇ yn + ∇ yn + ∇ yn + · · ·
dx2 x=xn h 12

Therefore,
 
00 1 11
y (5) = 23.474 + 10.164 + (2.904)
(1)2 12
= 23.474 + 10.164 + 2.662
= 36.3

3.2 Numerical integration

When the function is highly complex or given in the tabular form then we use numerical
integration to integrate the given function. The basic idea in numerical integration is to
replace the given function y = f (x) in tabular form by the interpolating polynomial F (x)
Rb Rb
and take f (x) dx = F (x) dx.
a a

3.2.1 Newton-Cotes’ quadrature formula

Let
Zb
I= f (x) dx
a

Let us divide the interval [a, b] into n sub-intervals of width h so that x0 = a, x1 =

x0 + h, x2 = x0 + 2h, · · · , xr = x0 + rh, · · · , xn = x0 + nh = b and let y = f (x) takes
the values y0 , y1 , · · · , yn for x = x0 , x1 , · · · , xn respectively. At this stage, we have the
function y = f (x) in the tabular form

x x0 x1 x2 ··· xn
y y0 y1 y2 ··· yn
34

Now replacing the given function by the interpolating polynomial(got by NFDIF) we have

Zb
I= f (x) dx
a
Zb 
1 1
= y0 + 4y0 (x − x0 ) + 42 y0 (x − x0 )(x − x1 )
h 2!h2
a

1 3
+ 3 4 y0 (x − x0 )(x − x1 )(x − x2 ) + · · · dx
3!h

x − x0 1
Put s = . This implies ds = dx i.e., dx = h ds
h h
When x = a (i.e., x = x0 ), s = 0
x0 + nh − x0
and when x = b (i.e., x = x0 + nh), s = = n. Therefore,
h
Zn 
42 y0 43 y0

I= y0 + 4y0 s + s(s − 1) + s(s − 1)(s − 2) + · · · h ds
2! 3!
0
Zn 
42 y0 2 43 y0 3

2
=h y0 + 4y0 s + (s − s) + (s − 3s + 2s) + · · · ds
2! 3!
0
n
s2 42 y0 s3 s2 43 y0 s4
    
3 2
= h y0 s + 4y0 + − + − s + s + ···
2 2! 3 2 3! 4
0
n2 42 y0 n3 n2 43 y0 n4
    
= h y0 n + 4y0 + − + − n3 + n2 + · · ·
2 2! 3 2 3! 4

. This is known as the Newton-Cotes’ quadrature formula.

Trapezoidal rule: Putting n = 1 in the Newton-Cotes’ quadrature formula, we get
 
1
I = h y0 + 4y0
2

i.e.,

Zb  
1
f (x) dx = h y0 + (y1 − y0 )
2
a
h
= (y0 + y1 )
2
i.e.,
 35

Rb
f (x) dx =(1/2)(length of the interval)[y at left end point+y at right end point]
a

This is known as the trapezoidal rule.

Note. If we do not divide the interval [a, b](i.e., when n = 1) then we get a very crude
Rb
value of the integral I = f (x) dx. To get more accurate value of I, we divide the interval
a
into n sub-intervals of width h, use trapezoidal rule to evaluate the integral of f (x) over
each of these sub-intervals and finally add all these values.

Dividing the interval [a, b] into n sub-intervals of width h and using trapezoidal rule
we get,
x
Zr+1
h
f (x) dx = (yr + yr+1 ) (0 ≤ r ≤ n − 1)
2
xr

Therefore,

Zb Zxn
f (x) dx = f (x) dx
a x0
Zx1 Zx2 Zx3 xZn−1 Zxn
= f (x) dx + f (x) dx + f (x) dx + · · · + f (x) dx + f (x) dx
x0 x1 x2 xn−2 xn−1
h h h h h
= (y0 + y1 ) + (y1 + y2 ) + (y2 + y3 ) + · · · + (yn−2 + yn−1 ) + (yn−1 + yn )
2 2 2 2 2
h
= [(y0 + y1 ) + (y1 + y2 ) + (y2 + y3 ) + · · · + (yn−2 + yn−1 ) + (yn−1 + yn )]
2
h
= [(y0 + yn ) + 2(y1 + y2 + · · · + yn−1 )]
2

This is known as the composite trapezoidal rule. Henceforth we drop the word “compos-
ite”and say trapezoidal rule for the composite trapezoidal rule.
Simpson’s one third rule:Putting n = 2 in the Newton-Cotes’ quadrature formula, we
get
 
1 2
I = h 2y0 + 24y0 + 4 y0
3
36

i.e.,

Zb  
1
f (x) dx = h 2y0 + 2(y1 − y0 ) + (y2 − 2y1 + y0 )
3
a
 
1
= h 2y0 + 2y1 − 2y0 + (y2 − 2y1 + y0 )
3
 
1
= h 2y1 + (y2 − 2y1 + y0 )
3
h
= (y0 + 4y1 + y2 )
3

Dividing the interval [a, b] into n sub-intervals of width h(n, an even number) and using
Simpson’s one third rule we get,
x
Zr+2
h
f (x) dx = (yr + 4yr+1 + yr+2 ) (0 ≤ r ≤ n − 2)
3
xr

Therefore,

Zb Zxn
f (x) dx = f (x) dx
a x0
Zx2 Zx4 Zx6 xZn−2 Zxn
= f (x) dx + f (x) dx + f (x) dx + · · · + f (x) dx + f (x) dx
x0 x2 x4 xn−4 xn−2
h h h
= (y0 + 4y1 + y2 ) + (y2 + 4y3 + y4 ) + (y4 + 4y5 + y6 )
3 3 3
h h
+ · · · + (yn−4 + 4yn−3 + yn−2 ) + (yn−2 + 4yn−1 + yn )
3 3
h
= [(y0 + 4y1 + y2 ) + (y2 + 4y3 + y4 ) + (y4 + 4y5 + y6 )
3
+ · · · + (yn−4 + 4yn−3 + yn−2 ) + (yn−2 + 4yn−1 + yn )]
h
= [(y0 + yn ) + 2(y2 + y4 + y6 · · · + yn−2 ) + 4(y1 + y3 + y5 · · · + yn−1 )]
3

This is known as the composite Simpson’s one third rule. Henceforth we drop the word
“composite”and say Simpson’s one third rule for the composite Simpson’s one third rule.

Simpson’s (3/8)th rule: Putting n = 3 in the Newton-Cotes’ quadrature formula, we

 37

get
 
9 9 2 3 3
I = h 3y0 + 4y0 + 4 y0 + 4 y0
2 4 8
3h
8y0 + 124y0 + 642 y0 + 43 y0
 
=
8
i.e.,
Zb
3h
f (x) dx = [8y0 + 12(y1 − y0 ) + 6(y2 − 2y1 + y0 ) + (y3 − 3y2 + 3y1 − y0 )]
8
a
3h
= [8y0 + 12y1 − 12y0 + 6y2 − 12y1 + 6y0 + y3 − 3y2 + 3y1 − y0 ]
8
3h
= [y0 + 3y1 + 3y2 + y3 ]
8
This is known as the Simpson’s (3/8)th rule.
Dividing the interval [a, b] into n sub-intervals of width h(n, a multiple of 3) and using
Simpson’s (3/8)th rule we get,
x
Zr+3
3h
f (x) dx = (yr + 3yr+1 + 3yr+2 + yr+3 ) (0 ≤ r ≤ n − 3)
8
xr

Therefore,
Zb Zxn
f (x) dx = f (x) dx
a x0
Zx3 Zx6 Zx9 xZn−3 Zxn
= f (x) dx + f (x) dx + f (x) dx + · · · + f (x) dx + f (x) dx
x0 x3 x6 xn−6 xn−3
3h 3h 3h
= (y0 + 3y1 + 3y2 + y3 ) + (y3 + 3y4 + 3y5 + y6 ) + (y6 + 3y7 + 3y8 + y9 )
8 8 8
3h 3h
+ · · · + (yn−6 + 3yn−5 + 3yn−4 + yn−3 ) + (yn−3 + 3yn−2 + 3yn−1 + yn )
8 8
3h
= [(y0 + 3y1 + 3y2 + y3 ) + (y3 + 3y4 + 3y5 + y6 ) + (y6 + 3y7 + 3y8 + y9 )
8
+ · · · + (yn−6 + 3yn−5 + 3yn−4 + yn−3 ) + (yn−3 + 3yn−2 + 3yn−1 + yn )]
i.e.,
Zb
3h
f (x) dx = [(y0 + yn ) + 2(y3 + y6 + y9 + · · · + yn−3 ) + 3(y1 + y2 + y4 + y5 + · · · + yn−2 + yn−1 )]
8
a
38

This is known as the composite Simpson’s (3/8)th rule. Henceforth we drop the word
“composite”and say Simpson’s (3/8)th rule for the composite Simpson’s (3/8)th rule.
Weddle’s rule: Putting n = 6 in the Newton-Cotes’ quadrature formula and replacing
41 3
140
46 y0 by 10 46 y0 we get

Zb
3h
y dx = (y0 + 5y1 + y2 + 6y3 + y4 + 5y5 + y6 )
10
a

R6 dx
Example 3.2.1. Evaluate 1+x2
, using trapezoidal rule by taking seven ordinates. Com-
0
pare your answer with the actual value of the integral.

1
Solution. Let y = . Since we need to take seven ordinates, we divide the interval
1 + x2
1
[0, 6] into six sub-intervals of width h = 1. The values of the function y = are as
1 + x2
given below.

x 0 1 2 3 4 5 6
y 1 0.5 0.2 0.1 0.0588 0.03846 0.027

According to trapezoidal rule,

Zxn
h
y dx = [(y0 + yn ) + 2(y1 + y2 + · · · + yn−1 )]
2
x0

Since n = 6 we get,

Zx6
h
y dx = [(y0 + y6 ) + 2(y1 + y2 + y3 + y4 + y5 )]
2
x0

1
Here, x0 = 0, x6 = 6, y = , h = 1, y0 = 1, y1 = 0.5, y2 = 0.2, y3 = 0.1, y4 = 0.0588,
1 + x2
y5 = 0.03846 and y6 = 0.027. Therefore,

Z6
dx 1
2
= [(1 + 0.027) + 2(0.5 + 0.2 + 0.1 + 0.0588 + 0.03846)]
1+x 2
0
1
= [1.027 + 2(0.89726)]
2
= 1.41076(= IT , say)
 39

Note that,

Z6
dx 6
2
= tan−1 (x) 0
1+x
0
= tan−1 (6) − tan−1 (0)
= tan−1 (6)
≈ 1.4056(= IA , say)

| IT − IA |= 0.00516. Therefore, the value of the given integral got by using trapezoidal
rule agrees with its actual value up to two decimal places.
R1 dx
Example 3.2.2. Evaluate 1+x
, using trapezoidal rule by taking h = 1/6 and hence find
0
an approximate value of log 2.
1 1
Solution. Let y = . The values of the function y = are as given below.
1+x 1+x

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

y 1 6/7 6/8 6/9 6/10 6/11 0.5

According to trapezoidal rule,

Zxn
h
y dx = [(y0 + yn ) + 2(y1 + y2 + · · · + yn−1 )]
2
x0

Since n = 6 we get,
Zx6
h
y dx = [(y0 + y6 ) + 2(y1 + y2 + y3 + y4 + y5 )]
2
x0

1
Here, x0 = 0, x6 = 1, y = , h = 1/6, y0 = 1, y1 = 6/7, y2 = 6/8, y3 = 6/9,
1+x
y4 = 6/10, y5 = 6/11 and y6 = 0.5. Therefore,

Z1   
dx 1/6 6 6 6 6 6
= (1 + 0.5) + 2 + + + +
1+x 2 7 8 9 10 11
0
1
≈ [1.5 + 2(3.41926)]
12
= 0.694877
40

Note that,

Z1
dx
= log(1 + x)|10
1+x
0
= log 2 − log 1
= log 2

Therefore, log 2 ≈ 0.694877

R2 2x+1
Example 3.2.3. Evaluate x2 +x
dx, using Simpson’s one third rule by taking h = 0.1
1
and hence find an approximate value of log 3.
2x+1 2x+1
Solution. Let y = x2 +x
. The values of the function y = x2 +x
are as given below.

x y
1 1.5
1.1 1.3853
1.2 1.2879
1.3 1.204
1.4 1.131
1.5 1.0667
1.6 1.0096
1.7 0.9586
1.8 0.9127
1.9 0.8711
2 0.8333

According to Simpson’s one third rule,

Zxn
h
y dx = [(y0 + yn ) + 2(y2 + y4 + · · · + yn−2 ) + 4(y1 + y3 + · · · + yn−1 )]
3
x0

Since n = 10 we get,
Zx10
h
y dx = [(y0 + y10 ) + 2(y2 + y4 + y6 + y8 ) + 4(y1 + y3 + y5 + y7 + y9 )]
3
x0

2x + 1
Here, x0 = 1, x10 = 2, y = 2 , h = 0.1, y0 = 1.5, y1 = 1.3853, y2 = 1.2879,
x +x
y3 = 1.204, y4 = 1.131, y5 = 1.0667, y6 = 1.0096, y7 = 0.9586, y8 = 0.9127, y9 = 0.8711
 41

and y10 = 0.8333. Therefore,

Z2
2x + 1 0.1
2
dx = [(1.5 + 0.8333) + 2 (1.2879 + 1.131 + 1.0096 + 0.9127)
x +x 3
1
+4(1.3853 + 1.204 + 1.0667 + 0.9586 + 0.8711)]
0.1
= [2.3333 + 2(4.3412) + 4(5.4857)]
3
≈ 1.098617

Note that,

Z2
2x + 1 2
2
dx = log(x2 + x) 1
x +x
1
 Z 0 
f (x)
because dx = log f (x) + C
f (x)
= log 6 − log 2
6
= log
2
= log 3

Therefore, log 3 ≈ 1.098617.

R2 2
Example 3.2.4. Evaluate ex dx, using Simpson’s one third rule by taking h = 0.2.
0

2 2
Solution. Let y = ex . The values of the function y = ex are as given below.

x y
0 1
0.2 1.0408
0.4 1.1735
0.6 1.4333
0.8 1.8965
1.0 2.7183
1.2 4.2207
1.4 7.0993
1.6 12.9358
1.8 25.5337
2 54.5982
42

According to Simpson’s one third rule,

Zxn
h
y dx = [(y0 + yn ) + 2(y2 + y4 + · · · + yn−2 ) + 4(y1 + y3 + · · · + yn−1 )]
3
x0

Since n = 10 we get,
Zx10
h
y dx = [(y0 + y10 ) + 2(y2 + y4 + y6 + y8 ) + 4(y1 + y3 + y5 + y7 + y9 )]
3
x0

2
Here, x0 = 0, x10 = 2, y = ex , h = 0.2, y0 = 1, y1 = 1.0408, y2 = 1.1735, y3 = 1.4333,
y4 = 1.8965, y5 = 2.7183, y6 = 4.2207, y7 = 7.0993, y8 = 12.9358, y9 = 25.5337 and
y10 = 54.5982. Therefore,

Z2
2 0.2
ex dx = [(1 + 54.5982) + 2 (1.1735 + 1.8965 + 4.2207 + 12.9358)
3
0
+4(1.0408 + 1.4333 + 2.7183 + 7.0993 + 25.5337)]
0.2
= [55.5982 + 2(20.2265) + 4(37.8254)]
3
≈ 16.49019

R1 dx
Example 3.2.5. Evaluate 1+x
, using Simpson’s (3/8)th rule by taking h = 1/6 and
0
hence find an approximate value of log 2.
1 1
Solution. Let y = . The values of the function y = are as given below.
1+x 1+x
x 0 1/6 2/6 3/6 4/6 5/6 1
y 1 6/7 6/8 6/9 6/10 6/11 0.5

According to Simpson’s (3/8)th rule,

Zxn
3h
y dx = [(y0 + yn ) + 2(y3 + y6 + · · · + yn−3 ) + 3(y1 + y2 + y4 + y5 + · · · + yn−2 + yn−1 )]
8
x0

Since n = 6 we get,
Zx6
3h
y dx = [(y0 + y6 ) + 2y3 + 3(y1 + y2 + y4 + y5 )]
8
x0
 43

1
Here, x0 = 0, x6 = 1, y = , h = 1/6, y0 = 1, y1 = 6/7, y2 = 6/8, y3 = 6/9,
1+x
y4 = 6/10, y5 = 6/11 and y6 = 0.5. Therefore,

Z1     
dx 3(1/6) 6 6 6 6 6
= (1 + 0.5) + 2 +3 + + +
1+x 8 9 7 8 10 11
0
1 4
≈ [1.5 + + 3(2.7526)]
16 3
≈ 0.6931958

Note that,

Z1
dx
= log(1 + x)|10
1+x
0
= log 2 − log 1
= log 2

Therefore, log 2 ≈ 0.6931958

Example 3.2.6. The velocity v(m/sec) of a particle at a distance s meters from a point
on its path is given in the following table.

s 0 10 20 30 40 50 60
v 47 58 64 65 61 52 38

Using Simpson’s one third rule, estimate the time taken by the particle to travel 60 meters.

ds 1
Solution. We know that v = . This implies dt = ds. The time taken by the particle
dt v
R60 1
to travel 60 meters is given by ds. To evaluate this integral we use Simpson’s one
0 v
third rule. According to Simpson’s one third rule,
Zxn
h
y dx = [(y0 + yn ) + 2(y2 + y4 + · · · + yn−2 ) + 4(y1 + y3 + · · · + yn−1 )]
3
x0

Since n = 6 we get,

Zx6
h
y ds = [(y0 + y6 ) + 2(y2 + y4 ) + 4(y1 + y3 + y5 )]
3
x0
44

1 1 1 1 1 1
Here, x0 = 0, x6 = 60, y = , h = 10, y0 = , y1 = , y2 = , y3 = , y4 = ,
v 47 58 64 65 61
1 1
y5 = , y6 = . Therefore,
52 38
Z60      
1 10 1 1 1 1 1 1 1
ds = + +2 + +4 + +
v 3 47 38 64 61 58 65 52
0
10
≈ [(0.04759) + 2(0.03202) + 4(0.05186)]
3
≈ 1.06357
Thus, the required time is 1.06357 seconds.
Example 3.2.7. A solid of revolution is formed by rotating about the x-axis, the area
between the lines x = 0, x = 1, the x-axis and a curve through the points with the
following coordinates.
x 0 0.25 0.5 0.75 1
y 1 0.9896 0.9589 0.9089 0.8415
Estimate the volume of the solid formed using Simpson’s rule.
R1 R1
Solution. The required volume is given by πy 2 dx i.e., π y 2 dx. To evaluate the integral
0 0
R1
y 2 dx we use Simpson’s one third rule. According to Simpson’s one third rule,
0

Zxn
h
y dx = [(y0 + yn ) + 2(y2 + y4 + · · · + yn−2 ) + 4(y1 + y3 + · · · + yn−1 )]
3
x0

Since n = 4 we get,
Zx4
h
y dx = [(y0 + y4 ) + 2y2 + 4(y1 + y3 )]
3
x0

Here, x0 = 0, x4 = 1, h = 0.25, y0 = 1, y1 = (0.9896)2 , y2 = (0.9589)2 , y3 = (0.9089)2 ,

y4 = (0.8415)2 . Therefore,
Z1
0.25 
y 2 dx = 1 + (0.8415)2 + 2 (0.9589)2 + 4 (0.9896)2 + (0.9089)2




3
0
0.25
≈ [1 + 0.7081 + 2(0.9195) + 4(0.9793 + 0.8261)]
3
≈ 0.8974
Thus, the required volume is π(0.8974) ≈ 2.8192.
 45

R1 1
Example 3.2.8. Evaluate 1+x2
dx using Weddle’s rule by taking seven ordinates.
0

1
Solution. Here y = . the function y in the tabular form is as shown below.
1 + x2
x 0 1/6 2/6 3/6 4/6 5/6 1
y 1 36/37 36/40 36/45 36/52 36/61 0.5

According to Weddle’s rule,

Zb
3h
y dx = (y0 + 5y1 + y2 + 6y3 + y4 + 5y5 + y6 )
10
a

1
Here a = 0, b = 1, y = , h = 1/6, y0 = 1, y1 = 36/37, y2 = 36/40, y3 = 36/45,
1 + x2
y4 = 36/52, y5 = 36/61 and y6 = 0.5
Therefore,

Z1            
1 3(1/6) 36 36 36 36 36
dx = 1+5 + +6 + +5 + 0.5
1 + x2 10 37 40 45 52 61
0
1
= [1 + 4.8649 + 0.9 + 4.8 + 0.6923 + 2.9508 + 0.5]
20
= 0.7854
Chapter 4

Numerical solution of equations

Many real world problems lead us to either single equation or a system of equations and
we are required solve such equation(s). When there is no analytical method for finding
the exact solution, we use numerical methods to get an approximate value of the required
solution.
In this chapter we shall discuss few numerical methods to solve single equation/system of
equations.

4.1 Solution of algebraic/transcendental equations

b
We know that the solution(or root) of the equation ax + b = 0(a 6= 0) is x = − and
√ a
−b ± b 2 − 4ac
the roots of the equation ax2 + bx + c = 0(a 6= 0) are x = . We do have
2a
formulae (in terms of coefficients and radicals) for cubic and biquadratic equations which
give exact value of the roots. But there is no such formula for an algebraic equation of
degree greater than four. There are several numerical methods which can be used to find
an approximate value of a root of the given algebraic/transcendental equation.
Definition 4.1.1. x = a is said to be a root of the equation f (x) = 0 if f (a) = 0
Example 4.1.2. 1. x = 0 is a root of the equation x2 = 0.
2. x = 0, π and 2π are roots of sin x = 0 (0 ≤ x ≤ 2π).
Note. A root of the equation f (x) = 0 can also be described as a point at which the
graph of y = f (x) either touches or crosses the x-axis.
The graphs of the functions are as shown in the Figures 4.1 and 4.2. The following
theorem gives a very useful information about the existence of a root of f (x) = 0 between
the two points.
Theorem 4.1.3. If f (x) is a continuous function such that f (a)f (b) < 0 then there exists
at least one root of f (x) = 0 between a and b.

46
 47

Figure 4.1: Curve in Example 4.1.2(1)

Figure 4.2: Curve in Example 4.1.2(2)

4.1.1 Bisection method

Let f (x) be a continuous function between x = a and x = b and let f (a)f (b) < 0.
Then there is at least one root of the equation f (x) = 0 between a and b. To find an
approximate value of a root of f (x) = 0, we divide the interval [a, b] into two equal parts
and take x1 := a+b 2
, the mid point of a and b, as the first approximation to the required
root. If f (x1 = 0 then x1 is the required root. Otherwise a root lies between either a and
x1 or x1 and b depending upon f (a)f (x1 ) < 0 or f (x1 )f (b) < 0. This process is repeated
till a root to the desired degree of accuracy is obtained.
Example 4.1.4. Find a root of the equation x3 − 2x − 5 = 0 in [2, 3] using the bisection
method correct to two decimal places.
Solution. Let f (x) = x3 − 2x − 5. Since f (2) = 23 − 2(2) − 5 = −1 < 0 and f (3) =
33 − 2(3) − 5 = 16 > 0, there is a root of x3 − 2x − 5 = 0 between 2 and 3.
2+3
The first approximation to the root is x1 = = 2.5. Note that f (x1 ) = f (2.5) =
2
3
(2.5) − 2(2.5) − 5 = 5.625 > 0. Therefore, a root lies between 2 and 2.5.
2 + 2.5
The second approximation to the root is x2 = = 2.25. Note that f (x2 ) =
2
f (2.25) = (2.25)3 − 2(2.25) − 5 = 1.890625 > 0. Therefore, a root lies between 2 and 2.25.
2 + 2.25
The third approximation to the root is x3 = = 2.125. Note that f (x3 ) =
2
f (2.125) = (2.125)3 − 2(2.125) − 5 ≈ 0.3457 > 0. Therefore, a root lies between 2 and
2.125.
The successive approximation to the root got by repeating this process are,
x4 = 2.0625 x5 = 2.09375 x6 = 2.159375
x7 = 2.1265625 x8 = 2.11015625 x9 = 2.101953125
Since | x9 − x8 |=| 2.101953125 − 2.11015625 |< 10−2 , a root of the given equation correct
to two decimal places is 2.101953125.
48

Figure 4.3: Successive approximations in bisection method

Example 4.1.5. Using the bisection method, find a real root of the equation cos x =
x exp(x) in [0, 1] correct to three decimal places.
Solution. Let f(x) := cos (x) − x exp(x). Since f(0) := cos (0) − 0 exp(0) = 1 > 0 and
f(1) := cos (1) − exp(1) ≈ −1.7184, there is a root of f(x) := cos (x) − x exp(x) between
0 and 1.
0+1
The first approximation to the root is x1 = = 0.5. Note that f(x1 ) = f(0.5) =
2
cos (0.5) − (0.5) exp(0.5) ≈ 0.05322 > 0. Therefore, a root lies between 0.5 and 1.
0.5 + 1
The second approximation to the root is x2 = = 0.75. Note that f(x2 ) = f(0.75) =
2
cos (0.75) − (0.75) exp(0.75) ≈ −0.85606 < 0. Therefore, a root lies between 0.5 and 0.75.
0.5 + 0.75
The third approximation to the root is x3 = = 0.625. Note that f(x3 ) =
2
f(0.625) = cos (0.625)−(0.625) exp(0.625) ≈ −0.35669 < 0. Therefore, a root lies between
0.5 and 0.625.
The successive approximation to the root got by repeating this process are,

x4 = 0.5625 x5 = 0.53125 x6 = 0.515625

x7 = 0.5234375 x8 = 0.51953125 x9 = 0.517578125

and x10 = 0.5185546875. Since | x10 − x9 |=| 0.5185546875 − 0.517578125 |< 10−3 , a root
of the given equation correct to three decimal places is 0.5185546875.
 49

Figure 4.4: Successive approximations in regula falsi method

4.1.2 Regula falsi method (or method of false position)

Let f (x) be a continuous function between x = a and x = b and let f (a)f (b) < 0. Then
there is at least one root of the equation f (x) = 0 between a and b. In this method, we
replace the curve y = f (x) (a ≤ x ≤ b) by the chord AB joining the two points (a, f (a))
and (b, f (b)) and take the point of intersection of the chord AB with the x-axis as an
approximation to the root. The x-coordinate of the point of intersection of the chord AB
with the x-axis is
b−a
c=a− f (a)
f (b) − f (a)
which is taken as the first approximation to the required root. The root lies between a
and c or c and b according as f (a)f (c) < 0 or f (c)f (b) < 0. This process is repeated till
the root is found to the desired degree of accuracy.

Example 4.1.6. Using the method of false position, find a real root of the equation
x3 − 3x + 1 in [0, 1]. Carry out three iterations.

Solution. Let f (x) = x3 − 3x + 1, a = 0 and let b = 1. Since f (a) = f (0) = 03 − 3(0) + 1 =

1 > 0 and f (b) = f (1) = 13 − 3(1) + 1 = −1 < 0, a root lies between 0 and 1.
In the method of false position, an approximate value of a root (lying between a and b)
50

of the equation f (x) = 0 is given by

b−a
c=a− f (a) (4.1.1)
f (b) − f (a)

Taking a = 0 and b = 1 in the RHS of Equation (4.1.1) we get

1−0
x1 := 0 − f (0)
f (1) − f (0)
Therefore,
1
x1 = − (1)
−1 − 1
= 0.5

Note that f (x1 ) = f (0.5) = −0.375 < 0. Therefore, a root lies between 0 and 0.5
Taking a = 0 and b = 0.5 in the RHS of Equation (4.1.1) we get
0.5 − 0
x2 := 0 − f (0)
f (0.5) − f (0)
Therefore,
0.5
x2 = − (1)
−0.375 − 1
≈ 0.36364

Note that f (x2 ) = f (0.36364) ≈ −0.0428 < 0. Therefore, a root lies between 0 and
0.36364
Taking a = 0 and b = 0.36364 in the RHS of Equation (4.1.1) we get
0.36364 − 0
x3 := 0 − f (0)
f (0.36364) − f (0)
Therefore,
0.36364
x3 = − (1)
−0.0428 − 1
≈ 0.3487

Example 4.1.7. Using the method of false position, find a real root of the equation
2x − log x = 7 in [4, 5]. Carry out three iterations.
Solution. Let f (x) = 2x − log x − 7, a = 4 and let b = 5. Since f (a) = f (4) = 2(4) −
log 4 − 7 ≈ −0.38629 < 0 and f (b) = f (5) = 2(5) − log 5 − 7 ≈ 1.39056 > 0, a root lies
between 4 and 5.
 51

In the method of false position, an approximate value of a root (lying between a and b)
of the equation f (x) = 0 is given by
b−a
c=a− f (a) (4.1.2)
f (b) − f (a)
Taking a = 4 and b = 5 in the RHS of Equation (4.1.2) we get
5−4
x1 := 4 − f (4)
f (5) − f (4)
Therefore,
1
x1 = 4 − (−0.38629)
1.39056 − (−0.38629)
≈ 4.2174
Note that f (x1 ) = f (4.2174) ≈ −0.0044 < 0. Therefore, a root lies between 4.2174 and 5
Taking a = 4.2174 and b = 5 in the RHS of Equation (4.1.2) we get
5 − 4.2174
x2 := 4.2174 − f (4.2174)
f (5) − f (4.2174)
Therefore,
0.7826
x2 = 4.2174 − (−0.0044)
1.39056 − (−0.0044)
≈ 4.2199
Note that f (x2 ) = f (4.2199) ≈ −0.00005 < 0. Therefore, a root lies between 4.2199 and
5
Taking a = 4.2199 and b = 5 in the RHS of Equation (4.1.2) we get
5 − 4.2199
x3 := 4.2199 − f (4.2199)
f (5) − f (4.2199)
Therefore,
0.7801
x3 = 4.2199 − (−0.00005)
1.39056 − (−0.00005)
≈ 4.22
Home Work

• Using regula falsi method, compute a real root of xex = 2 in [0, 1]. Carry out three
iterations.
• Find a root of cos x = 3x − 1 in [0, 1], using the method of false position. Carry out
three iterations.
52

4.1.3 Fixed point(or ordinary) iteration method

Definition 4.1.8. A point x = a is said to be a fixed point of f (x) if f (a) = a.

Example 4.1.9. 1. x = 0 is a fixed point of f (x) = sin x

2. x = 1 is a fixed point of f (x) = ex−1

There is a connection between the notion of a fixed point and that of a root. In fact,
x = a is a root of the equation f (x) = 0 if and only if it is a fixed point of g(x) := f (x)+x.
In the fixed point iteration method, we first rewrite the equation f (x) = 0 in the
form x = g(x) and then take x0 as an initial approximation to a root. We compute the
successive approximations of a root using xn+1 = g(xn ) (n = 0, 1, 2, · · · ).

Note. 1. The equation f (x) = 0 may be written in the form x = g(x) in several ways.
1
For example, x2 − x − 1 = 0 may be written as x = x2 − 1 or x = .
x−1
2. The fixed point iteration method does not always converge.

The following result gives a sufficient condition for the convergence of the fixed point
iteration method.

Theorem 4.1.10. Let x = a be a fixed point of g(x), g 0 (x) is continuous in some interval
I containing a and let | g 0 (x) |< 1 in I then lim xn = a where xn = g(xn−1 ) (n = 1, 2, · · · )
n→∞
and x0 is any point in I.

To find a root of the equation f (x) = 0 using fixed point iteration method, we proceed
as follows.

• Rewrite f (x) = 0 in the form x = g(x) in such a way that g 0 (x) is continuous in
some interval I and | g 0 (x) |< 1 in I.

• Take any point x0 ∈ I as an initial approximation.

• Use xn+1 = g(xn ) (n = 0, 1, 2, · · · ) to compute successive approximations of a root.

Note. | x |< a if and only if x ∈ (−a, a)

Example 4.1.11. Use the fixed point iteration method to find a root of x4 − 16x + 16 = 0
near x = 1. Carry out five iterations.
x4 + 16
Solution. Let f (x) = x4 − 16x + 16. We rewrite x4 − 16x + 16 = 0 as x = . Here
16
x4 + 16 x3 x3
g(x) = and | g 0 (x) |= . Clearly < 1 if and only if | x3 |< 4 if and
16 4 4
 53

√ √ √
only
√ | x |< 3 4. Therefore, | g 0 (x) |< 1 in the interval (− 3 4, 3 4). Clearly the interval
if √
(− 3 4, 3 4) contains 1. Now we take x0 = 1 as an initial approximation to a root and use

x4n + 16
xn+1 = (n = 0, 1, 2, · · · ) (4.1.3)
16
to compute successive approximations to a root.
Putting n = 0 in Equation (4.1.3), we get

x40 + 16
x1 =
16
1 + 16
=
16
= 1.0625
Putting n = 1 in Equation (4.1.3), we get

x41 + 16
x2 =
16
(1.0625)4 + 16
=
16
≈ 1.079652
Putting n = 2 in Equation (4.1.3), we get

x42 + 16
x3 =
16
(1.079652)4 + 16
=
16
≈ 1.0849
Putting n = 3 in Equation (4.1.3), we get

x43 + 16
x4 =
16
(1.0849)4 + 16
=
16
≈ 1.08658
Putting n = 4 in Equation (4.1.3), we get

x44 + 16
x5 =
16
(1.08658)4 + 16
=
16
≈ 1.087
54

Example 4.1.12. Use the fixed point iteration method to find a root of 5x = ex near
x = 0.2. Carry out five iterations.

ex ex
Solution. Let f (x) = 5x − ex . We rewrite 5x − ex = 0 as x = . Here g(x) = and
5 5
ex ex
| g 0 (x) |= . Clearly < 1 if and only if | ex |< 5 if and only if ex < 5 if and
5 5
only if x < log 5. Therefore, | g 0 (x) |< 1 in the interval (−∞, log 5). Clearly the interval
(−∞, log 5) contains 0.2. Now we take x0 = 0.2 as an initial approximation to a root and
use
exn
xn+1 = (n = 0, 1, 2, · · · ) (4.1.4)
5
to compute successive approximations to a root.
Putting n = 0 in Equation (4.1.4), we get

ex0
x1 =
5
e0.2
=
5
≈ 0.24428

Putting n = 1 in Equation (4.1.4), we get

ex1
x2 =
5
e0.24428
=
5
≈ 0.25534

Putting n = 2 in Equation (4.1.4), we get

ex2
x3 =
5
e0.25534
=
5
≈ 0.25818

Putting n = 3 in Equation (4.1.4), we get

ex3
x4 =
5
e0.25818
=
5
≈ 0.25891
 55

Putting n = 4 in Equation (4.1.4), we get

ex4
x5 =
5
e0.25891
=
5
≈ 0.2591

4.1.4 Newton-Raphson’s method

Let f (x) be a function such that f 0 (x) is continuous and let x0 be a point near a root α
of f (x) = 0. In Newton-Raphson’s method, we replace the curve y = f (x) by the tangent
T to the curve at x = x0 and take the point of intersection of the tangent T with the
x-axis as an approximation to the root. The point of intersection of the tangent T with
the x-axis is given by
f (x0 )
x1 = x0 − 0
f (x0 )
and is taken as the first approximation to the root. This process is repeated till the root
is found to the desired degree of accuracy.

Note. 1. Newton-Raphson’s method does not always converge.

2. In Newton-Raphson’s method, the successive approximations to a root of f (x) = 0

f (xn )
is computed using the formula xn+1 = xn − 0 (n = 0, 1, 2, · · · ) and this formula
f (xn )
is known as Newton-Raphson’s formula

Example 4.1.13. Using Newton-Raphson’s method, find a real root of x3 − 5x + 3 = 0

near 2. Carry out three iterations.

Solution. Let f (x) = x3 − 5x + 3 and let x0 = 2. Note that f 0 (x) = 3x2 − 5. Newton-
Raphson’s formula is
f (xn )
xn+1 = xn − 0 (n = 0, 1, 2, · · · ) (4.1.5)
f (xn )
Putting n = 0 in Equation (4.6), we get

f (x0 )
x1 = x0 −
f 0 (x0 )
f (2)
=2− 0
f (2)
1
= 2 − (as f (2) = 1 and f 0 (2) = 7)
7
≈ 1.8571
56

Figure 4.5: Successive approximations in Newton-Raphson method

Putting n = 1 in Equation (4.6), we get

f (x1 )
x2 = x1 −
f 0 (x1 )
f (1.8571)
= 1.8571 − 0
f (1.8571)
0.1193
= 1.8571 − (as f (1.8571) = 0.1193 and f 0 (1.8571) = 5.3465)
5.3465
≈ 1.8348

Putting n = 2 in Equation (4.6), we get

f (x2 )
x3 = x2 −
f 0 (x2 )
f (1.8348)
= 1.8348 − 0
f (1.8348)
0.0028
= 1.8348 − (as f (1.8348) = 0.0028 and f 0 (1.8348) = 5.0995)
5.0995
≈ 1.8343

Example 4.1.14. Using Newton-Raphson’s method, find a real root of 2x = cos x + 3

near π/2 correct to three decimal places.
 57

Solution. Let f (x) = 2x − cos x − 3 and let x0 = π/2. Note that f 0 (x) = 2 + sin x.
Newton-Raphson’s formula is
f (xn )
xn+1 = xn − (n = 0, 1, 2, · · · ) (4.1.6)
f 0 (xn )
Putting n = 0 in Equation (4.7), we get
f (x0 )
x1 = x0 −
f 0 (x0 )
π f (π/2)
= − 0
2 f (π/2)
π π−3
= − (as f (π/2) = π − 3 and f 0 (π/2) = 3)
2 3
≈ 1.5236
Putting n = 1 in Equation (4.7), we get
f (x1 )
x2 = x1 −
f 0 (x1 )
f (1.5236)
= 1.5236 − 0
f (1.5236)
0.00002
= 1.5236 − (as f (1.5236) = 0.00002 and f 0 (1.5236) = 2.9989)
2.9989
≈ 1.5236
Since | x2 − x1 |< 10−3 , a root of the given equation near π/2 correct to three decimal
places is 1.5236
Example 4.1.15. Use Newton-Raphson method to find a root of tan x = 1.5x near 1,
correct to three decimal places.
Solution. Let f (x) = tan x − 1.5x and let x0 = 1. Note that f 0 (x) = sec2 x − 1.5 =
(1 + tan2 x) − 1.5 = tan2 x − 0.5. Newton-Raphson’s formula is
f (xn )
xn+1 = xn − (n = 0, 1, 2, · · · ) (4.1.7)
f 0 (xn )
Putting n = 0 in Equation (4.1.7), we get
f (x0 )
x1 = x0 −
f 0 (x0 )
f (1)
=1− 0
f (1)
0.0574
=1− (as f (1) = 0.0574 and f 0 (1) = 1.9255)
1.9255
≈ 0.9702
58

Putting n = 1 in Equation (4.1.7), we get

f (x1 )
x2 = x1 −
f 0 (x1 )
f (0.9702)
= 0.9702 − 0
f (0.9702)
0.0045
= 0.9702 − (as f (0.9702) = 0.0045 and f 0 (0.9702) = 1.6311)
1.6311
≈ 0.9674

Putting n = 2 in Equation (4.1.7), we get

f (x2 )
x3 = x2 −
f 0 (x2 )
f (0.9674)
= 0.9674 − 0
f (0.9674)
0.000004
= 0.9674 − (as f (0.9674) = −0.000004 and f 0 (0.9674) = 1.6057)
1.6057
≈ 0.9674

Since | x3 −x2 |< 10−3 , a root of the given equation near 1, correct to three decimal places
is 0.9674

Convergence analysis

Let xn be an approximate value of a root α of f (x) = 0, got in the nth iteration of Newton-
Raphson’s method. Let xn differs from α by a small quantity n so that xn = α + n .
According to Newton-Raphson’s formula

f (xn )
xn+1 = xn −
f 0 (xn )

Therefore,

f (α + n )
α + n+1 = α + n −
f 0 (α + n )
 59

This implies,

f (α + n )
n+1 = n −
f 0 (α + n )
f (α) + n f 0 (α) + 12 2n f 00 (α) + · · ·
= n − 0
f (α) + n f 00 (α) + 12 2n f 000 (α) + · · ·
(by Taylor’s series expansion)
n f 0 (α) + 21 2n f 00 (α) + · · ·
= n −
f 0 (α) + n f 00 (α) + 12 2n f 000 (α) + · · ·
(as f (α) = 0)
n [f 0 (α) + n f 00 (α) + 12 2n f 000 (α) + · · · ] − [n f 0 (α) + 21 2n f 00 (α) + · · · ]
=
f 0 (α) + n f 00 (α) + 21 2n f 000 (α)
n f 0 (α) + 2n f 00 (α) + 12 3n f 000 (α) − n f 0 (α) − 12 2n f 00 (α) − · · ·
=
f 0 (α) + n f 00 (α) + 12 2n f 000 (α) + · · ·
n f 0 (α) + 2n f 00 (α) − n f 0 (α) − 12 2n f 00 (α)
≈
f 0 (α) + n f 00 (α) + 12 2n f 000 (α)
(neglecting third and higher powers of n )
1 2 00
 f (α)
2 n
=
f 0 (α) + n f 00 (α) + 12 2n f 000 (α)
2 f 00 (α)
≈ n 0
2f (α)
2 f 00 (α)
= n 0
2 f (α)

This shows that the subsequent error at each step is proportional to the square of the
previous error and hence the convergence is quadratic.

Note. Newton-Raphson’s method may some time diverge(see Figures 4.6 and 4.7)
60

Figure 4.6: The tangent is parallel to the x-axis

Figure 4.7: The successive approximations in N-R method go away from the root
 61

4.2 Solution of system of linear equations

In this section we shall study few methods to solve system of n linear equations in n
unknowns.

4.2.1 Cramer’s rule

Consider the following system of n linear equations in n unknowns.
a11 x1 + a12 x2 + · · · + a1n xn = b1
a21 x1 + a22 x2 + · · · + a2n xn = b2
(4.2.1)
··························· = ···
an1 x1 + an2 x2 + · · · + ann xn = bn
If D :=| A |6= 0 where A := (aij ) is the coefficient matrix of the system (4.2.1) then the
system has a unique solution and is given by
D1 D2 Dn
x1 = x2 = ··· xn =
D D D
where Dr (1 ≤ r ≤ n) is the determinant obtained from D by replacing in D the rth
column by the column with the entries b1 , b2 , · · · , bn .
As a special case, we shall consider now a system of three linear equations in three
unknowns.
a11 x1 + a12 x2 + a13 x3 = b1
a21 x1 + a22 x2 + a23 x3 = b2 (4.2.2)
a31 x1 + a32 x2 + a33 x3 = b3
 
a11 a12 a13
The coefficient matrix of the system (4.2.2) is A =  a21 a22 a23 . If D :=| A |6= 0
a31 a32 a33
then, by Cramer’s rule, the system (4.2.2) has a unique solution and is given by
D1 D2 D3
x1 = x2 = x3 =
D D D
where
b1 a12 a13 a11 b1 a13 a11 a12 b1
D1 = b2 a22 a23 D2 = a21 b2 a23 D3 = a21 a22 b2
b3 a32 a33 a31 b3 a33 a31 a32 b3
Example 4.2.1. Solve the following system of linear equations by Cramer’s rule.
3x − 2y − z = 3
2x + y + 2z = 5
x+y+z =2
62

 
3 −2 −1
Solution. The coefficient matrix of the given system of equations is A =  2 1 2 .
1 1 1
Let D =| A |. Note that
1 2 2 2 2 1
D =3· − (−2) · −1·
1 1 1 1 1 1
= 3(1 − 2) + 2(2 − 2) − (2 − 1)
= −4
6= 0
Therefore, by Cramer’s rule, the given system of linear equations has a unique solution
and is given by
D1 D2 D3
x= y= z=
D D D
where
3 −2 −1
D1 = 5 1 2
2 1 1
1 2 5 2 5 1
=3· − (−2) · −1·
1 1 2 1 2 1
= 3(1 − 2) + 2(5 − 4) − 1(5 − 2)
= −4

3 3 −1
D2 = 2 5 2
1 2 1
5 2 2 2 2 5
=3· −3· −1·
2 1 1 1 1 2
= 3(5 − 4) − 3(2 − 2) − 1(4 − 5)
=4
and
3 −2 3
D3 = 2 1 5
1 1 2
1 5 2 5 2 1
=3· − (−2) · +3·
1 2 1 2 1 1
= 3(2 − 5) + 2(4 − 5) + 3(2 − 1)
= −8
 63

−4 4 −8
Thus, x = = 1, y = = −1 and z = =2
−4 −4 −4
Example 4.2.2. Solve the following system of linear equations by Cramer’s rule.

x + 3y − 2z = 1
2x − 4y + 3z = 3
x + y + 2z = 9
 
1 3 −2
Solution. The coefficient matrix of the given system of equations is A =  2 −4 3 .
1 1 2
Let D =| A |. Clearly D = −26 6= 0. Therefore, by Cramer’s rule, the given system of
linear equations has a unique solution and is given by
D1 D2 D3
x= y= z=
D D D
where
1 3 −2 1 1 −2 1 3 1
D1 = 3 −4 3 D2 = 2 3 3 D3 = 2 −4 3
9 1 2 1 9 2 1 1 9

−26 −52 −78

Thus, x = = 1, y = = 2 and z = =3
−26 −26 −26

4.2.2 LU factorization method

When n is large, solving the system of n linear equations in n unknowns by Cramer’s rule
is impractical. In this section, we shall discuss one of the practical methods for solving
such systems.
An LU -factorization of a square matrix A is of the form A = LU where L is a lower
triangular matrix and U is an upper triangular matrix.
    
6 2 2 0 3 1
Example 4.2.3. 1. If A = then A = is an LU factor-
3 4 1 3 0 1
ization of A.
    
2 1 −1 0 −2 −1
2. If A = then A = is an LU factorization of A.
6 5 −3 2 0 1
    
10 2 4 2 0 0 5 1 2
3. If A =  15 4 3  then A =  3 −1 0   0 −1 3  is an LU factor-
0 −1 1 0 1 −2 0 0 1
ization of A.
64

 
0 2
Note. 1. Not all square matrices have an LU -factorization. For example A =
3 1
doesn’t have an LU -factorization. But it can be shown that, if a square matrix is
nonsingular then the matrix got by reordering the rows of the given matrix, has an
LU -factorization.
 Note
 that, by  interchanging
  therows of the matrix A we get the
3 1 1 0 3 1
matrix B = and B = .
0 2 0 1 0 2
2. Our aim is to solve the given system of linear equations. The solution of the given
system doesn’t change by reordering the equations of the system.
Consider the following system of equations.
a11 x1 + a12 x2 + · · · + a1n xn = b1
a21 x1 + a22 x2 + · · · + a2n xn = b2
(4.2.3)
··························· = ···
an1 x1 + an2 x2 + · · · + ann xn = bn
The system
 (4.2.3) can be written
 in a more convenient form as AX = B where
a11 a12 · · · a1n
 a21 a22 · · · a2n 
A=  · · · · · · · · · · · ·  (called the coefficient matrix of the system (4.2.3)),


a a · · · ann
 n1  n2  
x1 b1
 x2   b2 
X =  ..  and B =  .. . In LU -factorization method we write A = LU . Then
   
 .   . 
xn bn
AX = B gives LU X = B which implies L(U X) = B. This may be written as
LY = B (4.2.4)
where
UX = Y (4.2.5)
On solving Equation (4.2.4), we get Y . Using this Y , we solve the Equation (4.2.5) which
gives the solution of the given system also.

In LU -factorization method, if we take all the principal diagonal entries of the lower
triangular matrix L as 1 then the method is called Doolittle’s method and if we take all
the principal diagonal entries of the upper triangular matrix U as 1 then the method is
called Crout’s method.
Example 4.2.4. Solve the following system of equations by Doolittle’s method.
2x + 5y − 3z = 5
2x + 6y − 4z = 6
4x + 11y − 6z = 10
 65

 
2 5 −3
Solution. The coefficient matrix of the given system is A =  2 6 −4 .
4 11 −6
Let   
1 0 0 u11 u12 u13
A =  l21 1 0   0 u22 u23 
l31 l32 1 0 0 u33
Then    
2 5 −3 u11 u12 u13
 2 6 −4  =  l21 u11 l21 u12 + u22 l21 u13 + u23 
4 11 −6 l31 u11 l31 u12 + l32 u22 l31 u13 + l32 u23 + u33
Therefore,
u11 = 2 (4.2.6)
u12 = 5 (4.2.7)
u13 = −3 (4.2.8)
l21 u11 = 2 (4.2.9)
l21 u12 + u22 = 6 (4.2.10)
l21 u13 + u23 = −4 (4.2.11)
l31 u11 = 4 (4.2.12)
l31 u12 + l32 u22 = 11 (4.2.13)
l31 u13 + l32 u23 + u33 = −6 (4.2.14)
From the Equations (4.2.6) and (4.2.9)
2 2
l21 = = =1 (4.2.15)
u11 2
From the Equations (4.2.7), (4.2.10) and (4.2.15)

u22 = 6 − l21 u12 = 6 − 1(5) = 1 (4.2.16)

From the Equations (4.2.8), (4.2.11) and (4.2.15)

u23 = −4 − l21 u13 = −4 − 1(−3) = −1 (4.2.17)

From the Equations (4.2.12) and (4.2.6)

4 4
l31 = = =2 (4.2.18)
u11 2
From the Equations (4.2.13), (4.2.7), (4.2.16) and (4.2.18)
11 − l31 u12 11 − 2(5)
l32 = = =1 (4.2.19)
u22 1
66

From the Equations (4.2.14), (4.2.8), (4.2.19), (4.2.17) and (4.2.18)

u33 = −6 − l31 u13 − l32 u23 = −6 − 2(−3) − 1(−1) = 1 (4.2.20)

Thus,   
1 0 0 2 5 −3
A =  1 1 0   0 1 −1 
2 1 1 0 0 1
The given system of equations can be written as

AX = B
   
x 5
where X =  y  and B =  6 . Therefore,
z 10
      
1 0 0 2 5 −3 x 5
 1 1 0   0 1 −1   y  =  6 
2 1 1 0 0 1 z 10

i.e.,        
1 0 0 2 5 −3 x 5
 1 1 0   0 1 −1   y  =  6 
2 1 1 0 0 1 z 10
Writing     
2 5 −3 x y1
 0 1 −1   y  =  y2  (4.2.21)
0 0 1 z y3
we have,     
1 0 0 y1 5
 1 1 0   y2  =  6 
2 1 1 y3 10
The corresponding system of equations is

y1 = 5 (4.2.22)

y1 + y2 = 6 (4.2.23)
2y1 + y2 + y3 = 10 (4.2.24)
From the Equations (4.2.22) and (4.2.23),

y2 = 6 − y1 = 6 − 5 = 1 (4.2.25)
 67

From the Equations (4.2.22), (4.2.23) and (4.2.24),

y3 = 10 − 2y1 − y2 = 10 − 2(5) − 1 = −1 (4.2.26)

Now, from the Equations (4.2.21), (4.2.22), (4.2.25) and (4.2.26), we have
    
2 5 −3 x 5
 0 1 −1   y  =  1 
0 0 1 z −1

The corresponding system of equations is

2x + 5y − 3z = 5 (4.2.27)

y−z =1 (4.2.28)
z = −1 (4.2.29)
From the Equations (4.2.29) and (4.2.28),

y = 1 + z = 1 + (−1) = 0 (4.2.30)

From the Equations (4.2.29), (4.2.30) and (4.2.27),

5 − 5y + 3z 5 − 5(0) + 3(−1) 2
x= = = =1
2 2 2
Thus, the solution of the given system of equations is x = 1, y = 0 and z = −1

Example 4.2.5. Solve the following system of equations by Crout’s method.

3x − 6y + 3z = 9
x + 3y + z = −2
x − 3y + 3z = 4
 
3 −6 3
Solution. The coefficient matrix of the given system is A =  1 3 1 .
1 −3 3
Let   
l11 0 0 1 u12 u13
A=  l21 l22 0   0 1 u23 
l31 l32 l33 0 0 1
Then    
3 −6 3 l11 l11 u12 l11 u13
 1 3 1  =  l21 l21 u12 + l22 l21 u13 + l22 u23 
1 −3 3 l31 l31 u12 + l32 l31 u13 + l32 u23 + l33
68

Therefore,
l11 = 3 (4.2.31)
l21 = 1 (4.2.32)
l31 = 1 (4.2.33)
l11 u12 = −6 (4.2.34)
l21 u12 + l22 = 3 (4.2.35)
l31 u12 + l32 = −3 (4.2.36)
l11 u13 = 3 (4.2.37)
l21 u13 + l22 u23 = 1 (4.2.38)
l31 u13 + l32 u23 + l33 = 3 (4.2.39)
From the Equations (4.2.31) and (4.2.34)
−6 −6
u12 = = = −2 (4.2.40)
l11 3
From the Equations (4.2.32), (4.2.35) and (4.2.40)

l22 = 3 − l21 u12 = 3 − 1(−2) = 5 (4.2.41)

From the Equations (4.2.33), (4.2.36) and (4.2.40)

l32 = −3 − l31 u12 = −3 − 1(−2) = −1 (4.2.42)

From the Equations (4.2.37) and (4.2.31)

3 3
u13 = = =1 (4.2.43)
l11 3
From the Equations (4.2.38), (4.2.32), (4.2.41) and (4.2.43)

1 − l21 u13 1 − 1(1)

u23 = = =0 (4.2.44)
l22 5
From the Equations (4.2.39), (4.2.33), (4.2.44), (4.2.42) and (4.2.43)

l33 = 3 − l31 u13 − l32 u23 = 3 − 1(1) − 1(0) = 2 (4.2.45)

Thus,   
3 0 0 1 −2 1
A =  1 5 0  0 1 0 
1 −1 2 0 0 1
 69

The given system of equations can be written as

AX = B
   
x 9
where X =  y  and B =  −2 . Therefore,
z 4
      
3 0 0 1 −2 1 x 9
 1 5 0   0 1 0   y  =  −2 
1 −1 2 0 0 1 z 4

i.e.,        
3 0 0 1 −2 1 x 9
 1 5 0   0 1 0   y  =  −2 
1 −1 2 0 0 1 z 4
Writing     
1 −2 1 x y1
 0 1 0   y  =  y2  (4.2.46)
0 0 1 z y3
we have,     
3 0 0 y1 9
 1 5 0   y2  =  −2 
1 −1 2 y3 4
The corresponding system of equations is

3y1 = 9

i.e.,
y1 = 3 (4.2.47)
y1 + 5y2 = −2 (4.2.48)
y1 − y2 + 2y3 = 4 (4.2.49)
From the Equations (4.2.47) and (4.2.48),

−2 − y1 −2 − 3 −5
y2 = = = = −1 (4.2.50)
5 5 5
From the Equations (4.2.47), (4.2.48) and (4.2.49),

4 − y1 + y2 4−3−1 0
y3 = = = =0 (4.2.51)
2 2 2
70

Now, from the Equations (4.2.46), (4.2.47), (4.2.50) and (4.2.51), we have
    
1 −2 1 x 3
 0 1 0   y  =  −1 
0 0 1 z 0

The corresponding system of equations is

x − 2y + z = 3 (4.2.52)

y = −1 (4.2.53)
z=0 (4.2.54)
From the Equations (4.2.54), (4.2.53) and (4.2.52),

x = 3 + 2y − z = 3 + 2(−1) − 0 = 1

Thus, the solution of the given system of equations is x = 1, y = −1 and z = 0

Exercise

1. Solve the following system of linear equations by LU factorization method. Take

the diagonal entries of L as unity.
2x + y + 3z = 2, 4x + 7y + 5z = 8 and 2x − 14y + 9z = −7
             
2 1 3 1 0 0 2 1 3 y1 2 x −1
A =  4 7 5  =  2 1 0   0 5 −1  ,  y2  =  4  ,  y  =  1 
2 −14 9 1 −3 1 0 0 3 y3 3 z 1
 
2. Solve the following system of linear equations by LU factorization method. Take
the diagonal entries of L as unity.
3x + 2y + z = 6, 2x + y + 2z = 2 and 3x + 2y + 2z = 5
             
3 2 1 1 0 0 3 2 1 y1 6 x 1
A =  2 1 2  =  32 1 0   0 − 31 34  ,  y2  =  −2  ,  y  =  2 
3 2 2 1 0 1 0 0 1 y3 −5 z −1
 
3. Solve the following system of linear equations by LU factorization method. Take
the diagonal entries of U as unity.
x + 2y − z = 1, x + 4y − z = −1 and x + y + z = 5
             
1 2 −1 1 0 0 1 2 −1 y1 1 x 0
A=  1 4 −3  =  1 2 0   0 1 −1 , y2
   =  −1 , y
   =  2 
1 1 1 1 −1 1 0 0 1 y3 3 z 3
 
 71

4.2.3 Jacobi’s iteration method

Consider the following system of linear equations.
a11 x + a12 y + a13 z = b1
a21 x + a22 y + a23 z = b2
a31 x + a32 y + a33 z = b3
On solving for x, y and z we have
b1 − a12 y − a13 z
x=
a11
b2 − a21 x − a23 z
y=
a22
b3 − a31 x − a32 y
z=
a33
In Jacobi’s iteration method, we start with an initial approximation x = x(0) , y = y (0) ,
z = z (0) and get the successive approximations to the required solution by using

(n+1) b1 − a12 y (n) − a13 z (n)

x =
a11
b2 − a21 x(n) − a23 z (n)
y (n+1) = (n = 0, 1, 2, · · · )
a22
b3 − a31 x(n) − a32 y (n)
z (n+1) =
a33
The process is repeated till the difference between the successive approximations to the
solution is negligible.
Example 4.2.6. Solve the following system of linear equations by Jacobi’s iteration
method.
10x + y + z = 12
x + 15y − z = 15
x − y + 20z = 20
Take x(0) = 0, y (0) = 0 and z (0) = 0 as an initial approximation and carry out three
iterations.
Solution. We first rewrite the given system of equations as
12 − y − z
x=
10
15 − x + z
y=
15
20 − x + y
z=
20
72

Taking x(0) = 0, y (0) = 0 and z (0) = 0 as an initial approximation, we find the successive
approximations to the required solution by using

12 − y (n) − z (n)
x(n+1) =
10
15 − x(n) + z (n)
y (n+1) = (n = 0, 1, 2, · · · ) (4.2.55)
15
20 − x(n) + y (n)
z (n+1) =
20
Putting n = 0 in the System (4.2.55), we get

12 − y (0) − z (0)
x(1) =
10
12 − 0 − 0
= (because y (0) = z (0) = 0)
10
12
=
10
= 1.2

15 − x(0) + z (0)
y (1) =
15
15 − 0 + 0
= (because x(0) = z (0) = 0)
15
15
=
15
=1
20 − x(0) + y (0)
z (1) =
20
20 − 0 + 0
= (because x(0) = y (0) = 0)
20
20
=
20
=1
Putting n = 1 in the System (4.2.55), we get

12 − y (1) − z (1)
x(2) =
10
12 − 1 − 1
= (because y (1) = z (1) = 1)
10
10
=
10
=1
 73

15 − x(1) + z (1)
y (2) =
15
15 − 1.2 + 1
= (because x(1) = 1.2 and z (1) = 1)
15
14.8
=
15
= 0.9867
20 − x(1) + y (1)
z (2) =
20
20 − 1.2 + 1
= (because x(1) = 1.2 and y (1) = 1)
20
19.8
=
20
= 0.99
Putting n = 2 in the System (4.2.55), we get
12 − y (2) − z (2)
x(3) =
10
12 − 0.9867 − 0.99
= (because y (2) = 0.9867 and z (2) = 0.99)
10
10.0233
=
10
= 1.00233

15 − x(2) + z (2)
y (3) =
15
15 − 1 + 0.99
= (because x(2) = 1 and z (2) = 0.99)
15
14.99
=
15
= 0.9993
20 − x(2) + y (2)
z (3) =
20
20 − 1 + 0.9867
= (because x(2) = 1 and y (2) = 0.9867)
20
19.9867
=
20
= 0.999335
Excercise

1. Solve the following system of equations by Jacobi’s iteration method;

20x − y − z = 15, x + 10y − z = 18 and x − 2y + 15z = 42.
74

Take x(0) = y (0) = z (0) = 0 and carry out three iterations.



The exact solution is x = 1, y = 2 and z = 3 
2. Solve the following system of equations by Jacobi’s iteration method;
10x − y − z = 11, x − 12y + z = 13 and x − 2y + 20z = 3.
Take x(0) = y (0) = z (0) = 0 and carry out three iterations.
 

The exact solution is x = 1, y = −1 and z = 0 

4.2.4 Gauss-Seidel iteration method

Consider the following system of linear equations.
a11 x + a12 y + a13 z = b1
a21 x + a22 y + a23 z = b2
a31 x + a32 y + a33 z = b3
On solving for x, y and z we have
b1 − a12 y − a13 z
x=
a11
b2 − a21 x − a23 z
y=
a22
b3 − a31 x − a32 y
z=
a33
In Gauss-Seidel’s iteration method, we start with an initial approximation x = x(0) ,
y = y (0) , z = z (0) and get the successive approximations to the required solution by using
b1 − a12 y (n) − a13 z (n)
x(n+1) =
a11
b2 − a21 x(n+1) − a23 z (n)
y (n+1) = (n = 0, 1, 2, · · · )
a22
b3 − a31 x(n+1) − a32 y (n+1)
z (n+1) =
a33
The process is repeated till the difference between the successive approximations to the
solution is negligible.
Example 4.2.7. Solve the following system of linear equations by Gauss-seidel’s iteration
method.
20x + y − 2z = 17
3x + 20y − z = −18
2x − 3y + 20z = 25

Take x(0) = y (0) = z (0) = 0 as an initial approximation and carry out three iterations.
 75

Solution. We first rewrite the given system of equations as

17 − y + 2z
x=
20
−18 − 3x + z
y=
20
25 − 2x + 3y
z=
20
Taking x(0) = y (0) = z (0) = 0 as an initial approximation, we find the successive approxi-
mations to the required solution by using
17 − y (n) + 2z (n)
x(n+1) =
20
−18 − 3x(n+1) + z (n)
y (n+1) = (n = 0, 1, 2, · · · ) (4.2.56)
20
(n+1)
25 − 2x + 3y (n+1)
z (n+1) =
20
Putting n = 0 in the System (4.2.56), we get

17 − y (0) + 2z (0)
x(1) =
20
17 − 0 − 0
= (because y (0) = z (0) = 0)
20
= 0.85

−18 − 3x(1) + z (0)

y (1) =
20
−18 − 3(0.85) + 0
= (because x(1) = 0.85 and z (0) = 0)
20
= −1.0275
25 − 2x(1) + 3y (1)
z (1) =
20
25 − 2(0.85) + 3(−1.0275)
= (because x(1) = 0.85 and y (1) = −1.0275)
20
= 1.0109
Putting n = 1 in the System (4.2.56), we get

17 − y (1) + 2z (1)
x(2) =
20
17 − 1.0275 + 2(1.0109)
= (because y (1) = −1.0275 and z (1) = 1.0109)
20
= 1.0025
76

−18 − 3x(2) + z (1)

y (2) =
20
−18 − 3(1.0025) + 1.0109
= (because x(2) = 1.0025 and z (1) = 1.0109)
20
= −0.9998

25 − 2x(2) + 3y (2)
z (2) =
20
25 − 2(1.0025) + 3(−0.9998)
= (because x(2) = 1.0025 and y (2) = −0.9998)
20
= 0.9998
Putting n = 2 in the System (4.2.56), we get

17 − y (2) + 2z (2)
x(3) =
20
17 − (−0.9998) + 2(0.9998)
= (because y (2) = −0.9998 and z (2) = 0.9998)
20
=1

−18 − 3x(3) + z (2)

y (3) =
20
−18 − 3(1) + 0.9998
= (because x(3) = 1 and z (2) = 0.9998)
20
= −1

25 − 2x(3) + 3y (3)
z (3) =
20
25 − 2(1) + 3(−1)
= (because x(3) = 1 and y (3) = −1)
20
=1

Excercise

1. Solve the following system of equations by Gauss-Seidel’s iteration method;

20x − y − z = 19, x + 10y − z = 11 and x + 2y + 15z = 3.
Take x(0) = y (0) = z (0) = 0 and carry out three iterations.


The exact solution is x = 1, y = 1 and z = 0 

2. Solve the following system of equations by Gauss-Seidel’s iteration method;

25x − y − 2z = 1, x + 12y + z = 11 and x − 2y − 20z = 18.
Take x(0) = y (0) = z (0) = 0 and carry out three iterations.
 

The exact solution is x = 0, y = 1 and z = −1 
 77

4.3 Solution of system of non-linear equations

In this section we shall discuss Newton-Raphson’s method for solving system of non-linear
equations.
Consider the two equations
f (x, y) = 0
(4.3.1)
g(x, y) = 0

Let (x0 , y0 ) be an initial approximation to a solution of the System (4.3.1). A better

approximation (x1 , y1 ) is found as follows.
Let x1 = x0 + h and let y1 = y0 + k so that
f (x1 , y1 ) = 0
g(x1 , y1 ) = 0
i.e.,
f (x0 + h, y0 + k) = 0
(4.3.2)
g(x0 + h, y0 + k) = 0

Expanding each of the functions in the System (4.3.2) by Taylor’s series(up to first degree
terms), we get
f (x0 , y0 ) + hfx (x0 , y0 ) + kfy (x0 , y0 ) = 0 (4.3.3)
g(x0 , y0 ) + hgx (x0 , y0 ) + kgy (x0 , y0 ) = 0 (4.3.4)
Solving the Equations (4.3.3) and (4.3.4) for h and k we get a new approximation to the
required solution of the given system as x1 = x0 + h and y1 = y0 + k. This process is
repeated till we get the solution to the desired degree of accuracy.
Note. The Equations (4.3.3)and (4.3.4) are called Newton-Raphson’s equations.
Example 4.3.1. Solve the following system of non-linear equations: x2 +y = 3; y 2 +x = 5.
Taking x0 = y0 = 1 as an initial approximation, find x1 and y1 .
Solution. Let f (x, y) = x2 + y − 3 and let g(x, y) = x + y 2 − 5. Clearly fx = 2x,fy = 1,
gx = 1 and gy = 2y. We know that

f (x0 , y0 ) + hfx (x0 , y0 ) + kfy (x0 , y0 ) = 0

This implies
0 = f (1, 1) + hfx (1, 1) + kfy (1, 1) (as x0 = y0 = 1)
= [12 + 1 − 3] + h[2(1)] + k(1) (as f (x, y) = x2 + y − 3 fx (x, y) = 2x and fy (x, y) = 1)
Thus
2h + k = 1 (4.3.5)
78

g(x0 , y0 ) + hgx (x0 , y0 ) + kgy (x0 , y0 ) = 0

This implies

0 = g(1, 1) + hgx (1, 1) + kgy (1, 1) (as x0 = y0 = 1)

= [1 + 12 − 5] + h(1) + k[2(1)] (as g(x, y) = x + y 2 − 5 gx (x, y) = 1 and gy (x, y) = 2y)

Thus
h + 2k = 3 (4.3.6)
Solving the Equations (4.3.5) and (4.3.6) we get, h = − 31 ≈ −0.3333 and k = 53 ≈ 1.6667.
Therefore, x1 = x0 + h = 1 + (−0.3333) = 0.6667 and y1 = y0 + k = 1 + 1.6667 = 2.6667

Excercise

1. Solve the following system of non-linear equations: x2 + y = 4; y 2 + x = 10. Taking

x0 = y0 = 2 as an initial approximation, find x1 and y1 .

2. Solve the following system of non-linear equations: x2 + y = 7; y 2 + x = 11. Taking

x0 = 3.5 and y0 = −1.8 as an initial approximation, find x1 and y1 .
Chapter 5

Numerical solution of ODE and PDE

Many practical problems often lead to differential equations which can not solved by the
methods that we have learnt so far. In such cases we use numerical methods to solve the
differential equation under consideration.

5.1 Numerical solution of ODE

A number of numerical methods are available to solve the first order differential equation

dy
= f (x, y) given y(x0 ) = y0
dx
We shall now study few of them.

5.1.1 Taylor’s series method

Consider the differential equation

dy
= f (x, y) (5.1.1)
dx
together with the initial condition
y(x0 ) = y0 (5.1.2)
In this method we assume that y(x) can be expanded in powers of x − x0 , giving

(x − x0 )2 00 (x − x0 )3 000 (x − x0 )n (n)
y(x) = y(x0 )+y 0 (x0 )(x−x0 )+ y (x0 )+ y (x0 ) · · ·+ y (x0 )+· · ·
2! 3! n!
Taking h = x − x0 we get

h2 00 h3 hn
y(x0 + h) = y(x0 ) + hy 0 (x0 ) + y (x0 ) + y 000 (x0 ) · · · + y (n) (x0 ) + · · · (5.1.3)
2! 3! n!

79
80

y(x0 ) is known from the Equation (5.1.2) and y 0 (x0 ) is found from the Equation (5.1.1).
y 00 (x0 ), y 000 (x0 ) · · · are obtained by successively differentiating the Equation (5.1.1) and
evaluating at x = x0 . On substituting the values of h, y 0 (x0 ), y 00 (x0 ), · · · in the Equa-
tion (5.1.3), we get the required value of y.
Example 5.1.1. Employ Taylor’s series method to obtain approximate value of y at
dy
x = 0.1 for the differential equation − 2y = 3ex , y(0) = 0
dx
Solution. Given
dy
= 2y + 3ex (5.1.4)
dx
and
y(0) = 0 (5.1.5)
We know that
h2 00 h3 h4
y(x0 + h) = y(x0 ) + hy 0 (x0 ) +
y (x0 ) + y 000 (x0 ) + y (iv) (x0 ) + · · ·
2! 3! 4!
Here x0 = 0 and h = x − x0 = 0.1 − 0 = 0.1. Therefore, we have
(0.1)2 00 (0.1)3 000 (0.1)4 (iv)
y(0.1) = y(0) + (0.1)y 0 (0) + y (0) + y (0) + y (0) + · · · (5.1.6)
2! 3! 4!
From the Equation (5.1.4), we have
y 0 (0) = 2y(0) + 3e0
= 2(0) + 3(1) (as y(0) = 0 and e0 = 1) (5.1.7)
=3
Differentiating the Equation (5.1.4) w.r.t. x, we get
y 00 (x) = 2y 0 (x) + 3ex (5.1.8)
This implies
y 00 (0) = 2y 0 (0) + 3e0
= 2(3) + 3(1) (by the Equation (5.1.7)) (5.1.9)
=9
Differentiating the Equation (5.1.8) w.r.t. x, we get
y 000 (x) = 2y 00 (x) + 3ex (5.1.10)
This implies
y 000 (0) = 2y 00 (0) + 3e0
= 2(9) + 3(1) (by the Equation (5.1.9)) (5.1.11)
= 21
 81

Differentiating the Equation (5.1.10) w.r.t. x, we get

y (iv) (x) = 2y 000 (x) + 3ex (5.1.12)

This implies

y (iv) (0) = 2y 000 (0) + 3e0

= 2(21) + 3(1) (by the Equation (5.1.11)) (5.1.13)
= 45

From the Equations (5.1.5), (5.1.6), (5.1.7), (5.1.9), (5.1.11) and (5.1.13) we get

(0.1)2 (0.1)3 (0.1)4

y(0.1) = 0 + (0.1)3 + 9+ 21 + 45 + · · ·
2! 3! 4!
≈ 0.3 + 0.045 + 0.0035 + 0.0001875
= 0.3486875

Example 5.1.2. Using Taylor’s series method, find an approximate value of y(0.2) given
dy
that = 1 − 2xy, y(0) = 0 and h = 0.1.
dx
Solution. Given
dy
= 1 − 2xy (5.1.14)
dx
and
y(0) = 0 (5.1.15)
We know that
h2 00 h3 h4
y(x0 + h) = y(x0 ) + hy 0 (x0 ) + y (x0 ) + y 000 (x0 ) + y (iv) (x0 ) + · · · (5.1.16)
2! 3! 4!
Here x0 = 0 and h = 0.1. Therefore, we have

(0.1)2 00 (0.1)3 000 (0.1)4 (iv)

y(0.1) = y(0) + (0.1)y 0 (0) + y (0) + y (0) + y (0) + · · · (5.1.17)
2! 3! 4!
From the Equation (5.1.14), we have

y 0 (0) = 1 − 2(0)y(0)
(5.1.18)
=1

Differentiating the Equation (5.1.14) w.r.t. x, we get

y 00 (x) = −2[xy 0 (x) + y(x)] (5.1.19)

82

This implies
y 00 (0) = −2[(0)y 0 (0) + y(0)]
= −2[0(1) + 0] (by the Equations (5.1.15) and (5.1.18)) (5.1.20)
=0
Differentiating the Equation (5.1.19) w.r.t. x, we get
y 000 (x) = −2[xy 00 (x) + 2y 0 (x)] (5.1.21)
This implies
y 000 (0) = −2[(0)y 00 (0) + 2y 0 (0)]
= −2[(0)(0) + 2(1)] (by the Equations (5.1.18) and (5.1.20)) (5.1.22)
= −4
Differentiating the Equation (5.1.21) w.r.t. x, we get
y (iv) (x) = −2[xy 000 (x) + 3y 00 (x)] (5.1.23)
This implies
y (iv) (0) = −2[(0)y 000 (0) + 3y 00 (0)]
= −2[(0)(−4) + 3(0)] (by the Equations (5.1.20) and (5.1.22)) (5.1.24)
=0
From the Equations (5.1.15), (5.1.17), (5.1.18), (5.1.20), (5.1.22) and (5.1.24) we get
(0.1)2 (0.1)3 (0.1)4
y(0.1) = 0 + (0.1)1 + (0) + (−4) + (0) + · · ·
2! 3! 4!
(5.1.25)
≈ 0.1 − 0.0007
= 0.0993
To find y(0.2), we replace x0 by x1 := x0 + h in the Equation (5.1.16).
h2 00 h3 h4
y(x1 + h) = y(x1 ) + hy 0 (x1 ) + y (x1 ) + y 000 (x1 ) + y (iv) (x1 ) + · · ·
2! 3! 4!
Here h = 0.1 and x1 = x0 + h = 0 + 0.1 = 0.1. Therefore, we have
(0.1)2 00 (0.1)3 000 (0.1)4 (iv)
y(0.2) = y(0.1) + (0.1)y 0 (0.1) + y (0.1) + y (0.1) + y (0.1) + · · ·
2! 3! 4!
(5.1.26)
From the Equation (5.1.14), we have
y 0 (0.1) = 1 − 2(0.1)y(0.1)
= 1 − 2(0.1)(0.0993) (from the Equation (5.1.25)) (5.1.27)
= 0.9801
 83

From the Equation (5.1.19), we have

y 00 (0.1) = −2[(0.1)y 0 (0.1) + y(0.1)]

= −2[(0.1)(0.9801) + 0.0993] (by the Equations (5.1.25) and (5.1.27)) (5.1.28)
= −0.3946

From the Equation (5.1.21), we have

y 000 (0.1) = −2[(0.1)y 00 (0.1) + 2y 0 (0.1)]

= −2[(0.1)(−0.3946) + 2(0.9801)] (by the Equations (5.1.27) and (5.1.28))
= −3.8415
(5.1.29)

From the Equation (5.1.23), we have

y (iv) (0) = −2[(0.1)y 000 (0.1) + 3y 00 (0.1)]

= −2[(0.1)(−3.8415) + 3(−0.3946)] (by the Equations (5.1.28) and (5.1.29))
= 3.1359
(5.1.30)

From the Equations (5.1.26), (5.1.27), (5.1.28), (5.1.29) and (5.1.30) we get

(0.1)2 (0.1)3 (0.1)4

y(0.2) = 0.0993 + (0.1)(0.9801) + (−0.3946) + (−3.8415) + (3.1359) + · · ·
2! 3! 4!
≈ 0.09933 + 0.09801 − 0.001973 − 0.00064025 + 0.00001306625
= 0.1947

Exercise

1. Using Taylor’s series method, find an approximate value of y(0.1) given that
y 0 = x − y 2 , y(0) = 1 and h = 0.1.

2. Using Taylor’s series method, find an approximate value of y(0.2) given that
dy
= x2 + y 2 , y(0) = 1 and h = 0.1.
dx

5.1.2 Modified Euler’s method

Let y = f (x) be a differentiable curve. The curve y = f (x) (x0 −  ≤ x ≤ x0 + ) can
be approximated by the tangent line about the point x0 . Using this, we can find an
approximate value of y at any point near x0 .
Consider the differential equation
dy
= f (x, y) (5.1.31)
dx
84

Figure 5.1: Approximation of a curve by the tangent

together with the initial condition

y(x0 ) = y0 (5.1.32)
Let x1 = x0 +h. An approximate value of y at x1 can be found using the tangent at x = x0
to the solution curve of the given differential equation. The equation of the tangent at
x = x0 to the solution curve of the given differential equation is given by(see Figure 5.1)

y − y0 = m(x − x0 ) (5.1.33)

where m is the slope of the tangent at x = x0 . If y1 = y(x1 ) then from the Equa-
tion (5.1.33), we get

y1 = m(x1 − x0 ) + y0
= mh + y0 (as x1 − x0 = h)

Note that from the Equations (5.1.31) and (5.1.32), m = f (x0 , y0 ). Therefore, y1 =
y0 + hf (x0 , y0 ). Similarly, we can show that y2 = y1 + hf (x1 , y1 ), y3 = y2 + hf (x2 , y2 ) and,
in general,
yn+1 = yn + hf (xn , yn ) (5.1.34)
where yn+1 = y(xn+1 ) and xn+1 = xn + h (n ≥ 0). When the step size h is quite large
then by using the Equation (5.1.34) we get a very crude value of y. To minimize the error
we use
h ∗
yn+1 = yn + [f (xn , yn ) + f (xn+1 , yn+1 )] (5.1.35)
2
 85

where
∗
yn+1 = yn + hf (xn , yn ) (5.1.36)
Example 5.1.3. Using Modified Euler’s method find an approximate value of y(0.1)
given that y 0 = −xy 2 , h = 0.1 and y(0) = 2.
Solution. We know that
h ∗
yn+1 = yn + [f (xn , yn ) + f (xn+1 , yn+1 )] (n ≥ 0) (5.1.37)
2
where
∗
yn+1 = yn + hf (xn , yn ) (n ≥ 0) (5.1.38)
Here f (x, y) = −xy 2 , h = 0.1, x0 = 0 and y0 = 2. Now, putting n = 0 in the Equa-
tion (5.1.38) we get

y1∗ = y0 + hf (x0 , y0 )
= y0 + h[−x0 (y0 )2 ]
(5.1.39)
= 2 + (0.1)[−0(22 )]
=2

Putting n = 0 in the Equation (5.1.37), we get

(1) h
y1 = y0 + [f (x0 , y0 ) + f (x1 , y1∗ )]
2
h
= y0 + [−x0 (y0 )2 − x1 (y1∗ )2 ]
2
0.1 (5.1.40)
=2+ [−0(22 ) − (0.1)(22 )]
2
= 2 − 0.02
= 1.98

(2) h (1)
y1 = y0 + [f (x0 , y0 ) + f (x1 , y1 )]
2
h (1)
= y0 + [−x0 (y0 )2 − x1 (y1 )2 ]
2
0.1 (5.1.41)
=2+ [−0(22 ) − (0.1)((1.98)2 )]
2
= 2 − 0.019602
= 1.980398
Therefore, the value of y at x = 0.1 correct to three decimal places is 1.98
Example 5.1.4. Using Modified Euler’s method find an approximate value of y(0.1)
given that y 0 = x2 + y, h = 0.05 and y(0) = 1. Carry out two iterations at each step.
86

Solution. We know that

h ∗
yn+1 = yn + [f (xn , yn ) + f (xn+1 , yn+1 )] (n ≥ 0) (5.1.42)
2

where
∗
yn+1 = yn + hf (xn , yn ) (n ≥ 0) (5.1.43)

Here f (x, y) = x2 + y, h = 0.05, x0 = 0 and y0 = 1. Now, putting n = 0 in the

Equation (5.1.43) we get

y1∗ = y0 + hf (x0 , y0 )
= y0 + h[(x0 )2 + y0 ]
(5.1.44)
= 1 + (0.05)[02 + 1]
= 1.05

Putting n = 0 in the Equation (5.1.42), we get

(1) h
y1 = y0 + [f (x0 , y0 ) + f (x1 , y1∗ )]
2
h
= y0 + [((x0 )2 + y0 ) + ((x1 )2 + y1∗ )]
2
0.05 2 (5.1.45)
=1+ [(0 + 1) + ((0.05)2 + 1.05)]
2
= 1 + 0.051313
= 1.051313

(2) h (1)
y1 = y0 + [f (x0 , y0 ) + f (x1 , y1 )]
2
h (1)
= y0 + [((x0 )2 + y0 ) + ((x1 )2 + y1 )]
2
0.05 2 (5.1.46)
=1+ [(0 + 1) + ((0.05)2 + 1.051313)]
2
= 1 + 0.051345
= 1.051345

Therefore, y(0.05) = 1.0513. To find y(0.1), we put n = 1 in the Equation (5.1.43).

y2∗ = y1 + hf (x1 , y1 )
= y1 + h[(x1 )2 + y1 ]
(5.1.47)
= 1.0513 + (0.05)[(0.05)2 + 1.0513]
= 1.10399
 87

Putting n = 1 in the Equation (5.1.42), we get

(1) h
y2 = y1 + [f (x1 , y1 ) + f (x2 , y2∗ )]
2
h
= y1 + [((x1 )2 + y1 ) + ((x2 )2 + y2∗ )]
2
0.05 (5.1.48)
= 1.0513 + [((0.05)2 + 1.0513) + ((0.1)2 + 1.10399)]
2
= 1.0513 + 0.05419475
= 1.10549475

(2) h (1)
y2 = y1 + [f (x1 , y1 ) + f (x2 , y2 )]
2
h (1)
= y1 + [((x1 )2 + y1 ) + ((x2 )2 + y2 )]
2
0.05 (5.1.49)
= 1.0513 + [((0.05)2 + 1.0513) + ((0.1)2 + 1.10549475)]
2
= 1.0513 + 0.05423
= 1.10553
Therefore, the value of y at x = 0.1 correct to four decimal places is 1.1055

Exercise

1. Using Modified Euler’s method find an approximate value of y(0.1) given that y 0 =
x + y, h = 0.1 and y(0) = 1.

2. Using Modified Euler’s method find an approximate value of y(1.5) given that y 0 −
√
xy = 2, h = 0.25 and y(1) = 1. Carry out two iterations at each step.

5.1.3 Runge - Kutta method of order four

The Ruge - Kutta method of order four is of great practical importance. This method
gives us more accurate solution than the Modified Euler’s method.
Consider the differential equation
dy
= f (x, y)
dx
together with the initial condition
y(x0 ) = y0
In this method, to find an approximate value of y at x = a, we first divide the interval
[x0 , a] into n subintervals of width h and then for each i = 0, 1, · · · , n − 1 we calculate

yi+1 = yi + k
88

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

k1 = hf (xi , yi )
 
h k1
k2 = hf x i + , yi +
2 2
 
h k2
k3 = hf x i + , yi +
2 2
and
k4 = hf (xi + h, yi + k3 )

Example 5.1.5. Apply Runge - Kutta method of order four, to find an approximate
value of y(0.1) given that y 0 = x + y 2 , h = 0.1 and y(0) = 1.

Solution. Given y 0 = x + y 2 , h = 0.1 and y(0) = 1. Here, f (x, y) = x + y 2 , x0 = 0 and

y0 = 1. We know that
k1 = hf (xi , yi )
 
h k1
k2 = hf x i + , yi +
2 2
 
h k2
k3 = hf x i + , yi +
2 2
and
k4 = hf (xi + h, yi + k3 )

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

Putting i = 0 we get

k1 = hf (x0 , y0 )
= h(x0 + y02 )
= (0.1)(0 + 1)
= 0.1
 89
 
h k1
k2 = hf x0 + , y0 +
2 2
   2 !
h k1
= h x0 + + y0 +
2 2
   2 !
0.1 0.1
= (0.1) 0 + + 1+
2 2
= (0.1) (0.05 + 1.1025)
≈ 0.1153
 
h k2
k3 = hf x0 + , y0 +
2 2
   2 !
h k2
= h x0 + + y0 +
2 2
   2 !
0.1 0.1153
= (0.1) 0 + + 1+
2 2
≈ (0.1) (0.05 + 1.1186)
≈ 0.1169
k4 = hf (x0 + h, y0 + k3 )
= (0.1) [x0 + h] + [y0 + k3 ]2

= (0.1) [0 + 0.1] + [1 + 0.1169]2

≈ (0.1) (0.1 + 1.2475)

≈ 0.1348
1
k = (k1 + 2k2 + 2k3 + k4 )
6
1
= (0.1 + 2(0.1153) + 2(0.1169) + 0.1348)
6
= 0.1165
and
y1 = y0 + k
= 1 + 0.1165
= 1.1165
 
y0 = y(x0 ), y1 = y(x1 ) = y(x0 + h), and here x0 = 0 and h = 0.1. Thus, y1 = y(0 + 0.1) = y(0.1)
 
Example 5.1.6. Apply Runge - Kutta method of order four, to find an approximate
value of y(0.4) given that y 0 = x + y, h = 0.2 and y(0) = 1.
90

Solution. Given y 0 = x + y, h = 0.2 and y(0) = 1. Here, f (x, y) = x + y, x0 = 0 and

y0 = 1. We know that
k1 = hf (xi , yi )
 
h k1
k2 = hf xi + , yi +
2 2
 
h k2
k3 = hf xi + , yi +
2 2
and
k4 = hf (xi + h, yi + k3 )
1
k = (k1 + 2k2 + 2k3 + k4 )
6
and
yi+1 = yi + k
Putting i = 0 we get

k1 = hf (x0 , y0 )
= h(x0 + y0 )
= (0.2)(0 + 1)
= 0.2
 
h k1
k2 = hf x0 + , y0 +
2 2
   
h k1
= h x0 + + y0 +
2 2
   
0.2 0.2
= (0.2) 0 + + 1+
2 2
= (0.2) (0.1 + 1.1)
= 0.24
 
h k2
k3 = hf x0 + , y0 +
2 2
   
h k2
= h x0 + + y0 +
2 2
   
0.2 0.24
= (0.2) 0 + + 1+
2 2
= (0.2) (0.1 + 1.12)
= 0.244
 91

k4 = hf (x0 + h, y0 + k3 )
= h ([x0 + h] + [y0 + k3 ])
= (0.2) ([0 + 0.2] + [1 + 0.244])
= (0.2) (0.2 + 1.244)
= 0.2888

1
k = (k1 + 2k2 + 2k3 + k4 )
6
1
= (0.2 + 2(0.24) + 2(0.244) + 0.2888)
6
= 0.2428

and

y1 = y0 + k
= 1 + 0.2428
= 1.2428

 
y0 = y(x0 ), y1 = y(x1 ) = y(x0 + h), and here x0 = 0 and h = 0.2. Thus, y1 = y(0 + 0.2) = y(0.2)

  
What is asked?...y(0.4) =?
 
To find y(0.4), put i = 1.

k1 = hf (x1 , y1 )
= h(x1 + y1 )
= (0.2)(0.2 + 1.2428)
≈ 0.2886

 
h k1
k2 = hf x1 + , y1 +
2 2
   
h k1
= h x1 + + y1 +
2 2
   
0.2 0.2886
= (0.2) 0.2 + + 1.2428 +
2 2
= (0.2) (0.3 + 1.3871)
≈ 0.3374
92
 
h k2
k3 = hf x1 + , y1 +
2 2
   
h k2
= h x1 + + y1 +
2 2
   
0.2 0.3374
= (0.2) 0.2 + + 1.2428 +
2 2
= (0.2) (0.3 + 1.4115)
= 0.3423

k4 = hf (x1 + h, y1 + k3 )
= h ([x1 + h] + [y1 + k3 ])
= (0.2) ([0.2 + 0.2] + [1.2428 + 0.3423])
= (0.2) (0.4 + 1.5851)
≈ 0.397

1
k = (k1 + 2k2 + 2k3 + k4 )
6
1
= (0.2886 + 2(0.3374) + 2(0.3423) + 0.397)
6
≈ 0.3408

and

y2 = y1 + k
= 1.2428 + 0.3408
= 1.5836
 
y2 = y(x2 ) = y(x1 + h), and here x1 = x0 + h, x0 = 0 and h = 0.2. Thus, x1 = 0.2

  
and y2 = y(0.2 + 0.2) = y(0.4)
 

Exercise

1. Apply Runge - Kutta method of order four, to find an approximate value of y(0.2)
given that y 0 = 3x + 12 y, h = 0.1 and y(0) = 1.

2. Apply Runge - Kutta method of order four, to find an approximate value of y(0.2)
given that 10y 0 = x2 + y 2 , h = 0.1 and y(0) = 1.

 
1. y1 = y(0.1) = 1.0666 and y2 = y(0.2) = 1.1673 2. y1 = y(0.1) = 1.0101 and y2 = y(0.2) = 1.0206
 
 93

5.1.4 Milne’s method

Consider
dy
= f (x, y)
dx
together with the initial condition
y(x0 ) = y0
In Milne’s method, to find an approximate value of y at x = a , we first divide the interval
[x0 , a] into four subintervals of width h and then compute y1 := y(x0 +h), y2 := y(x0 +2h)
and y3 := y(x0 + 3h) using suitable methods. Next we find f0 := f (x0 , y0 ), f1 := f (x1 , y1 ),
f2 := f (x2 , y2 ) and f3 := f (x3 , y3 ). Now a crude value of y at x = a is got by using

(p) 4h
y4 = y0 + (2f1 − f2 + 2f3 )
3
which is called a predictor. A better value of y4 is found using

(c) h
y4 = y2 + (f2 + 4f3 + f4 )
3
(p)
which is called a corrector where f4 = f (x4 , y4 ). Using this new value of y4 , f4 is
computed and applying corrector formula again we get an improved value of y4 . We
repeat this until the absolute difference between the successive values of y4 is negligible.

Example 5.1.7. Using Milne’s method find an approximate value of y(0.4) given that
y 0 = e−x − y and

x 0 0.1 0.2 0.3

y 0 0.0905 0.1637 0.2222

Use corrector formula twice.

Solution. Here f (x, y) = e−x − y, h = 0.1, x0 = 0, x1 = 0.1, x2 = 0.2, x3 = 0.3, y0 = 0,

y1 = 0.0905, y2 = 0.1637 and y3 = 0.2222. Now we compute

f0 = f (x0 , y0 )
= e−x0 − y0
= e−0 − 0
=1

f1 = f (x1 , y1 )
= e−x1 − y1
= e−0.1 − 0.0905
≈ 0.8143
94

f2 = f (x2 , y2 )
= e−x2 − y2
= e−0.2 − 0.1637
≈ 0.6550
f3 = f (x3 , y3 )
= e−x3 − y3
= e−0.3 − 0.2222
≈ 0.5186
Now using the predictor formula we have

(p) 4h
y4 = y0 + (2f1 − f2 + 2f3 )
3
4(0.1)
=0+ (2(0.8143) − 0.6550 + 2(0.5186))
3
= 0.2681
(p)
f4 = f (x4 , y4 )
(p)
= e−x4 − y4
= e−0.4 − 0.2681
≈ 0.4022
Now using corrector formula we have

(c) h
y4 = y2 + (f2 + 4f3 + f4 )
3
0.1
= 0.1637 + (0.655 + 4(0.5186) + 0.4022)
3
≈ 0.268087
(c)
f4 = f (x4 , y4 )
(c)
= e−x4 − y4
= e−0.4 − 0.268087
≈ 0.40223
Now using corrector formula again we have

(c) h
y4 = y2 + (f2 + 4f3 + f4 )
3
0.1
= 0.1637 + (0.655 + 4(0.5186) + 0.40223)
3
≈ 0.268088
 95

Exercise

1. Using Milne’s method find an approximate value of y(0.8) given that y 0 = x − y 2

and

x 0 0.2 0.4 0.6

y 0 0.02 0.0795 0.1762

Use corrector formula twice.

2. Using Milne’s method find an approximate value of y(0.8) given that y 0 = 1 + y 2

and

x 0 0.2 0.4 0.6

y 0 0.2027 0.4228 0.6841

Use corrector formula twice.

5.2 Numerical solution of PDE

The partial differential equations arise quite naturally in science and engineering. In most
of the cases we use numerical methods to solve PDE, as only few of the PDE’s can be
solved by analytical methods. The method of finite differences is a most commonly used
numerical method to solve the given PDE. In this method, the derivatives in the given
PDE and the boundary conditions are replaced by their finite difference approximations.
As a result, the given PDE gets converted to the difference equation which is solved by
the iterative methods.

5.2.1 Classification of second order linear PDE

The general linear PDE of second order in two independent variables is of the form

∂2u ∂2u ∂2u

A + B + C +D =0 (5.2.1)
∂x2 ∂x ∂y ∂y 2
∂u ∂u
where A, B, C are functions of x and y and D is a function of x, y, u, and . The
∂x ∂y
PDE (5.2.1) is said to be

• elliptic if B 2 − 4AC < 0

• parabolic if B 2 − 4AC = 0

• hyperbolic if B 2 − 4AC > 0

96

Example 5.2.1. Classify the following PDE’s:

∂2u ∂2u ∂ 2 u ∂u ∂u
1. 2
+ 4 + 4 2
− +2 =0
∂x ∂x ∂y ∂y ∂x ∂y
2 2
2∂
u 2 ∂ u
2. x + (1 − y ) 2 = 0 (x 6= 0, −1 < y < 1)
∂x2 ∂y
∂2u ∂2u 2
2 ∂ u
3. (1 + x2 ) + (5 + 2x 2
) + (4 + x ) =0
∂x2 ∂x ∂y ∂y 2
Solution. The general linear PDE of second order in two independent variables is of the
form
∂2u ∂2u ∂2u
A 2 +B +C 2 +D =0 (5.2.2)
∂x ∂x ∂y ∂y
1. Given
∂2u ∂2u ∂ 2 u ∂u ∂u
2
+ 4 + 4 2
− +2 =0 (5.2.3)
∂x ∂x ∂y ∂y ∂x ∂y
On comparing the given equation with (5.2.2) we get A = 1, B = 4 and C = 4.
Therefore, B 2 − 4AC = 0 and hence the given PDE is parabolic.

2. Given
∂2u 2
2 ∂ u
x2 + (1 − y ) =0 (5.2.4)
∂x2 ∂y 2
(x 6= 0, −1 < y < 1). On comparing the given equation with (5.2.2) we get A = x2 ,
B = 0 and C = 1 − y 2 . Therefore,

B 2 − 4AC = −4x2 (1 − y 2 )
= 4x2 (y 2 − 1)

Since x 6= 0, x2 > 0 and as −1 < y < 1, y 2 − 1 < 0. Thus, B 2 − 4AC < 0 and hence
the given PDE is elliptic.

3. Given
∂2u ∂2u 2
2 ∂ u
(1 + x2 ) + (5 + 2x 2
) + (4 + x ) =0 (5.2.5)
∂x2 ∂x ∂y ∂y 2
On comparing the given equation with (5.2.2) we get A = 1 + x2 , B = 5 + 2x2 and
C = 4 + x2 . Therefore,

B 2 − 4AC = (5 + 2x2 )2 − 4(1 + x2 )(4 + x2 )

= 25 + 20x2 + 4x4 − 4(4 + 5x2 + x4 )
=9

Thus, B 2 − 4AC > 0 and hence the given PDE is hyperbolic.

 97

Figure 5.2: The region divided into a network of square meshes

Exercise

Classify the following PDE’s:

∂2u ∂2u ∂2u ∂u ∂u
1. 2
+ 6 + 9 2
+2 − =0
∂x ∂x ∂y ∂y ∂x ∂y
∂2u 2
2∂ u
2. (1 − x2 ) + y = 0 (y 6= 0, −1 < x < 1)
∂x2 ∂y 2
∂2u ∂2u ∂2u
3. x(1 + x) 2 + x(5 + 2x) + x(4 + x) 2 = 0 (x 6= 0)
∂x ∂x ∂y ∂y

5.2.2 Solution of Laplace’s equation

Consider
∂2u ∂2u
+ =0 (5.2.6)
∂x2 ∂y 2
98

Let R be a rectangular region for which u(x, y) is known at the boundary. Divide this
region into a network of square meshes of side h as shown in the Figure 5.2. Note that
∂2u
the finite difference approximation for is given by
∂x2
∂2u u(x − h, y) − 2u(x, y) + u(x + h, y)
2
≈
∂x h2
Writing u(x, y) = u(ih, jh) as simply ui,j we get

∂2u ui−1,j − 2ui,j + ui+1,j

≈ (5.2.7)
∂x2 h2
∂2u
The finite difference approximation for is given by
∂y 2

∂2u u(x, y − h) − 2u(x, y) + u(x, y + h)

2
≈
∂y h2
i.e.,
∂2u ui,j−1 − 2ui,j + ui,j+1
2
≈ (5.2.8)
∂y h2
Replacing the derivatives in the equation (5.2.6) by their finite difference approxima-
tions(as given in equations (5.2.7) and (5.2.8) ) we get
ui−1,j − 2ui,j + ui+1,j ui,j−1 − 2ui,j + ui,j+1
+ =0
h2 h2
⇒ ui−1,j − 2ui,j + ui+1,j + ui,j−1 − 2ui,j + ui,j+1 = 0
1
⇒ ui,j = [ui−1,j + ui+1,j + ui,j−1 + ui,j+1 ]
4
This is known as the standard five point formula.
∂2u ∂2u
Example 5.2.2. Using the standard five point formula find the solution of 2 + 2 = 0
∂x ∂y
at the interior nodes of the square region shown in the Figure 5.4, taking the values
provided at the boundary.

Solution. According to the standard five point formula

1
ui,j = [ui−1,j + ui+1,j + ui,j−1 + ui,j+1 ]
4
If u1 , u2 , u3 and u4 are the values of u at the interior mesh points then
1 1
u1 = [2 + 4 + u2 + u3 ] = [6 + u2 + u3 ] (5.2.9)
4 4
 99

Figure 5.3: The schematic form of standard 5 point formula

Figure 5.4: The square mesh in Example 5.2.2

100

Figure 5.5: The square mesh in Example 5.2.3

1 1
u2 = [5 + 5 + u1 + u4 ] = [10 + u1 + u4 ] (5.2.10)
4 4
1 1
u3 = [1 + 1 + u1 + u4 ] = [2 + u1 + u4 ] (5.2.11)
4 4
1 1
u4 = [2 + 4 + u2 + u3 ] = [6 + u2 + u3 ] (5.2.12)
4 4
From equations (5.2.9) and (5.2.12),

u1 = u4 (5.2.13)

From equations (5.2.10) and (5.2.13) we have

1
u2 = [10 + 2u1 ] (5.2.14)
4
From equations (5.2.11) and (5.2.13) we have
1
u3 = [2 + 2u1 ] (5.2.15)
4
On solving equations (5.2.9), (5.2.14) and (5.2.15) we get u1 = 3, u2 = 4 and u3 = 2.
From equation (5.2.13) we have u4 = u1 = 3.

Example 5.2.3. Solve uxx + uyy = 0 for the square mesh with the boundary values as
shown in the Figure 5.5.
 101

Solution. According to the standard five point formula

1
ui,j = [ui−1,j + ui+1,j + ui,j−1 + ui,j+1 ]
4
If u1 , u2 , u3 and u4 are the values of u at the interior mesh points then
1 1
u1 = [0 + 0 + u2 + u3 ] = [u2 + u3 ] (5.2.16)
4 4
1 1
u2 = [2 + 2 + u1 + u4 ] = [4 + u1 + u4 ] (5.2.17)
4 4
1 1
u3 = [2 + 2 + u1 + u4 ] = [4 + u1 + u4 ] (5.2.18)
4 4
1 1
u4 = [3 + 5 + u2 + u3 ] = [8 + u2 + u3 ] (5.2.19)
4 4
From equations (5.2.17) and (5.2.18),

u2 = u3 (5.2.20)

From equations (5.2.16) and (5.2.20) we have

1
u1 = [2u2 ] (5.2.21)
4
From equations (5.2.19) and (5.2.20) we have
1
u4 = [8 + 2u2 ] (5.2.22)
4
On solving equations (5.2.17), (5.2.21) and (5.2.22) we get u1 = 1, u2 = 2 and u4 = 3.
From equation (5.2.20) we have u3 = u2 = 2.

5.2.3 Solution of Poisson’s equation

Consider
∂2u ∂2u
+ = f (x, y) (5.2.23)
∂x2 ∂y 2
Dividing the region R into a network of square meshes with spacing h and replacing the
partial derivatives in (5.2.23), by their finite difference approximations we get
1
ui,j = [ui−1,j + ui+1,j + ui,j−1 + ui,j+1 − h2 fi,j ] (5.2.24)
4
Example 5.2.4. Solve the equation ∇2 u = −10(x2 + y 2 + 10) over the square mesh with
sides x = 0 = y, x = 3 = y with u = 0 at the boundary and mesh length h = 1.
102

Figure 5.6: The square mesh in Example 5.2.4

Solution. Let u1 , u2 , u3 and u4 be the values of u at the interior mesh points (1, 1), (1, 2), (2, 1)
and (2, 2) respectively(see Figure 5.6). We know that
1
ui,j = [ui−1,j + ui+1,j + ui,j−1 + ui,j+1 − h2 fi,j ] (5.2.25)
4
Therefore,
1
u1 = [0 + u2 + 0 + u3 − (12 )(−120)]
4
(as h = 1, f (x, y) = −10(x2 + y 2 + 10) and f1,1 = f (1, 1) = −120)
i.e.,
1
u1 = [u2 + u3 + 120] (5.2.26)
4
1
u2 = [u1 + 0 + 0 + u4 − (12 )(−150)]
4
(as h = 1, f (x, y) = −10(x2 + y 2 + 10) and f2,1 = f (2, 1) = −150)
i.e.,
1
u2 = [u1 + u4 + 150] (5.2.27)
4
1
u3 = [0 + u4 + u1 + 0 − (12 )(−150)]
4
(as h = 1, f (x, y) = −10(x2 + y 2 + 10) and f1,2 = f (1, 2) = −150)
 103

i.e.,
1
u3 = [u1 + u4 + 150] (5.2.28)
4
1
u4 = [u3 + 0 + u2 + 0 − (12 )(−180)]
4
(as h = 1, f (x, y) = −10(x2 + y 2 + 10) and f2,2 = f (2, 2) = −180)
i.e.,
1
u4 = [u2 + u3 + 180] (5.2.29)
4
From equations (5.2.27)and (5.2.28) we get

u2 = u3 (5.2.30)

From equations (5.2.26) and (5.2.30)

1
u1 = [2u2 + 120] (5.2.31)
4
From equations (5.2.29) and (5.2.30)
1
u4 = [2u2 + 180] (5.2.32)
4
On solving equations (5.2.27), (5.2.31) and (5.2.32) we get u1 = 67.5, u4 = 82.5 and
u2 = 75 = u3 .

Exercise

Solve the PDE ∇2 u = −81xy, 0 < x, y < 1, h = 1/3, u(0, y) = 0 = u(x, 0) and
u(1, y) = 100 = u(x, 1)

5.2.4 Solution of heat equation

The one dimensional heat equation is given by
∂u ∂2u
= c2 2 (5.2.33)
∂t ∂x
Consider a rectangular mesh in the x − t plane with spacing h along x-direction and
k along the time t-direction. Replacing the partial derivatives by their finite difference
approximations we get

ui,j+1 = αui−1,j + (1 − 2α)ui,j + αui+1,j (5.2.34)

kc2 1
where α = 2
. This is known as Schmidt explicit formula and is valid only for 0 < α ≤ .
h 2
The schematic form of the equation (5.2.34) is as shown in the Figure 5.7.
104

Figure 5.7: The schematic form of the equation (5.2.34)

Note. If α = 0.5 then the equation (5.2.34) reduces to

1
ui,j+1 = [ui−1,j + ui+1,j ] (5.2.35)
2
This relation is known as Bendre-Schmidt recurrence relation.

Example 5.2.5. Use Schmidt explicit method to obtain numerical solution of ut = uxx ,
u(0, t) = u(1, t) = 0, u(x, 0) = sin πx with h = 1/3, k = 1/36 for two time levels.
kc2 9 1
Solution. Here c2 = 1, h = 1/3 and k = 1/36. Therefore, α = 2 = = . The
h 36 4
Schmidt explicit formula for the heat equation is given by

ui,j+1 = αui−1,j + (1 − 2α)ui,j + αui+1,j

1
Substituting α = , we have
4
  
1 1 1
ui,j+1 = ui−1,j + 1 − 2 ui,j + ui+1,j
4 4 4
1 1 1
= ui−1,j + ui,j + ui+1,j
4 2 4
Therefore,
1
ui,j+1 = [ui−1,j + 2ui,j + ui+1,j ] (5.2.36)
4
Since u(x, 0) = sin πx we get
√
1 3
u( , 0) = sin(π/3) = ≈ 0.866
3 2
√
2 3
u( , 0) = sin(2π/3) = ≈ 0.866
3 2
 105

Figure 5.8: The rectangular mesh in Example 5.2.5

106

j \i 0 1 2 3
0 0 0.866 0.866 0
1 0 0.64952 0.64952 0
2 0 0.48714 0.48714 0

Table 5.1: Values of u in Example 5.2.5 at the boundary and interior mesh points
j \i 0 1 2 3 4 5
0 0 0.5878 0.9511 0.9511 0.5878 0
1 0 0.47555 0.76945 0.76945 0.47555 0
2 0 0.384725 0.6225 0.6225 0.384725 0
3 0 0.31125 0.5036125 0.5036125 0.31125 0
4 0 0.25180625 0.40743125 0.40743125 0.25180625 0

Table 5.2: Values of u in Example 5.2.6 at the boundary and interior mesh points

Thus, by the equation (5.2.36), we get

√ ! √
1 3 √ 3 3
u1 = + 3 = ≈ 0.64952
4 2 8
√ ! √
1 3 √ 3 3
u2 = + 3 = ≈ 0.64952
4 2 8
√ √ ! √
1 6 3 3 3 9 3
u3 = + = ≈ 0.48714
4 8 8 32
√ √ ! √
1 3 3 6 3 9 3
u4 = + = ≈ 0.48714
4 8 8 32

Example 5.2.6. Solve the boundary value problem ut = uxx under the conditions
u(0, t) = u(1, t) = 0, u(x, 0) = sin πx (0 ≤ x ≤ 1) with h = 0.2 and k = 0.02 for
four time levels.
kc2 0.02 1
Solution. Here c2 = 1, h = 0.2 and k = 0.02. Therefore, α = 2 = = . The
h 0.04 2
Schmidt explicit formula for the heat equation is given by
ui,j+1 = αui−1,j + (1 − 2α)ui,j + αui+1,j
1
Substituting α = , we have
2
  
1 1 1
ui,j+1 = ui−1,j + 1 − 2 ui,j + ui+1,j
2 2 2
 107

Figure 5.9: The rectangular mesh in Example 5.2.6

108

Therefore,
1
ui,j+1 = [ui−1,j + ui+1,j ] (5.2.37)
2
Since u(x, 0) = sin πx we get
u(0.2, 0) = sin(0.2)π ≈ 0.5878
u(0.4, 0) = sin(0.4)π ≈ 0.9511
u(0.6, 0) = sin(0.6)π ≈ 0.9511
u(0.8, 0) = sin(0.8)π ≈ 0.5878
Thus, by the equation (5.2.37), we get
1
u1 = (0 + 0.9511) = 0.47555
2
1
u2 = (0.5878 + 0.9511) = 0.76945
2
1
u3 = (0.5878 + 0.9511) = 0.76945
2
1
u4 = (0 + 0.9511) = 0.47555
2
1 1
u5 = (0 + u2 ) = (0 + 0.76945) = 0.384725
2 2
1 1
u6 = (u1 + u3 ) = (0.47555 + 0.76945) = 0.6225
2 2
1 1
u7 = (u2 + u4 ) = (0.76945 + 0.47555) = 0.6225
2 2
1 1
u8 = (u3 + 0) = (0.76945 + 0) = 0.384725
2 2
1 1
u9 = (0 + u6 ) = (0 + 0.6225) = 0.31125
2 2
1 1
u10 = (u5 + u7 ) = (0.384725 + 0.6225) = 0.5036125
2 2
1 1
u11 = (u6 + u8 ) = (0.6225 + 0.384725) = 0.5036125
2 2
1 1
u12 = (u7 + 0) = (0.6225 + 0) = 0.31125
2 2
1 1
u13 = (0 + u10 ) = (0 + 0.5036125) = 0.25180625
2 2
1 1
u14 = (u9 + u11 ) = (0.31125 + 0.5036125) = 0.40743125
2 2
1 1
u15 = (u10 + u12 ) = (0.5036125 + 0.31125) = 0.40743125
2 2
1 1
u16 = (u11 + 0) = (0.5036125 + 0) = 0.25180625
2 2
 109

Figure 5.10: The schematic form of the equation (5.2.40)

5.2.5 Solution of wave equation

Consider the one dimensional wave equation
∂2u 2
2∂ u
=c (5.2.38)
∂t2 ∂x2
subject to the conditions u(x, 0) = f (x), ut (x, 0) = 0 (a ≤ x ≤ b), u(a, t) = φ(t) and
u(b, t) = ψ(t). Now consider a rectangular mesh in the x − t plane with spacing h along
x-direction and k along the time t-direction. Replacing the partial derivatives by their
finite difference approximations we get
ui,j+1 = 2(1 − α2 c2 )ui,j + α2 c2 (ui−1,j + ui+1,j − ui,j−1 ) (5.2.39)
k
where α = . If α = 1/c then the equation (5.2.39) reduces to
h
ui,j+1 = ui−1,j + ui+1,j − ui,j−1 (5.2.40)
The schematic form of the equation (5.2.40) is as shown in the Figure 5.10.
Note. The equation (5.2.40) gives the value of u at the (j + 1)th level in terms of the
values of u at the j th and (j − 1)th levels. So, in order to use the formula (5.2.40) to solve
the given wave equation, we need to know the values of u at the two consecutive levels.
The value of u when j = 0(i.e., at 0th level) is got by using the given data u(x, 0) = f (x).
To find u at the first level, we replace the partial derivative in the given data ut (x, 0) = 0
(a ≤ x ≤ b) by its finite difference approximation and get ui,1 = ui,−1 . Now, by putting
j = 0 in the equation (5.2.40), we get
ui,1 = ui−1,0 + ui+1,0 − ui,−1
= ui−1,0 + ui+1,0 − ui,1 (as ui,1 = ui,−1 )
110

Figure 5.11: The schematic form of the equation (5.2.41)

This implies
ui−1,0 + ui+1,0
ui,1 = (5.2.41)
2
which gives the values of u at the first level. The schematic form of the equation (5.2.41)
is as shown in the Figure 5.11.

Example 5.2.7. Evaluate the pivotal values of the equation utt = 16uxx up to t = 1
given that u(0, t) = u(5, t) = 0, ut (x, 0) = 0, u(x, 0) = x2 (5 − x), k = 0.25 and h = 1.

Solution. The difference equation for solving the given wave equation is

ui,j+1 = 2(1 − α2 c2 )ui,j + α2 c2 (ui−1,j + ui+1,j − ui,j−1 ) (5.2.42)

k
Here c2 = 16 which implies c = 4. Given, k = 0.25 and h = 1 which implies α = =
h
1 1
0.25 = = . Thus, the equation (5.2.42) reduces to
4 c
ui,j+1 = ui−1,j + ui+1,j − ui,j−1

Since u(x, 0) = x2 (5 − x) we have

u(1, 0) = 4
u(2, 0) = 12
u(3, 0) = 18
u(4, 0) = 16
The values of u at the first level :
1
u1 = (0 + 12) = 6
2
1
u2 = (4 + 18) = 11
2
 111

Figure 5.12: The rectangular mesh in Example 5.2.7

1
u3 = (12 + 16) = 14
2
1
u4 = (18 + 0) = 9
2
The values of u at the second level :

u5 = 0 + 11 − 4 = 7
u6 = 6 + 14 − 12 = 8
u7 = 11 + 9 − 18 = 2
u8 = 14 + 0 − 16 = −2
The values of u at the third level :

u9 = 0 + 8 − 6 = 2
u10 = 7 + 2 − 11 = −2
u11 = 8 + (−2) − 14 = −8
u12 = 2 + 0 − 9 = −7
112

j\i 0 1 2 3 4 5
0 0 4 12 18 16 0
1 0 6 11 14 9 0
2 0 7 8 2 -2 0
3 0 2 -2 -8 -7 0
4 0 -9 -14 -11 -6 0

Table 5.3: Values of u in Example 5.2.7 at the boundary and interior mesh points

The values of u at the fourth level :

u13 = 0 + (−2) − 7 = −9
u14 = 2 + (−8) − 8 = −14
u15 = −2 + (−7) − 2 = −11
u16 = −8 + 0 − (−2) = −6