Page 35 - Numerical methods for chemical engineering

Elimination methods for solving linear systems 21

zero,
 
a 11 a 12 a 13 a 14 ... a 1N b 1
(2, 1) (2, 1) (2, 1) (2, 1) (2, 1) 
 0 a a a ... a b

22 23 24 2N 2 
 
(3, 2) (3, 2) (3, 2)
 0 0 a 33 a 34 ... a 3N b (3, 2) 

3

 (4, 3) 
(A, b) (N, 3) =  0 0 0 0 ... a 4N b (4, 3)  (1.107)

4

 (5, 3) (5, 3) 
 0 0 0 0 ... a 5N b 5 
 
 . . . . . . 
 . . . . . . 
 . . . . . . 
(N, 3) (N, 3)
0 0 0 0 ... a b
NN N
Thus, we can write the fourth column as a linear combination of the first three columns, for
some c 1 , c 2 , c 3 ,
a j4 = c 1 a j1 + c 2 a j2 + c 3 a j3 ∀ j = 1, 2,..., N (1.108)
so that any equation j in this system can be written as
a j1 x 1 + a j2 x 2 + a j3 x 3 + (c 1 a j1 + c 2 a j2 + c 3 a j3 )x 4 + ··· + a jN x N = b j
a j1 (x 1 + c 1 x 4 ) + a j2 (x 2 + c 2 x 4 ) + a j3 (x 3 + c 3 x 4 ) +· · · + a jN x N = b j (1.109)
This means that we can make any change to x of the form
x 1 → x 1 − c 1 x 4 x 2 → x 2 − c 2 x 4
(1.110)
x 3 → x 3 − c 3 x 4 x 4 → x 4 + x 4
without changing the value of Ax. Obviously, there can be no unique solution of Ax = b,
and we thus stop the elimination procedure at this point.
If a unique solution exists, Gaussian elimination with partial pivoting will find it, even
if there are zeros along the principal diagonal of the original augmented matrix. It is called
partial pivoting because we only swap rows; if we swap columns as well (which necessi-
tates more complex book-keeping, but results in better propagation of round-off error), we
perform Gaussian elimination with full pivoting.
The algorithm for Gaussian elimination with partial pivoting is

for i = 1, 2,..., N − 1; iterate over columns
select row j ≥ i such that |a ji |= max j≥i {|a ii |, |a i+1, i |,..., |a N1 |}
if a ji = 0, no unique solution exists, STOP
if j = i, interchange rows i and j of augmented matrix
for j = i + 1, i + 2,..., N; iterate over rows j > i
λ ← a ji /a ii
for k = i, i + 1,..., N; iterate over elements in row j from left
a jk ← a jk − λa ik
end for k = i, i + 1,..., N
b j ← b j − λb i
end for j = i + 1, i + 2,..., N
end for i = 1, 2,..., N − 1