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

Назначение

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

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

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

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

    int ass2r_c (integer *ia, integer *ja, real *an, real *ad,
            integer *n, real *b, real *x, real *q, real *eps, integer *itmax, 
            integer *iflag)

   Параметры 

   Версии 

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

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

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

struct {
    int iter;
} ass2rr_;

#define ass2rr_1 ass2rr_

int main(void)
{
    /* Initialized data */
    static int ia[6] = { 1,2,3,5,6,8 };
    static int ja[7] = { 5,1,2,1,2,3,1 };
    static float an[7] = { 1.f,1.f,1.f,1.f,1.f,1.f,2.f };
    static float ad[5] = { 4.f,2.f,2.f,8.f,16.f };
    static float b[5] = { 1.f,1.f,1.f,1.f,1.f };

    /* System generated locals */
    int i__1;

    /* Local variables */
    extern int ass2r_c(int *, int *, float *, float *, int *, float *,
                       float *, float *, float *, int *, int *);
    static int i__, n;
    static float q;
    static int iflag;
    static float x[5];
    static int itmax;
    static float eps;

    n = 5;
    q = 1.5f;
    i__1 = n;
    for (i__ = 1; i__ <= i__1; ++i__) {
/* l1: */
        x[i__ - 1] = b[i__ - 1] / ad[i__ - 1];
    }
    eps = .001f;
    itmax = 500;
    ass2r_c(ia, ja, an, ad, &n, b, x, &q, &eps, &itmax, &iflag);

    printf("\n  %5i \n", iflag);
    printf("\n  %5i \n", ass2rr_1.iter);
    printf("\n  %15.6e %15.6e %15.6e %15.6e %15.6e \n",
           x[0], x[1], x[2], x[3], x[4]);
    return 0;
} /* main */


Результаты: 

       iflag = 0 ,   iter = 7 ,

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