Below is the output of a MATLAB session on the example we discussed in class.
>> % In this example we work out the LU decomposition to solve the system
>> % on page 124 of the text book.
>> A = [4 -2 1; -3 -1 4; 1 -1 3]
A =
4 -2 1
-3 -1 4
1 -1 3
>> b = [15; 8; 13]
b =
15
8
13
>> M21 = eye(3); % multiplier matrix to make A(2,1) zero
>> M21(2,1) = -A(2,1)/A(1,1)
M21 =
1.0000 0 0
0.7500 1.0000 0
0 0 1.0000
>> U21 = M21*A
U21 =
4.0000 -2.0000 1.0000
0 -2.5000 4.7500
1.0000 -1.0000 3.0000
>> M31 = eye(3); % multiplier matrix to make A(3,1) zero
>> M31(3,1) = -A(3,1)/A(1,1)
M31 =
1.0000 0 0
0 1.0000 0
-0.2500 0 1.0000
>> U31 = M31*U21
U31 =
4.0000 -2.0000 1.0000
0 -2.5000 4.7500
0 -0.5000 2.7500
>> M32 = eye(3); % multiplier to make U31(3,2) zero
>> M32(3,2) = -U31(3,2)/U31(3,3)
M32 =
1.0000 0 0
0 1.0000 0
0 0.1818 1.0000
>> U = M32*U31
U =
4.0000 -2.0000 1.0000
0 -2.5000 4.7500
0 -0.9545 3.6136
>> M32(3,2) = -U31(3,2)/U31(2,2)
M32 =
1.0000 0 0
0 1.0000 0
0 -0.2000 1.0000
>> U = M32*U31
U =
4.0000 -2.0000 1.0000
0 -2.5000 4.7500
0 0 1.8000
>> L = eye(3); % collect all multipliers with inverse sign
>> L(2,1) = -M21(2,1);
>> L(3,1) = -M31(3,1);
>> L(3,2) = -M32(3,2);
>> L
L =
1.0000 0 0
-0.7500 1.0000 0
0.2500 0.2000 1.0000
>> L*U
ans =
4 -2 1
-3 -1 4
1 -1 3
>> det(A)
ans =
-18
>> det(U)
ans =
-18
In practice, we can store the coefficients of L and U in the matrix A, as is done by the .m file below:
function A = lufac(A)
% LU factorization without pivoting
n = size(A,2);
for j = 1:n-1
for i = j+1:n % store multipliers
A(i,j) = A(i,j)/A(j,j);
end;
for i = j+1:n % eliminate
for k = j+1:n
A(i,k) = A(i,k) - A(i,j)*A(j,k);
end;
end;
end;
Then we continue our MATLAB session:
>> % Now we test our little implementation on this example:
>> A
A =
4 -2 1
-3 -1 4
1 -1 3
>> lufac(A)
ans =
4.0000 -2.0000 1.0000
-0.7500 -2.5000 4.7500
0.2500 0.2000 1.8000
>> % Observe that the matrix on return contains both L and U.
The original motivation for the LU factorization was of course to solve the linear system A*x = b. This problem is now reduced to the solving of two triangular systems:
>> % Now we solve the system A*x = b with backsubstitution
>> % First consider L*y = b (y = U*x):
>> y(1) = b(1)
y =
15
>> y(2) = b(2) - L(2,1)*y(1)
y =
15.0000 19.2500
>> y(3) = b(3) - L(3,1)*y(1) - L(3,2)*y(2)
y =
15.0000 19.2500 5.4000
>> L*y' % check
ans =
15
8
13
>> x(3) = y(3)/U(3,3)
x =
0 0 3.0000
>> x(2) = (y(2) - U(2,3)*x(3))/U(2,2)
x =
0 -2.0000 3.0000
>> x(1) = (y(1) - U(1,2)*x(2) - U(1,3)*x(3))/U(1,1)
x =
2.0000 -2.0000 3.0000
>> A*x'
ans =
15.0000
8.0000
13.0000
>> b
b =
15
8
13
Formally, we define the forward and backward back substitution by the .m files "forsubs" and "backsubs":
function y = forsubs(A,b)
%
% Solves A*y = b with forward back substitution using the lower
% triangular matrix L in the LU factorization of A.
%
n = size(A,2);
for i = 1:n
yy(i) = b(i);
for j = 1:i-1
yy(i) = yy(i) - A(i,j)*yy(j);
end;
end;
y = yy';
function x = backsubs(A,y)
%
% Solves U*x = y with backward back substitution using the upper
% triangular matrix U in the LU decomposition of A.
%
n = size(A,2);
for i = n:-1:1
xx(i) = y(i);
for j = i+1:n
xx(i) = xx(i) - A(i,j)*xx(j);
end;
xx(i) = xx(i)/A(i,i);
end;
x = xx';
We use our three .m files on our running example in the MATLAB session below:
>> % we solve the system with forward and backward back substitution:
>> A = [4 -2 1; -3 -1 4; 1 -1 3] % coefficient matrix
A =
4 -2 1
-3 -1 4
1 -1 3
>> b = [15;8;13] % right hand side vector
b =
15
8
13
>> LU = lufac(A) % LU contains LU factorization of A
LU =
4.0000 -2.0000 1.0000
-0.7500 -2.5000 4.7500
0.2500 0.2000 1.8000
>> y = forsubs(LU,b) % L*y = b forward back substitution
y =
15.0000
19.2500
5.4000
>> x = backsubs(LU,y) % U*x = y backward back substitution
x =
2.0000
-2.0000
3.0000
>> A*x % check whether A*x equals b
ans =
15.0000
8.0000
13.0000