Systems of Linear Algebraic Equations - Numerical Methods in Engineering with MATLAB - page 71

%--------------Elimination pass---------------
fori=k+1:n
if A(i,k) ˜= 0
lambda = A(i,k)/A(k,k);
A(i,k+1:n) = A(i,k+1:n) - lambda*A(k,k+1:n);
b(i) = b(i) - lambda*b(k);

end

end

end
%------------Back substitution phase----------
fork=n:-1:1
b(k) = (b(k) - A(k,k+1:n)*b(k+1:n))/A(k,k);

end
x=b;

LUdecPiv

The Gauss eliminationalgorithm can bechanged to Doolittle's decompositionwith
minor changes. The most important of these is keeping a record of the rowinter-
changes during the decompositionphase. In LUdecPiv this record iskept in the
permutation array perm , initially set to [1 n ] T .Whenever two rows are inter-

changed, the corresponding interchange is also carried out in perm . Thus perm shows
how the original rows were permuted. This informationis thenpassed to the function
LUsolPiv , which rearranges the elements of the constant vectorinthe sameorder
beforecarrying outforward and back substitutions.

, 2 ,...,

function [A,perm] = LUdecPiv(A)
%LUdecompositionofmatrixA;returnsA=[L\U]
% and the row permutation vector 'perm'.
% USAGE: [A,perm] = LUdecPiv(A)

n = size(A,1); s = zeros(n,1);
perm = (1:n)';
%----------Set up scale factor array----------
fori=1:n;s(i)=max(abs(A(i,1:n)));end
%---------Exchange rows if necessary----------
fork=1:n-1
[Amax,p] = max(abs(A(k:n,k))./s(k:n));
p=p+k-1;
if Amax < eps
error('Matrix is singular')

end