CC BY
9УДК 517.6:517.958
ВЕКТОРНЫЙ аналог метода прогонки для решения
ТРЕХ- И ПЯТИДИАГОНАЛЬНЫХ МАТРИЧНЫХ УРАВНЕНИЙ
Н.К. ВОЛОСОВА
(Московский государственный технический университет им. Н.Э. Баумана); д-р физ.-мат. наук, проф. К.А. ВОЛОСОВ (Российский университет транспорта, Москва); канд. физ.-мат. наук А.К. ВОЛОСОВА (ООО «Трамплин», Москва);
канд. физ.-мат. наук, доц. Д.Ф. ПАСТУХОВ, канд. физ.-мат. наук, доц. Ю.Ф. ПАСТУХОВ (Полоцкий государственный университет)
Предложен алгоритм векторного аналога прогонки для решения произвольных матричных уравнений с квадратными трех- и пятидиагональными матрицами за конечное число арифметических вычислений. Доказаны достаточные условия корректности векторных формул прогонки для произвольных трех-диагональных матриц (теорема 1) и достаточные условия для пятидиагональных симметрических матриц Теплица (теорема 2). Приведенные программа и два примера показывают, что данные алгоритмы являются точными. Предложен численный алгоритм поиска предельных значений для коэффициентов прогонки вперед (теорема 3), показано, что полученные численные предельные значения не противоречат теореме 2.
Ключевые слова: векторный аналог метода прогонки, трех- и пятидиагональные матрицы, матрица Теплица, выпуклые множества, численные методы математической физики, параллельные вычисления.
Введение. Матрицы и матричные уравнения специального типа применяются во многих разделах прикладной математики. В квантовой механике динамика частиц со спином определяется матрицами кватернионов (полукватернионов) [1, 2]. Другой пример: одним из методов решения эллиптических уравнений математической физики численными методами является метод прогонки [3, 4, 5, 12, 13, 21, 22]. Здесь используются матрицы диагонального вида. Алгебраический метод прогонки, используемый построчно на прямоугольной сетке совместно с формулой простой итерации [5] является приближенным методом, так как число итераций не ограничено, но имея формулы с аппроксимацией дифференциальных операторов с высоким порядком погрешности можно значительно снизить число и время вычислений [5]. В данной работе рассмотрен векторный аналог метода прогонки для решения матричных уравнений с квадратными матрицами трехи пятидиагонального типа за конечное число арифметических действий. Если диагональная матрица, соответствующая разностным уравнениям прогонки, имеет постоянные коэффициенты на главной диагонали и на двух (четырех) диагоналях параллельным главной, то матрица коэффициентов называется матрицей Теплица. В данной работе доказаны необходимые условия корректности формул прогонки для произвольных трехдиагональных матриц и для пятидиагональных симметрических матриц Теплица, решаемых векторным аналогом метода прогонки. Сегодня необходимо рассматривать также численные задачи с параллельными вычислениями [3, 4, 7, 11, 14]. Поэтому для решения пятидиагональных матричных уравнений в работе рассмотрены два алгоритма последовательного и параллельного вычисления.
Уравнения математической физики питают своими идеями не только традиционные разделы численной математики, такие как матричные вычисления, которым посвящена данная работа, но и новые ветви прикладной математики, такие как стеганография (впервые эту идею применила Н.К. Волосова [17-20]).
Постановка задачи. Рассмотрим матричное уравнение, в котором неизвестная матрица X, а также заданные матрицы А левой части и Е правой части уравнения (1) являются квадратными порядка п
Кроме того, в матричном уравнении (1) рассмотрим матрицу А трехдиагонального или пятидиагонального типа соответственно, у которой коэффициенты удовлетворяют условиям (2)
АХ = ^.
(1)
С учетом условия (2), уравнение (1) запишем подробно для трехдиагональных матриц
а11 а12 0 0
а21 а22 а23 0.... 0
0 а32 азз аз4 0... 0
0 0... ап_х 0 0...
п-2 "п-1,п-1 "п-1,п
0
/11 /и/21 /22"
/1п-1 /1.п Л,п-1 /2,п
/ / ... / /
л п1 л п2 '" л п,п-1 л п,
(3)
Транспонируя уравнение (1) и его подробную запись (3), получим соответственно формулы (4), (5):
(4)
АХ = ^ < Хт ■ Ат = Ет,
а11 а21 0 0... 0 "
а12 а22 а32 0.... 0 ' /11 /21... /п-1,1 /п,1
0 а 23 «33 а43 0... 0 = /22... /п-1,2, ,2
0 0... ап-2,п-1 ап-1,п-1 ап,п-1 _ /1,п /2,п ... / п-1,п Г ■У п,п
0 0... 0 а п-1,п а п,п
(5)
Последнее матричное уравнение (5) с учетом условий (2) равносильно системе векторных уравнений (6):
*11а11 + "^21 а12 = /11 , -^-12 а11 + Х22а12 = /12'...Х1 ¡а11 + Х2 ¡а12 = /1 ; ' = 1 п << ацХ + а12Х = /
к+1, а к+1 = /к, ;, = 1, пЛ = 2, п -1<
Г*-1,1а*,*-1 + Хк,1акк + Хк+1,1ак,к+1 = /к,1, Хк-1,;ак,к-1 + Хк,;акк + Хк^ >а<'
(6)
к-1 к < ак,к-1Х + аккХ + ак,к+1Х
к+1 = /к, Ук = 2, п - 1
Хп-1,1ап,п-1 + Хп,1ап,п = /п,1,...Хп-1,;ап,п-1 + Хп;ап,п = /„,; >= 1 п < ап,п-1Х"-1 + ап,пХ" = Г
В векторной системе уравнений (6) в к-е уравнение входят строки с номерами к - 1, к, к + 1 решения Х-матрицы с коэффициентами из к-й строки матрицы А и из к-й строкой матрицы Е или в векторном виде
1 2 _ г1
а ^ х I а 2 Х — /
ак,к-Х-1 + ашхк + акк+1 хк+1 = /к,Ук = 2,п -1
(7)
а ,хп-1 + а хп = /п
п,п-1 п,п ^
Будем искать решение рекуррентно заданной системы векторных уравнений (7) в виде
х = Хкх+ +ук, к = 1,п-1 .
(8)
Из первого уравнения системы (7) имеем
х1= х2<"^=_/.
/
Поскольку из (8) хк 1 =\_ 1хк , то преобразуем среднее уравнение системы (7):
(1 к \ к к+1 гк к ( 1 \ к+1 гк
к-1Х +vk-l)+ аккХ + акл+1х = / , Х (акк-А-1 + акк)=-ак,к+1Х + / - ак,к-1ук-1 <
< X =
гк
/ - ак,к-1ук-1 (акк-Лк-1 + акк ) (акк-Лк-1 + акк )
-х +
<1к =
гк
-ак,к+1_ _ / - ак,к-1Ук-1
(акк-Лк-1 + акк ) (акк-Лк-1 + акк )
,к = 2, п -1.
(9)
11
11
Анализ размерности [6, 9] показывает, что в формулах (8), (9) величины Xk, a,
'к' /t,к-1' кк' "■к,к+1
являются
числами, а /к, \к - векторами. Кроме того, последнее уравнение системы (7) имеет на одно слагаемое меньше, чем среднее уравнение, поэтому и решение для последнего уравнения (7) следует искать не в виде (8), но в виде хп = \п. Используя (8) при к = п - 1 и подставляя уравнение хп—1 = Хп_ 1хп + в последнее уравнение системы (7), получим
п—1 п / Л п \
а ,х + а х = а ,1л ,х , ) + а
п,п—1 п,п п,п— 1 у п—1 п—1 у
= fn « х" (a„,„_iX„_i + an n ) = fn - a
^ х =
f" - an,n-1Vn-1
и-Ли-i + an
■ = v.
(10)
Уравнения (9) называются формулами прогонки вперед, а уравнения (10), (8) - формулами прогонки назад. Рассмотрим тестовый пример (11), в котором вычисления проверяются напрямую перемножением матриц:
4 -1 0 0 0 0 0 1 2 1 2 1 2 1 2 7 2 7 2 7 2
-1 4 -1 0 0 0 0 2 1 2 1 2 1 2 6 0 6 0 6 0 6
0 -1 4 -1 0 0 0 1 2 1 2 1 2 1 0 6 0 6 0 6 0
0 0 -1 4 -1 0 0 2 1 2 1 2 1 2 = 6 0 6 0 6 0 6
0 0 0 -1 4 -1 0 1 2 1 2 1 2 1 0 6 0 6 0 6 0
0 0 0 0 -1 4 -1 2 1 2 1 2 1 2 6 0 6 0 6 0 6
0 0 0 0 0 -1 4 1 2 1 2 1 2 1 2 7 2 7 2 7 2
(11)
Трехдиагональная матрица из тестового примера (11) применяется для решения уравнения Пуассона на прямоугольнике с шаблоном «крест» [3, с. 584]. Программа, написанная нами на языке FORTRAN [7, 8] с использованием алгоритма (8)-(11), возвращает решение (неизвестную матрицу X в примере (11) по заданным матрицам A, F (таблица 1).
Таблица 1. - Решение, полученное программой с использованием алгоритма (8)—(10)
(a
i/j 1 2 3 4 5 6 7
1 1.000000000 2.000000000 1.000000000 2.000000000 1.000000000 2.000000000 1.000000000
00000 00000 00000 00000 00000 00000 00000
2 2.000000000 1.000000000 2.000000000 1.000000000 2.000000000 1.000000000 2.000000000
00000 00000 00000 00000 00000 00000 00000
3 1.000000000 2.000000000 1.000000000 2.000000000 1.000000000 2.000000000 1.000000000
00000 00000 00000 00000 00000 00000 00000
4 2.000000000 1.000000000 2.000000000 1.000000000 2.000000000 1.000000000 2.000000000
00000 00000 00000 00000 00000 00000 00000
5 1.000000000 2.000000000 1.000000000 2.000000000 1.000000000 2.000000000 1.000000000
00000 00000 00000 00000 00000 00000 00000
6 2.000000000 1.000000000 2.000000000 1.000000000 2.000000000 1.000000000 2.000000000
00000 00000 00000 00000 00000 00000 00000
7 1.000000000 2.000000000 1.000000000 2.000000000 1.000000000 2.000000000 1.000000000
00000 00000 00000 00000 00000 00000 00000
Замечание 1. Сравнение значений таблицы 1 и второй матрицы X из примера (11) показывает их полное совпадение с двойной точностью. Таким образом, алгоритм (8)-(10) является точным методом решения трехдиагональных уравнений (1) за конечное число арифметических операций [3, 4]. Оценим число арифметических операций. Для вычисления хп по формуле (10) необходимо 3п + 2 операций, для вычисления Л и ук по формуле (9) необходимо 3(п - 2) и (3п + 2)(п - 2) операций соответственно. Для вычисления хк по формуле (8) число операций составит 2п(п - 1). Общее число арифметических операций N = 3п + 2 + 3п — 6 + 3п2 — 4п — 4 + 2п2 — 2п = 5п2 — 8~5я2.
Теорема 1 (достаточные условия корректности алгоритма (8)-(10). Пусть выполнены условия:
1) |а (| > |а|Н | + |а ,+11 > \аи+1| > 0, V? = 2, п — 1 трехдиагональная матрица А уравнения (1) с нестрогим
диагональным преобладанием элементов, а для первой и последней строк Ц11 > \а12 \ > 0, \апп| > |ая | > 0 ;
2) И Ц/II <«>, V? = й=1.
Тогда:
1) |^|< 1, Vk = 1,n-1;
2) формулы прогонки (9), (10) корректны, т.е.
Vj. < да, k = 1, n — 1, х < да, k = 1, n.
Доказательство проведем по индукции. 1) Для базы индукции при к = 1 имеем
I I I I I I "12
KI - Ы ^М =W <1N =
( a2 + &22)
|— |Я2Д| N
■=1.
База индукции проверена. Преобразуем неравенство:
Ки+1 - K+U + К+и+2 > ai+U+1\ — Км - К+V+2
Индуктивный переход. Пусть верно
,(9)
\кк\< 1, Vk = 1,i ^+J =
ai+u ^i + ai+1i+1
i+1i+1 "i+1.
i m
^ I i+1,i+ 2 | ^
i+1i+1 i+1,i
Первая часть теоремы 1 доказана Vk = 1, n — 1.
2) Обозначим IIf'\ = max ||f |I, IIf || = max ||f'\ , ||v' || = max ||v,. 11J|v|| = max ||v1| . Так как
II II j = 1,n II J II 11 11 ;=!,„ II II II II j = 1,n II J II 11 11 i=1,„—1 И И
| < 1, Vk = 1, n — 1, |akд—Ak—1 + ak I - Ы — |ak,k—1 | К—11 - К I — |ak,k—J - |ak,k+1 | > 0(ak,k+1 * 0)
1 1
\ak ,k—1K k—1 + a
■ < да, Vk = 1, n — 1.
kk "k ,k+1
База индукции a1,11 - |a121 > 0,||v J = -—j- < -—- < да проверена. Индуктивный переход. Пусть
vk_J <да, vj = J
\\/k — ak,k—1Vk—J \\/k — ak,k—1VkJ \\/k\\ + 1 ak,k—J| |Vk—11|
+ a.
k,k—1 vk—1
\ak ,k—A k—1 + akk
< да, Vk = 2, n — 1:
f — a„,„—1V„—J /1 + an,n—1\\\Vn—1¡ f + К „—1 V 1)
"-" < -—p—n—¡——¡— < п—г—¡——г <да.
a ,Х , +a
n,n—1 n—1 n
По формуле (8)
|^|< 1, к = 1, п -1;
Цх1 || <|х2|| +|| < ||х2 || +<|х3 || + 2|VI <|х4|| + 3|<|хп\ + (и-1)||у|| < да, ||хк\ <да,Ук = . Теорема 1 доказана.
Рассмотрим матричное уравнение (1) с пятидиагональной матрицей, то есть со вторым условием на коэффициенты (2). Повторяя рассуждения, аналогичные (3)-(7), получим систему векторных уравнений (12):
1 + 2 3 4 у 2
ak,k —2xk —2 + ak,k—1 xk—1 + akkXk + ak¿+,хк+1 + ak¿+2xk+2 = fk,Vk = 3,n — 2 .
a^,^ х"—3 + a„—1^2 x"—2 + aи-l,и-l*и_1 + a^x" = f"—1
a 1x"—2 + a .х^1 + a xn = fn
a
a
a
—a
2,3
2,3
2,3
<
<
<
a
a
1.1
2,3
ai+1,i+1 «7+1,1 \ai+1,i+2
—a
a
a
<
<
a
a
a
a
k ,k+1
k kk+1
k,k+1
xn\\ =
a — a ,
n,n n,n—1
a — a ,
n,n n,n—1
В системе уравнений (12), кроме известных элементов аг . пятидиагональной матрицы, заданы вектор-строки /',' = 1, п правой части уравнения (1), х',' = 1,п - неизвестные вектор-строки уравнения (1). Аналогично (7), (8) будем искать решение третьей строки системы (12) в виде
хк = Л хк+1 + Л2 кхк +2 + ,к = 1, п — 2 .
(13)
Теория размерностей [6, 9] показывает, что в (13) Х1к, Х2к являются числами, а Ук, как и хк, -векторами.
Выразим из первого уравнения (12) х1 = ——х2 — — х3 + —. Сравнивая последнее выражение х1 с уравнением (13) при к = 1, получим
Л1,1 =— XЛ2, =— XУ1,, = ^,] = 1,п .
а,,
(14)
Подставив х1 во второе уравнение (12), из которого выразим х2, получим
1 2 3 4 _
1 х I г, х I ^ х I л х — а-^ •
12 2 13 3
—— х — —х +
/
1
2 3 4 _
-*11 у
х2 +
Л1 у
х3 + а24 х4 = /2 —
а21/1
11 у
х2 =
^ а13а21 а23а11 1 „3
х —
А ( х4 +
а11/ а21/
1
Сравнивая последнее выражение и формулу (13) при к = 2, получим коэффициенты
Л1,2 = -
Л2,2 = —"
а11/] а21/1
] = 1, п .
(15)
В работе [5, с. 69] получены формулы для решения скалярного разностного уравнения
Лкхк—2 + 4Л—1 —Скхк + В1кхк+1 + В2кхк +2 = / , Vk = 2, п — 2
(16)
с коэффициентами прогонки
Л1, к =
В1к + А2к Л2к— 1 + А1к Л1к—2Л2к—1
С — А Л Л — А Л —А Л Ск А1к Л1к— 2Л1к—1 А1к Л2к—2 А2к Л1к—1
Л 2,к =
В,.
С — А Л Л —А Л — А Л Ск А1к Л1к— 2Л1к—1 А1к Л2к—2 А2к Л1к—1
А1к Л1к— 2Ук—1 + А1к У к— 2 + А2к Ук— 1 — Рк
С — А Л Л — А Л — А Л Ск А1кЛ1к—2Л1к— 1 А1к Л2к— 2 А2к Л1к—1
к = 3, п — 2 .
(17)
Сравнивая уравнения (16) с третьим уравнением системы (12), получим формулы для векторных формул метода прогонки в соответствии с (17):
Л1,к = —
_ак ,к+1 + ак,к—1Л2к— 1 + ак ,к—2Л1к— 2Л2к— 1
\ак ,к + ак ,к—2Л1к— 2 Л1к—1 + ак ,к— 2Л2к— 2 + ак ,к—1Л1к—1 )
Л 2,к =
ак ,к + ак ,к— 2Л1к— 2Л1к—1 + ак ,к—2Л2к—2 + ак ,к—1Л1к— 1
Ук, ] = —
Ч ,к—2У 2Ук—1, ] + ак ,к—2Ук—2, ] + ак ,к—1Ук—1, ]—*к, , ^
V ак,к + ак,к—2Л1к—2Л1к— 1 + ак,к— 2Л2к—2 + ак,к— 1Л1к— 1 У
к = 3, п — 2, ] = 1, п .
(18)
Формулы (14), (15) совместно с (18) называются формулами прогонки вперед. В настоящее время в численных методах актуально рассматривать задачи с параллельными вычислениями, когда несколько ядер процессора выполняют однотипные операции для сокращения времени работы программы. Например, задачу (1) могут параллельно решать два ядра, если известно решение одной
а
а
11
1
а 2 а 1
а22 —
а23 —
а
11
V
V.
а24 а11
а13а21 а23а11
У 2, ] =
к
строки матрицы X. Для простоты будем считать известной последнюю строку решения хп и укажем формулы получения остальных строк. Выражая из последнего уравнения (12) хп—2, получим
а а п
п- 2 _ _ и,и—1 п-1 _ п,п п , / х — х х +
а_ „о а „о _
(19)
Используя уравнение (13) можем записать хп 3 = Л1л_3хп 2 + Л2л_3хп 1 + ул_3, к = п — 3 , которое подставим в четвертое уравнение системы (12):
п—1,п—3 п—1,п— 2 п—1,п—1 п—1,п
= р 1 « а , ,(Л, х" 2 +Л х" 1 +у ,) + а , ,хп 2 + а , ,хп 1 + а , хп =
л п—1,п—3 у 1,п—3 2,п—3 п—3 у п—1,п—2 п—1,п—1 п—1,п
= ( а , + а , хп—2 + (а , + а , Л хп—1 + а , хп = /п—1 — а , ,у. ..
у п—1 ,п—3 1,п—3 п—1,п—2 у у п—1,п—3 2,п—3 п—1,п—1 у п—1,п л п—1,п—3 к—3
В последнее уравнение подставим хп 2 - правую часть формулы (19):
(ап— 1,п—3Л1,п—3 + ап—1,п— 2 )
^ а , , а ^
п,п— 1 п—1 _ п,п У
V ап,п— 2 ап,п— 2 ап,п— 2 У
+ (ап—1,п—3Л2,п—3 + ап—1,п—1 ) 1 + ¿V 1,пхп = Г~1 - ап- 1 , п-3 Ук-3 »
» х
ап—1,п— 3Л2,п—3 + ап— 1,п—1 '
1 (а , Д, ,+ а , ,,) = х" —а , +--^^ (а , А, ,+ а , ,,)
\ п—1,п—3 1,п—3 п—1,п—2 / п—1,п \ п—1,п—3 1,п—3 п—1,п—2 /
+/п—1 — ап—1,п—3Ук—3--(ап—1,п—3Л1,п—3 + ап—1,п—2 ) '
- 4- ( 5, + )
п \ ' ' /
п
+ У-"1 — ап—1,п—3Ук—3,] — (ап—1,п—3Л1, п—3 + ап—1,п—2 )
« х,. = -
У
а
а
л ап,п—1 / л \
an-1,n-3Л2,n-3 + an-1,n-1 (an-1,n-3Л1,n-3 + an-1,n-2 )
, ] = 1, п. (20)
V
а
У
Сравнивая формулу (20) с решением второй строки системы (23) видим, что
Г а 1
ап п / \
—ап—1,п +-'—( ап—1,п—3Л1,п—3 + ап—1,п—2 )
V ап, п—2 У
Л1,п—1 = Т-^
ап п—1 / \
ап—1,п—3Л2,п—3 + ап—1,п—1--'-( ап—1,п—3Л1,п—3 + ап—1,п—2 )
ап, п—2
] = 1, п (21)
— ап—1,п—3Ук—3, ]----( ап—1,п—3Л1,п—3 + ап—1, п—2 )
ап, п—2
Уп—1 = Т-^
ап п—1 / \
ап—1,п—3Л2,п—3 + ап—1,п—1--'-( ап—1,п—3Л1,п—3 + ап—1,п—2 )
ап, п—2
Таким образом, получен алгоритм параллельного вычисления. По этому алгоритму сначала вычисляем коэффициенты прогонки вперед (14), (15), (18), (21). Далее по известной строке хп по формуле (20) получим хп—1, а по формулам прогонки назад (20), (13) - строки решения хк, к = п — 2,1. Имея уравнение (1) с матрицами порядка 2п + 1 с известной строкой хп+1, первый процессор вычисляет строки с 1 по п сверху вниз (для него последней является строка с номером п + 1). Второй процессор вычисляет строки с 2п + 1 по п + 2 снизу вверх (для него последней является строка с номером п + 1).
Рассмотрим тестовый пример (22), в котором коэффициенты матрицы взяты из работы [5, с. 73, формула (34)].
х
п
Здесь первая и последняя строка пятидиагональной матрицы системы содержат по 3 ненулевых элемента, вторая и предпоследняя строки - по 4 ненулевых элемента, остальные строки - по 5 ненулевых элементов.
10 3 1 6 2 3 0 0 0 0
1 6 10 3 1 6 2 3 0 0 0
2 1 10 1 2 0 0
3 6 3 6 3
0 2 1 10 1 2 0
3 6 3 6 3
0 0 2 1 10 1 2
3 6 3 6 3
0 0 0 2 3 1 6 10 3 1 6
0 0 0 0 2 3 1 6 10 3 _
-7 -31 2 -7 -31 2 -7 -31 2 -7
-15 -6 -15 -6 -15 -6 -15
-4 -11 -4 -11 -4 -11 -4
-11 -4 -11 -4 -11 -4 -11
-4 6 -4 6 -4 6 -4
-15 -6 -15 -6 -15 -6 -15
-7 -31 2 -7 -31 2 -7 -31 2 -7
(22)
Программа на FORTRAN с учетом алгоритма (14), (15), (18), (21), (20), (13) с известной последней строкой решения xn = (3,6,3,6,3,6,3), возвращает остальные строки решения, записанные в таблице 2.
Таблица 2. - Решение, полученное программой с использованием алгоритма (14), (15), (18), (21), (20), (13), xn = (3,6,3,6,3,6,3)
i/j 1 2 3 4 5 6 7
1 3.000000000 6.000000000 3.000000000 6.000000000 3.000000000 6.000000000 3.000000000
00000 00000 00000 00000 00000 00000 00000
2 6.000000000 3.000000000 6.000000000 3.000000000 6.000000000 3.000000000 6.000000000
00000 00000 00000 00000 00000 00000 00000
3 3.000000000 6.000000000 3.000000000 6.000000000 3.000000000 6.000000000 3.000000000
00000 00000 00000 00000 00000 00000 00000
4 6.000000000 3.000000000 6.000000000 3.000000000 6.000000000 3.000000000 6.000000000
00000 00000 00000 00000 00000 00000 00000
5 3.000000000 6.000000000 3.000000000 6.000000000 3.000000000 6.000000000 3.000000000
00000 00000 00000 00000 00000 00000 00000
6 6.000000000 3.000000000 6.000000000 3.000000000 6.000000000 3.000000000 6.000000000
00000 00000 00000 00000 00000 00000 00000
7 3.000000000 6.000000000 3.000000000 6.000000000 3.000000000 6.000000000 3.000000000
00000 00000 00000 00000 00000 00000 00000
Сравнение таблицы 2 и решения примера (22) показывает, что алгоритм (14), (15), (18), (21), (20), (13) с одной известной строкой хп = (3,6,3,6,3,6,3) является точным методом [3], решаемым за конечное число арифметических действий (совпадают 15 значащих цифр у всех элементов неизвестной матрицы).
Рассмотрим алгоритм решения задачи (1) с пятидиагональной матрицей одним ядром процессора в случае, если все строки решения неизвестны. В системе уравнений (12) третье, четвертое, пятое разностные уравнения содержат соответственно 5, 4, 3 разностных слагаемых. Поскольку решение третьего уравнения (12) имеет вид (13) и содержит два разностных слагаемых и одно постоянное слагаемое, то четвертое и пятое уравнения имеют решение на одно, на два разностных слагаемых меньше соответственно:
п-2 л п-1 , л п ,
х = Х1,п-2х +Х2,п-2х +Уп-2
хп-1 =\п-1 хп +Уп-1 (23)
х .
Подставив в последнее уравнение системы (12) первых две формулы из системы (23), получим а Х-2 +а .хп-1 + а X = Г о а , (X ,хп-1 + Х, X + v ,) + а ,хп-1 + а X = Г о
п,п-2 п,п-1 п,п л п,п-2 у 1,п-2 2,п-2 п-2/ п,п-1 п,п ^
о( а -X, , + а ,) хп-1 +(а -X, , + а ) хп = Тп - а -V , о
у п,п-2 1,п-2 п,п-1 / у п,п-2 2,п-2 п,п у •■> п,п-2 п-2
«(a А, , + a ,)(L ,х" +v ,) + (a А, , + a )х" = f" -
У n,n-2 1,n-2 n,n-1 / у 1,n-1 n-1 j у n,n-2 2,n-2 n,n / J
«((an,n-2^1,n-2 + an,n-1 )4n-1 + an,n-2^2,n-2 + an,n ) X" = f" - an,n-2Vn-2 - Vn-1 (an,n-2^1,n-2 + an,n-1 )
n— 2 i.n-z п,п-1 / i.n-i n,n-2/^2,n-2 1 an,n )х f an,n-2Vn-2 Vn-1 (an,n-2 1,n-2 + an,n- ' '' n ^ - an,n-2vn-2,j - Vn-1,j ( an,n-2^1,n-2 + an,n-1 ) . :-
« Xj =T-T---Чт, j = 1, n ; (24)
(( an,n—2\n—2 + an,n-1 ) 1 + an,n-2^2,n-2 + an,n ) n ./j - an,n-2Vn-2,j - Vn-1,j (an,n-2^1,и-2 + an,n-1 ) . :--,„
v" = 77-;-ь-n-\' j =1 n. (25)
((an,n-2Kl,n-2 + an,n-1 )Л1,И-1 + an,n-2K2,n-2 + an,n )
Решение, полученное программой на FORTRAN с использованием алгоритма (14), (15), (18), (21), (24), (20), (13), приведено в таблице 3.
Таблица 3. — Решение, полученное программой с использованием алгоритма (14), (15), (18), (21), (24), (20), (13)
ij 1 2 3 4 5 6 7
1 3.000000000 6.000000000 3.000000000 6.000000000 3.000000000 6.000000000 3.000000000
00000 00000 00000 00000 00000 00000 00000
2 6.000000000 3.000000000 6.000000000 3.000000000 6.000000000 3.000000000 6.000000000
00000 00000 00000 00000 00000 00000 00000
3 3.000000000 6.000000000 3.000000000 6.000000000 3.000000000 6.000000000 3.000000000
00000 00000 00000 00000 00000 00000 00000
4 6.000000000 3.000000000 6.000000000 3.000000000 6.000000000 3.000000000 6.000000000
00000 00000 00000 00000 00000 00000 00000
5 3.000000000 6.000000000 3.000000000 6.000000000 3.000000000 6.000000000 3.000000000
00000 00000 00000 00000 00000 00000 00000
6 6.000000000 3.000000000 6.000000000 3.000000000 6.000000000 3.000000000 6.000000000
00000 00000 00000 00000 00000 00000 00000
7 3.000000000 6.000000000 3.000000000 6.000000000 3.000000000 6.000000000 3.000000000
00000 00000 00000 00000 00000 00000 00000
Замечание 2. При решении матричного уравнения (22) с квадратной матрицей порядка п = 151 решение, как и в таблице 3 (при п = 7), алгоритмом (14), (15), (18), (21), (24), (20), (13) возвращается программой с двойной точностью за конечное число арифметических операций.
Последние уравнения (24), (25) согласуются с последним уравнением системы (23). В двух приведенных примерах (11), (22) шаблоны трех- и пятидиагональных матриц используются для аппроксимации дифференциального оператора Пуассона, из-за чего сумма весовых коэффициентов шаблона равна нулю (т.к. производная константы есть ноль) [3]. Разностные схемы для лапласиана на шаблоне «крест» и девятиточечном шаблоне имеют вид
Аик к = ТГ ["к к-1 + "к к+1 + "к-1,к + ик+1,к - 4ик ,к ] :
А", , = —
к,к h2
1/ ч 2, ,10
~ ("к ,к-1 + "к,к+1 + "к-1,к + "к+1,к ) + ~ ("к-1,к-1 + "к+1,к-1 + "к-1,к+1 + "к+1,к+1) ~ и<-
[5].
Поэтому центральный (диагональный) коэффициент имеет знак, противоположный знакам других коэффициентов шаблона (недиагональным коэффициентам строки матрицы). Из приведенных примеров видно, что диагональный элемент имеет максимальный модуль. Наименьший модуль 1/6 коэффициента расположен в узлах, удаленных от центра на шаг по одной координатной прямой. Промежуточное значение 2/3 находится в узлах, удаленных на шаг по двум координатным прямым. Выразим все сказанное в виде условий. Для удобства введем обозначения:
ft =
a12 > q2 = a13 a1,3
, z =
a11 a11 a1,2
Для пятидиагональной матрицы Теплица потребуем нестрогое двойное диагональное преобладание ее элементов
2(1 a
> 4 (| a
2 I) < I , = 3 n - 2 « 2 (2| ak,*+11 + 2 |a*,*+2 I)
<\a,r,r «
M
1 111
< - « q1 + q2 < -, q2 > q1 « q1 <- > q2 > - .
4 i2 ^^ ^848
a
Теорема 2 (о корректности алгоритма прогонки (14), (15), (18), (21), (24), (20), (13). Пусть на пяти-диагональную матрицу Теплица А уравнения (1) наложены условия:
1) А - симметрическая ц . = а.,,,, ] = 1,п, аи~\ = а,-1, , = 2, п — 1, ц ,._2 = аи+2, , = 3, и — 2;
2) элементы матрицы А имеют нестрогое двойное диагональное преобладание:
0 < 21аЬЛ —2 | < 2 (|акк—2 I + \ак,к—11 + |аА,к+11 + |аА,А +2 |) < |аА,А | ,Vk = 3П — 2 , 0 < 2|а13| < 2(|ц 2| + Ц31) < 1ц , 0 < 2|а24| < 2(|а21| +|а2 3| + |а241)< |а2 21,
0 < 2| ап—1,п—з\ < 2 (| ап—1,п—з| + \ап—1,п—2\ + |ап—1,и |) < |ап—1,п—^ 0 < 2| ап,п—2 | < 2 (| ап,п—2 | + |ап,п—] ) < \ап,п\ ■
3) акМ1 • ак,к < 0 ак,к +2 • ак,к < 0 ■ а,
Тогда У г =
•*1,3
1,2
[14]:
21
1) 0<Хи<-91, г = 1,п — 1,0<А,2_, <— 92, г = 1,п — 2;
2) формулы прогонки (24), (25), (21), (20), (18), (15), (14), (13) - корректны.
Доказательство проведем по индукции. Левые части условий 2 Теоремы 2 обеспечивают корректность формул (14), (19) и ненулевые элементы крайних диагоналей матрицы Теплица А, а следовательно, ненулевые диагональные элементы матрицы А.
,(14)
1) База индукции Хц =
4 (14)
= 91 <" 91, |Х2Д| =
21 3) 3) = 92 < — 92, Хц > 0, Х2Д > 0 проверена.
Далее числитель и знаменатель формул (15) делим тождественно на число ак
ы=
а13а21 а23а11
а22 а11 — а12а21
1+
1
< ап а23 + а1з а21 = 9 + 99 < 4 о Щ и ! 69 < 4 = ! 3), ........ 1 — 912 3У1 1 — а2 , 1 3 ^
1 —
64
I I С') Х2,2 =
92 21 21/ , , , 1 63 20
—^Т < — 92 о —(1 — д12 )> 1 о 1 — 91 > 1--= —> — .
1 — 92 20 20 у ' 64 64 21
21
Индуктивный переход. Пусть выполнены условия Х1 . < — 9, , = 1, к — 1, Х2< — 9, , = 1, к — 1.
Х1> 0, Х2> 0. Тогда
1Хи| < I
(18) 91 + 91 Iх 2к—1 | + 92 Iх 1к—2 Х2к—1|
21 4•21
91 + 91— 92 + 9192
— 92 |Х1к—2Х1к—1 | — 92 |Х2к—2 | — 91 |Х1к—11
91 7 . А
20
3 • 20
16 2 21 2 4 2
1— — 9291 — 92 — 91
< 91
,21 7 2
1 + — 92 + —92
20 2 5 2
16 1 2 21 2 4 2
1----9.--92 — 9,
V 9 4 20 3
4
< 3 91 о
о
, 21 7 2
1 + — 92 + - 92
20 2 5 2
Л
, 16 1 2 21 2 4 2
1----91--92 — 91
9 4 20 2 3 1
4 21 7 2 4 16 2 16 2 7 2 14 2 21 64 2 1
< —о 1 +--9 + — 9. <---9.--9--9. о— 9. +--9 +--9, < —
3 20 5 3 27 9 5 5 20 27 3
14
12 Л
2 21 5 3
92 +---92 +—-
ч2 20 14 162 ,
64 2 1 14 9 ^ 14 ( 3 64 2 1658 829
+ — 912 <- +---, Е ■— I 92 ^ — 912 <-=-. (26)
27 3 5 256 5 Г2 161 2у 1 3840 1920
£
а
а
12
13
2 .
а24а11
5
Неравенство (26) определяет внутреннюю область эллипса E с центром (0, -3/16):
1"
1)q = 0, q -I27 = 0.375 > I q e [0,0.375] з
^ I 1,1
3840 5 256) 64
,21 7 -
1 +--q2 +— q 2
20 2 5 2
0,
Vk = 3, n - 2 ,
q, 1
.20'5
21 1 7 1
1 +---+--- .
-20 8 5 64-* 1.202 < 4 = 1,3(3),
16 1 21 1 4 1 3
, 16 2 21 2 4 2 , 16 1
1---q2q.--q2--q1 1----r--
9 21 20 3 1 9 83 20 64 3 64
,21 7 2
1 +--q 2 + — q 2
20 2 5 2
q1
l+21-1 + 7^-1 . 16 120 5215 125 4 1 * 1 327353 < 4 = 1,3(3) •
'16 2 21 2 4 2 , 16
1---q2q1--q2--q1 1--
9 21 20 2 3 1 9 5 - 400 20 25 3 400
Эллипс (26) представляет выпуклое множество, поэтому весь отрезок прямой ^ + д2 = — между точ-
ками | 1,11 z = 1, \ Iz = 4 целиком расположен внутри эллипса [10, с. 33] Vz =
: [1,4]. Если
К-2, hk-1, 4-1, 4k-2 > 0 ^ SiSn ) = sign
I I (18)
4k =
( q1 + q1^2k-1 + qTk\k-2^2k-1 ^
I- q2^1k-2^1k-1 - q2^2k-2 - q1^1k-1
>;2) Чг
= +1, т.е. 0<A,U <-qui = 1,n-2,
ak,k + ak ,k-2^1k-2^1k-1 + ak ,k-2^2t-2 + ak,k-Att-1
, 16 2 21 2 4 2
1-q2 — q1--q2 — q1
9 20 3
21
< — q2 О 20 2
1
-<-
21
16 2 21 2 4 2 20
1-q2— q1--q2 — q1
2 9 1 20 2 3 1
, 1 16 2 21 2 4 2 20 ^ 16 2 21 2 1
1---qj2--q22--q2 >—, E : — q2 + — q22 < —
4 9 20 3 21 9 20 21
1 9 3
1)q = 0,q{ <---,q <-./—* 0.1637q e
21 16 4
1
8
L 3 /г] 4 1 ]
0,-J— 4\21 з 0,8
(27)
1 1
1
Ч2
1
20,5
'16 2 21 2 4 2 ,
1---q2 q1--q2--a 1 - ,
9 21 20 2 3 1 9 83 20 64 3 64
——-r^-—- * 1.042 < — = 1,05(0) .
16 1 21 1 4 1 20
1
Ч2
—------— * 1.048 < — = 1,05(0) •
16 1 21 1 4 1 20
'16 2 21 2 4 2 ,
1---qq--q2--q1 1 -
9 21 20 3 1 9 5 - 400 20 25 3 400
Поскольку эллипс (27) представляет выпуклое множество [10, с. 33], то весь отрезок прямой ^ + д2 =
1
i 1 11 „ | 11,,
между точками | —, — | z = 1, | — , — | z = 4 целиком расположен внутри эллипса Vz =
е [1,4]. Если
4-2, 4k-1, 4*-1, 4k-2 > 0 ^ Si8n (4k ) = Si8n
q2
1 q2^1k-2 4k-1 ' 2k-2 - q^ 1 k -1
= +1, т.е.
21
0<^, . < — q, i = 1, n - 2 . u 20
Первая часть теоремы 2 доказана Vz =
[1,4] •
a
1.3
a
1,2
-a
2,k
Л,
2,k
a
1,3
a
1,2
a
1,3
2) Доказательство второй части (корректность формул прогонки).
Найдем условие, при котором знаменатель формул (21) сохраняет знак, что обеспечит корректность формул (21). С учетом условий Теоремы 2 получим
91 + 92 < ^ =
=е[м]," 91 < 8л 92 < 8,91 < 20л 92 < 5,
1п —1,п —3Х2,п —3 + ап—1,п—1 (ап—1,п—3Х1,п— 3 + °п—1,п—2 )
2п—1,п—3Х2,п—3 + ап— 1,п— 1 (ап—1,п—3Х1,п— 3 + °п—1,п—2 )
> 0 о
п—3 2,п—3 п—1,п—1 п—1,п—3 п—1,п—2 п—1,п—3 3 п—1,п—2
■1,«—3Х1,п—3 — ап—1,п—22\ > 0 о |922Х2,И—3 — 92 — 9г91Х1,И—3 — > 0 о
2 2 л 21 3 Г 21 1 1 479 л ^
о 92 + 92 91Х1,п—3 + 91 — 92 Х2,п—3 > 92 — — 92 > 921 1 —^^1 = > 0, У92 £
1 1
85
(21)
IV1| =
—ап—\,п + ( ЯИ—1,И—3Х1,И—3 + аи—1,И—2 )
+ 1,и—1 — ^ (Щ—3Х^ 3 + «„—1,—2)
aи-1.и-2 + (aи-1.и-3Х1.и-3 + aи-1.и-2 )
п х '
—пА-3 + " ^^(а^Х^ + 0-^2)
(-а , 2а , ,+ а , ,а , А, ,+ а , ,а , 2)
( ап-1,И-32Х2
-1.и-3Х1.и-3 ап-1,п-2 )
-92Х1п-3 - 9192 - 91
92 Х2,И-3 ^2 3 -<?1
929^ + 9192 + 91
9291- + 91
91
, 7
1 + 3 92
92 + 9192Х1,-3 + 912 — 922Х2,-3 " 92 + 9192Х1„—3 + ^ — 923 Ц 92 | 1 + ^ — ^ 21
, 20 1 9 41 20
< 992
1 7 5 ^ 3
+1
1 +
21 1
202 (1/5) 20 25
22
9-—¿Чт и1.511* < 23 * = 1.53(3) 91.
92 1 + -___21 92 1 5 92 92
80 500
(28)
Таким образом, формулы (21) корректны, т.е. существуют конечное значение Х1и-1 и конечная норма вектора уп_х. Рассмотрим корректность формул (24), (25):
*/ =
1] ] УИ-1,] (aи.и-2Х1.и-2 + а„,И-1 )
((aИ.И-2Х1.И-2 + а„,„-1 )Х1,п-1 + а„,„-2Х2,„-2 + ап,п )
В последней формуле разделим знаменатель на диагональный элемент \апп | и оценим дробь по модулю:
|( ап,п-2Х1,,
'•и.и-2Х1.и-2 + ап,п-1 )Х,И-1 + ап,п-2Х 2,п-2 +
> 0 о |(-92Х1,И-2 - 91 )Х1,и-1 - 92Х2,п-2 + 1
(28)
= 1 — (92Х1,И—2 + 91 ) Х1,и—1 — 9-Х-.и—2 >
(28\ ( 4 123 а 2 21 92 2 23 9,2 21 1 92 1 23 1 21 Л
> 1 -I 9291 — + 91 I— -ф, — = 1--912--^-----> 1---------е-т--и 0,734 > 0 .
| 1 3 ) 15 92 20 45 15 92 20 25 45 64 15 64(1/8) 500
а
О
<
а
Таким образом, корректны формулы (24), (25), то есть ограничено значение \х"\.] = 1,п . А также
корректны формулы обратной прогонки, указанные ниже, т.к. конечность величин хч , Х.
1,п-1 ' II уп-1||
Х,
Х
к = 1, п - 2 была показана нами ранее
= хихк+1 + Х2Дхк +2 + V*, к = 1, п - 2 хп-1 = Х1 и_ хп + vи_1
Теорема 2 доказана. Оказывается, что при больших порядках (п > 50) матричного уравнения (1) коэффициенты прямой прогонки (18) имеют предельные значения, на что указывает распечатка коэффи-
циентов. Обозначим предельные значения Х1 дельными значениями, имеем:
к — Х, Х2,к ,
к —да к—
у. Заменяя все коэффициенты в (18) их пре-
(
Х1,к = -
_ак ,к+1 + ак,к-1Х2к-1 + ак ,к-2Х1к-2Х2к-1
ак ,к + ак ,к-2Х1к-2 Х1к-1 + ак ,к-2Х2к-2 + ак,к-1Х1к-1 )
Х 2,к =
ак ,к + ак ,к-2Х1к-2Х1к-1 + ак ,к-2Х2к-2 + ак ,к-1Х1к-1
ак,к+1 + ак,к-1 У + ак,к-2 ХУ ак,к + ак,к-2 Х + ак,к-2 У + ак,к-1 Х
У =-
ак,к + ак,к-2 Х + ак,к-2 У + ак,к-1 Х
Для симметрической матрицы Теплица имеем:
-?2
^ 41 + <?1 У + 42 ХУ Л
1 + ^Х +4 У +
■ У = -
1 + д2Х + д2У + 4 х
=■
. 42 =■
(29)
Теорема 3. Нелинейная система уравнений (29) численно разрешима методом Ньютона - Зейделя по итерационным формулам с диагональными элементами матрицы Якоби, полученными в работе [15]:
У'+1 = У'
х1 + 4 (х! ) + 2д2Х!У* + 4 (х11 ) + 4 + 4У 1 + 34 (х1) + + 2д1х1
У * + 42 У * ( X1+1 )2 + Ч2 (У * )2 + 4Х+1У * + 42 1 + 4 (х1+1 )2 + 2д2У1 + д,х1+1 '
* = 0,1,2,
^ = 0,1,2,
(30)
Доказательство. В работе [15, с. 14] показано, что для решения нелинейной системы уравнений (29) достаточно переписать нелинейные уравнения (29) в каноническом виде с нулевой правой частью, т.е. получить функции четырех переменных (4,4 - параметры, х*у - независимые переменные). Имеем:
/ (х1, У1,41,42) = х1 + 42 (х1 )3 + 242Х У + 4х (х1 )2 + 4Х + 4ХУ = 0 ,
/2 (х1, У1,41,42) = У1 + 42У1 (х1 )2 + 42 (У* )2 + 4ХУ + 42 = 0 . Их частные производные - диагональные элементы матрицы Якоби, согласно [15, с. 14] получим
/х' (х1, У1,41,42) = 1 + 342 (х1 )2 + 242У + 24Х ,
г I / 1+1 1 \ 1 . / 1+1 \2 . о 1 I 1+1
ЛУ (х , У, 42 ) =1 + 42 (х ) + 242У + 41х . Согласно [15, с. 14] имеем систему итерационных уравнений:
У- , я
[(х , У , 41,42
) _ , х + 42 (х* )3 + 242Х''У + 41 (х* )2 + 41 + 41У
у1+1 = У1
/1х*' (х*, У, 41,42) 1 + 34 (х1 )2 + 24У1 + 24Х1
/2 (х1+1, У1,41,42) , У* + 42У* (х1+1 )2 + 42 (У* )2 + 41 х1+У + 42
/2 у' ( х1+1, У1
41,42)
= У
1 + 4 (х11+1) + 242У1 + 4Х1
1' = 0,1,2,
* = 0,1,2,
Теорема 3 доказана, поскольку формула (31) совпадает с формулой (30).
v
к
хп = V
п
х = -
а
а
а
а
к .к+1
к .к-1
к .к+2
к .к-2
к .к
к .к
к .к
к .к
Замечание 3. Для матричного уравнения (22) второго численного примера q1 = 1/6 3 = —~,
2/3 1
q2 =-= — программа с начальным значением x° = y° = 0.1, щ = и2 = 101 по формулам (18) дает ите-
—10/3 5
рационные значения коэффициентов прогонки Х126 = x = 6,6324473666998-10—2, ^226 = y = 0,2096722801763936 уже на 26 шаге, и далее коэффициенты прогонки имеют стационарные значения. Итерационные формулы (30) дают x = 6,6324473666998 • 10—2, y = 0,209672280176393 , x° = y0 = 0.1, т.е. разность между текущими значениями коэффициентов прогонки и их предельными значениями падает к нулю с двойной точностью по геометрической прогрессии за несколько десятков шагов. Кроме того,
h I = X = 6,6324473666998• 10—22 <-\q, I = — = 6,6(6)• 10—22, I i.«l 3^н 60
h I = y = 0,209672280176393 < — |q7| = — = 0,21(0), I 201 1 100
что не противоречит Теореме 2.
Рассмотрим другой параметрический случай q1 = — 1, q2 = — 1.
Итерационные формулы (29) дают X = 0,1492864354 y = 0,1298948902, x0 = y0 = 0.1. Кроме того,
— 4 4 — 21 21
h I = x = 0,1492864354 < -|q, I = — = 0,166(6),I = y = 0,1298948902 < — IqJ =-= 0,13125(0),
I 31^1 24 1 2-"l ^ 201 1 160
(28) 23 a f 1 1 Л (28) 23 q f 1 1 A
M = о,376148815 < -f = wx [ q1 = q2 = 5 )|Ч J =105 < vq 1 [ q1 == 8 * = 8 J,
что также не противоречит Теореме 2.
Ниже приведена программа для примера (22), написанная на FORTRAN, в которой переменные и массивы заданы с двойной точностью [7, 15].
Замечание 3. Для использования программы необходимо скопировать ее текст из файла *.pdf в файл формата *.docx, сохранить, удалить строки с информацией колонтитула и номером страницы, перенести его в новый проект FORTRAN *.f90.
program matprogonka;
integer(8),parameter::n1=7,n2=7;integer(8)::k,i,j; real(8)::a(n1,n2),f(n1,n2),nu(n1,n2),x(n1,n2),lamda(n1); real(8)::lamda1(n1+1),lamda2(n1+1),c1 ,c2,c3,c7,c4,c5,c6;
a(1,1)=-1d1/3d0;a(1,2)=1d0/6d0;a(1,3)=2d0/3d0;a(1,4)=0d0;a(1,5)=0d0;a(1,6)=0d0;a(1,7)=0d0; a(2,1)=1d0/6d0;a(2,2)=-1d1/3d0; a(2,3)=1d0/6d0;a(2,4)=2d0/3d0;a(2,5)=0d0;a(2,6)=0d0; a(2,7)=0d0;a(3,1)=2d0/3d0; a(3,2)=1d0/6d0;a(3,3)=-1d1/3d0; a(3,4)=1d0/6d0;a(3,5)=2d0/3d0;a(3,6)=0d0; a(3,7)=0d0;a(4,1)=0d0;a(4,2)=2d0/3d0; a(4,3)=1d0/6d0;a(4,4)=-1d1/3d0;a(4,5)=1d0/6d0;a(4,6)=2d0/3d0;a(4,7)=0d0; a(5,1)=0d0;a(5,2)=0d0;a(5,3)=2d0/3d0;a(5,4)=1d0/6d0;
a(5,5)=-1d1/3d0;a(5,6)=1d0/6d0;a(5,7)=2d0/3d0;a(6,1)=0d0;a(6,2)=0d0; a(6,3)=0d0;a(6,4)=2d0/3d0; a(6,5)=1d0/6d0; a(6,6)=-1d1/3d0;a(6,7)=1d0/6d0;
a(7,1)=0d0;a(7,2)=0d0;a(7,3)=0d0;a(7,4)=0d0;a(7,5)=2d0/3d0;a(7,6)=1d0/6d0;a(7,7)=-1d1/3d0; f(1,1)=-7d0;f(1,2)=-31d0/2d0;f(1,3)=-7d0;f(1,4)=-31d0/2d0;
f(1,5)=-7d0;f(1,6)=-31d0/2d0;f(1,7)=-7d0; f(2,1)=-15d0;f(2,2)=-6d0;f(2,3)=-15d0;f(2,4)=-6d0;
f(2,5)=-15d0;f(2,6)=-6d0;f(2,7)=-15d0;f(3,1)=-4d0;f(3,2)=-11d0;f(3,3)=-4d0;f(3,4)=-11d0;
f(3,5)=-4d0;f(3,6)=-11d0;f(3,7)=-4d0; f(4,1)=-11d0;f(4,2)=-4d0;f(4,3)=-11d0;f(4,4)=-4d0;
f(4,5)=-11d0;f(4,6)=-4d0;f(4,7)=-11d0;f(5,1)=-4d0;f(5,2)=-11d0;f(5,3)=-4d0;f(5,4)=-11d0;
f(5,5)=-4d0;f(5,6)=-11d0;f(5,7)=-4d0;f(6,1)=-15d0;f(6,2)=-6d0;f(6,3)=-15d0;f(6,4)=-6d0;
f(6,5)=-15d0;f(6,6)=-6d0;f(6,7)=-15d0;f(7,2)=-31d0/2d0;f(7,1)=-7d0;f(7,4)=-31d0/2d0;
f(7,3)=-7d0; f(7,6)=-31d0/2d0;f(7,5)=-7d0;f(7,7)=-7d0;
lamda1( 1)=-a(1,2)/a(1,1);lamda2(1)=-a(1,3)/a(1,1);
lamda1(2)=(a(2,1)*a(1,3)-a(2,3)*a(1,1))/(a(2,2)*a(1,1)-a(1,2)*a(2,1));
lamda2(2)=-a(2,4)*a(1,1)/(a(2,2)*a(1,1)-a(1,2)*a(2,1));do j=1,n2;
nu(1,j)=f(1,j)/a(1,1);nu(2,j)=(a(1,1)*f(2,j)-a(2,1)*f(1,j))/(a(2,2)*a(1,1)-a(1,2)*a(2,1));enddo; do i=3,n1-2;do j=1,n2;c1=a(i,i+1)+a(i,i-1)*lamda2(i-1)+a(i,i-2)*lamda2(i-1)*lamda1(i-2); c2=-a(i,i)-a(i,i-2)*lamda1(i-2)*lamda1(i-1)-a(i,i-2)*lamda2(i-2)-a(i,i-1)*lamda1(i-1);
c3=a(i,i-2)*lamda1(i-2)*nu(i-1 ,j)+a(i,i-2)*nu(i-2,j)+a(i,i- 1)*nu(i-1 ,j)-f(i,j); lamda1(i)=c1/c2;lamda2(i)=a(i,i+2)/c2;nu(ij)=c3/c2;enddo;enddo;
c4=(a(n1-1,n2-1)+a(n1-1,n2-3)*lamda2(n1-3)-(a(n1,n2-1)/a(n1,n2-2))*(lamda1(n1-3)*a(n1-1,n2-3)+a(n1-1,n2-2))); c6=a(n1-1,n1-3)*lamda1(n1-3)+a(n1-1,n1-2);
c5=(a(n1,n1)/a(n1,n1-2))*(a(n1-1,n1-3)*lamda1(n1-3)+a(n1-1,n1-2))-a(n1-1,n1);lamda1(n1-1)=c5/c4; do j=1,n2;nu(n1-1,j)=(f(n1-1,j)-a(n1-1,n1-3)*nu(n1-3,j)-f(n1,j)*c6/a(n1,n1-2))/c4;enddo; do j=1,n2;c4=f(n1,j)-a(n1,n1-2)*nu(n1-2,j)-nu(n1-1,j)*(a(n1,n1-2)*lamda1(n1-2)+a(n1,n1-1)); c5=a(n1,n1-2)*lamda1(n1-2)*lamda1(n1-1)+a(n1,n1-1)*lamda1(n1-1)+a(n1,n1-2)*lamda2(n1-2)+a(n1,n1); nu(n1j)=c4/c5;x(n1j)=nu(n1j);enddo;do j=1,n2;x(n1-1,j)=lamda1(n1-1)*x(n1,j)+nu(n1-1,j); enddo;do i=n1-2,1,-1;do j=1,n2;x(i,j)=lamda1(i)*x(i+1j)+lamda2(i)*x(i+2j)+nu(ij);enddo;enddo; do i=1,n1;do j=1 ,n2;print*,ij,x(ij);enddo;enddo; end program matprogonka;
( Au, u )
В работе С.К. Годунова [16, с. 233] показано, что множество отношений Релея —-является
(u, u )
хаусдорфовым множеством для любой матрицы А и выпукло. В примере [16, с. 233] матрица содержит три
"-1 2q
а множество Релея имеет вид эллипса. Отметим, что пятидиа-
различных ненулевых элемента
0 1
тональная матрица Теплица из примера (22) также имеет три различных ненулевых элемента и ее множества корректности формул прямой прогонки (26), (27) представляют собой эллипсы.
В статье «Обратная теорема Гамильтона»1 (теорема 9) и работе [8] возникает блочно-диагональная матрица - матрица Гессе от функции Лагранжа либо функции Гамильтона. Достаточные условия в теореме 9 численно возможно проверять с помощью векторного метода прогонки для диагональных матриц, описанного в данной работе.
В работе получены результаты:
1. Предложен векторный аналог решения матричных уравнений с трехдиагональной матрицей формулы (8), (9), (10).
2. В теореме 1 получены достаточные условия корректности алгоритма (8)-(10).
3. Рассмотрен алгоритм (14), (15), (18), (21), (20), (13), решения задачи (1) с пятидиагональной матрицей для параллельного вычисления 2 процессорами и алгоритм (14), (15), (18), (21), (24), (20), (13) решения этой же задачи одним процессором.
4. В теореме 2 получены достаточные условия корректности алгоритма (14), (15), (18), (21), (24), (20), (13).
5. Приведенные тестовые примеры и программа на FORTRAN показывают, что алгоритмы (8)-(10), (14), (15), (18), (21), (24), (20), (13) являются точными, так как дают программой значение неизвестной матрицы за конечное число операций, совпадающей поэлементно с точной матрицей в 15 значащих цифрах.
6. В теореме 3 получен численный алгоритм для предельных значений коэффициентов прогонки вперед, показано, что для двух параметрических случаев предельные значения подчиняются теореме 2.
ЛИТЕРАТУРА
1. Козлов, А.А. Преобразование подобия на множестве полукватернионов / А.А. Козлов, К.С. Сурав-нева, И.Л. Жалейко // Вестн. Полоц. гос. ун-та. Сер. С, Фундам. науки. - 2019. - № 4. - С. 115-123.
2. Козлов, А.А. Множество полуоктав / А.А. Козлов // Вестн. Полоц. гос. ун-та. Сер. С, Фундам. науки. -2016. - № 12. - С. 75-85.
3. Бахвалов, Н.С. Численные методы / Н.С. Бахвалов, Н.П. Жидков, Г.М. Кобельков. - 7-е изд. - М. : БИНОМ. Лаборатория знаний, 2011. - 636 с.
4. Бахвалов, Н.С. Численные методы в задачах и упражнениях / Н.С. Бахвалов, А.В. Лапин, Е.В. Чижон-ков. - М. : БИНОМ, 2010. - 240 с.
5. Пастухов, Д.Ф. Аппроксимация уравнения Пуассона на прямоугольнике повышенной точности / Д.Ф. Пастухов, Ю.Ф. Пастухов // Вестн. Полоц. гос. ун-та. Сер. С, Фундам. науки. - 2017. -№ 12. - С. 62-77.
6. Колмогоров, А.Н. Элементы теории функций и функционального анализа / А.Н. Колмогоров, С.В. Фомин. - М. : Наука, 1976. - 543 с.
7. Бартеньев, О.В. Фортран для профессионалов. Математическая библиотека IMSL: Ч. 1. - М. : ДИАЛОГ : МИФИ, 2001. - 437 с.
1 Пастухов, Ю.Ф. Обратная теорема Гамильтона / Ю.Ф. Пастухов, Д.Ф. Пастухов // Вестн. Полоц. гос. ун-та. Сер. С, Фундам. науки. - 2019. - № 12. - С. 86-100.
8. Пастухов, Ю.Ф. Свойства функции Гамильтона в вариационных задачах со старшими производными / Ю.Ф. Пастухов, Д.Ф. Пастухов // Вестн. Полоц. гос. ун-та. Сер. С, Фундам. науки. - 2019. -№ 4. - С. 137-153.
9. Александров, П.С. Введение в теорию размерностей / П.С. Александров, Б.А. Пасынков. - М. : Наука, 1973. - 577 с.
10. Галеев, Э.М., Краткий курс теории экстремальных задач / Э.М. Галеев, В.М. Тихомиров - М. : Изд-во Моск. ун-та, 1989. - 204 с.
11. Пикулин, В.П. Практический курс по уравнениям математической физики : учеб. пособие / В.П. Пи-кулин, С.И. Похожаев. - М. : Наука, 1995. - 224 с.
12. Пастухов, Д.Ф. Оптимальный порядок аппроксимации разностной схемы волнового уравнения на отрезке / Д.Ф. Пастухов, Ю.Ф. Пастухов, Н.К. Волосова // Вестн. Полоц. гос. ун-та. Сер. С, Фундам. науки. - 2018. - № 12. - С. 60-74.
13. Пастухов, Д.Ф. К вопросу о редукции неоднородной краевой задачи Дирихле для волнового уравнения на отрезке / Д.Ф. Пастухов, Ю.Ф. Пастухов, Н.К. Волосова // Вестн. Полоц. гос. ун-та. Сер. С, Фундам. науки. - 2018. - № 4. - С. 167-186.
14. Пастухов, Ю.Ф. Необходимые условия в обратной вариационной задаче / Ю.Ф. Пастухов // Фундаментальная и прикладная математика. - 2001. - Т. 7, вып. 1. - С. 285-288.
15. Пастухов, Д.Ф. Численные методы. Лекции. Численный практикум [Электронный ресурс] / Д.Ф. Пастухов, Ю.Ф. Пастухов - Новополоцк : ПГУ, 2019. - 227 с. - Режим доступа: http:/elib.psu.by:8080/ handle/123456789/21502. - Дата доступа: 15.04.2019.
16. Годунов, С.К. Современные аспекты линейной алгебры / С.К. Годунов. - Новосибирск : Науч. книга, 1997. - 407 с.
17. Вакуленко, С.П. Способы передачи QR кода в стеганографии / С.П. Вакуленко, Н.К. Волосова, Д.Ф. Пастухов // Мир транспорта. - 2018. - Т. 16, № 5 (78). - С. 14-25.
18. Пастухов, Д.Ф. Некоторые методы передачи QR кода в стеганографии / Д.Ф. Пастухов, Н.К. Волосова, А.К. Волосова // Мир транспорта. - 2019. - Т. 17, № 3 (82). - С. 16-39.
19. Волосова, Н.К. Применение преобразования Радона в стеганографии // Некоторые актуальные проблемы современной математики и математического образования : сб. материалов науч. конф. , Герце-новские чтения - 2018, СПб., 9-13 апр. 2018 г. / Рос. гос. пед. ун-т им. А.И. Герцена. - СПб., 2018. -С. 234-238.
20. Волосова, Н.К. Преобразование Радона и уравнение Пуассона в компьютерной стеганографии // Между-нар. конф. по дифференциальным уравнениям и динамическим системам : сб. ст. - Суздаль, 2018. - С. 61.
21. Пастухов, Д.Ф. Минимальная разностная схема для уравнения Пуассона на параллелепипеде с шестым порядком погрешности / Д.Ф. Пастухов, Ю.Ф. Пастухов, Н.К. Волосова // Вестн. Полоц. гос. ун-та. Сер. С, Фундам. науки. - 2019. - № 4. - С. 154-173.
22. Модифицированное разностное уравнение К.Н. Волкова для уравнения Пуассона на прямоугольнике с четвертым порядком погрешности / Н.К. Волосова [и др.] // Евразийское Научное Объединение. - 2019. -№ 6-1 (52). - С. 4-11.
Поступила 17.09.2019
VECTOR ANALOGUE OF THE METHOD PROGONKI FOR DECISION THREE AND FIVE DIAGONAL MATRIX EQUATIONS
N. VOLOSOVA, K. VOLOSOV, A. VOLOSOVA, D. PASTUHOV, Y. PASTUHOV
An algorithm is proposedfor a vector analogue of sweep for solving arbitrary matrix equations with square three- andfive-diagonal matrices for a finite number of arithmetic calculations. We prove sufficient conditions for the correctness of vector sweep formulas for arbitrary three-diagonal matrices (Theorem 1) and sufficient conditions for five-diagonal symmetric Toeplitz matrices (Theorem 2). The above program and two examples show that these algorithms are accurate. A numerical algorithm is proposed for finding limit values for forward sweep coefficients (Theorem 3), and it is shown that the obtained numerical limit values do not contradict Theorem 2.
Keywords: vector analogue of the method of the racing, three and five diagonal matrixes, Toeplitz matrix, convex sets, numerical methods of mathematical physics, parallel calculations.
