LU Факторизация с частичным поворотом Matlab - PullRequest
Новая
       7

LU Факторизация с частичным поворотом Matlab

0 голосов
/ 28 февраля 2020

Я пытаюсь реализовать факторизацию LU с частичным поворотом на PA (P - матрица перестановок, nxn) без явного обмена строками или формирования P. Я никогда не создавал код факторизации LU без явного обмена строками, и это доказывает быть трудным для меня. Кто-нибудь сможет мне помочь? 

function [L,U,P] = LUFact(A)
    %%LU Decomposition
    Ab = [A,b];
    n = length(A);
    L = eye(n);
    P = eye(n);

    %A(1,1) pivot
    %Row exchange 
    col1=Ab(:,1);
    [temp,idx] = max(col1);
    temp = Ab(1,:);
    Ab(1,:) = Ab(idx,:);
    Ab(idx,:) = temp;
    %Computation in the pivot column
    for i=2:n
        alpha = Ab(i,1) / Ab(1,1);
        L(i,1) = alpha;
        Ab(i,:) = Ab(i,:) - alpha*Ab(1,:);
    end
    %A(2,2) pivot
    %Row exchange
    col2 = Ab(2:end,2);
    [temp,idx] = max(col2);
    temp = Ab(2,:);
    Ab(2,:) = Ab(idx,:);
    Ab(idx,:) = temp;
    for i=3
        alpha = Ab(i,2) / Ab(2,2);
        L(i,2) = alpha;
        Ab(i,:) = Ab(i,:) - alpha*Ab(2,:);
    end
    U=Ab(1:n,1:n);

РЕДАКТИРОВАТЬ: 

A = [2 1 1 0; 4 3 3 1; 8 7 9 5; 6 7 9 8];
b = [-1; 3; 2; 1];
n = size(A, 1);
L = eye(n);
U = eye(n);
Ab = [A b];

%Permutation Vector
p = zeros(n, 1);
for i = 1: 1: n
    p(i) = i;
end

%Finding L
for k = 1 : i-1
    L(i,k) = A(i,k);
    for j = 1 : k-1
        L(i,k) = L(i,k) - L(i,j)*U(j,k);
    end
    L(i,k) = L(i,k) / U(k,k);
end

%Finding U
for k = i : n
    U(i,k) = A(i,k);
    for j = 1 : i-1
        U(i,k) = U(i,k) - L(i,j)*U(j,k);
    end
end

for i = 1 : n
    L(i,i) = 1;
end

%Singularity
for z = 1 : n
if A(z,z) <= 0.000000001
    disp('Matrix A is singular.');
    disp('L =');
    disp(L);
    disp('U = ');
    disp(U);
else
%%Gaussian Elimination
x = zeros(n,1);
for k = 1 : 1 : n
    max = abs(A( p(k), p(k) ) );
    max_pos = k;
    if abs(A( p(1), p(k) ) ) > max
        max = abs(A(p(1), p(k)));
        max_pos = 1;
    end 
end

%Replace p_k element with max p element
temp_p = p;
p(k) = temp_p(max_pos);
p(max_pos) = temp_p(k);

%Elimination
for i = 1 : 1 : n
    if i~= k
        alpha = A(p(i), k) / A(p(k), k);
        for j = k : 1 : n
            A(p(i), j) = A(p(i), j) - alpha*A(p(k), j);
        end
    end
end

%Back-Substitution
for i = n: -1: 1
    x(i) = ( Ab(i, end) - Ab(i, i+1 : n)*x(i+1 : n) / Ab(i, i));
end

    disp('x = ');
    disp(x);
    disp('L = ');
    disp(L);
    disp('U = ');
    disp(U);
    disp(A);

end
end
...