How to get the correct answer for LU Decomposition with partial pivoting?

256 views (last 30 days)
I am trying to solve a matrix using LU decomposition with crout's method with partial pivoting. I found this code online at this website.
I added few lines to get the my x from Ax=b equation but it seems that the code is not giving me the correct matrices at the end. Sometimes it even gives me error. I can't figure out where in the code is giving a problem.
A = [0 1 -1 1 0 0 0 ;
1 1 0 0 0 0 0 ;
0 0 0 -1 0 1 0 ;
0 1 0 0 -1 -0.3 0 ;
0 0 0 0 0 0 1 ;
0 1 0 0 0 0 -1 ;
0 1 0 0 1 0 0 ] % given matrix to solve
b = [0; 120; 0; 100; 0; 0; 0]
[n,n]=size(A);
L=eye(n); P=L; U=A;
for k=1:n
[pivot m]=max(abs(U(k:n,k)));
m=m+k-1;
if m~=k
% interchange rows m and k in U
temp=U(k,:);
U(k,:)=U(m,:);
U(m,:)=temp;
% interchange rows m and k in P
temp=P(k,:);
P(k,:)=P(m,:);
P(m,:)=temp;
if k >= 2
temp=L(k,1:k-1);
L(k,1:k-1)=L(m,1:k-1);
L(m,1:k-1)=temp;
end
end
for j=k+1:n
L(j,k)=U(j,k)/U(k,k);
U(j,:)=U(j,:)-L(j,k)*U(k,:);
end
end
L,U
Here is the code I added to get x.
Y=zeros(n,1);
Y(1)=b(1)/L(1,1);
for j=2:n
Y(j)=(b(j)-L(j,1:j-1)*Y(1:j-1))/L(j,j);
end
Y
X=zeros(n,1);
X(n)=Y(n)/U(n,n);
for j=n-1:-1:1
X(j)=(Y(j)-U(j,j+1:n)*X(j+1:n))/U(j,j);
end
X
When I solve it with x=A\b, i am getting completely different answers. What should I do?

Accepted Answer

Jan
Jan on 25 Feb 2022
Edited: Jan on 25 Feb 2022
You forgot to apply the permutation matrix P to b. The results are correct:
P * L * U - A % == zeros, fine
b = P * b;
Then A\b replies the same as your X.
A simplification of your code:
% Replace:
temp=U(k,:);
U(k,:)=U(m,:);
U(m,:)=temp;
% by:
U([m, k], :) = U([k, m], :);

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by