部分旋转的LU分解是如何工作的？

Question

我使用的方法是，最初将L的主对角线上的元素设置为1(我认为这是杜利特尔的方法，但不确定，因为我看到它的名称不同)。我知道有大量的文献、论文和书籍，但我找不到不使用高斯消元法来寻找L的简单例子。

Answer 1

与完全透视相比，部分透视仅使用行交换，而完全透视也可透视列。如下图和代码所示的部分旋转的主要目的是交换行以找到最大的u，以避免在for循环中被一个非常小的1除以，这将导致一个很大的条件数。

如果你尝试LU分解和一些病态矩阵的天真实现，比如一些任意的对角占优矩阵，它应该会爆炸。

​

​

function [L,U,P] = my_lu_piv(A)
n = size(A,1);
I = eye(n);
O = zeros(n);
L = I;
U = O;
P = I;
    function change_rows(k,p)
        x = P(k,:); P(k,:) = P(p,:); P(p,:) = x;
        x = A(k,:); A(k,:) = A(p,:); A(p,:) = x;
        x = v(k); v(k) = v(p); v(p) = x;
    end

    function change_L(k,p)
        x = L(k,1:k-1); L(k,1:k-1) = L(p,1:k-1);
        L(p,1:k-1) = x;
    end
for k = 1:n
    if k == 1, v(k:n) = A(k:n,k);
    else
        z = L(1:k-1,1:k -1)\ A(1:k-1,k);
        U(1:k-1,k) = z;
        v(k:n) = A(k:n,k)-L(k:n,1:k-1)*z;
    end
    if k<n
        x = v(k:n); p = (k-1)+find(abs(x) == max(abs(x))); % find index p
        change_rows(k,p);
        L(k+1:n,k) = v(k+1:n)/v(k);
        if k > 1, change_L(k,p); end
    end
    U(k,k) = v(k);
end
end