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

27 次查看（过去 30 天）
Kenny Gee
Kenny Gee 2022-2-25
评论： Kenny Gee 2022-2-26
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.
https://gist.github.com/02015678/353ea94c27e4c5f89c9d
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?

请先登录，再回答此问题。

采纳的回答

Jan
Jan 2022-2-25
编辑：Jan 2022-2-25
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], :);

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Linear Algebra 的更多信息

标签

产品


版本

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by