我想实现我自己的LU分解与部分pivoting。我的代码在下面,显然工作正常,但对于某些矩阵,当与内置的[L, U, P] = lu(A)
函数进行比较时,它会给出不同的结果LU分解与部分旋转Matlab
任何人都可以发现它错在哪里?
function [L, U, P] = lu_decomposition_pivot(A)
n = size(A,1);
Ak = A;
L = zeros(n);
U = zeros(n);
P = eye(n);
for k = 1:n-1
for i = k+1:n
[~,r] = max(abs(Ak(:,k)));
Ak([k r],:) = Ak([r k],:);
P([k r],:) = P([r k],:);
L(i,k) = Ak(i,k)/Ak(k,k);
for j = k+1:n
U(k,j-1) = Ak(k,j-1);
Ak(i,j) = Ak(i,j) - L(i,k)*Ak(k,j);
end
end
end
L(1:n+1:end) = 1;
U(:,end) = Ak(:,end);
return
这是我测试过的两个矩阵。第一个是正确的,而第二个是倒转的。
A = [1 2 0; 2 4 8; 3 -1 2];
A = [0.8443 0.1707 0.3111;
0.1948 0.2277 0.9234;
0.2259 0.4357 0.4302];
UPDATE
我检查了我的代码,并纠正了一些错误,但还是有一些用的部分回转下落不明。在第一列的最后两行总是倒
function [L, U, P] = lu_decomposition_pivot(A)
n = size(A,1);
Ak = A;
L = eye(n);
U = zeros(n);
P = eye(n);
for k = 1:n-1
[~,r] = max(abs(Ak(k:end,k)));
r = n-(n-k+1)+r;
Ak([k r],:) = Ak([r k],:);
P([k r],:) = P([r k],:);
for i = k+1:n
L(i,k) = Ak(i,k)/Ak(k,k);
for j = 1:n
U(k,j) = Ak(k,j);
Ak(i,j) = Ak(i,j) - L(i,k)*Ak(k,j);
end
end
end
U(:,end) = Ak(:,end);
return
你能确切地指定你放这条线的位置吗?我测试了你的代码,并且L矩阵有问题 – CTZStef 2014-01-10 13:52:59
应该在'P([k r],:) = P([r k],:); – BRabbit27 2014-01-10 15:10:22