> V1 = [-1;2;3];
> V2 = [3;5;2];
> V3 = [8;1;9];
> A = [V1,V2,V3]
A =
-1 3 8
2 5 1
3 2 9
Here is the LU decomposition code in MatLab
function [L,U] = GE(A)
%
% A is nxn matrix
% L is nxn lower triangular
% U is nxn upper triangular
%
% We compute the LU decomposition of A using
% Gaussian Elimination
%
[n,n] = size(A);
for k=1:n-1
% find multiplier
if (A(k,k) ~= 0)
A(k+1:n,k) = A(k+1:n,k)/A(k,k);
% zero out column
A(k+1:n,k+1:n) = A(k+1:n,k+1:n) - A(k+1:n,k)*A(k,k+1:n);
else
break;
end
end
L = eye(n,n) + tril(A,-1);
U = triu(A);
end
We are going to explain all of this code.
For our example here, the matrix is 3x3
and do find the LU decomposition we start in column 1.
Step One:
this is k = 1 in the code above
Is A(1,1) != 0 Yes so go ahead an
zero out the rest of column one
Look at A in row 2 and row 3, col 1
> A(2:3,1)
ans =
2
3
Store the multipliers in this part of A
This rewrites A
> A(2:3,1) = A(2:3,1)/A(1,1)
A =
-1 3 8
-2 5 1
-3 2 9
Zero out column 1 below A(1,1)
R2new = R2 + 2 *R1 ==> multiplier is -1 stored in old A(2,1)
A(2,2:3) = A(2,2:3) - A(2,1)*A(1,2:3)
R3new = R3 + 3 *R1 ==> multiplier is -3 stored in old A(3,1)
A(3,2:3) = A(3,2:3) - A(3,1)*A(1,2:3)
The above two steps can be written more
succinctly using MatLab as
> A(2:3,2:3) = A(2:3,2:3) - A(2:3,1)*A(1,2:3)
A =
-1 3 8
-2 11 17
-3 11 33
Step Two:
Now work on col 2
This is k = 2 in the code above
A(2,2) is not zero so ok to divide with it
Store multiplier
> A(3:3,2) = A(3:3,2)/A(2,2)
A =
-1 3 8
-2 11 17
-3 1 33
Zero out column 2 below A(2,2)
R3new = R3 -1 *R2 ==> multipier is 1 stored in old A(3,2)
A(3:3,3:3) = A(3:3,3:3) - A(3:3,2)*A(2.3:3)
The above step can be written in MatLab like this
> A(3:3,3:3) = A(3:3,3:3) - A(3:3,2)*A(2,3:3)
A =
-1 3 8
-2 11 17
-3 1 16
so
Now eye(3,3) in Matlab gives the 3x3 identity matrix
and tril(A,-1) extracts the lower triangular part of A
skipping the main diagonal -- that is why we have a -1
in the tril(A,-1)
So the L we want is
> L = eye(3,3) + tril(A,-1)
The U we want is
> U = triu(A)
and we don't use the -1 option because we want the
main diagonal now.
So we have
> L
L =
1 0 0
-2 1 0
-3 1 1
> U
U =
-1 3 8
0 11 17
0 0 16
We could repackage what we did like this
for k = 1:2
% find multiplier
if A(k,k) ~= 0
% store multipliers
A(k+1:3,k) = A(k+1:3,k)/A(k,k);
% zero out column
A(k+1:3,k+1:3) = A(k+1:3,k+1:3) - A(k+1:3,k)*A(k,k+1:3);
else
break;
end
end
% eye(3,3) gives the 3x3 identity matrix
% tril(A,-1) finds the lower triangular part of A below
% the main diagonal: the L we want is the sum of these
L = eye(n,n) + tril(A,-1);
% triu(A) finds the upper triangular part of A including the
% main diagonal.
U = triu(A);
There is no reason we couldn't do this for n > 3.
We would have
[n,n] = size(A);
for k=1:n-1
% find multiplier
if (A(k,k) ~= 0)
A(k+1:n,k) = A(k+1:n,k)/A(k,k);
% zero out column
A(k+1:n,k+1:n) = A(k+1:n,k+1:n) - A(k+1:n,k)*A(k,k+1:n);
else
break;
end
end
L = eye(n,n) + tril(A,-1);
U = triu(A);
If we use this as a MatLab function, we have
> A = [V1,V2,V3];
> [L2,U2] = GE(A);
> L2
L2 =
1 0 0
-2 1 0
-3 1 1
> U2
U2 =
-1 3 8
0 11 17
0 0 16
which is what we had before.