Решение набора систем уравнений с общими левыми частями

Этот алгоритм является развитием алгоритма решения системы линейных уравнений методом QR-разложения. Задача, стоящая перед нами, такова: есть набор систем N уравнений с N неизвестными, общей левой частью и разными правыми частями (всего M систем):

Можно объединить набор систем уравнений относительно неизвестных-векторов и получить одно матричное уравнение, в котором неизвестным, левой и правой частями являются матрицы: AX = B. Как и в предыдущем алгоритме, мы приводим левую часть к верхнетреугольной форме преобразованиями вращений: (QA)X=QB. Отсюда легко получается набор неизвестных X.

Поскольку матрица вращений одна для всех уравнений, то можно ускорить вычисления в M раз по сравнению с решением систем уравнений по одному.

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

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

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

unit QRSolveSystemUnit;

interface QRSolveSystem;

implementation

(*******************************************************************
Функция одновременного решения набора систем линейных уравнений
с общей правой частью методом QR-разложения.

Набор систем имеет вид:
CX = B, где C - матрица NxN, X - матрица NxM, B - матрица NxM.

Входные данные:
    A - в столбцах с 1 по N хранится матрица С, в следующих M столбцах
        хранится матрица B. Вся матрица имеет размер N*(N+M).
    N - число уравнений и неизвестных.
    M - число правых частей.
    X - матрица N*M, столбцы с 1 по M хранят решения.

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

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

Чтобы проводить вычисления независимо от того, насколько
матрица близка к вырожденной, следует задать epsilon
равным 0.
*******************************************************************)
function QRSolveSystem(
    A:          array of array of Real;
    const N:    Integer;
    const M:    Integer ;
    out X:      array of array of Real;
    Epsilon:    Real):Boolean;
var
    Col     :   Integer;
    Row     :   Integer;
    I       :   Integer;
    J       :   Integer;
    X1      :   Real;
    X2      :   Real;
    T       :   Real;
    T1      :   Real;
    T2      :   Real;
    C       :   Real;
    S       :   Real;
    Denomin :   Real;
begin
    SetBounds(X, [1,N], [1,M]);
    
    //вращение
    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
            begin
                T:=X2/X1;
                Denomin:=AbsReal(X1)*Sqrt(1+T*T)
            end
            else
            begin
                T:=X1/X2;
                Denomin:=AbsReal(X2)*Sqrt(1+T*T);
            end;
            C:=X1/Denomin;
            S:=-X2/Denomin;

            //поворот матрицы коэффициентов и правых частей
            for I:=Col to N+M 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
            for J:=1 to M do
            begin
                T1:=A[Row,N+J];
                for I:=Row+1 to N do
                    T1:=T1-A[Row,I]*X[I,J];
                X[Row,J]:=T1/A[Row,Row];
            end;
        end
        else
        begin
            Result:=False;
            Exit;
        end;
end;

end.

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