Решение методом вращений

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

здесь L - нижняя треугольная матрица (т.е. элемент l_ij =0, если j > i), а R_i,j  - матрица вращения. Матрица вращения это "почти единичная" матрица у которой r_ii =r_jj =cos(u), и r_ij =-r_ji =sin(u). Действительно нетрудно проверить, что умножая, например, матрицу 2-го порядка

на некоторый вектор x, мы осуществим поворот этого вектора на угол u относительно начала координат. Заметим, что определитель матрицы вращения равен 1, и R_ij *R ^T_ij  = E, т.е. обратная матрица совпадает с транспонированной. Итак представив исходную матрицу в виде произведения нижней треугольной матрицы на последовательность матриц вращения, мы можем решить вначале систему Ly = b, а затем систему R_n-1,n * ... *R_1,2 *x=y. Первая система легко решается благодаря простому виду матрицы L, а вторую систему нетрудно разрешить, благодаря свойствам матрицы вращения.

Несколько слов непосредственно о реализации алгоритма. Вначале последовательно для каждого ненулевого элемента расположенного выше главной диагонали матрицы A находится преобразование вращения, которое переводит его в нуль. Для хранения преобразование вращения R_ij  достаточно хранить соответствующий угол u, но проще хранить t = tg(2u), тогда

В результате последовательности преобразований, в матрице A будет, непосредственно матрица L (l_ij  = a_ij , i >= j; l_ij  = 0, i < j) и соответствующие преобразование R_ij  (t_ij  = a_ij , i < j, здесь t_ij  - тангес удвоенного угла вращения для преобразования R_ij ). Теперь решаем вначале систему Ly = b, а затем из вектора y, получаем решения исходной системы, обращая преобразования вращения.

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

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

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

Блоксхемы:

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

unit SolveLUUnit;

interface LUSolve;

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

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

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

Чтобы проводить вычисления независимо от того, насколько
матрица близка к вырожденной, следует задать epsilon
равным 0.
*******************************************************************)
function LUSolve(
            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;
    r1: Real;
    r2: Real;
    p:  Real;
    r:  Real;
    t:  Real;
    c:  Real;
    s:  Real;
begin
    SetBounds(X, [1,N]);
    
    //приводим к нужной форме
    i:=1;
    while i<=n do
    begin
        x[i]:=a[i,n+1];
        i:=i+1
    end;
    k:=1;
    while k<=n do
    begin
        j:=k+1;
        while j<=n do
        begin
            r2:=a[k,j];
            if r2<>0 then
            begin
                r1:=a[k,k];
                if r1<>0 then
                begin
                    if AbsReal(r1)<=AbsReal(r2) then
                    begin
                        p:=r1/r2;
                        r:=sqrt(1+p*p);
                        t:=1/(r+AbsReal(p));
                        if p<0 then
                        begin
                            t:=-t;
                        end;
                        if r1<0 then
                        begin
                            a[k,k]:=-r*AbsReal(r2)
                        end
                        else
                        begin
                            a[k,k]:=r*AbsReal(r2)
                        end;
                    end
                    else
                    begin
                        p:=r2/r1;
                        r:=sqrt(1+p*p);
                        t:=p/(r+1);
                        a[k,k]:=r*r1
                    end;
                end
                else
                begin
                    t:=1;
                    a[k,k]:=a[k,j]
                end;
                a[k,j]:=t;
                r:=t*t;
                r1:=1/(1+r);
                c:=(1-r)*r1;
                s:=2*t*r1;
                i:=k+1;
                while i<=n do
                begin
                    r1:=a[i,k];
                    r2:=a[i,j];
                    a[i,k]:=r1*c+r2*s;
                    a[i,j]:=r2*c-r1*s;
                    i:=i+1
                end;
            end;
            j:=j+1
        end;
        k:=k+1
    end;
    
    //решение промежуточной системы Ly=b - ошибки возможны именно здесь
    Result:=True;
    if AbsReal(a[1,1])<Epsilon then
    begin
        Result:=False;
        Exit;
    end;
    x[1]:=x[1]/a[1,1];
    k:=2;
    while k<=n do
    begin
        s:=0;
        j:=1;
        while j<=k-1 do
        begin
            s:=s+a[k,j]*x[j];
            j:=j+1
        end;
        if AbsReal(a[k,k])<Epsilon then
        begin
            Result:=False;
            Exit;
        end;
        x[k]:=(x[k]-s)/a[k,k];
        k:=k+1
    end;


    //решение Qx=y
    for k:=n-1 downto 1 do
        for j:=n downto k+1 do
        begin
            t:=a[k,j];
            c:=(1-t*t)/(1+t*t);
            s:=-2*t/(1+t*t);
            r1:=c*x[k]+s*x[j];
            r2:=-s*x[k]+c*x[j];
            x[k]:=r1;
            x[j]:=r2;
        end;
end;

end.

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