Метод вращений решения линейных систем

МЕТОД ВРАЩЕНИЙ РЕШЕНИЯ  ЛИНЕЙНЫХ СИСТЕМ

( QR – разложение )

Рассмотрим  линейную систему общего вида . Предположим, что методом Гаусса решается система Ах =b с матрицей коэффициентов

 
 

Приведение  такой системы к треугольному виду прямым ходом метода Гаусса равносильно следующей последовательности эквивалентных преобразований матрицы А: 

Очевидно, в случае п х п -матрицы такого типа прямой ход метода Гаусса

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

    Рассмотренный пример оправдывает поиск других подходов к построению прямых методов решения линейных систем (2.1), возможно, более сложных, чем метод Гаусса, но не допускающих большого роста элементов в процессе преобразований и как следствие численно более устойчивых.

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

    Пусть и - некоторые отличные от нуля числа. Умножим первое уравнение системы

    

      на  , второе — на и сложим их; полученным уравнением заменим первое уравнение системы. Затем первое уравнение исходной системы умножаем на -    второе - на   и результатом их сложения заменяем второе уравнение. Таким образом, первые два уравнения системы заменяются уравнениями

      
 
 

    

    На  введенные два параметра  , наложим два условия:

      

    —условие  обнуления (т.е. исключение из второго уравнения) и

                            

• условие нормировки. Легко проверить, что за и удовлетворяющие этим условиям, можно принять

           (1)

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

    После фиксирования и способом (1) система принимает вид

         (2)

    где

     ,  ( ),

     , ( ),  

    Далее первое уравнение системы (2) заменяется новым, полученным сложением результатов  умножения первого и третьего уравнений (2) соответственно на 
 

   

    а третье - уравнением, полученным сложением  результатов умножения тех же уравнений соответственно на - и у. Получаем систему 
 
 

    

где 

, ( ),

, ( ),

Проделав  такие преобразования п -1 раз, придем к системе

   (3)

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

    Оно показывает сохранение величины суммы квадратов изменяемых на данном промежуточном шаге преобразований пары элементов любого столбца; так как остальные элементы столбцов при этом остаются неизменными, то, значит, на любом этапе преобразований длина столбца будет одной и той же, т.е. не будет наблюдаться роста элементов (но какой ценой это достигается! коэффициенты первого уравнения пересчитываются n-1 раз!).

    Дальше  точно так же за n-2 промежуточных шага преобразуем подсистему

    системы (3), создавая нули под элементом  и т.д.

    В результате п - 1 таких этапов прямого хода исходная система будет приведена к треугольному виду: 
 
 
 
 

                             

Нахождение отсюда неизвестных х ,х ,...,x не составляет труда. 

 

     Результаты работы программы

    Решение СЛАУ методом вращения (QR-разложение)

    =Исходная  система========================================

    0.00   1.00   2.00   3.00   4.00   5.00   6.00   7.00   8.00   9.00  | 1.000

    1.00   2.00   3.00   4.00   5.00   6.00   7.00   8.00   9.00  10.00  | 1.000

    2.00   3.00   4.00   5.00   6.00   7.00   8.00   9.00  10.00  11.00  | 1.000

    3.00   4.00   5.00   6.00   7.00   8.00   9.00  10.00  11.00  12.00  | 1.000

    4.00   5.00   6.00   7.00   8.00   9.00 10.00  11.00  12.00  13.00  | 1.000

    5.00   6.00   7.00   8.00   9.00  10.00  11.00  12.00  13.00  14.00  | 1.000

    6.00   7.00   8.00   9.00  10.00  11.00  12.00  13.00  14.00  15.00  | 1.000

    7.00   8.00   9.00  10.00  11.00  12.00  13.00  14.00  15.00  16.00  | 1.000

    8.00   9.00  10.00  11.00  12.00  13.00  14.00  15.00  16.00  17.00  | 1.000

    9.00  10.00  11.00  12.00  13.00  14.00  15.00  16.00  17.00  18.00  | 1.000 

    =Решение  системы=========================================

    x[0]=0.046  x[1]=-0.174  x[2]=-0.081  x[3]=-0.031  x[4]=-0.032  x[5]=0.043

    x[6]=0.052  x[7]=0.064  x[8]=0.041  x[9]=0.027

    Макс. невязка: 0.000000

 

    Листинг программы

    #include <stdio.h>

    #include <conio.h>

    #include <math.h>

 

    #define _N 10   // число уравнений

    #define _M 10   // число неизвестных

 

    class matrix {

      public:

       matrix(int n, int m);

       ~matrix();

       void print_slau();

       void print_x();

       void qr();

       void seta(int i, int j, float value);

       void setb(int i, float value);

      private:

       float **a;

       float *b;

       float *x;

       int n,m;

    };

    // ------------------------------------------------------

    // Конструктор

    // ------------------------------------------------------

    matrix::matrix(int N,int M) {

     n = N; m = M;

     a = new float* [n];

     for (int i=0; i<n; i++) {

          a[i] = new float [m];

          for (int j=0; j<m; j++) a[i][j]=0; }

     b = new float [n];

     x = new float [n];

     for (i=0; i<n; i++) { b[i]=0; x[i]=0; }

    }

 

    // ------------------------------------------------------

    // Деструктор

    // ------------------------------------------------------

    matrix::~matrix() {

     for (int i=0; i<n; i++) delete a[i];

     delete[] a;

     delete b; delete x;

    }

 

    // ------------------------------------------------------

    // Вывод СЛАУ на экран

    // ------------------------------------------------------

    void matrix::print_slau() {

     for (int i=0; i<n; i++) {

        for (int j=0; j<m; j++) printf("%.2f  ",a[i][j]);

        printf("| %.3f  ",b[i]);

        printf("\n");

     }

    }

 

    // ------------------------------------------------------

    // Печать решения и невязки решения

    // ------------------------------------------------------

    void matrix::print_x() {

     float max=0,h;

     for (int i=0; i<n; i++) {

          h=0;

          for (int j=0; j<n; j++) h=h+x[j]*a[i][j];

          if (max<fabs(b[i]-h)) max=fabs(b[i]-h);

          printf("x[%i]=%.3f  ",i,x[i]);

          if (i==5) printf("\n");

     }

     printf("\nМакс. невязка: %f\n",max);

    }

    // ------------------------------------------------------

    // Установить значение матрицы  a[][]

    // ------------------------------------------------------

    void matrix::seta(int i, int j, float value) {

      a[i][j] = value;

    }

    // ------------------------------------------------------

    // Установить значение вектора  b[]

    // ------------------------------------------------------

    void matrix::setb(int i, float value) {

      b[i] = value;

    }

 

    // ------------------------------------------------------

    // QR-алгоритм

    // ------------------------------------------------------

    void matrix::qr () {

      int l,k;

      float c[_N][_M];

      float s[_N][_M];

      float akk,akl,alk,all,bk,bl;

 
 
 
 

      // Прямой ход

      for (k=0; k<n-1; k++) {   // "Большой"  шаг (исключение переменных)

       for (l=k+1; l<n; l++) {     // "Малый"  шаг

        c[k][l] = a[k][k] / (sqrt( a[k][k]*a[k][k] + a[l][k]*a[l][k] ));