Page 106 - Applied Numerical Methods Using MATLAB

DECOMPOSITION (FACTORIZATION) 95
(cf) The number of floating-point multiplications required in this routine lu_dcmp() is

NA−1 NA−1

2
(NA − k)(NA − k + 1) = {NA(NA + 1) − (2NA + 1)k + k }
k=1 k=1
1 1
= (NA − 1)NA(NA + 1) − (2NA + 1)(NA − 1)NA + (NA − 1)NA(2NA − 1)
2 6
1 1 3
= (NA − 1)NA(NA + 1) ≈ NA (2.4.7)
3 3
with NA: the size of matrix A

(0)
0. Initialize A (0) = A, or equivalently, a mn = a mn for m, n = 1: NA.
1. Let k = 1.
2. If a (k−1) = 0, do an appropriate row switching operation so that
kk
(k−1)
a = 0.
kk
When it is not possible, then declare the case of singularity and stop.
3. a (k) = a (k−1) for n = k : NA (Just leave the kth row as it is.)
kn kn = u kn
(2.4.8a)
(k) (k−1) (k−1)
a = a /a = l mk for m = k + 1: NA (2.4.8b)
mk mk kk
(k) (k)
4. a (k) = a (k−1) − a a for m, n = k + 1: NA (2.4.9)
mn mn mk kn
5. Increment k by 1 and if k< NA − 1, go to step 1; otherwise, go to step 6.
6. Set the part of the matrix A (NA−1) below the diagonal to L (lower tri-
angular matrix with the diagonal of 1’s) and the part on and above the
diagonal to U (upper triangular matrix).

>>A = [1 2 5;0.2 1.6 7.4; 0.5 4 8.5];
>>[L,U,P] = lu_dcmp(A) %LU decomposition
L = 1.0 0 0 U = 1 2 5 P = 1 0 0
0.5 1.0 0 0 3 6 0 0 1
0.2 0.4 1.0 0 0 4 0 1 0
>>P’*L*U - A %check the validity of the result (P’ = P^-1)
ans = 0 0 0
0 0 0
0 0 0
>>[L,U,P] = lu(A) %for comparison with the MATLAB built-in function

What is the LU decomposition for? It can be used for solving a system of
linear equations as
Ax = b (2.4.10)

T
Once we have the LU decomposition of the coefficient matrix A = P LU,it is
more efficient to use the lower/upper triangular matrices for solving Eq. (2.4.10)