Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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!

share|improve this question
2  
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. –  fbloni Jun 2 at 20:18

2 Answers 2

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$.

share|improve this answer

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!

share|improve this answer
    
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

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.