How to get the correct answer for LU Decomposition with partial pivoting?
65 views (last 30 days)
Show older comments
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?
0 Comments
Accepted Answer
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)
See Also
Categories
Find more on Linear Algebra in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!