5

部分ピボットを使用して独自のLU分解を実装しようとしています。[L, U, P] = lu(A)私のコードは以下のとおりで、明らかに正常に機能していますが、一部の行列では、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

これが私がテストした2つのマトリックスです。最初のものは正しいですが、2番目のものはいくつかの要素が反転しています。

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];

アップデート

コードをチェックしていくつかのバグを修正しましたが、部分的なピボットにはまだ何かが欠けています。最初の列では、最後の2行が常に反転しています（matlabのlu（）の結果と比較して）

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
4

3 に答える 3

3

マトリックスPIにスワップがあった場合、マトリックスLもスワップする必要があることを忘れました。したがって、Pをスワップした後に次の行を追加するだけで、すべてが正常に機能します。

L([k r],:) = L([r k],:);
于 2013-03-09T10:34:45.190 に答える
3

両方の機能が正しくありません。これが正しいものです。

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
于 2014-08-06T16:41:30.427 に答える
1

私の答えはここにあります：

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
于 2014-03-04T20:55:05.157 に答える