2013-03-08 55 views
3

我想实现我自己的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 

回答

2

我忘了,如果有在矩阵PI交换还曾因此交换矩阵L.(用MATLAB路()的结果进行比较)只需在交换P之后添加下一行,并且一切都会很好。

L([k r],:) = L([r k],:); 
+0

你能确切地指定你放这条线的位置吗?我测试了你的代码,并且L矩阵有问题 – CTZStef 2014-01-10 13:52:59

+0

应该在'P([k r],:) = P([r k],:); – BRabbit27 2014-01-10 15:10:22

0

我的答案就在这里:

function [L, U, P] = LU_pivot(A) 
[n, n1] = size(A); L=eye(n); P=eye(n); U=A; 
for j = 1:n 
    [pivot m] = max(abs(U(j:n, j)));  
    m = m+j-1; 
    if m ~= j 
    U([m,j],:) = U([j,m], :); % interchange rows m and j in U 
    P([m,j],:) = P([j,m], :); % interchange rows m and j in P 
    if j >= 2; % very_important_point 
     L([m,j],1:j-1) = L([j,m], 1:j-1); % interchange rows m and j in columns 1:j-1 of L 
    end; 
    end 
    for i = j+1:n  
    L(i, j) = U(i, j)/U(j, j); 
    U(i, :) = U(i, :) - L(i, j)*U(j, :); 
    end 
end 
1

这两个功能是不正确的。 这是正确的。

function [L, U, P] = LU_pivot(A) 
    [m, n] = size(A); L=eye(n); P=eye(n); U=A; 
    for k=1:m-1 
     pivot=max(abs(U(k:m,k))) 
     for j=k:m 
      if(abs(U(j,k))==pivot) 
       ind=j 
       break; 
      end 
     end 
     U([k,ind],k:m)=U([ind,k],k:m) 
     L([k,ind],1:k-1)=L([ind,k],1:k-1) 
     P([k,ind],:)=P([ind,k],:) 
     for j=k+1:m 
      L(j,k)=U(j,k)/U(k,k) 
      U(j,k:m)=U(j,k:m)-L(j,k)*U(k,k:m) 
     end 
     pause; 
    end 
end