Numerical Linear Algebra:
Solving
Ax=b using Elimination Methods
Outline:
I.
Gaussian
Elimination
II.
Gaussian
Elimination with Pivoting
III.
Gauss-Jordan
Elimination
IV.
LU Decomposition
V.
Matrix Inverses, Norms
and Condition Numbers
I. Naïve Gaussian elimination
Eliminate 1st to (n-1)st
equation from suceeding equations.
Subtract multiple of equation multiplied by ratio of coefficients 
for 
Problems:
1.
Division by zero
2.
Round off error
3.
II. Pivoting –
Switch equations so that at any step, 
is largest
Partial Pivoting:
interchange rows (equations)
Complete Pivoting: interchanges rows (equations) and order of variables (columns)
Matlab script for Gaussian elimination – Gauss.m
%GAUSS
Solves a set of linear algebraic equations by the Gauss
% elimination method.
%
% GAUSS(A,C) finds
unknowns of a set of linear algebraic
% equations. A is the matrix
of coefficients and C is the
% vector of constants.
%
% See also
%
(c)
%
c
= (c(:).')';
% Make sure it's a column vector
n
= length(c);
[nr nc] = size(A);
%
Check coefficient matrix and vector of constants
if nr ~= nc
error('Coefficient
matrix is not square.')
end
if nr ~= n
error('Coefficient
matrix and vector of constants do not have the same length.')
end
%
Check if the coefficient matrix is singular
if det(A) == 0
fprintf('\n Rank = %7.3g\n',rank(A))
error('The
coefficient matrix is singular.')
end
unit = diag(ones( 1 , n)); % Unit matrix
order = [1 : n]; % Order
of unknowns
aug = [A c]; % Augmented matrix
%
Gauss elimination
for k = 1 : n - 1
pivot = abs(aug(k , k));
prow = k;
pcol
= k;
% Locating the
maximum pivot element
for row = k : n
for col = k : n
if abs(aug(row , col)) > pivot
pivot =
abs(aug(row , col));
prow =
row;
pcol = col;
end
end
end
% Interchanging the
rows
pr = unit;
tmp
= pr(k , :);
pr(k , : ) =
pr(prow , : );
pr(prow , : ) = tmp;
aug
= pr * aug;
% Interchanging the
columns
pc = unit;
tmp = pc(k , : );
pc(k , : ) = pc(pcol , : );
pc(pcol , : ) = tmp;
aug(1 : n , 1 : n)
= aug(1 : n , 1 : n) * pc;
order = order * pc; % Keep track of the column interchanges
% Reducing the
elements below diagonal to zero in the column k
lk
= unit;
for m = k + 1 : n
lk(m , k) = - aug(m , k) / aug(k , k);
end
aug
= lk * aug;
end
x
= zeros(n , 1);
%
Back substitution
t(n) = aug(n , n + 1) / aug(n , n);
x(order(n)) = t(n);
for k = n - 1 : -1 : 1
t(k) = (aug(k,n+1)
- sum(aug(k,k+1:n) .* t(k+1:n))) / aug(k,k);
x (order(k)) =
t(k);
end
III. Gauss
%
% Gauss-Jordan method.
%
% JORDAN(A,C) finds
unknowns of a set of linear algebraic
% equations. A is the matrix of coefficients and C is the
% vector of constants.
%
% See also GAUSS, SEIDEL
% (c)
%
c = (c(:).')'; % Make sure it's a column vector
n = length(c);
[nr
nc] = size(A);
% Check coefficient matrix
and vector of constants
if nr ~= nc
error('Coefficient
matrix is not square.')
end
if nr ~= n
error('Coefficient
matrix and vector of constants do not have the same length.')
end
% Check if the coefficient
matrix is singular
if det(A) == 0
fprintf('\n Rank = %7.3g\n',rank(A))
error('The
coefficient matrix is singular.')
end
unit = diag(ones( 1 , n)); % Unit matrix
order = [1 : n]; %
Order of unknowns
aug = [A c]; % Augmented matrix
% Gauss -
for k = 1 : n
pivot = abs(aug(k , k));
prow = k;
pcol
= k;
% Locating the
maximum pivot element
for row = k : n
for col = k : n
if abs(aug(row , col)) > pivot
pivot =
abs(aug(row , col));
prow =
row;
pcol = col;
end
end
end
% Interchanging the
rows
pr = unit;
tmp
= pr(k , :);
pr(k , : ) =
pr(prow , : );
pr(prow , : ) = tmp;
aug
= pr * aug;
% Interchanging the
columns
pc = unit;
tmp
= pc(k , : );
pc(k , : ) = pc(pcol , : );
pc(pcol , : ) = tmp;
aug(1 : n , 1 : n)
= aug(1 : n , 1 : n) * pc;
order = order * pc; % Keep track of the column interchanges
% Reducing the
elements above and below diagonal to zero
lk
= unit;
for m = 1 : n
if m == k
lk(m , k) = 1 / aug(k
, k);
else
lk(m , k) = - aug(m , k) / aug(k , k);
end
end
aug
= lk * aug;
end
x = zeros(n
, 1);
% Solution
for k = 1 : n
x(order(k)) = aug(k , n + 1);
end
Example Problem 9.11
A=[1, 1, -1; 6, 2, 2; -3, 4, 1];
b=[3, 2, 1];
Gauss(A,b)
ans =
0.5000
1.0000
-1.5000
0.5000
1.0000
-1.5000
IV. LU Decomposition
Separates elimination from manipulation of the right hand side.
Strategy:
1.
Decompose A=LU (L is unit lower, U is upper triangular)
2.
Substitution steps
a.
Solve Ld=b (forward substitution)
b.
Solve Ux=d (back
substitution)
Advantage:
If you have many right hand sides, you only have to compute A=LU once, so you
can save computation.
The
Matrices L and U can be stored in the space occupied by the matrix A. The matrix L is the matrix of factors used
for reduction and U is the result of Gaussian elimination. These factors are put in the place of the
zeroes that result from elimination.
Nomenclature:
Doolittle
decomposition L is unit lower
triangular
Crout decomposition U is unit
upper triangular
[L,U,P]=lu(A)
1.0000 0 0
-0.5000 1.0000 0
0.1667 0.1333 1.0000
U =
6.0000 2.0000 2.0000
0 5.0000 2.0000
0 0 -1.6000
P =
0 1
0
0 0 1
1 0
0
[L,U]=lu(A)
0.1667 0.1333 1.0000
1.0000 0 0
-0.5000 1.0000 0
U =
6.0000 2.0000 2.0000
0 5.0000 2.0000
0 0 -1.6000
d=L\b'
x=U\d
3
2 1
d =
2.0000
2.0000
2.4000
x =
0.5000
1.0000
-1.5000
V. Matrix Inverses, Norms and Condition Numbers
Efficient
ways to compute inverses
can be written as a series of matrix
vector products whose right hand sides are column vectors of I, e.g,:
which can be solved by Gaussian elimination, Gauss Jordan elimination or LU decomposition.
Matrix Norms
Measure
of the size; magnitude of the transformation
2-norm,
1-norm or the 
-norm
row sum norm
column sum norm
Euclidean norm
norm(A,1)
norm(A,2)
10
ans =
10
ans =
6.9099
Condition number
cond(A) = 
always greater than or equal to 1
magnification of relative error in solving Ax=b
Example:
If
the elements of A are known to t digits of precision and the condition
number is 10c, then the solution is only valid to t-c digits.
4.7882