Подпрограмма:  ASS2R

Назначение

Решение невырожденной разреженной линейной системы итерационным методом Гаусса - Зейделя (матрица системы представлена в формате RR (LU) U).

Математическое описание

Сокращенное название формата RR (LU) U происходит от английского словосочетания "Row - wise Representation, Lower - Upper, Unordered" (строчное представление, нижний треугольник - верхний треугольник, неупорядоченное).

Данный формат используется для прямоугольных (квадратных) матриц, у которых все или большинство диагональных элементов не равны нулю. Матрица  A в этом формате представляется в виде суммы L + D + U, где  L - нижняя треугольная,  D - диагональная,  U - верхняя треугольная матрицы. Диагональные элементы хранятся в отдельном одномерном массиве AD, а элементы матриц  L и  U содержатся в одномерном массиве AN и связаны между собой списками (портретами) IA и JA.

Например, пусть задана матрица

               |    4    0    1    0    2    |
               |    1    2    0    0    0    |
   A   =    |    0    0    2    1    0    |
               |    1    1    0    1    1    |
               |    0    0    0    0    16  | 

Тогда в формате RR (LU) U данная матрица с точностью до упорядоченности ее ненулевых внедиагональных элементов имеет вид:

                IA  =  ( 1; 3; 4; 5; 8; 8 )
                JA  =  ( 5, 3; 1; 4; 5, 1, 2 )
               AN  =  ( 2, 1, 1, 1, 1, 1, 1 )
               AD  =  ( 4, 2, 2, 1, 16 ) 

Пусть квадратная матрица  A порядка  n линейной системы Ax = b задана в формате RR (LU) U и пусть известно, что все диагональные элементы матрицы  A не равны нулю.

Метод Гаусса - Зейделя заключается в следующем. Перепишем  i - е уравнение системы в виде:

         i -1                                  n
          ∑   ai j x j  +  ai i x i  +   ∑     ai j x j  =  b i 
         j =1                               j= i +1 

Поскольку по предположению все  a_i i ≠ 0, то это уравнение можно записать в виде:

                             i -1                  n
          x  =  ( b i -   ∑   ai j x j -    ∑     ai j x j ) / ai i
                            j =1               j= i +1 

Это уравнение решается итерационно, если задать начальное приближение к решению x_j ⁽¹⁾ ( j = 1, 2, ..., n). В случае, когда начальное приближение x_j ⁽¹⁾ заранее неизвестно, то можно положить

        xj(1)  =  b j / aj j ,   j = 1, 2, ..., n . 

Следует отметить, что скорость сходимости итерационного процесса Гаусса - Зейделя может быть значительно увеличена, если начальное приближение к решению выбрано удачно.

Далее на  m - м итерационном шаге процесса Гаусса - Зейделя имеем (m = 2, 3, 4, ...)

                                i -1
     xi(m)  =  ( bi  -     ∑    ai j xj(m) -
                               j =1
                                                     n
                                                 -  ∑     ai j xj(m -1) ) / ai i
                                                   j= i +1 

Заметим, что здесь вновь вычисленные компоненты решения на  m - м шаге сразу же используются для вычисления следующих компонент.

В некоторых случаях сходимость метода Гаусса - Зейделя можно ускорить, если использовать прием верхней релаксации:

               xi(m)  =  xi(m -1)  +  Q( xi(m)  -  xi(m -1) ) 

Релаксационный множитель  Q обычно выбирается в пределах от 1 до 2. Если матрица положительно определенная, то  Q выбирается в пределах от 0 < Q < 1. Иногда в качестве  Q используется диагональная матрица, позволяющая подбирать релаксационные множители индивидуально для каждого уравнения системы.

Метод Гаусса - Зейделя хорошо сходится, если матрица  A является нижней или "почти нижней" треугольной. Для его сходимости не требуется диагонального преобладания в матрице  A. Необходимым и достаточным условием этого метода является следующее: все корни уравнения

               |   a1 1λ    a1 2      a1 3    ...  a1 n    |   
               |   a2 1λ    a2 2λ    a2 3    ...  a2 n    |
       det   |   .  .  .  .  .  .  .  .  .  .  .  .  .  .  .         =   0
               |   an 1λ    an 2λ    an 3λ  ... an nλ  | 

должны быть по модулю меньше 1. Если матрица  A является вещественной симметричной положительно определенной матрицей, то метод Гаусса - Зейделя всегда сходится. Скорость сходимости будет быстрее, если матрица  A близка к диагональной.

Описанный метод Гаусса - Зейделя реализован в виде подпрограммы ASS2R. В этой подпрограмме не требуется задавать на входе начальное приближение к решению: оно вычисляется в ASS2R по формуле x_i ⁽¹⁾ = b_i /a_i i. Подпрограмма предполагает задание максимально допустимого количества итераций. Контроль точности ведется по абсолютной погрешности EPS: текущая итерация считается приемлемым приближением к решению, если | x_i ^(m) - x_i ^(m - 1) | < EPS для всех  i.

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

Н.С.Бахвалов. Численные методы. Изд - во "Наука", 1973.

Использование

    SUBROUTINE  ASS2R ( IA, JA, AN, AD, N, B, X, Q, EPS, ITMAX,
                                             IFLAG) 

   Параметры 

   Версии 

   Вызываемые подпрограммы  нет 

   Замечания по использованию  

Пример использования

       DIMENSION  IA(6), JA(7), AN(7), AD(5), B(5), X(5) 
       COMMON /ASS2RR/ ITER  
       DATA  IA /1, 2, 3, 5, 6, 8/,
      *           JA /5, 1, 2, 1, 2, 3, 1/,
      *           AN /1., 1., 1., 1., 1., 1., 2./,
      *           AD /4., 2., 2., 8., 16./,
      *           B /1., 1., 1., 1., 1./
       N = 5 
       Q = 1.5 
       DO 1  I = 1, N
    1 X(I) = B(I) /AD(I)
       EPS = 0.001
       ITMAX = 500 
       CALL  ASS2R (IA, JA, AN, AD, N, B, X, Q, EPS, ITMAX, IFLAG)  

Результаты: 

       IFLAG = 0 ,   ITER = 7 ,

       X = (0.245396, 0.377041, 0.188364, 0.0778308, 0.0203379)