Решение методом вращений (другой вариант, QR-разложение)

Решаем систему n линейных уравнений Ax = b с квадратной матрицей. Вначале приведем матрицу A к верхне-треугольному виду преобразованиями вращения:

здесь R - верхняя треугольная матрица (т.е. элемент r_ij =0, если j < i), а Q_i,j  - матрица вращения. Матрица вращения это "почти единичная" матрица у которой q_ii =q_jj =cos(u), и q_ij =-r_ji =sin(u). Действительно нетрудно проверить, что умножая, например, матрицу 2-го порядка

на некоторый вектор x, мы осуществим поворот этого вектора на угол u относительно начала координат.

Существует теорема, согласно которой любую матрицу A можно привести к верхне-треугольному виду умножением на матрицы вращений. Несколько слов непосредственно о реализации алгоритма. Во время работы алгоритма мы обходим матрицу по строчкам сверху вниз, начиная с первой. Каждая строка обходится слева направо. К элементам, расположенным ниже диагонали, мы применяем преобразование вращения, которое обнуляет их. При этом мы будем умножать как матрицу коэффициентов, так и правую часть уравнения. Полученная система с треугольной матрицей решается тривиально.

Достоинством этого алгоритма в сравнении с методом Гаусса является бОльшая точность в случае решения системы с плохо обусловленной матрицей, хотя алгоритм несколько проигрывает по скорости работы.

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

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

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

unit QRSolveUnit;

interface QRSolve;

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]
методом QR-разложения.

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

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

Чтобы проводить вычисления независимо от того, насколько
матрица близка к вырожденной, следует задать epsilon
равным 0.
*******************************************************************)
function QRSolve(
    A: array of array of Real;
    const N:Integer;
    out X:array of Real;
    Epsilon:Real):Boolean;
var
    Col     :   Integer;
    Row     :   Integer;
    I       :   Integer;
    X1      :   Real;
    X2      :   Real;
    T1      :   Real;
    T2      :   Real;
    C       :   Real;
    S       :   Real;
    Denomin :   Real;
begin
    SetBounds(X, [1,N]);
    
    //вращение
    for Col:=1 to N do
        for Row:=Col+1 to N do
        begin
            //выбор угла вращения
            X1:=A[Col,Col];
            X2:=A[Row,Col];
            if (X1=0) and (X2=0) then
                Continue;
            if AbsReal(X1)>AbsReal(X2) then
                Denomin:=AbsReal(X1)*Sqrt(1+Sqr(X2/X1))
            else
                Denomin:=AbsReal(X2)*Sqrt(1+Sqr(X1/X2));
            C:=X1/Denomin;
            S:=-X2/Denomin;

            //поворот матрицы коэффициентов и правой части
            for I:=Col to N+1 do
            begin
                T1:=A[Col,I];
                T2:=A[Row,I];
                A[Col,I]:=C*T1-S*T2;
                A[Row,I]:=S*T1+C*T2;
            end;
        end;

    //обратный ход
    Result:=True;
    for Row:=N downto 1 do
        if AbsReal(A[Row,Row])>=Epsilon then
        begin
            X[Row]:=A[Row,N+1];
            for I:=Row+1 to N do
                X[Row]:=X[Row]-A[Row,I]*X[I];
            X[Row]:=X[Row]/A[Row,Row];
        end
        else
        begin
            Result:=False;
            Exit;
        end;
end;

end.

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