Page 105 - Applied Numerical Methods Using MATLAB

94 SYSTEM OF LINEAR EQUATIONS
 
a 11 a 12 a 13
step 1:  a 21 a 22 a 23 
a 31 a 32 a 33
 
u 11 = a 11 u 12 = a 12 u 13 = a 13
(1) (1)
 
→  l 21 = a 21 /u 11 a = a 22 − l 21 u 12 a (2.4.6a)
22 23 = a 23 − l 21 u 13 
(1) (1)
l 31 = a 31 /u 11 a = a 32 − l 31 u 12 a = a 33 − l 31 u 13
32 33
 
u 11 u 12 u 13
(1) (1)
 
step 2: →  l 21 u 22 = a u 23 = a (2.4.6b)
22 23 
(1) (2) (1)
l 31 l 32 = a /u 22 a = a − l 32 u 23
32 33 33
This leads to an LU decomposition algorithm generalized for an NA × NA
nonsingular matrix as described in the following box. The MATLAB routine
“lu_dcmp()” implements this algorithm to find not only the lower/upper
triangular matrix L and U, but also the permutation matrix P .Werun it for
a3 × 3matrixtoget L, U,and P and then reconstruct the matrix P −1 LU = A
from L, U,and P to ascertain whether the result is right.
function [L,U,P] = lu_dcmp(A)
%This gives LU decomposition of A with the permutation matrix P
% denoting the row switch(exchange) during factorization
NA = size(A,1);
AP = [A eye(NA)]; %augment with the permutation matrix.
fork=1:NA-1
%Partial Pivoting at AP(k,k)
[akx, kx] = max(abs(AP(k:NA,k)));
if akx < eps
error(’Singular matrix and No LU decomposition’)
end
mx = k+kx-1;
if kx>1%Row change if necessary
tmp_row = AP(k,:);
AP(k,:) = AP(mx,:);
AP(mx,:) = tmp_row;
end
% LU decomposition
form=k+1:NA
AP(m,k) = AP(m,k)/AP(k,k); %Eq.(2.4.8.2)
AP(m,k+1:NA) = AP(m,k + 1:NA)-AP(m,k)*AP(k,k + 1:NA); %Eq.(2.4.9)
end
end
P = AP(1:NA, NA + 1:NA + NA); %Permutation matrix
for m = 1:NA
for n = 1:NA
if m == n, L(m,m) = 1.; U(m,m) = AP(m,m);
elseifm>n, L(m,n) = AP(m,n); U(m,n) = 0.;
else L(m,n) = 0.; U(m,n) = AP(m,n);
end
end
end
if nargout == 0, disp(’L*U = P*A with’); L,U,P, end
%You can check if P’*L*U = A?