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
add comment

1 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
add comment

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.