Решение СЛАУ методом Гаусса с частичным выбором главного элемента

Процедура решает неоднородную систему n линейных алгебраических уравнений с n неизвестными:

В отличие от предыдущего алгоритма на каждом шаге мы ищем не просто отличный от нуля коэффициент при x_k , а максимальный из них по абсолютной величине, в остальном практически та же схема приведения системы к треугольному виду:

Если при поиске отличного от нуля коэффициента такого не окажется, то матрица системы вырождена и алгоритм неприменим. Для сравнения с нолем в алгоритм передается малое число epsilon, и любое число, по модулю меньшее epsilon, считается нолем. В случае вырожденой матрицы функция возвращает False. Если матрица невырождена, то функция возвращает True, а переменная X содержит решение системы.

Если нашли ошибку в алгоритме - сообщите!

Реализации алгоритма на различных языках:

Блоксхемы:

Реализация алгоритма на AlgoPascal:

unit LESGaussMUnit;

interface LESGaussSolveM;

implementation

(*******************************************************************
Функция решения системы линейных уравнений вида:
A[1,1]*X[1] + .. + A[1,N]*X[N] = A[1,N+1]
        ....................
A[N,1]*X[1] + .. + A[N,N]*X[N] = A[N,N+1]
Гауссом с частичным выбором максимального элемента.

Функция   возвращает   True,   если    det(A') <>0
(Х хранит решение), и False если det(A')=0,     в
таком случае X не хранит решение.
A' - матрица из первых N столбцов матрицы A.

Epsilon определяет точность сравнения чисел. Если
число по модулю меньше или равно Epsilon, то оно
считается равным 0. Это число позволяет указать
алгоритму, насколько именно переданная в него матрица
должна быть близка к вырожденной, чтобы вызвать выход
со значением True.

Чтобы проводить вычисления независимо от того, насколько
матрица близка к вырожденной, следует задать epsilon
равным 0.
*******************************************************************)
function LESGaussSolveM(
            A   :   array of array of Real;
    const   N   :   Integer;
    out     X   :   array of Real;
    const   Epsilon : Real):Boolean;
    
var
    I   :   Integer;
    J   :   Integer;
    K   :   Integer;
    
    M   :   Real;
    T   :   Real;
begin
    SetBounds(X, [1,N]);
    Result:=True;
    for i:=1 to n do
    begin
        k:=i;
        m:=AbsReal(a[i,i]);
        for j:=i+1 to n do
        begin
            if m<AbsReal(a[j,i]) then
            begin
                m:=AbsReal(a[j,i]);
                k:=j;
            end;
        end;
        if AbsReal(m)>Epsilon then
        begin
            for j:=i to n+1 do
            begin
                t:=a[i,j];
                a[i,j]:=a[k,j];
                a[k,j]:=t;
            end;
            for k:=i+1 to n do
            begin
                t:=a[k,i]/a[i,i];
                a[k,i]:=0;
                for j:=i+1 to n+1 do
                begin
                    a[k,j]:=a[k,j]-t*a[i,j];
                end;
            end;
        end
        else
        begin
            Result:=False;
            Break;
        end;
    end;
    if Result then
    begin
        i:=n;
        repeat
            x[i]:=a[i,n+1];
            j:=i+1;
            while j<=n do
            begin
                x[i]:=x[i]-a[i,j]*x[j];
                j:=j+1;
            end;
            x[i]:=x[i]/a[i,i];
            i:=i-1
        until not(i>=1);
    end;
end;

end.

Бочканов Сергей, Быстрицкий Владимир
Copyright © 1999-2004
При поддержке проекта MANUAL.RU