CC BY
52
АЛГОРИТМЫ И ВЫЧИСЛИТЕЛЬНЫЕ МЕТОДЫ
УДК 519.612
А. О. Синюхин, А. А. Буздин
МЕТОД ПОСТРОЕНИЯ АЛГОРИТМОВ ПОЛНОЙ РЕДУКЦИИ ДЛЯ РЕШЕНИЯ СИСТЕМ УРАВНЕНИЙ
Рассмотрен подход к реализации метода полной редукции для системы уравнений с блочной трехдиагональной матрицей произвольной размерности. Выведены формулы, позволяющие производить исключение неизвестных в произвольном порядке. Изложены алгоритмы для первой и второй краевой задачи.
An approach to the method of cyclic reduction for the system of equations with block-tridiagonal matrix of arbitrary dimension is presented. Formulas for elimination of unknowns in arbitrary order are obtained. Algorithms for the first and second boundary problems are described.
Ключевые слова: система линейных уравнений, блочная трехдиагональ-ная матрица, метод полной редукции.
Key words: system of linear equations, block-tridiagonal matrix, method of cyclic reduction.
Метод редукции является одним из эффективных прямых методов решения СЛАУ с блочной трехдиагональной матрицей
где С — линейные функции от одной матрицы. Такие СЛАУ возникают, например, в результате конечно-разностной дискретизации краевых задач для уравнения вида
Ранние алгоритмы метода полной редукции имели следующие недостатки: они были неустойчивы и количество блоков системы не могло быть произвольным. От недостатков удавалось избавиться путем использования искусственных приемов, затрудняющих понимание и реализацию алгоритмов [1; 2]. В данной работе предложен подход к методу редукции, позволяющий с единой точки зрения естественным образом получить легко реализуемый и устойчивый алгоритм для произвольного числа узлов и различных краевых задач.
С БЛОЧНОЙ ТРЕХДИАГОНАЛЬНОЙ МАТРИЦЕЙ
43
Введение
A = blocktridiag{ -1, Ct, -1},
(1)
© А. О. Синюхин, А. А. Буздин, 2015 Вестник Балтийского федерального университета им. И. Канта. 2015. Вып. 10. С. 43-51.
1. Вычисление обратной матрицы
44
Определим полином А р (x) следующим образом:
АР (х) =
С, (х) -1 0
-1
-,+11
С„+п (х) -1
-1 С 0
Р-1(х) -1
0
0 -1 СР (х)
(2)
Тогда через Ар (С) обозначим соответствующий матричный полином, где С е Ш:х: . При этом будем полагать Ар+1 (С) = I.
Матричные полиномы Ар (С) обладают свойствами, которые легко получить из разложения определителя (2) по строке с номером и
1) А ? = С А |+1 -А £2,
2) А N = -А 1- 2 А N1 + С, А1-1А N1 -А1-1 А* 2.
Сформулируем и докажем лемму, которая будет использоваться для дальнейшего вывода расчетных формул. Для сокращения количество индексов в лемме и в последующих рассуждениях обозначим через
А / а: =Ак (А: )-1 =(А:)-1 А .
Лемма 1. Пусть А - матрица вида (1). Тогда элементы обратной матрицы равны
А. = (АN)-1 А1-1А^!,, ^ А. = (АN)-1 А1-1АГ+ 1, , > ;. (3)
Доказательство. Покажем, что ААШ — единичная матрица. Используя свойство 2), получаем: (АА-1).. = -Д_1,, + СА,i - А+1,, = I. Для случая , > применяя свойство 1), получим
(ДА-1) = -А . . + СА . . - А+ . = (-А1-1 А14 + С. А1-1А1+1 - А1- 1А1+, ) =
\ Д . .-1,. . .1+1,. Д^ 1 I I 1 г+1 1 1+2 у
А1
.А-
АN
•(-А А N1 -А N2 ) = 0.
Для оставшегося случая г < . равенство нулю проверяется совершенно аналогично. □
Отметим, что формулы (3) являются обобщением выражений, которые приведены в [3, с. 61 — 62] и [4] для первой краевой задачи.
2. Метод исключения неизвестных в произвольном порядке
Выведем формулы, позволяющие вычислять коэффициенты уравнений и правую часть после исключения произвольной группы неизвестных (м;+1 ... иг-1) из системы (1). Перепишем ее преобразуемую часть: -н,-1 + Си-(I 0 ... о)ис = ^, _
-(I 0 ... 0)) + АДС-(0 ... 0 I) = £, (4)
-(0 ... 0 I)ис + Сгиг -иг+1 = /Т.
В этих формулах Uc = (u;+1 ... ur), Fc = (f;+1 ... fr). Выразим блок Uc из 2-го уравнения системы (4) и поставим в 1-е и 3-е, получим
_и,_1 + с(\ _= ¿(2), _+ C[2)ur _ur+1 = f \ (5)
где коэффициенты уравнений (5) рассчитываются по следующим формулам:
R(2)=_(I 0 ... 0))(0 ... 0 I),42)=_(0 ... 0 i))(I 0 ... 0)T, (6) с(2)= с_(I 0 ... 0) ( I 0 ... 0)T,cf = Cr_(0 ... 0 i))(0 ... 0 I), )
а правые части — по формулам
f>= f + (i 0 ... 0) A_1 Fc, fr(2>= fr + (0 ... 0 i )) 11. (7)
В случае, когда исключаемый блок расположен в начале или конце, можно формально полагать l = 0 или r = N + 1 соответственно.
Лемма 2. (Общий вид коэффициентов уравнения и правой части при исключении произвольной группы неизвестных). Пусть матрица СЛАУ имеет вид (1). Тогда в результате исключения неизвестных с номерами l +1, ..., r _1 коэффициенты в (5) выражаются формулами
R( 2)= , с( 2>=^_1, cr2>=#1-, (8)
; r дr ' ; дr ' r Дr _1 д;+1 д;+1 д;+1
а правые части в (5) — формулами
¿(2)=/+х Дй /, /(2)=/+х /. (9)
I=+1 д+1 I=+1
При этом, поскольку результат не зависит от способа исключения, (8) и (9) не зависят от того, в каком порядке это исключение производилось.
Доказательство. Выражения для 1(2), К,(2) а также правой части можно получить, если подставить необходимые элементы матрицы А-1 в
формулы (6), (7). Для расчета С2', С^ требуется кроме этого применить свойство 1) после приведения к общему знаменателю:
с(2' = с _ дг_1 / дг_1 = дr_1 / дr_1 с(2) = с _ Дr_2 / Дr_1 = Дг / Дг_1
И Н д;+2 / д;+1 д; / д;+1' r r д;+1 / д;+1 д;+1 / д;+1 •
□
Теорема. Уравнение, связывающее неизвестные и;, ис, и (1), имеет вид
1 щ Uc ^ur = f^. (9.1)
Дс_1 Ac_: Ar_: c Ar_:
-1;+1 +1" c+1 "c+1
В случае, когда исключается блоки (щ ... u_) и (u+1 ... u_), имеем
Ч 1 1 Ар)
—L_и--и — tyt '
Дг 1
Tc^ uc _Д_Т ur = , (9.2)
Д1 Дc+1 Дc+1
а когда блоки (u;+ ... u_ и (u+1 ... un), то
1 ДN
дс -1 и + Д-гД^ис = /Г. (9.3)
д+1 д;+1 дс+1
Доказательство. Для вывода формул (9.1) —(9.3) достаточно дважды воспользоваться леммой 2, исключив блоки неизвестных (иг+1 ... ис-1) и (Мс+1 ... иг-1), а затем применить свойство 1). □
Вычисление правой части по формуле (9) нецелесообразно, поскольку требует большое количество арифметических действий. Вме-
45
сто этого построим процесс преобразования правых частей следующим образом. Пусть к определенному моменту оказались исключенными
блоки (мг+1... Uc-i) и (uc+i ... Mr_i), а правые части f/p), , f^' известны. Исключим теперь неизвестную uc, выразив ее из формулы (9.1) и подставив в формулы для соседних уравнений, что повлечет преобразование правой части этих уравнений по формулам
f(+1] = f( +§-[, fp+1) = Р +¿г) (10) Ai+1 Ai+1
Выражения (10) позволяют описать эволюцию правой части при ис-
- ключении одной произвольной неизвестной. Таким образом, прямой
46 ход сводится к заданию порядка исключения и применению для пересчета правых частей формул (10). Процесс исключения повторяется до тех пор, пока в СЛАУ не останется одно уравнение с одной неизвестной
an
A1 U = f(n) A1-1Af+1 Uc " fc ,
решая которое, получаем uc.
Обратный ход метода реализуется по формулам (9.1) , (9.2) и (9.3) как
рекуррентный процесс. Если известны ui, ur и правая часть f p), то из выражения (9.1) неизвестная uc может быть вычислена следующим образом:
Ac-1 A r-1 Ac-1 A'-1
u =Ai±1A±L f( r) +Ai+2_ +A^ . (11)
c АГ-1 -1 c АГ-1 r АГ-1 1 4 '
Ai + 1 Ai + 1 Ai + 1
Для реализации алгоритма по (10), (11) необходимо находить произведения отношений полиномов на некоторый вектор: Ak / A™f . Как и в случае скалярного аргумента, отношение матричных полиномов Ak / Am можно разложить на элементарные дроби, если полином в знаменателе не имеет кратных корней и степень числителя меньше степени знаменателя. Разложение определяется следующими формулами:
Ai
Am = £а,(С -V ) (12)
n i
= <(X,) , (13,
i <A-->'W п,,,(X,-X,)
где X, — корни полинома A^1 (x). Значение Ak (X,) можно посчитать по рекуррентной формуле, которая вытекает из свойства 1): Ak(x ) = Ct(x )A" (x )-A'k2 (x),
а корни X, могут быть эффективно посчитаны методом бисекции.
Выражения (12), (13) могут быть рассмотрены как следующий вычислительный процесс: для каждого корня X, считаем а, и находим (C -X,I) 1 a,f методом прогонки. Если числитель на некотором корне X,
оказывается равным нулю, то X, = 0 и соответствующее слагаемое в (12) обнуляется.
В ряде случаев корни полинома известны заранее, например для первой и второй краевой задачи, что упрощает алгоритм.
3. Первая краевая задача для эллиптического уравнения
Матрица СЛАУ для первой краевой задачи имеет следующий вид: А = blocktridiag {-7, С, - 7}.
В этом случае легко проверить, что Ар = и + 1 (С /2), где и — полином Чебышева второго рода. Тогда уравнение, содержащее неизвестные иг, и, Иу, будет иметь вид
1 . И г-1 1 г(р)
--1! Л--у г 1_1!--1! —
u--— u = Г'. (14)
U , ' U , ,U c U 1 r Jc y '
c-l-1 c-l-1 r-c-1 r-c-1
При этом правая часть при исключении произвольной неизвестной преобразуется по формулам
/l(P+1) = + Ur^ /^ , /r(P+1) = /( + U--1 /c(P). (15)
Ur-l-1 Ur-l-1
Для пересчета правой части требуется вычислять Uk(C/2)/Un(C/2) / k < n, а для применения формулы (14) — Uk (C/2) Ui (C/2)/Un(C/2) /. Корни полинома Un(x) известны: Хв = cos (вк / (n +1)). Обозначим
6s,n =кв / (n +1). Значение производной U'n (x) на этих корнях известно:
Un (Хв ) = (-1 Г1 (n +1 )/sin2 (0в,n ).
Известно и значение полинома
U k (\): Uk (Х ) = sin ((k + 1)9в,„)/sin (9в,„).
Таким образом, согласно выражениям (12), (13) можно записать коэффициенты разложения в явном виде:
«в = (- 1)в-1 sin ((k +1) 6в ,„) sin (0в ,„) / (n +1) .
Так же на элементарные дроби раскладывается отношение
Uk (x / 2)Ul (x / 2)(Un (x /2))1, при этом коэффициенты разложения
«в = (-1 )в-1 sin ((k + 1)0в) sin (( + 1)0в) / (n + 1) .
Отметим, что если (k + 1)в mod (n + 1) = 0 или (l + 1)в mod (n + 1) = 0, то в разложении соответствующая элементарная дробь отсутствует.
Алгоритм прямого хода можно построить следующим образом: сначала исключим неизвестные с нечетными номерами, затем с номерами, кратными двум и не кратными четырем, затем не кратными восьми и т. д., пока не останется одно уравнение с одной неизвестной. Далее обратный ход метода реализуется согласно формуле (14).
Алгоритм 1. В массиве P задается правая часть СЛАу. В результате работы алгоритма полученное решение записывается в этот же массив. 1. Прямой ход. for (k = 0; k < [log2(N)j; k+=1 ) { for (i = 2k; i < N; i = i + 2k + 1) { l := i - 2k; r := min(i + 2k, N + 1); if (l > 1) P[l] := P[l] + UoverU(r - i - 1, r - l - 1, P[i]); if (r < N) P[r] := P[r] + UoverU(i - l - 1, r - l - 1, P[i]);
}
}
47
48
2. Обратный ход.
for (k = [log2(N)]; k > 0; k-=1) { for (i = 2k; i < N; i = i + 2k + 1) { l := i - 2k; r := min(i + 2k, N + 1); P[i] := UUoverU(i - l - 1, r - i - 1, r - l - 1, P[i]); if (l > 1) P[i] := p[i] + UoverU(r - i - 1, r - l - 1, P[l]); if (r < N) P[i] := P[i] + UoverU(i - l - 1, r - l - 1, P[r]);
}
}
function UoverU(k, n, f) R := 0;
for (s = 1; s < n; s += 1) {
if (k + 1)s mod (n + 1) == 0 continue;
as = 2( -1)s-1 sin((k + 1)%s /(n +1)) sin(rcs /(n +1)) /(n +1);
R := R + Sweeping(C - 2cos(ns /(n +1)) • I, asf);
}
return R; function UUoverU(k, l, n, f) R := 0;
for (s = 1; s < n; s += 1) {
if ((k + 1)s mod (n + 1) == 0) or ((l + 1)s mod (n + 1) == 0) continue;
as = 2 (-1)"-1 sin((k + 1)rcs /(n + 1))sin((l + 1)rcs /(n +1)) /(n +1);
R := R + Sweeping(C - 2cos(rcs /(n +1)) • I, asf);
}
return R;
Здесь Sweeping(M, f) — процедура решения СЛАУ с трехдиагональ-ной матрицей M и правой частью f методом прогонки.
4. Вторая краевая задача д^я эллиптического уравнения
Матрица СЛАУ для второй краевой задачи для эллиптического уравнения с постоянными коэффициентами имеет вид
(C/2 -I 0 ••• 0 л -I с -I ■■. ; а = 0 ■•. ■ ■ 0
: ■ ■ C -I 0 • 0 -I с/2, Отметим, что Arq = U + г (C / 2) для q > 1 и p < N.
В силу одного из свойств, связывающих полиномы Чебышева первого и второго рода, А1 = C /2 • Up-! (C /2)- Up-2 (C /2) = Tp (C /2), где
T — полином Чебышева первого рода. Аналогично AN = TN- +1 (C /2).
В результате разложения Af по первой и последней строке определителя получим
АN = С / 2АN - AN = C /2 (C /2A2n-1 - AN-2)- C / 2АN-1 + AN-2 = = C2 /4 Un-2 - (C UN-3 - UN-4 ) = (C2 /4 -1) U-2 (C /2).
Таким образом, все Ар представляют собой полиномы Чебышева первого и второго рода от одной матрицы.
Построим процесс следующим образом: будем исключать только «внутренние» неизвестные, то есть на первом шаге имеющие номера 2, 4, 6, 8, 10, 12, ..., на втором — 3, 7, 11 и т. д. При таком порядке исключения на прямом ходу необходимо вычислять только Ик/Ип f. В результате останется два уравнения, содержащие неизвестные иг и им
T
а?
N 1 М1 тт UN = /w
U 1 U
N-2 N-2
ui +т TN1 un = fN' .
И 1 и
^-'N-2 N-2
Произведем еще одно исключение любой из оставшихся неизвестных, например и^ что повлечет пересчет
((р+1 _ ((Р) , 1 ((Р)
11 11 Тт ] N '
^-1
TN -1_/('+1)
Из единственного уравнения с одной неизвестной получим u1
T
U = 7-lNЛ
1 (С2/4 -1)
и 1
Затем вычислим uN =—/N' +--u
TN -1 TN-1
Далее аналогично обратному ходу в первой краевой задаче по формуле (14) вычисляем uc, используя найденные значения u, ur и • Таким образом, кроме отношений Uk/Un и UjcU;/Un необходимо по одному разу найти Ut/T„, Tk/ ((х2 / 4-1) Un) и два раза 1/T„. Корни полинома T„ известны: xs = cos ((2s -1) ж / (2n)).
Обозначим isn = (2s - 1)ж /(2n). Выпишем значения числителя и производной знаменателя на этих корнях:
Uk (xs ) = sin ((k + 1) ^ls,n ) /sin (лs,n ) ,
Tn (xs )= n sin (nls,n ) / sin (l,n ) = (-1)S-1 n / sin (ls,n ) •
Для отношения Uk/Tn получаем коэффициенты разложения:
«s =(-l)S- 1sin ((k + l)ls,n ) / n •
Обратимся к проблеме разложения на элементарные дроби выражения
Tn+! (x/2)/((x2 /4-1)Un (x/2)).
Производная знаменателя на его корнях принимает следующие значения:
((х2 /4-1)Un)'=fx2/4"1)Un (х/2), х = 21-, U ' 1 х/2Un (х/2), х = ±2,
где = cos (rcs / (n +1)) — корни полинома ^(х). Значение U'n (Xs) приводилось ранее. В силу одного из свойств полиномов Чебышева
49
50
Tk (cos8) = cos(k8), имеем: Tn+1 (Xs ) = cos (ns) = (-1)".
Итак, получаем следующие коэффициенты разложения на элементарные дроби:
|2/( n +1) ,x = 2Xs, 1/(n +1), x = +2.
Алгоритм 2.
1. Прямой ход.
for (k = 0; k < [log2(N)j; k+=1 ) { for (i = 2k + 1; i < N; i = i + 2k+1) { l := i - 2k; r: = min(i + 2k, N); P[l] := P[l] + UoverU(r - i - 1, r - l - 1, P[i]); P[r] := P[r] + UoverU(i - l - 1, r - l - 1, P[i]);
}
}
P[1] := P[1] + UoverT(0, N - 1, P[N]);
2. Обратный ход.
P[1] := SolveFirstLine(N - 2, P[1]);
P[N] := UoverT(N - 2, N - 1, P[N]) + UoverT(0, N - 1, P[1]); for (k = [log2(N)]; k >= 0; k -= 1) { for (i = 2k + 1; i < N; i = i + 2k+1) { l := i - 2k; r: = min(i + 2k, N); P[i] := UUoverU(r - i - 1, i - l - 1, r - l - 1, P[i]); P[i] := P[i] + UoverU(r - i - 1, r - l - 1, P[l]); P[i] := P[i] + UoverU(i - l - 1, r - l - 1, P[r]);
}
}
function UoverT(k, n, f) R := 0;
for (s = 1; s < n; s += 1) {
=(2s -1)71 / (2n); as = 2 (-1)s1 sin ((k + 1)л,) / n;
R := R + Sweeping (C - 2 cos ) • I, asf);
}
return R; function SolveFirstLine(n, f) R := 0; a = 2/(n +1); for (s = 1; s < n; s += 1) {
= sn / (n +1); R := R + Sweeping (C - 2 cos ) • I, af);
}
a = 1/( n +1);
R := R + Sweeping (C + 21, af) + Sweeping (C - 2I, af); return R;
5. Третья краевая задача для эллиптического уравнения
Граничные условия третьего рода приводят к следующей СЛАУ: (С + 2а1) ч1 - 2и2 = /1,
-Ч-1 + СЧ - Ч+1 = , -1 +(С + 2Р/) = /ы.
Чтобы привести матрицу данной системы к виду (1), можно разделить первое и последнее уравнения на 2. При ц > 1, р < N полиномы Ар по-прежнему равны И +1(С / 2).
Для решения третьей краевой задачи можно использовать алгоритм, аналогичный тому, что изложен для второй краевой задачи, однако необходимо вычислить корни полиномов
А1 (х), АN (х), АN (х), например, методом бисекции.
Заключение
Метод исключения, исследованный в данной работе, позволяет достаточно легко реализовать алгоритмах для решения краевых задач для уравнения с постоянными коэффициентами. Алгоритмы являются устойчивыми, когда матрица С - 27 положительно определена. Исключение неизвестных возможно в любом удобном порядке по единообразным формулам (10), (11).
Если количество блоков системы равно 2п - 1, то рассмотренный алгоритм совпадает с алгоритмом, который описан в [3, с. 136].
Изложенный метод допускает обобщение на случай, когда матрица СЛАУ имеет вид А = blocktridiag {-Р1!,С1,-у,Ь}.
51
Список литературы
1. Buzbee B. L., Golub G. H., Nielson C. W. On direct methods for solving Poisson's equation // SIAM Journal of Numerical Analysis. 1970. № 7. P. 627-656.
2. Sweet R. A Cyclic Reduction Algorithm for Solving Block Tridiagona! Systems of Arbitrary Dimensions // SIAM J. Number. Anal. 1977. Vol. 14, № 4. P. 706-720.
3. Самарский А. А, Николаев Е. С. Методы решения сеточных уравнений. М., 1978.
4. Bank R. E., Rose D. J. Marching Algorithms For Elliptic Boundary Value Problems. The Constant Coefficient Case // SIAM J. Number. Anal. 1977. Vol. 14, № 5.
Об авторах
Александр Олегович Синюхин — асп., Балтийский федеральный университет им. И. Канта, Калининград. E-mail: [email protected]
Алексей Алексеевич Буздин — канд. физ.-мат. наук, доц., Балтийский федеральный университет им. И. Канта, Калининград. E-mail: [email protected]
About the authors
Alexander Sinyukhin — PhD student, I. Kant Baltic Federal University, Kaliningrad. E-mail: [email protected]
Dr Alexey Buzdin — Ass. Prof., I. Kant Baltic Federal University, Kaliningrad. E-mail: [email protected]
52
