用部分旋转Matlab进行LU分解 [英] LU decomposition with partial pivoting Matlab

问题描述

我正在尝试通过部分数据透视实现我自己的LU分解.我的代码在下面，显然可以正常工作，但是对于某些矩阵，与matlab中内置的[L, U, P] = lu(A)函数进行比较时，得出的结果不同

I am trying to implement my own LU decomposition with partial pivoting. My code is below and apparently is working fine, but for some matrices it gives different results when comparing with the built-in [L, U, P] = lu(A) function in matlab

谁能发现哪里出问题了?

Can anyone spot where is it wrong?

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

这是我测试过的两个矩阵.第一个是正确的，而第二个是一些倒置的元素.

Here are the two matrices I've tested with. The first one is correct, whereas the second has some elements inverted.

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

更新

我已经检查了我的代码并更正了一些错误，但是部分枢轴仍然缺少一些东西.在第一列中，最后两行始终是反转的(与matlab中lu()的结果相比)

I have checked my code and corrected some bugs, but still there's something missing with the partial pivoting. In the first column the last two rows are always inverted (compared with the result of lu() in matlab)

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

推荐答案

我忘记了，如果矩阵PI中存在交换，则PI也必须交换矩阵L.因此，只需在交换P之后添加下一行，一切就可以了很棒.

I forgot that If there was a swap in matrix P I had to swap also the matrix L. So just add the next line after after swapping P and everything will work excellent.

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

这篇关于用部分旋转Matlab进行LU分解的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！