MATLAB LU Разложение Частичное вращение

Question

Я пытаюсь работать с моей декомпозицией lu, в значительной степени основанной на декомпозиции LU с частичным поворотом Matlab

function [L,U,P] = lup(A)
n = length(A);
L = eye(n);
U = zeros(n);
P = eye(n);

for k=1:n-1
% find the entry in the left column with the largest abs value (pivot)
[~,r] = max(abs(A(k:end,k)));
r = n-(n-k+1)+r;

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

% from the pivot down divide by the pivot
L(k+1:n,k) = A(k+1:n,k) / A(k,k);

U(k,1:n) = A(k,1:n);
A(k+1:n,1:n) = A(k+1:n,1:n) - L(k+1:n,k)*A(k,1:n);

end
U(:,end) = A(:,end);

end

Похоже, что он работает для большинства матриц (равен функции matlab lu), однако следующая матрица, похоже, дает другие результаты:

A = [
 3    -7    -2     2
-3     5     1     0
 6    -4     0    -5
-9     5    -5    12
];

Я просто не могу понять, что происходит не так. Кажется, работает нормально на матрицах, упомянутых в связанном посте

Answer 1

Ты был довольно близко Я изменил всего три строки

for k=1:n-1 стал for k=1:n мы не делаем -1, потому что мы также хотим получить L(n,n)=u(n,n)/u(n,n)=1 с вашим методом мы оставляли это

L(k+1:n,k) = A(k+1:n,k) / A(k,k); стал L(k:n,k) = A(k:n,k) / A(k,k); потому что ты уходил L(k,k)=A(k,k)/A(k,k)=1

поскольку k+1 изменение нам не нужно начинать с единичной матрицы для L, так как мы теперь воспроизводим 1 на диагонали так L=eyes(n); стал L=zeros(n);

и завершенный код

function [L,U,P] = lup(A)
% lup factorization with partial pivoting
% [L,U,P] = lup(A) returns unit lower triangular matrix L, upper
% triangular matrix U, and permutation matrix P so that P*A = L*U.
    n = length(A);
    L = zeros(n);
    U = zeros(n);
    P = eye(n);


    for k=1:n
        % find the entry in the left column with the largest abs value (pivot)
        [~,r] = max(abs(A(k:end,k)));
        r = n-(n-k+1)+r;    

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

        % from the pivot down divide by the pivot
        L(k:n,k) = A(k:n,k) / A(k,k);

        U(k,1:n) = A(k,1:n);
        A(k+1:n,1:n) = A(k+1:n,1:n) - L(k+1:n,k)*A(k,1:n);

    end
    U(:,end) = A(:,end);

end