# Implement a program in Matlab for LU decomposition with pivoting

I need to write a program to solve matrix equations Ax=b where A is an nxn matrix, and b is a vector with n entries using LU decomposition. Unfortunately I'm not allowed to use any prewritten codes in Matlab. I am having problems with the first part of my code where i decompose the matrix in to an upper and lower matrix.

I have written the following code, and cannot work out why it is giving me all zeros on the diagonal in my lower matrix. Any improvements would be greatly appreciated.

function[L R]=LR2(A)

%Decomposition of Matrix AA: A = L R

z=size(A,1);

L=zeros(z,z);

R=zeros(z,z);

for i=1:z

% Finding L

for k=1:i-1

L(i,k)=A(i,k);

for j=1:k-1

L(i,k)= L(i,k)-L(i,j)*R(j,k);

end

L(i,k) = L(i,k)/R(k,k);

end

% Finding R

for k=i:z

R(i,k) = A(i,k);

for j=1:i-1

R(i,k)= R(i,k)-L(i,j)*R(j,k);

end

end

end

R

L

end


I know that i could simply assign all diagonal components of L to be 1, but would like to understand what the problem is with my program!

I am also wondering how to change this program to include pivoting. I understand I need to say that if a diagonal element is equal to zero something needs to be changed. How would I go about this? Thanks in advance!

-
You have assigned $L$ to be a zero matrix and the for loop for $i^{th}$ row of $L$ runs from the first column to the $(i-1)^{th}$ column. The entry $L(i,i)$ is still left unchanged. So the diagonal entries still remain zero. – user17762 Jan 8 '11 at 22:29
@Sivaram so would it be "cheating" to simply assign all of my L(i,i)s to be 1? Thanks alot :) – Lucy Marshall Jan 8 '11 at 22:32
No it wouldn't. The best way is to decompose $L$ and $U$ directly in $A$. Since you know there will be ones on the diagonal, for huge matrix dimensions saving the previous known zeros and ones is a waste of memory space. – TheWaveLad Jun 2 '14 at 20:18

For backward and forward elimination I used

% Now use a vector y to solve 'Ly=b'
for j=1:z

for k=1:j-1
b(j)=b(j)-L(j,k)*b(k);
end;
b(j) =b(j)/L(j,j);
end;

% Now we use this y to solve Rx = y


x = zeros( z, 1 );

for i=z:-1:1
x(i) = ( b(i) - R(i, :)*x )/R(i, i);
end


I put that into your code, and it works ;) For helping you with pivot strategies it would be helpful to know what strategie you want/have to use as there are various versions. One simple version would be to just swap rows such that the diagonal element $a_{ii} \neq 0$ for all $i = 1, \dots, z$.

-

function[L R x]=LR2(A,b)

% This program will find a solution to Ax=b using first giving the decomposition of the matrix into L and R and then solving.

% Part 1 - Is this matrix square and nonsingular? give an error if not

[z y]=size(A);

if (z ~= y )

disp ( 'LR2 error: Matrix must be square' );

return;

end;

if det(A)==0

disp('L singular error');

return;


end

% Part 2 : Decomposition of matrix into L and R

L=zeros(z,z);

R=zeros(z,z);

for i=1:z

% Finding L

for k=1:i-1

L(i,k)=A(i,k);

for j=1:k-1

L(i,k)= L(i,k)-L(i,j)*R(j,k);

end

L(i,k) = L(i,k)/R(k,k);

end

% Finding R

for k=i:z

R(i,k) = A(i,k);

for j=1:i-1

R(i,k)= R(i,k)-L(i,j)*R(j,k);

end

end

end

for i=1:z

L(i,i)=1;


end

% Program shows R and L

R

L

% Now use a vector y to solve 'Ly=b'

y=zeros(z,1);

y(1)=b(1)/L(1,1);

for i=2:z

y(i)=-L(i,1)*y(1);

for k=2:i-1

y(i)=y(i)-L(i,k)*y(k);

y(i)=(b(i)+y(i))/L(i,i);

end;


end;

% Now we use this y to solve Rx = y

x=zeros(z,1);

x(1)=y(1)/R(1,1);

for i=2:z

x(i)=-R(i,1)*x(1);

for k=i:z

x(i)=x(i)-R(i,k)*x(k);

x(i)=(y(i)+x(i))/R(i,i);

end;

end


%print x x end

I have solved the original question. However there are still some problems with my solving of Ly=b and Rx=y. I also need to include pivoting. Any improvements would be greatly appreciated again!

-
It is easy to solve triangular systems by backward or forward substitutions. For pivoting you need some way to keep track of the permutations. – timur Jan 9 '11 at 4:06