Page 49 - Numerical Methods for Chemical Engineering

38 1 Linear algebra

that every linear system A˜a (k) = e [k] to be solved has the same matrix A. Therefore, dur-
ing Gaussian elimination, the sequence of row operations and pivoting is the same for
each system. It would be very helpful if we could do this elimination only once, and
store all of the details of how Gaussian elimination is performed for A, so that for sub-
sequent solution of any system Ax = b, we need only perform the backward substitution
3
step for the new b. Then, the work of computing the inverse scales only as N rather
4
than N .

Matrix factorization

We have seen that the column vectors of A −1 may be computed one-by-one by solving the
N linear systems

A˜a (k) = e [k] k = 1, 2,..., N (1.187)

More generally, let us say that we wish to solve some set of M linear systems that have
[2]
[1]
the same matrix A but different vectors b , b ,..., b [M] . In this section, we show that
it is possible in this case to do the work of Gaussian elimination only once, and to factor
A into a product of two triangular matrices. Thus, each system with a new vector b [k] can
2
be solved with two successive substitutions. As these substitutions require only N FLOPs
3
each, a much smaller number than the 2N /3 FLOPs required for Gaussian elimination,
this factorization can save much effort.

LU decomposition
Let us say that we want to solve M linear systems that have the same matrix,

Ax [k] = b [k] k = 1, 2,..., M (1.188)

Let us assume that it is possible to decompose A into the product of a lower triangular matrix
L and an upper triangular matrix U, A = LU. Then,

Ax [k] = LUx [k] = b [k] (1.189)

We obtain x [k] quickly by solving in succession two triangular problems, the first by forward
substitution and the second by backward substitution,

Lc [k] = b [k]
Ux [k] = c [k] (1.190)

We now show that with some additional book-keeping, Gaussian elimination without par-
tial pivoting returns just such an LU factorization, A = LU. We now perform Gaussian