44 // Numerical Methods //
The process of computing L and U for a given A is known as LU Decomposition or LU Factorisation. LU
decomposition is not unique (the combinations of L and U for a prescribed A are endless), unless certain
constraints are placed on L or U. These constraints distinguish one type of decomposition from another.
Two commonly used decompositions are given below:
1. Cholesky’s decomposition: Constraints are L = UT
2. Crout’s decomposition: Constrains are Uii = 1, i = 1, 2, ..., n.
After decomposing the matrix A, it is easier to solve the equations Ax = b.
We can rewrite the equations as
LUx = b
or denoting Ux = y, the above equation becomes
Ly = b
This equation Ly = b can be solved for y by forward substitution. Then Ux = y will yield x by the
backward substitution process. The advantage of LU decomposition method over the Gauss elimination
method is that once A is decomposed, we can solve Ax = b for as many constant vectors b as we please.
Also, the forward and backward substitutions operations are much less time consuming than the
decomposition process.
2.7 CHOLESKY’S TRIANGULARISATION METHOD
Cholesky’s decomposition method is faster than the LU decomposition. There is no need for pivoting. If the
decomposition fails, the matrix is not positive definite.
Consider the system of linear equations:
a11x1 + a12x2 + a13x3 = b 1
a21x1 + a22x2 + a23x3 = b 2
a31x1 + a32x2 + a33x3 = b 3 (2.21)
The above system can be written as (2.22)
Ax = b
a11 a12 a13 x1 b1
where A = a21 a22
a23 , x = x2 , b = b2
a31 a32 a33 x3 b3
Let A = LU… (2.23)
1 0 0 u11 u12 u13
L = l21 1 0 and U = 0 u22
u23
l31 l32 1 0 0 u33
Equation (2.21) can be written as
LUX = b (2.24)
If we write UX = V (2.25)
// Linear System of Equations // 45
Equation (2.24) becomes
LV = b (2.26)
Equation (2.26) is equivalent to the system
v1 = b 1
l21v1 + v2 = b 2
l31v1 + l32v2 + v3 = b 3 (2.27)
The above system can be solved to find the values of v1, v2 and v3 which give us the matrix V.
UX = V
then becomes
u11x1 + u12x2 + u13x3 = v 1
u22x2 + u23x3 = v 2
u33x 3 = v 3 (2.28)
which can be solved for x3, x2 and x1 by the backward substitution process.
In order to compute the matrices L and U, we write Eq. (2.23) as
1 0 0 u11 u12 u13 a11 a12 a13
l 0 0 u22 u23 = a21 a22 a23
21 1 (2.29)
l31 l32 1 0 0 u33 a31 a32 a33
Multiplying the matrices on the left and equating the corresponding elements of both sides, we obtain
u11 = a11, u12 = a12, u13 = a13 (2.30)
a21
l21u11 = a21 ⇒ l21 =
a11
a (2.31)
l31u11 = a31 ⇒ l31 = 31
a11
a21
l21u12 + u22 = a22 ⇒ u22 = a22 − a12
a11
a21 (2.32)
l21u13 + u23 = a23 ⇒ u23 = a23 − a13
a11
1 a31
l31u12 + l32 u22 = a32 ⇒ l32 = a32 − a a12 (2.33)
u22 11
and l31u13 + l32u23 + u33 = a33 (2.34)
The value of u33 can be computed from Eq. (2.34).
To obtain the elements of L and U, we first find the first row of U and the first column of L. Then, we
determine the second row of U and the second column of L. Finally, we compute the third row of U.
46 // Numerical Methods //
Cholesky’s triangularisation method is also known as Crout’s triangularisation method or method of
factorisation.
Example E2.11
Solve the following equations by Cholesky’s triangularisation method.
2x + y + 4z = 12
8x – 3y + 2z = 20
4x + 11y – z = 33
Solution:
2 1 4 x 12
We have A = 8 −3 2 , X = y , B = 20
4 11 −1 z 33
1 0 0 u11 u12 u13 2 1 4
l 0 0 u22
u23 = 8 −3 2
Let 21 1
l31 l32 1 0 0 u33 4 11 −1
Multiplying and equating we get:
l × u11 = 2 ⇒ u11 = 2
l × u12 = 1 ⇒ u12 = 1
l × u13 = 4 ⇒ u13 = 4
8 8
l21 × u11 = 8 ⇒ l21 = = =4
u11 2
l21 × u12 + u22 = –3 ⇒ u22 = –3 – l21 × u12 = –3 – 4 × 1 = –7
l21 × u13 + u23 = 2 ⇒ u23 = 2 – l21 × u13 = 2 – 4 × 4 = –14
4 4
l31 × u11 = 4 ⇒ l31 = = =2
u11 2
11 − l31 × u12 11 − 2 × 1 9
l31 × u12 + l32 × u22 = 11 ⇒ l32 = = =−
u22 −7 7
9
l31 × u13 + l32 × u23 + l × u33 = –1 ⇒ u33 = –1 – l31 × u13 – l32 × u23 = –1 – 2 × 4 – − ( −14) = –27
7
1 0 0
4 2 1 4
1 0
We get: A= 0 −7 − 14
9
2 − 1 0 0 −27
7
// Linear System of Equations // 47
and the given system can be written as:
1 0 0
4 2 1 4 x 12
1 0
0 −7 −14 y = 20
9
2 − 1 0 0 −27 z 33
7
Writing: LV = B, we get
1 0 0
4 V1 12
1 0 V2 = 20
9
2 − 1 V3 33
7
which gives V1 = 12 ⇒ V2
4V1 + V2 = 20 ⇒ V2 = 20 – 4 × 12 = –28
9 9
2V1 – V2 + V3 = 33 ⇒ V3 = 33 + (–28) – 2 × 12 = –27
7 7
The solution to the original system is given by:
UX = V
2 1 4 x 12
0 −7 −14 y = −28
0 0 −27 z −27
2x + y + 4z = 12
–7y – 14z = –28
–27z = –27 ⇒ z =1
14
7y = 28 – 14 × 1 ⇒ y = ⇒ y=2
7
6
2x = 12 – y – 4z = 12 – 2 – 4 × 1 ⇒ x = ⇒ x = 3
2
Example E2.12
Solve the system of equations using Cholesky’s factorisations.
x1 + x2 + x3 – x4 = 2
x1 – x2 – x3 + 2x4 = 0
4x1 + 4x2 + x3 + x4 = 11
2x1 + x2 + 2x3 – 2x4 = 2
48 // Numerical Methods //
Solution:
The set of equations can be written in the matrix form [A]{x} = {b}
1 1 1 −1 x1 2
1 −1 −1 2 x 0
2 =
4 4 1 1 x3 11
2 1 2 −2 x4 2
Let us decompose [A] in the form
[A] = [L] [U]
1 0 0 0 u11 u12 u13 u14
l 1 0 0 0 u u23 u24
[ L] = and [U ] =
21 22
where
l31 l32 1 0 0 0 u33 u34
l41 l42 l43 1 0 0 0 u44
The product of [L][U] gives
u11 u12 u13 u14
l u l21u12 + u22 l21u13 + u23 l21u14 + u24
[ L][U ] =
21 11
l31u11 l31u12 + l32 u22 l31u13 + l32 u23 + u33 l31u14 + l32 u24 + u34
l41u11 l41u12 + l42 u22 l41u13 + l42 u23 + l43u33 l41u14 + l42 u24 + l43u34 + u44
Equating the elements of this matrix to the [A] matrix yields the following equations
u11 = 1 l21u11 = 1 l31u11 = 4 l41u11 = 2
u12 = 1 l21u12 + u22 = –1 l31u12 + l32u22 = 4 l41u12 + l42u22 = 1
u13 = 1 l21u13 + u23 = –1 l31u13 + l32u23 + u33 = 1 l41u13 + l42u23 + l23u33 = 2
u14 = –1 l21u14 + u24 = 2 l31u14 + l32u24 + u34 = 1 l41u14 + l42u24 + l43u34 + u44 = –2
By solving these sixteen equations we get
1 0 0 0 1 1 1 −1
1 1 0 0
[ L] = and [U ] = 0 −2 −2 3
4 0 1 0 0 0 −3 5
1
2 2 − 3 1
1 1
0 0 0 6
To solve [A]{x} = {b} we have to solve the two systems
[L]{Y} = {b}
[U]{x} = {Y}
// Linear System of Equations // 49
1 0 0 0 y1 2
1 1 0 0
y2 = 0
4 0 1 0 y3 11
i.e.,
2 2 3 1 y4 2
1 1
which gives by forward substitution
y1 = 2, y2 = –2, y3 = 3, y4 = 0
and hence [U]{x} = {y} becomes
1 1 1 −1 x1 2
0 −2 −2 3
x2 = −2
0 0 −3 5 x3 3
1
0 0 0 x4 0
6
Then by back substitution we obtain
x4 = 0, x3 = –1, x2 = 2, x1 = 1.
Example E2.13
Solve the system of linear equations using Cholesky’s factorisation method.
2x – 6y + 8z = 24
5x + 4y – 3z = 2
3x + y + 2z = 16
Solution:
1 0 0 u11 u12 u13 2 −6 8
l 0 0 u22 u23 = 5 4 −3
21 1
l31 l32 1 0 0 u33 3 1 2
u11 u12 u13 2 −6 8
l u l21u12 + u22 l21u13 + u23 = 5 4 −3
21 11
l31u11 l31u12 + l32 u22 l31u13 + l32 u23 + u33 3 1 2
u11 = 2, u12 = –6, u13 = 8
5
l21 = = 2.5
u11
3
l31 = = 1.5
u11
u22 = 4 – l21u12 = 19
u23 = –3 – l21u13 = –23
50 // Numerical Methods //
1 − l31u12 10
l32 = =
u22 19
40
l33 = 2 – l31u13 – l32u23 =
19
1 0 0 2 −6 8
0 , U = 0 19 −23
L = 2.5 1
1.5 10
19
1 0 0 40
19
1 0 0 v1 24
LV = B ⇒ 2.5 1 0 v2 = 2
1.5 10
19
1 v3 16
⇒ v1 = 24
v2 = 2 – 2.5 × 24 = –58
10 200
v3 = 16 – 1.5 × 24 – ( −58) =
19 19
24
2 −6 8 x
−58
UX = V ⇒ 0 19 −23 y =
200
0 0 40
19 z
19
2x – 6y + 8z = 24 (E.1)
19y – 23z = –58 (E.2)
40 200
z= ⇒ z=5 (E.3)
19 19
From Eqs.(E.2) and (E.3), we have
y=3
From Eqs.(E.1), (E.2) and (E.3), we get
x =1
2.8 CROUT’S METHOD
This method is based on the fact that every square matrix A can be expressed as the product of a lower
triangular matrix and an upper triangular matrix, provided all the principle minors of A are non-singular. Also,
such a factorisation, if exists, is unique.
This method is also called triangularisation or factorisation method. Here, we factorise the given matrix
as A = LU, where L is a lower triangular matrix with unit diagonal elements and U is an upper triangular
matrix. Then,
A–1 = (LU)–1 = U–1L–1
// Linear System of Equations // 51
Consider the system
a11x1 + a12x2 + a13x3 = b 1
a21x1 + a22x2 + a23x3 = b 2
a31x1 + a32x2 + a33x3 = b 3 (2.35)
The above system can be written as
Ax = b
Let A = LU (2.36)
l11 0 0 1 u12 u13
where L = l21 l22
0 and U = 0 1 u23 (2.37)
l31 l32 l33 0 0 1
Here, L is a lower triangular matrix and U is an upper triangular matrix with diagonal elements equal to unity.
A = LU ⇒ A–1 = U–1L–1 (2.38)
a11 a12 a13 l11 0 0 1 u12 u13
Now A = LU ⇒ a21 a22 a23 = l21 l22 0 0 1 u23
a31 a32 a33 l31 l32 l33 0 0 1
a11 a12 a13 l11 l11u12 l11u13
a
a23 = l21 l21u12 + l22 l21u13 + l22u23
or 21 a22
a31 a32 a33 l31 l31u12 +l32 l31u13 +l32u23 +l33
Equating the corresponding elements, we obtain
l11 = a11 l21 = a21 l31 = a31 (2.39)
l11u12 = a12 l11u13 = a13 (2.40)
l21u12 + l22 = a22 l31u12 + l32 = a32 (2.41)
l21u13 + l22u23 = a23 (2.42)
and l31u13 + l32u23 + l33 = a33 (2.43)
from (2.40) we find
u12 = a12 /l11 = a12/a11
from (2.41) we obtain
l 22 = a22 – l21u12 (2.44)
l 32 = a32 – l31u12 (2.45)
Equation (2.42) gives
u 23 = (a23 – l21u23)/l22 (2.46)
from the relation (2.43) we get
l 33 = a33 – l31u13 – l32u23 (2.47)
Thus, we have determined all the elements of L and U.
From Eqs.(2.36) and (2.37) we have
LUx = b (2.48)
52 // Numerical Methods //
Let UX = V
v1
v
V=
2
where
vn
From Eq. (2.48) we have LV = b, which on forward substitution yields V.
From UX = V, we find x (by backward substitution).
Example E2.14
Solve the following set of equations by Crout’s method:
2x + y + 4z = 12
8x – 3y + 2z = 20
4x + 11y – z = 33
Solution:
2 1 4 x 12
We have A = 8 −3 2 , X = y , B = 20
4 11 −1 z 33
AX = B
Let A = LU
l11 0 0 1 u12 u13
L = l21 l22 0 U = 0 1 u23
l31 l32 l33 0 0 1
2 1 4 l11 0 0 1 u12 u13
8 −3 2 = l 0 0 1 u23
21 l22
4 11 −1 l31 l32 l33 0 0 1
2 1 4 l11 l11u12 l11u13
8 −3 2 = l l21u13 + l22u23
21 l21u12 + l22
4 11 −1 l31 l31u12 + l32 l31u13 + l32 u23 + l33
1
l 11u 12 = 1 ⇒ u12 =
2
4
l 11u 12 = 4 ⇒ u13 = =2
2
1
l22 + l21u12 = –3 ⇒ l22 = −3 − 8 = −7
2
1
l32 + l21u12 = –3 ⇒ l32 = 11 − 4 = 9
2
// Linear System of Equations // 53
2−8× 2
l21u13 + l22u23 = 2 ⇒ u23 = =2
−7
l31u13 + l32u23 + l33 = –1 ⇒ l33 = −1 − 4 × 2 − 9 × 2 = −27
2 0 0 1 12 2
L = 8 −7 0 and U = 0 1 2
4 9 −27 0 0 1
LV = B
2 0 0 v1 12
8 −7
0 v2 = 20
4 9 −27 v3 33
2v1 = 12 ⇒ v1 = 6
−20 + 8 × 6
8v1 – 7v2 = 20 ⇒ v2 = =4
7
−33 + 4 × 6 + 9 × 4
4v1 + 9v2 – 27v3 = 33 ⇒ v3 = =1
27
V1 6
V = V2 = 4 ; Ux = V
V3 1
1 12 2 x 6
0 1 2 y = 4
0 0 1 z 1
1
x+ y + 2z = 6
2
y + 2z = 4
z =1
y = 4 – 2× 1
⇒ y=2
1
x=6– × 2–2× 1 ⇒ x = 3
2
Example E2.15
Solve the following set of equations by using the Crout’s method:
2x1 + x2 + x3 = 7
x1 + 2x2 + x3 = 8
x1 + x2 + 2x3 = 9
54 // Numerical Methods //
Solution:
2 1 1 x 7
A = 1 2 1 , x = y , B = 8
1 1 2 z 9
Let A = LU
l11 0 0 1 u12 u13
L = l21 l22 0 U = 0 1 u23
l31 l32 l33 0 0 1
2 1 1 l11 l11u12 l11u13
1 2 1 = l l21u13 + l22u23
21 l21u12 + l22
1 1 2 l31 l31u12 + l32 l31u13 + l32 u23 + l33
l11 = 2, l21 = 1, l31 = 1
1 1
u12 = , u13 =
2 2
1 3
l22 = 2 – l21u12 = 2 – 1 × =
2 2
1 1
l32 = 1 – l31u12 = 1 – 1 × =
2 2
1 − l21u13 1
u23 = =
l22 3
1 1 1 4
l33 = 2 – l31u13 – l32u23 = 2 – − × =
2 2 3 3
2 0 0 1 1/ 2 1/ 2
L = 1 3/ 2 0 , U = 0
1 1/ 3
1 1/ 2 4 / 3 0 0 1
Ax = B, LU.x = B, Ux = V
2 0 0 v1 7
LV = B ⇒ 1 3 / 2 0 v2 = 8
1 1/ 2 4 / 3 v3 9
2v1 = 7 ⇒ v1 = 3.5
3 3.5
v1 + v2 = 8 ⇒ v2 = 3 3
2 ⇒ V =
1 4 3
v1 + v2 + v3 = 9 ⇒ v3 = 3
2 3
// Linear System of Equations // 55
1 1/ 2 1/ 2 x1 3.5
Ux = V ⇒ 0 1 1/ 3 x2 = 3
0 0 1 x3 3
1 1
x1 + x2 + x3 = 3.5 (E.1)
2 2
1
x2 + x3 = 3 (E.2)
3
x3 = 3 (E.3)
From Eqs.(E.2) and (E.3), we have
x2 = 2
From Eq.(E.1), we get
x1 = 1
2.9 THOMAS ALGORITHM FOR TRIDIAGONAL SYSTEM
Consider the system of linear simultaneous algebraic equations given by
Ax = b
where A is a tridiagonal matrix, x = [x1, x2, …, xn]T and b = [b1, b2, …, bn]T. Hence, we consider a 4 × 4
tridiagonal system of equations given by
a12 a13 0 0 x1 b1
a b
21 a22 a23 0 x2 2
=
0 a31 a32 a33 x3 b3 (2.48a)
0 0 a41 a42 x4 b4
Equation (2.48a) can be written as
a12x1 + a13x2 = b 1
a21x1 + a22x2 + a23x3 = b 2
a31x2 + a32x3 + a33x4 = b 3
a41x3 + a42x4 = b 4 (2.48b)
The system of equations given by Eq.(2.48b) is solved using Thomas Algorithm which is described in three
steps as shown below:
Step 1: Set y1 = a12 and compute
ai1a(i −1)3
yi = ai 2 − i = 2, 3, …, n
yi −1