Chapter 2

Root finding

Dr. Gokul K. C.
Department of Mathematics,
School of Science,
Kathmandu University, Nepal
Email: gokul.kc@ku.edu.np

Numerical Methods (MCSC 202)

1 / 47
Objective
To find the roots of equation

f (x) = 0 (1)

Where, f (x) may be algebraic or transcendental equations.

For eg: f (x) = x 3 − x − 1, f (x) = sin(x) + 3x, f (x) = x − 1 + e x etc.

Generally, root finding methods broadly classified into two parts:

Classification of root finding methods

Bracketing methods
1 Bisection method
2 False position method
Open methods
1 Secant method
2 Newton Raphson method
3 Iteration method

2 / 47
Rate of Convergence
If a sequence x1 , x2 , x3 , · · · , xn converges to a value α and if there exist real
numbers C > 0 and p > 0 such that

|xn+1 − α|
lim =C (2)
n→∞ |xn − α|p

then we say that p is the rate of convergence of the sequence.

If α is a real root of the equation f (x) = 0 and x1 , x2 , x3 , · · · , xn be the
sequence of approximations then the error in nth approximation is defined
by |n | = |xn − α|. Now, equation (2) becomes

|n+1 |
lim =C (3)
n→∞ |n |p

When n sufficiently large, (3) becomes

|n+1 | = C |n |p (4)

When p = 1, we say the sequence converges linearly and when p = 2 we

say the sequence converges quadratically. If 1 < p < 2 then the sequence
exhibits superlinear convergence.
3 / 47
Bisection Method
Let a and b be chosen such that f (a) × f (b) < 0. i.e. either
f (a) > 0, f (b) < 0 or f (a) < 0, f (b) > 0. The root lies on [a, b].
If x1 be initial approximation then compute x1 = (a + b)/2.
If f (x1 ) = 0, x1 is the desired root. If not then f (x1 ) < 0 or f (x1 ) > 0.
If f (x1 ) × f (a) < 0, then the root lies on [x1 , a] otherwise on [x1 , b].
Let the root lies on [x1 , b] and this new interval is denoted by [a1 , b1 ] such
that |b1 − a1 | = |b − a|/2.
Now, the next approximation is x2 = (x1 + b)/2.
If f (x2 ) = 0 thenx2 is the required root. If not, the root lies either on
[x1 , x2 ] or on [x2 , b]
Let the root lies on [x1 , x2 ] such that x3 = (x1 + x2 )/2
If this new interval is denoted by [a2 , b2 ], then |b2 − a2 | = |b1 − a1 |/2 = |b−a|
22
After nth approximation, the interval [an , bn ] be such that |bn − an | = |b−a|2n .
The computation can be terminated if the error |bn − an | becomes less than
prescribed tolerance(or error()). That is,
|bn − an | = |b−a|
2n ≤ . On simplification gives, n ≥
ln(|b−a|/)
ln2

Now for rate of convergence, we have |bn+1 − an+1 | = |bn −a

2
n|
⇒ n+1 = n
2
⇒ n+1 = C n 1 , C = 1/2. ∴ p = 1. The method linearly converges.
4 / 47
5 / 47
Bisection Method

Computatioal Steps
1 Input: Function(f (x)), tolerance() and initial approximations (a) and
(b), such that f (a) × f (b) < 0.
2 Compute x = (a + b)/2.
3 If f (x) × f (a) < 0, the root lies on [x, a], set b = x and go to step 2.
4 If f (x) × f (b) < 0, the root lies on [x, b], set a = x and go to step 2.
5 Repeat steps 2-4 until |b − a| <  or |f (x)| < 

6 / 47
Example
Solve using bisection method the equation x 3 − x − 1 = 0, correct to 4
decimal places.

Given: f (x) = x 3 − x − 1,
Tolerance=4 decimal places
 = 1/2 × 10−4 = 0.00005
Let a = 1 such that f (a) = f (1) = −1 < 0
and b = 2 such that f (b) = f (2) = 5 > 0.
Total steps required: (n) ≥ ln(|2 − 1|/0.00005)/ln(2) = 14.28
ThereforeTotal steps required: (n) ≈ 15
Our first approximation is
x = (b + a)/2 = (2 + 1)/2 = 1.5.
Further calculations are tabulated below:

7 / 47
 f (a) < 0(−) f (b) > 0(+)
n a b x f(x) Remarks
1 1 2 1.5 0.87 f (x) > 0, replace b by x
2 1 1.5 1.25 -0.29 f (x) < 0, replace a by x
3 1.25 1.5 1.375 0.22 f (x) > 0, replace b by x
4 1.25 1.375 1.3125 -0.05 f (x) < 0, replace a by x
5 1.3125 1.3750 1.3438 0.08 f (x) > 0, replace b by x
6 1.3125 1.3438 1.3281 0.01 f (x) > 0, replace b by x
7 1.3125 1.3281 1.3203 -0.01 f (x) < 0, replace a by x
8 1.3203 1.3281 1.3242 -0.00 f (x) < 0, replace a by x
9 1.3242 1.3281 1.3262 0.00 f (x) > 0, replace b by x
10 1.3242 1.3262 1.3252 0.00 f (x) > 0, replace b by x
11 1.3242 1.3252 1.3247 -0.00 f (x) < 0, replace a by x
12 1.3247 1.3252 1.3250 0.00 f (x) > 0, replace b by x
13 1.3247 1.3250 1.3248 0.00 f (x) > 0, replace b by x
14 1.3247 1.3248 1.3248 0.00 f (x) > 0, replace b by x
Exercise: Find a root, correct to three decimal places of the function
xe x − 1 = 0 Ans : 0.5671

8 / 47
False position method

Let a and b be chosen such that f (a) × f (b) < 0. i.e. either
f (a) > 0, f (b) < 0 or f (a) < 0, f (b) > 0. The root lies on [a, b].
Now, the equation of chord joining (a, f (a)) and (b, f (b)) is given by
y − f (a) = f (b)−f
b−a
(a)
(x − a)
If this chord intersects X-axis at (x1 , 0), we have
0 − f (a) = f (b)−f
b−a
(a)
(x1 − a)
Simplifying, We get first approximation,
f (a) af (b)−bf (a)
x1 = a − f (b)−f (a) (b − a) = f (b)−f (a)
If f (x1 ) = 0, x1 is the desired root. If not then f (x1 ) < 0 or f (x1 ) > 0.
If f (x1 ) × f (a) < 0, then the root lies on [x1 , a] otherwise on [x1 , b].
Let the root lies on [x1 , b], the next approximation(x2 ) is
x2 = x1ff(b)−f
(b)−bf (x1 )
(x1 )
and so on.

9 / 47
10 / 47
False position method

Computatioal Steps
1 Input: Function(f (x)), tolerance() and initial approximations (a) and
(b), such that f (a) × f (b) < 0.
af (b)−bf (a)
2 Compute x = f (b)−f (a)
3 If f (x) × f (a) < 0, the root lies on [x, a], set b = x and go to step 2.
4 If f (x) × f (b) < 0, the root lies on [x, b], set a = x and go to step 2.
5 Repeat steps 2-4 until |xn − xn−1 | <  or |f (x)| < 

11 / 47
Example
Solve using false position method the equation x 3 − x − 1 = 0, correct to 4
decimal places.

Given: f (x) = x 3 − x − 1, Tolerance=4 decimal places,

 = 1/2 × 10−4 = 0.00005
Let a = 1 such that f (a) = f (1) = −1 < 0 and b = 2 such that
f (b) = f (2) = 5 > 0.
Our first approximation is x = aff (b)−bf (a)
(b)−f (a) =
1×5−2×(−1)
5−(−1) = 1.1667. Further
calculations are tabulated below:
n a(-) b(+) x f(x)
1 1 2 1.1667 -0.5787
2 1.1667 2 1.2531 -0.2854
3 1.2531 2 1.2934 -0.1295
Exercise: Find a root,
4 1.2934 2 1.3113 -0.0566
correct to three decimal
5 1.3113 2 1.3190 -0.0243
places of the function
6 1.3190 2 1.3223 -0.0104
2x = log10 x + 7
7 1.3223 2 1.3237 -0.0044
Ans:3.789
8 1.3237 2 1.3243 -0.0019
9 1.3243 2 1.3245 -0.0008
10 1.3245 2 1.3246 -0.0003
11 1.3246 2 1.3247 -0.0001 12 / 47
Secant method

Let a and b be two initial approximations.

The equation of secant line passing through (a, f (a)) and (b, f (b)) is
given by y − f (a) = f (b)−f
b−a
(a)
(x − a)
If this secant line intersects X-axis at (x1 , 0), we have
0 − f (a) = f (b)−f
b−a
(a)
(x1 − a)
Simplifying, We get first approximation,
f (a) af (b)−bf (a)
x1 = a − f (b)−f (a) (b − a) = f (b)−f (a)
If f (x1 ) = 0, x1 is the desired root.
x1 f (b)−bf (x1 )
If not,the next approximation(x2 ) is x2 = f (b)−f (x1 )
and so on.

13 / 47
14 / 47
Secant method

Computatioal Steps
1 Input: Function(f (x)), tolerance() and initial approximations (a) and
(b).
af (b)−bf (a)
2 Compute x = f (b)−f (a)
3 Set a = b and b = x.
4 Repeat steps 2-3 until |xn − xn−1 | <  or |f (x)| < 

15 / 47
Example
Solve using secant method the equation x 3 − x − 1 = 0, correct to 4
decimal places.

Given: f (x) = x 3 − x − 1, Tolerance=4 decimal places,

 = 1/2 × 10−4 = 0.00005
Let a = 1 and b = 2 such that f (a) = f (1) = −1 and f (b) = f (2) = 5. Our
first approximation is x = aff (b)−bf (a)
(b)−f (a) =
1×5−2×(−1)
5−(−1) = 1.1667. Further
calculations are tabulated below:
n a b x
1 1 2 1.1667
2 2 1.1667 1.2531
3 1.1667 1.2531 1.3372
4 1.2531 1.3372 1.3239
5 1.3372 1.3239 1.3247
Exercise: Find a root, correct to three decimal places of the function
4e −x sin x = 1. Ans : 0.370697

16 / 47
Rate of convergence (Secant method)
If xk+1 , xk , xk−1 are three consecutive approximations for the root α of
equation f (x) = 0, then the formula for secant method is

(xk − xk−1 )f (xk )

xk+1 = xk −
f (xk ) − f (xk−1 )

Now k = xk − α, substituting xk = α + k , we obtain

(k − k−1 )f (α + k )
k+1 = k −
f (α + k ) − f (α + k−1 )

Expanding f (α + k ) and f (α + k−1 in Taylor series about the point α and

noting that f (α) = 0, we get

(k − k−1 )[f (α) + k f 0 (α) + 12 2k f 00 (α) + · · · ]

k+1 = k −
[f (α) + k f 0 (α) + 12 2k f 00 (α) + · · · ] − [f (α) + k−1 f 0 (α) + 21 2k−1 f 00 (α) +

(k − k−1 )[k f 0 (α) + 12 2k f 00 (α) + · · · ]

k+1 = k −
(k − k−1 )f 0 (α) + 12 (2k − 2k−1 )f 00 (α) + · · ·
17 / 47
Taking common (k − k−1 )f 0 (α) and cancel out, we get
00
h i
k + 12 2k ff 0 (α)
(α)
+ ···
k+1 = k − h 00 (α)
i
1 + 21 (k + k−1 ) ff 0 (α) + ···

1
Again expanding 1+x = 1 − x + x 2 − x 3 + · · · assuming
00
x = 12 (k + k−1 ) ff 0 (α)
(α)
+ · · · , we get

1 2 f 00 (α) f 00 (α)
    
1
k+1 = k − k + k 0 + ··· 1 − (k + k−1 ) 0 + ··· + ···
2 f (α) 2 f (α)

1 f 00 (α)
k+1 = k k−1 + · · · ⇒ k+1 = Ak k−1
2 f 0 (α)
1/p
Since k+1 = C pk ⇒ k = C pk−1 so that k−1 = C −1/p k , we get

1+1/p
pk = AC −(1+1/p) k

Equating the powers of k , p = 1 + 1/p ⇒ p 2 − p − 1 = 0. Solving we get

p = 1.618. Hence Secant method has superlinear convergence.
18 / 47
Convergence (False position method)

The formula for both Secant method and False position method are same.
The difference is on algorithm while updating the values. In secant
method both ends [a, b] varies, but in false position method one end of [a,
b] is always fixed. Let us suppose ’a’ is fixed, so that the error due to this
end is constant (k−1 = C1 ).
Now, the error equation according to secant method

k+1 = Ak k−1

becomes
k+1 = A k C1 ⇒ k+1 = C2 k C2 = AC1
Substituting k+1 = C pk and equating the powers of k ,

C pk = C2 k

We get p = 1. The false position method converges linearly.

19 / 47
Newton-Raphson method
Let x0 be an approximate root of f (x) = 0 and let x1 = x0 + h be the exact
root such that f (x1 ) = f (x0 + h) = 0.
Expanding f (x0 + h) = 0 by Taylor’s series, we obtain
h2 00
f (x0 + h) = f (x0 ) + hf 0 (x0 ) +
f (x0 ) + · · · = 0 (5)
2!
Neglecting the second and higher order derivatives, we have
f (x0 )
f (x0 ) + hf 0 (x0 ) = 0, h = − (6)
f 0 (x0 )
Substituting the value of h in x1 , we get first approximation
f (x0 )
x1 = x0 − (7)
f 0 (x0 )
After nth approximation, we have
f (xn )
xn+1 = xn − (8)
f 0 (xn )
Which is the Newton-Raphson formula
20 / 47
Method of tangent

21 / 47
Newton-Raphson method

Computatioal Steps
1 Input: Function(f (x)), Derivative of function(f 0 (x)), tolerance() and
initial approximation(a)
f (a)
2 Compute x = a − f 0 (a)
3 Set a = x.
4 Repeat steps 2-3 until |xn − xn−1 | <  or |f (x)| < 

Generalized Newton’s method

If x = a is a root of f (x) = 0 with multiplicity p (Repeated root p times),
then the iteration formula for Newton-Raphson method is given by the
following equation:
f (xn )
xn+1 = xn − p 0 (9)
f (xn )

22 / 47
Example
Solve using Newton-Raphson method the equation x 3 − x − 1 = 0, correct
to 4 decimal places.

Given: f (x) = x 3 − x − 1 and f 0 (x) = 3x 2 − 1

Tolerance = 4 decimal places,  = 1/2 × 10−4 = 0.00005
Let a = 1 such that f (a) = f (1) = −1 and f 0 (a) = f 0 (1) = 2.
Then the first approximation is x = a − ff0(a) (−1)
(a) = 1 − 2 = 1.5.
Further calculations are tabulated below:

n a x
1 1 1.5
2 1.5 1.3478
3 1.3478 1.3252
4 1.3252 1.3247

Exercise: Find a root, correct to three decimal places of the function

2x − 3 = cos x Ans : 1.524
Exercise: Find a double root correct to 4 decimal places of the function
f (x) = x 3 − x 2 − x + 1 = 0
23 / 47
Rate of convergence
Substituting xk = α + k and expanding f (α + k ) and f 0 (α + k ) in Taylor
series about the point α, we get
[k f 0 (α) + 12 2k f 00 (α) + · · · ]
k+1 = k −
[f 0 (α) + k f 00 (α) + · · · ]
Taking common f 0 (α) and cancel out, we get
00
h i
k + 12 2k ff 0 (α)
(α)
+ ···
k+1 = k − h 00 (α)
i
1 + k ff 0 (α) + ···
00
Again expanding 1
1+x = 1 − x + x 2 − x 3 + · · · assuming x = k ff 0 (α)
(α)
+ ···,

1 f 00 (α) f 00 (α)
   
k+1 = k − k + 2k 0 + · · · 1 + k 0 + ···
2 f (α) f (α)
1 f 00 (α) 2
k+1 =  + · · · ⇒ k+1 = A2k
2 f 0 (α) k
⇒ C pk = A2k
Equating the powers of k , we get p = 2. Hence Newton’s method
converges quadratically. 24 / 47
Iteration method

Given equation,
f (x) = 0 (10)
Rewriting equation (10) in the form

x = g (x) (11)

If x0 be initial approximation then substituting in (11), we get

x1 = g (x0 ) (12)

Substituting updated value successively in (11), after nth approximation, we

get the general formula:
xn+1 = g (xn ) (13)
Convergent Criteria
|g 0 (x)| < 1 (14)

25 / 47
Fixed point iteration method

Convergent function (g (x)) Divergent function (g (x))

26 / 47
Fixed point iteration method

Computatioal Steps
1 Input: Function(f (x)), Function(g (x)), tolerance() and initial
approximation(a)
2 Compute x = g (a)
3 Set a = x.
4 Repeat steps 2-3 until |xn − xn−1 | <  or |f (x)| < 

27 / 47
Example
Solve using Iteration method the equation x 3 − x − 1 = 0, correct to 4
decimal places.

Given: f (x) = x 3 − x − 1,Tolerance = 4 decimal places,

 = 1/2 × 10−4 = 0.00005
The possible choices for g(x) are
√ q
g (x) = x 3 + 1 OR g (x) = 3 x + 1 OR g (x) = x 21−1 OR g (x) = x+1 x
√ q
Let us choose g (x) = 3 x + 1 Again if g (x) = x+1 x
and a = 1 such that√ and a = 1 such thatq
x1 = g (a) = g (1) = 3 2 = 1.2599 x1 = g (a) = g (1) = 21 = 1.4142
x2 = g (1.2599) = 1.3123 x2 = g (1.4142) = 1.3066
x3 = g (1.3123) = 1.3224 x3 = g (1.3066) = 1.3287
x4 = g (1.3224) = 1.3243 x4 = g (1.3287) = 1.3239
x5 = g (1.3243) = 1.3246 x5 = g (1.3239) = 1.3249
x6 = g (1.3246) = 1.3247 x6 = g (1.3249) = 1.3247
If we choose g (x) = x + 1 OR g (x) = x 21−1 , then what happens?
3

Exercise: Find a root, correct to three decimal places of the function

sin x = 10(x − 1) Ans : 1.089
28 / 47
Rate of Convergence
Since α is a real root of f (x) = 0, substituting x = α in (11), we get

α = g (α) (15)

Subtracting (15) from (13)

xn+1 − α = g (xn ) − g (α)

Using mean value theorem on RHS,

n+1 = (xn − α)g 0 (η) α < η < xn

Taking modulus
|n+1 | = |xn − α||g 0 (η)| = |n ||g 0 (η)|
|n+1 | = k|n | ∵ |g 0 (η)| = k < 1
|C pn | = k|n |
Equating the powers of k , we get p = 1. Hence iteration method
converges linearly.
29 / 47
Solution to systems of non-linear equations
Iteration Method
Consider the system of two non-linear equations:

f (x, y ) = 0
g (x, y ) = 0
Convergent Criteria
Rewriting the equation in the form:
∂F ∂F
x = F (x, y ) y = G (x, y ) + < 1
∂x ∂y
OR
∂G ∂G
x = G (x, y ) y = F (x, y ) + < 1
∂x ∂y
If (x0 , y0 ) be an initial approximation to the root (α, β), then we get,
x1 = F (x0 , y0 ) y1 = G (x0 , y0 )
x2 = F (x1 , y1 ) y2 = G (x1 , y1 )
.. ..
. .
xn+1 = F (xn , yn ) yn+1 = G (xn , yn )
30 / 47
Exampale
Exmaple: Find a real root of the equations y 2 − 5y + 4 = 0 and
3yx 2 − 10x + 7 = 0 correct to 4 decimal places using initial approximation
(0, 0)

Rewriting the equations in the form:

s
3yx 2 + 7 y2 + 4 10x − 7 p
x= ,y = OR x= , y = 5y − 4
10 5 3y

n x y n x y
1 0.70000 0.80000 8 0.99001 0.99972
2 0.81760 0.92800 9 0.99395 0.99989
3 0.88610 0.97224 10 0.99635 0.99996
4 0.92901 0.98905 11 0.99780 0.99998
5 0.95608 0.99564 12 0.99868 0.99999
6 0.97303 0.99826 13 0.99920 1.00000
7 0.98354 0.99931 14 0.99952 1.00000

Exercise: solve the system: x 2 + y = 11, x + y 2 = 7.

31 / 47
Newton-Raphson Method
Consider the system of two non-linear equations:

f (x, y ) = 0
g (x, y ) = 0

Let (x0 , y0 ) be initial approximation and let (x1 = x0 + h, y1 − y0 + h) be the

exact root such that
f (x1 , y1 ) = f (x0 + h, y0 + k) = 0, g (x1 , y1 ) = g (x0 + h, y0 + k) = 0
Expanding r.h.s using Taylor’s series, we obtain
∂f ∂f
f (x0 + h, y0 + k) = f (x0 , y0 ) + h +k + ··· = 0
∂x ∂y
∂g ∂g
g (x0 + h, y0 + k) = g (x0 , y0 ) + h +k + ··· = 0
∂x ∂y
Neglecting the second and higher order derivatives, we have
∂f ∂f
h +k = −f (x0 , y0 )
∂x ∂y
∂g ∂g
h +k = −g (x0 , y0 )
∂x ∂y
32 / 47
Solving using Cramer’s rule, we get
∂f ∂f ∂f ∂f
∂x0 ∂y0 −f0 ∂y0 ∂x0 −f0
D= ∂g ∂g 6= 0, D1 = ∂g , D2 = ∂g
∂x0 ∂y0 −g0 ∂y0 ∂x0 −g0

Now,
D1 D2
,
h= k=
D D
The new approximations are, therefore,

x1 = x0 + h, y1 = y0 + k

After nth approximation, we get the general formula

xn+1 = xn + h, yn+1 = yn + k

33 / 47
Exampale
Exmaple: Find a real root of the equations 3yx 2 − 10x + 7 = 0 and
y 2 − 5y + 4 = 0 correct to 4 decimal places using initial approximation (0, 0)

Given, f (x, y ) = 3yx 2 − 10x + 7 = 0, g (x, y ) = y 2 − 5y + 4 = 0

∂f ∂f 2 ∂g ∂g
Then, ∂x = 6yx − 10, ∂y = 3x , ∂x = 0, ∂y = 2y − 5
∂f ∂f ∂g ∂g
x y f g D D1 D2 h k
∂x ∂y ∂x ∂y
D1 D2
a b c d ad − bc −fd + gb −ga + fc
D D
0.000 0.000 7.000 4.000 -10.000 0.000 0.000 -5.000 50.000 35.000 40.000 0.700 0.800
0.700 0.800 1.176 0.640 -6.640 1.470 0.000 -3.400 22.576 4.939 4.250 0.219 0.188
0.919 0.988 0.315 0.035 -4.552 2.532 0.000 -3.024 13.764 1.042 0.161 0.076 0.012
0.994 1.000 0.022 0.000 -4.033 2.967 0.000 -3.000 12.101 0.067 0.001 0.006 0.000
1.000 1.000 0.000 0.000 -4.000 3.000 0.000 -3.000 12.000 0.000 0.000 0.000 0.000

Exercise: solve the system: x 2 − y 2 = 4, x 2 + y 2 = 16.

34 / 47
System of linear equations
Consider a system of linear equations:

a11 x1 + a12 x2 + a13 x3 = b1

a21 x1 + a22 x2 + a23 x3 = b2
a31 x1 + a32 x2 + a33 x3 = b3 (16)

This system(16) can be represented by matrix equation

AX = B (17)
     
a11 a12 a13 x1 b1
Where, A =  a21 a22 a23 , X =  x2 , B =  b2 
a31 a32 a33 x3 b3
In this chapter, we will study the following 3 methods to solve linear
system.
LU decomposition method
Tri-diagonal system
Iterative methods
35 / 47
LU decomposition method
 
a11 a12 a13
Let A =  a21 a22 a23  be a nonsingular square matrix. Then A can be
a31 a32 a33
factorized in the form LU such that A = LU, where L is lower triangular
matrix
 and U is upper
 triangular
 matrix represented
 by:
1 0 0 u11 u12 u13
L =  l21 1 0  and U =  0 u22 u23 .
l31 l32 1 0 0 u33
Substituting A = LU in (17), we get

LUX = B (18)

If we set
UX = Y (19)
then,
LY = B (20)
Solving (20) by forward substitution, we obtain vector Y . Substituting this
Y in (19) and solving by backward substitution, we get the required
solution vector X .
36 / 47
Solve the system by LU  
decomposition 2 3 1
∴ U =  0 1/2 5/2 .
2x + 3y + z = 9
0 0 18
x + 2y + 3z = 6
3x + y + 2z = 8 Now for L, we collect all pivot
  columns from pivot element in
2 3 1 downward direction.
Given A =  1 2 3 .  
2 0 0
3 1 2  1 1/2 0 
Since 2 is a pivot element, using
3 −7/2 18
Gaussian elimination, perform
R2 = R2 − 1/2R1 , R3 = R3 − 3/2R1 Dividing every column by corre-
  sponding pivot element, we get
2 3 1  
 0 1/2 5/2  1 0 0
0 −7/2 1/2 ∴ L =  1/2 1 0 
Again, 1/2 is second pivot element, 3/2 −7 1
perform R3 = R3 + 7R2
  Never transform pivot ele-
2 3 1 ment to identity 1
 0 1/2 5/2  = U
0 0 18
37 / 47
  
9
Given, B =  6 . Now, UX = Y is:
8 
2 3 1

x
 
9


Setting UX = Y =⇒ LY = B.  0 1/2 5/2   y  =  3/2 

  0 0 18 z 5
y1
Let Y =  y2  so that LY = B is: solving using backward substitution:
y3
18z = 5 =⇒ z = 5/18
    
1 0 0 y1 9 1/2y + 5/2z = 3/2 =⇒ y = 29/18
 1/2 1 0   y2  =  6  2x + 3y + z = 9 =⇒ x = 35/18
3/2 −7 1 y3 8 
35/18


Solving using forward substitution: ∴ X =  29/18 

5/18
y1 = 9
1/2y1 + y2 = 6 =⇒ y2 = 3/2 Exercise: Find LU and solve AX = B
3/2y1 − 7y2 + y3 = 8 =⇒ y3 = 5 
3 −7 −2
 
−7


A =  −3 5
 
9 1 ,B =  5 
∴ Y =  3/2  6 −4 0 2
5

38 / 47
Tridiagonal system
Consider the system of equations defined by:

b1 x1 + c1 x2 = d1 Thomas algorithm
a2 x1 + b2 x2 + c2 x3 = d2 Set α1 = b1 and
a3 x2 + b3 x3 + c3 x4 = d3 compute
.. αi = bi − aαi ci−1
i−1
,
. i = 2, 3, · · · , n.
an−1 xn−2 + bn−1 xn−1 + cn−1 xn = dn−1 Set β1 = d1 /b1 and
an xn−1 + bn xn = dn compute
βi = di −aαi iβi−1 ,
The matrix of coefficient is i = 2, 3, · · · , n.
 
b1 c1 0 0 ··· ··· ··· 0 Set xn = βn and

 a2 b2 c2 0 ··· ··· ··· 0 
 compute
0 a3 b3 c3 ··· ··· ··· 0
xi = βi − ci αxi+1
 
  ,
 ··· ··· ··· ··· ··· ··· ··· ···  i

 0 0 0 ··· 0 an−1 bn−1 cn−1

 i = n − 1, n − 2, · · · , 1.
0 0 0 ··· 0 0 an bn

39 / 47
 2. Setting β1 = d1 /b1
Solve using Thomas algorithm:
=⇒ β1 = 0/2 = 0 and compute
2x − y = 0 d2 −a2 β1 0−(−1)(0)
−x + 2y − z = 0 β2 = α2 = 3/2 =0
d3 −a3 β2 0−(−1)(0)
−y + 2z − u = 0 β3 = α3 = 4/3 =0
−z + 2u = 1 d4 −a4 β3 1−(−1)(0)
β4 = α4 = 5/4 = 4/5
The augmented matrix is:

2 −1 0 0 : 0
 3. Setting x4 = β4
 −1 2 −1 0 : 0  =⇒ x4 = 4/5 and compute
 
 0 −1 2 −1 : 0  x3 = β 3 − c3 x4
=0− (−1)(4/5)
= 3/5
α3 4/3
0 0 −1 2 : 1 c2 x3 (−1)(3/5)
x2 = β 2 − α2 =0− 3/2 = 2/5
Given i = 4, c1 x2 (−1)(2/5)
x1 = β 1 − α1 =0− 2 = 1/5
1. Setting α1 = b1
=⇒ α1 = 2 and compute
Exercise: Solve the system using
a2 c1 (−1)(−1)
α2 = b2 − α1 =2− 2 = 3/2 Thomas algorithm
a3 c2 (−1)(−1)
α3 = b3 − α2 =2− 3/2 = 4/3 3x1 − x2 = −1
α4 = b4 − a4 c3
=2− (−1)(−1)
= 5/4 −x1 + 3x2 − x3 = 7
α3 4/3
−x2 + 3x3 = 7

40 / 47
Vector Norms
A norm is a way to measure the size of a vector. In other words, norms are
a class of functions that enable us to quantify the magnitude of a vector.
Generally, the norm of a vector is a measure of its length from origin.
A function k.k : Rn → R is called a vector norm if it has the following
properties:
kxk ≥ 0 for any vector x ∈ Rn and kxk = 0 if and only if x = 0
kαxk = |α|kxk for any vector x ∈ Rn and any scalar α ∈ R
kx + y k ≤ kxk + ky k for any vectors x, y ∈ Rn
The most commonly used vector norms belong to the family of p − norms
1/p
or lp − norms defined by kxkp = (Σni=1 |xi |p ) The following p − norms are of
particular interest:
p = 1 : The l1 or Taxicab norm kxk1 = |x1 | + |x2 | + · · · + |xn |
p
p = 2 : The l2 or Euclidean norm kxk2 = x12 + x22 + · · · + xn2
p = ∞ : The l∞ or Maximum norm kxk∞ = max |xi |
i

41 / 47
Matrix Norms
The norm of a matrix is a measure of how large its elements are. It is a
way of determining the size” of a matrix that is not necessarily related to
how many rows or columns the matrix has.
The norm of a square matrix A is a non-negative real number denoted kAk
satisfying the following properties:
kAk ≥ 0 and kAk = 0 if and only if A = 0
kαAk = |α|kAk (any scalar α)
kA + Bk ≤ kAk + kBk
kABk ≤ kAkkBk
The most commonly used matrix norms are:
The l1 − norm or the column norm kAk1 = max Σi |aij |
j
 1/2
The Euclidean norm kAke = Σi,j |aij |2
The l∞ − norm or the row norm kAk∞ = max Σj |aij |
i
The l2 − norm or spectral norm
kAk2 = (Maximum eigenvalue of AT A)1/2

42 / 47
Examples

Example 1: Find l1 , l2 , l∞ norms of ~x = [1, −13, 5, 3, −4]

l1 − norm : k~x k1 = |1| + | − 13| + |5| + |3| + | − 4| = 26
p
l2 − norm : k~x k2 = 12 + (−13)2 + 52 + 32 + (−4)2 = 14.8324
l∞ − norm : max|1|, | − 13|, |5|, |3|, | − 4| = 13
 
1 2 3
Example 2: Find column, Euclidean and row norm of A =  4 5 6 
7 8 9
kAk1 = max[1 + 4 + 7, 2 + 5 + 8, 3 + 6 + 9] = max[12, 15, 18] = 18
√
kAke = 12 + 22 + 32 + 42 + 52 + 62 + 72 + 82 + 92 = 16.88
kAk∞ = max[1 + 2 + 3, 4 + 5 + 6, 7 + 8 + 9] = max[6, 15, 24] = 24

43 / 47
Iteration methods
Consider a system:
a11 x1 + a12 x2 + a13 x3 + · · · + a1n xn = b1
a21 x1 + a22 x2 + a23 x3 + · · · + a2n = b2
a31 x1 + a32 x2 + a33 x3 + · · · + a3n = b3 · · · (∗)
..
.
an1 x1 + an2 x2 + an3 x3 + · · · + ann = bn
If the coefficient matrix of (∗) is diagonally dominant, i.e.
a11 > a12 + a13 + · · · + a1n , a22 > a21 + a23 + · · · + a2n ,
a33 > a31 + a32 + · · · + a3n , · · · , ann > an1 + an2 + · · · + an,n−1 then we can
solve using iterative methods. Rewriting (∗) in the form:
b1 a12 a13 a1n
x1 = a11 − a11 x2 − a11 x3 − ··· − a11 xn
b2 a21 a23 a2n
x2 = a22 − a22 x1 − a22 x3 − ··· − a22 xn

x3 = ab333 − aa31 x1 − aa32 x2 − · · · − aa3n xn · · · (∗∗)

33 33 33
..
.
a
xn = abnnn − aann
n1
x1 − aann
n2
x2 − · · · − n,n−1
ann xn−1
Rewriting (∗∗) in the form
X = BX + C 44 / 47
Jacobi method (Method of simultaneous displacements)
Let (x10 , x20 , x30 , · · · , xn0 ) be initial approximation. Substituting initial
approximations in (∗∗) we get first update (x11 , x21 , x31 , · · · , xn1 ). Again
substituting first update set in (∗∗) we get (x12 , x22 , x32 , · · · , xn2 ) and so on.
Same initial set is used in (∗∗) to update next value.
Convergent Criteria: A sufficient condition for convergence kBk < 1

Gauss Seidel method (Method of successive displacements)

Let (x10 , x20 , x30 , · · · , xn0 ) be initial approximation. Substituting initial
approximations in first equation of (∗∗) to get x11 . Again, substituting
(x11 , x20 , x30 , · · · , xn0 ) in second equation to get x21 . Again, substituting
(x11 , x21 , x30 , · · · , xn0 ) in third equation to get x31 and so on.
Updated values of previous steps are used in (∗∗) to update next value.
Convergent criteria:
Pn aij
j=1,j6=i aii ≤ 1, i = 1, 2, · · · , n,
where the 0 <0 sign should be valid on at least one equation.
Gauss seidel method is 2 times faster than Jacobi method.

45 / 47
 Jacobi method
Solve using iteration methods:
2x − y = 1
−x + 3y − z = 8
−y + 2z = −5

Rewriting,
1+y 8+x+z y −5
x= 2
, y= 3
, z= 2
Using initial approximation (0, 0, 0), the
results are tabulated below:
Gauss-Seidel method

46 / 47
Exercise
Solve using iterative methods:

10x − 2y − z − u = 3
−2x + 10y − z − u = 15
−x − y + 10z − 2u = 27
−x − y − 2z + 10u = −9

47 / 47