KIL1005: Numerical Methods for Engineering
Semester 2, Session 2019/2020
Lecture 4: LU Decomposition and Matrix Inversion
7 May 2020
Dr. Mahar Diana Hamid
Department of Chemical Engineering
University Malaya
LU Decomposition and Matrix Inversion
› Provides an efficient way to compute matrix inverse
by separating the time consuming elimination of the
Matrix [A] from manipulations of the right-hand side
{B}.
› Gauss elimination, in which the forward elimination
comprises the bulk of the computational effort, can be
implemented as an LU decomposition.
2
If
L- lower triangular matrix
U- upper triangular matrix
Then,
[A]{X}={B} can be decomposed into two matrices [L] and
[U] such that
[L][U]=[A]
[L][U]{X}={B}
Similar to first phase of Gauss elimination, consider
[U]{X}={D}
[L]{D}={B}
– [L]{D}={B} is used to generate an intermediate vector {D}
by forward substitution
– Then, [U]{X}={D} is used to get {X} by back substitution.
3
Fig.10.1
4
LU decomposition
› requires the same total FLOPS (floating point operations per
second) as for Gauss elimination.
› Saves computing time by separating time-consuming
elimination step from the manipulations of the right hand side.
› Provides efficient means to compute the matrix inverse
› Note: In computing, floating point operations per
second(FLOPS, flops or flop/s) is a measure of computer
performance, useful in fields of scientific computations that
require floating-point calculations. For such cases it is a more
accurate measure than measuring instructions per second.
5
Error Analysis and System Condition
› Inverse of a matrix provides a means to test whether
systems are ill-conditioned.
Vector and Matrix Norms
• Norm is a real-valued function that provides a
measure of size or “length” of vectors and
matrices. Norms are useful in studying the error
behavior of algorithms.
6
A simple example is a vector in three - dimensiona l Euclidean space
that can be represente d as
F a b c
where a, b, and c are the distances along x, y, and z axes, repectivel y.
7
Figure 10.6
8
› The length of this vector can be simply computed as
F e a b c 2 2 2
Length or Euclidean norm of [F]
• For an n dimensional vector
X x1 x2 xn
a Euclidean norm is computed as
n
Xe i
x 2
i 1
For a matrix [A]
n n Frobenius norm
Ae
i 1 j 1
ai2, j
9
› Frobenius norm provides a single value to quantify the
“size” of [A].
Matrix Condition Number
› Defined as
Cond A A A 1
•For a matrix [A], this number will be greater
than or equal to 1.
10
X A
Cond A
X A
• That is, the relative error of the norm of the
computed solution can be as large as the relative
error of the norm of the coefficients of [A]
multiplied by the condition number.
• For example, if the coefficients of [A] are
known to t-digit precision (rounding errors~10 -t)
and Cond [A]=10c, the solution [X] may be valid
to only t-c digits (rounding errors~10c-t).
11
MATLAB code for Gauss Elimination
Initialization and input descriptions
clc
close all
disp('Gauss Elmination method to solve "n" equations for in the format of:
"[A][x]=[b]" ');
A=input ('Enter the coefficients in the equations in matrix form i.e. Input
[A] matrix: ');
b=input ('Enter the right hand side column vector [b]: ');
n=length(b);
tic;
Creating an Upper Triangular Matrix
for j=1:n-1
if A(j,j)==0 Pivoting Operation because we do not want zero denominator
m=j; Pivot the row
for m=m+1:n For the subsequent row of m+1
if A(m,j)==0
continue
end
break
end
B=A(j,:);
C=b(j);
A(j,:)= A(m,:);
b(j)=b(m);
A(m,:)= B;
b(m)=C;
end
Forward elimination
for i=1+ss:n-1
fac = A(i+1,j)/A(j,j);
A(i+1,:) = A(i+1,:) - fac*A(j,:);
b(i+1) = b(i+1) - fac*b(j);
end
ss=ss+1;
end
Backward substitution
x(n)=b(n)/A(n,n);
for i=n-1:-1:1 syntax for count down from n-1 until 1
coefsum=0;
for j=i+1:n
coefsum=coefsum+A(i,j)*x(j);
end
x(i)=(1/A(i,i))*(b(i)-coefsum);
Displaying the output
toc
disp('Upper Triangular Form of Matrix [A] =');
disp(A)
disp('Solution of the system of equations is :');
disp(x')