БЧА НИВЦ МГУ. ASS2R_P. Решение систем с матрицами специального вида (невырожденными)

Текст подпрограммы и версий
ass2r_p.zip , ass2e_p.zip

Тексты тестовых примеров
tass2r_p.zip , tass2e_p.zip

Подпрограмма: ASS2R (модуль ASS2R_p)

Назначение

Решение невырожденной разреженной линейной системы итерационным методом Гаусса - Зейделя (матрица системы представлена в формате 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
          ∑   a_{i j} x_j  +  a_{i i} x_i  +   ∑     a_{i j} x_j  =  b_i 
         ^{j =1}                               ^{j= i +1}

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

                             _{i -1}                  _n
          x  =  ( b_i -   ∑   a_{i j} x_j -    ∑     a_{i j} x_j ) / a_{i i}
                            ^{j =1}               ^{j= i +1}

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

        x_j⁽¹⁾  =  b_j / a_{j j} ,   j = 1, 2, ..., n .

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

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

                                _{i -1}
     x_i^(m)  =  ( b_i  -     ∑    a_{i j} x_j^(m) -
                               ^{j =1}
                                                     _n
                                                 -  ∑     a_{i j} x_j^{(m -1)} ) / a_{i i}
                                                   ^{j= i +1}

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

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

               x_i^(m)  =  x_i^{(m -1)}  +  Q( x_i^(m)  -  x_i^{(m -1)} )

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

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

               |   a_{1 1}λ    a_{1 2}      a_{1 3}    ...  a_{1 n}    |   
               |   a_{2 1}λ    a_{2 2}λ    a_{2 3}    ...  a_{2 n}    |
       det   |   .  .  .  .  .  .  .  .  .  .  .  .  .  .  .         =   0
               |   a_{n 1}λ    a_{n 2}λ    a_{n 3}λ  ... a_{n n}λ  |

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

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

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

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

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

procedure ASS2R(var IA :Array of Integer; var JA :Array of Integer;
                var AN :Array of Real; var AD :Array of Real;
                N :Integer; var B :Array of Real;
                var X :Array of Real; Q :Real; EPS :Real;
                ITMAX :Integer; var IFLAG :Integer);

   Параметры

IA, JA, - AN	заданные портрет и ненулевые внедиагональные элементы матрицы A в формате RR (LU) U;
AD -	вещественный одномерный массив длины N, содержащий диагональные элементы матрицы A; требуется, чтобы ни один диагональный элемент не был равен нулю;
N -	заданный порядок матрицы A (тип: целый);
B -	вещественный одномерный массив длины N, содержащий вектор правой части системы;
X -	вещественный одномерный массив длины N, на выходе из подпрограммы содержащий вычисленные компоненты вектора решения системы;
Q -	заданный релаксационный множитель (тип: вещественный);
EPS -	заданная абсолютная погрешность, с которой требуется вычислить решение;
ITMAX -	заданное максимальное количество итераций по методу Гаусса - Зейделя;
IFLAG -	целая переменная, служащая для сообщения о режиме окончания работы подпрограммы:

IFLAG=0 -	решение системы вычислено с заданной точностью;
IFLAG=1 -	решение не получено с заданной точностью за заданное максимальное количество итераций.

   Версии

ASS2E -

решение невырожденной разреженной линейной системы итерационным методом Гаусса - Зейделя в режиме расширенной (Extended) точности; при этом параметры AN, AD, B, X, Q и EPS должны иметь тип Extended.

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

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

Подпрограммы ASS2R и ASS2E имеют глобальную запись (структуру данных) с именем _ASS2RR и элементом ITER. Значение целой переменной ITER полагается равным количеству в действительности выполненных итераций, если IFLAG = 0. Если IFLAG = 1, то ITER = ITMAX .

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

Unit TASS2R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, ASS2R_p;

function TASS2R: String;

implementation

function TASS2R: String;
var
N,I,ITMAX,_i,IFLAG :Integer;
Q,EPS :Real;
X :Array [0..4] of Real;
const
IA :Array [0..5] of Integer = ( 1,2,3,5,6,8 );
JA :Array [0..6] of Integer = ( 5,1,2,1,2,3,1 );
AN :Array [0..6] of Real = ( 1.0,1.0,1.0,1.0,1.0,1.0,2.0 );
AD :Array [0..4] of Real = ( 4.0,2.0,2.0,8.0,16.0 );
B :Array [0..4] of Real = ( 1.0,1.0,1.0,1.0,1.0 );
label
_1;
begin
Result := '';  { результат функции }

{      ТЕСТ ДЛЯ ПОДПРОГРАММЫ ASS2R }

N := 5;
Q := 1.5;
for I:=1 to N do
 begin
_1:
  X[I-1] := B[I-1]/AD[I-1];
 end;
EPS := 0.001;
IТМАХ := 500;
ASS2R(IA,JA,AN,AD,N,B,X,Q,EPS,ITMAX,IFLAG);
Result := Result + Format('%s',[' IFLAG=']);
Result := Result + Format('%5d ',[IFLAG]);
Result := Result + Format('%s',[' _ASS2RR.elm1=']);
Result := Result + Format('%5d ',[_ASS2RR.elm1]) + #$0D#$0A;
Result := Result + Format('%s',[' X=']);
Result := Result + #$0D#$0A;
for _i:=0 to 4 do
 begin
  Result := Result + Format('%20.16f ',[X[_i]]);
  if ( ((_i+1) mod 4)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
UtRes('TASS2R',Result);  { вывод результатов в файл TASS2R.res }
exit;
end;

end.


Результаты: 

       IFLAG = 0 ,   ITER = 7 ,

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