LU разложение квадратной матрицы устранения Matlab Гаусса

Question

Я пытаюсь создать программу, которая принимает квадратную (n-на-n) матрицу в качестве входных данных, и, если она обратима, LU будет разлагать матрицу с использованием гауссовского исключения.

Вот моя проблема: в классе мы узнали, что лучше менять строки так, чтобы ваш столбец всегда был наибольшим числом (в абсолютном значении) в своем столбце. Например, если матрица была A = [1,2;3,4] затем переключение строк это [3,4;1,2] и тогда мы можем приступить к устранению Гаусса.

Мой код работает правильно для матриц, которые не требуют изменения строки, но для тех, которые делают, это не так. Это мой код:

function newgauss(A)
    [rows,columns]=size(A);
    P=eye(rows,columns); %P is permutation matrix
    if(det(A)==0) %% determinante is 0 means no single solution
        disp('No solutions or infinite number of solutions')
        return;
    end
    U=A;
    L=eye(rows,columns);
    pivot=1;
    while(pivot<rows)
        max=abs(U(pivot,pivot));
        maxi=0;%%find maximum abs value in column pivot
        for i=pivot+1:rows
            if(abs(U(i,pivot))>max)
                max=abs(U(i,pivot));
                maxi=i;
            end
        end %%if needed then switch
        if(maxi~=0)
            temp=U(pivot,:);
            U(pivot,:)=U(maxi,:);
            U(maxi,:)=temp;
            temp=P(pivot,:);
            P(pivot,:)=P(maxi,:);
            P(maxi,:)=temp;
        end %%Grade the column pivot using gauss elimination
        for i=pivot+1:rows
            num=U(i,pivot)/U(pivot,pivot);
            U(i,:)=U(i,:)-num*U(pivot,:);
            L(i,pivot)=num;
        end
        pivot=pivot+1;
    end
    disp('PA is:');
    disp(P*A);
    disp('LU is:');
    disp(L*U);
end

Пояснение: так как мы переключаем строки, мы хотим разложить P (матрица перестановок) раз Aа не оригинал A что мы имели в качестве входных данных.

Пояснение к коду:

1. Сначала я проверяю, является ли матрица обратимой, если нет, остановитесь. Если это так, то разворот (1,1)
2. Я нахожу наибольшее число в столбце 1 и переключаю строки
3. Оцените столбец 1 с использованием исключения Гаусса, поворачивая все, кроме точки (1,1), к нулю
4. Pivot теперь (2,2), найдите наибольшее число в столбце 2... Промыть, повторить

Answer 1

Ваш код, кажется, работает хорошо из того, что я могу сказать, по крайней мере, для основных примеров A=[1,2;3,4] или же A=[3,4;1,2], Измените определение вашей функции на:

function [L,U,P] = newgauss(A)

так что вы можете вывести свои рассчитанные значения (гораздо лучше, чем при использовании disp, но это показывает правильные результаты тоже). Тогда вы увидите, что P*A = L*U, Может быть, вы ожидали L*U равному A напрямую? Вы также можете подтвердить, что вы правы с помощью Matlab's lu функция:

[L,U,P] = lu(A);
L*U
P*A

Матрицы перестановок являются ортогональными матрицами, поэтому P ⁻¹ = PT. Если вы хотите вернуться A в своем коде вы можете сделать:

P'*L*U

Точно так же, используя Matlab's lu с выводом матрицы перестановок вы можете сделать:

[L,U,P] = lu(A);
P'*L*U

(Вы также должны использовать error или же warning а не как вы используете disp в проверке определителя, но они, вероятно, не учат этому.)

Answer 2

Обратите внимание, что det Функция реализована с использованием самой декомпозиции LU для вычисления детерминанта... рекурсивный любой:)

Кроме того, в конце страницы есть напоминание, в котором предлагается использовать cond вместо det проверить на матричную особенность:

Тестирование сингулярности с использованием abs(det(X)) <= tolerance не рекомендуется, так как трудно выбрать правильный допуск. Функция cond(X) можно проверить на единичные и почти особые матрицы.

COND использует разложение по единственному значению (см. Его реализацию: edit cond.m)

Answer 3

Для тех, кто ищет это в будущем и нуждается в рабочем решении:

Код ОП не содержит логики для переключения элементов в L при создании матрицы перестановок P, Скорректированный код, который дает тот же вывод, что и Matlab lu(A) функция:

function [L,U,P] = newgauss(A)
    [rows,columns]=size(A);
    P=eye(rows,columns); %P is permutation matrix
    tol = 1E-16; % I believe this is what matlab uses as a warning level
    if( rcond(A) <= tol) %% bad condition number
        error('Matrix is nearly singular')
    end
    U=A;
    L=eye(rows,columns);
    pivot=1;
    while(pivot<rows)
        max=abs(U(pivot,pivot));
        maxi=0;%%find maximum abs value in column pivot
        for i=pivot+1:rows
            if(abs(U(i,pivot))>max)
                max=abs(U(i,pivot));
                maxi=i;
            end
        end %%if needed then switch
        if(maxi~=0)
            temp=U(pivot,:);
            U(pivot,:)=U(maxi,:);
            U(maxi,:)=temp;
            temp=P(pivot,:);
            P(pivot,:)=P(maxi,:);
            P(maxi,:)=temp;

            % change elements in L-----
            if pivot >= 2
                temp=L(pivot,1:pivot-1);
                L(pivot,1:pivot-1)=L(maxi,1:pivot-1);
                L(maxi,1:pivot-1)=temp;
            end
        end %%Grade the column pivot using gauss elimination
        for i=pivot+1:rows
            num=U(i,pivot)/U(pivot,pivot);
            U(i,:)=U(i,:)-num*U(pivot,:);
            L(i,pivot)=num;
        end
        pivot=pivot+1;
    end
end

Надеюсь, что это поможет кому-то наткнуться на это в будущем.