# L U decomposition

984 views (last 30 days)
John on 15 Feb 2011
Answered: Mahesh Prajapati on 21 Sep 2020
Below I have a code written for solving the L U decomposition of a system of equations however I need my code to just output the answers with this format it outputs the variables in the matrix for example i need the function to output x [1;2;3;4] any suggestions?
function[L,U,X]=LU_Parker(A,B)
[m n]=size(A);
if (m ~= n )
disp ( 'LR2 error: Matrix must be square' );
return;
end;
% Part 2 : Decomposition of matrix into L and U
L=zeros(m,m);
U=zeros(m,m);
for i=1:m
% 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)*U(j,k);
end
L(i,k) = L(i,k)/U(k,k);
end
% Finding U
for k=i:m
U(i,k) = A(i,k);
for j=1:i-1
U(i,k)= U(i,k)-L(i,j)*U(j,k);
end
end
end
for i=1:m
L(i,i)=1;
end
% Program shows U and L
U
L
% Now use a vector y to solve 'Ly=b'
y=zeros(m,1);
y(1)=B(1)/L(1,1);
for i=2:m
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 Ux = y
x=zeros(m,1);
x(1)=y(1)/U(1,1);
for i=2:m
x(i)=-U(i,1)*x(1);
for k=i:m
x(i)=x(i)-U(i,k)*x(k);
x(i)=(y(i)+x(i))/U(i,i);
end;
end

#### 1 Comment

Daz on 3 Feb 2015
Aren't you going to get a divide by 0 error? At the very end of what I quoted, you have L(i,k) = L(i,k)/U(k,k);
But the first time through, U is a zero matrix.
L=zeros(m,m); U=zeros(m,m); for i=1:m % 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)*U(j,k); end L(i,k) = L(i,k)/U(k,k); end

Oleg Komarov on 15 Feb 2011
Matlab is case-sensitive, if you want to store the output of x then in the first line change X to lowercase.
Oleg

John on 15 Feb 2011
I tried this but it still outputs my answer the same way, I originally had it as a lowercase x but I changed it to upper case after I realized it didn't change anything.

#### 1 Comment

Oleg Komarov on 15 Feb 2011
Then can you post the undesired result and the desired one? It's not very clear from your first description.

Tan Edwin on 15 Feb 2011
Maybe u can try adding X=x to allow it to ouput the values of x?
not sure if this is what u want.
edwin

John on 15 Feb 2011
Yes, redefining the x like you said allowed the function to output what I was needing, however I must have an error in my coding because I inputed the following matrices and got the following answer but I am getting a 0 for one of the answers which should not be there. Any possible solutions?
INPUT
A=[ 6 0 0 0 0; 0 1 0 -2 0; 1 0 -3 0 0; 0 8 -4 -3 -2; 0 2 0 0 -1];
B=[1;0;0;1;0];
LU_Parker(A,B)
Output:
X =
0.1667
0
0.0432
0.1841
1.7778
ans =
1.0000 0 0 0 0
0 1.0000 0 0 0
0.1667 0 1.0000 0 0
0 8.0000 1.3333 1.0000 0
0 2.0000 0 0.3077 1.0000

Mohamed Said Attia on 4 Jun 2011
*there is a problem with the way you are solving the equation to get y & x try*
% Now use a vector y to solve 'Ly=b'
y=zeros(m,1); % initiation for y
y(1)=B(1)/L(1,1);
for i=2:m
%y(i)=B(i)-L(i,1)*y(1)-L(i,2)*y(2)-L(i,3)*y(3);
y(i)=-L(i,1)*y(1);
for k=2:i-1
y(i)=y(i)-L(i,k)*y(k);
end;
y(i)=(B(i)+y(i))/L(i,i);
end;
y
% Now we use this y to solve Ux = y
x=zeros(m,1);
x(m)=y(m)/U(m,m);
i=m-1;
q=0;
while (i~= 0)
x(i)=-U(i,m)*x(m);
q=i+1;
while (q~=m)
x(i)=x(i)-U(i,q)*x(q);
q=q+1;
end;
x(i)=(y(i)+x(i))/U(i,i);
i=i-1;
end;
x

Mohamed Said Attia on 4 Jun 2011
and when you call the function from matlab use
[L,U,X]=LU_Parker(A,B) not LU_Parker(A,B)

#### 1 Comment

Walter Roberson on 4 Jun 2011
Not really relevant: if you do not specify output variables and do not put a semi-colon at the end of the line, you will get
ans =
for each of the output variables, in left-to-right order.

Mahesh Prajapati on 21 Sep 2020
any suggestions? function[L,U,X]=LU_Parker(A,B) [m n]=size(A); if (m ~= n ) disp ( 'LR2 error: Matrix must be square' ); return; end; % Part 2 : Decomposition of matrix into L and U L=zeros(m,m); U=zeros(m,m); for i=1:m % 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)*U(j,k); end L(i,k) = L(i,k)/U(k,k); end % Finding U for k=i:m U(i,k) = A(i,k); for j=1:i-1 U(i,k)= U(i,k)-L(i,j)*U(j,k); end end end for i=1:m L(i,i)=1; end % Program shows U and L U L % Now use a vector y to solve 'Ly=b' y=zeros(m,1); y(1)=B(1)/L(1,1); for i=2:m 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 Ux = y x=zeros(m,1); x(1)=y(1)/U(1,1); for i=2:m x(i)=-U(i,1)*x(1); for k=i:m x(i)=x(i)-U(i,k)*x(k); x(i)=(y(i)+x(i))/U(i,i); end; end