Редукция неивестных при решении систем линейных алгебраических уравнений с помощью прогонки Яненко

В статье пойдет речь об организации параллельных (одновременных) вычислений при решении линейных алгебраических уравнений, возникающих при конечно-разностной аппроксимации задач для дифференциальных уравнений. В 1978 году в работе [1] был предложен способ обращения трехдиагональной матрицы, при котором оказалось возможным организовать параллельные (одновременные)

вычисления. Основная идея способа состояла в выделении части искомых неизвестных в качестве параметрических (обозначаемых далее через { а к }), формировании системы уравнений относительно { а к }, называемой далее « модулем обмена », определении параметрических неизвестных { а к } и нахождении по ним остальных искомых неизвестных. При численной реализации такого подхода оказывается возможным организовать параллельные (одновременные) вычисления. Для построения модуля обмена необходимо найти зависимость искомых неизвестных от параметров { а к }. Найденные зависимости будем называть « предрешением задачи », так как знание их явного вида и известные значения { а к } решают задачу. Узлы сетки, в которых значения искомых неизвестных принимаются за параметрические, будем называть далее « узлами распараллеливания ». За более подробным описанием этого подхода можно обратиться к [2]. В среде вычислительных математиков, знакомых с этим подходом, за ним закрепился термин «параметрическая параллельная прогонка Яненко» по имени одного из создателей этого приема вычислений, [3]. Ниже используем этот подход для нахождения решений систем линейных алгебраических уравнений, возникающих при конечно-разностой аппроксимации задач для дифференциальных уравнений. Отмечается определенная аналогия с циклической редукцией или четным-нечетным исключением искомых неизвестных, [4].

Реализация параметрического способа организации параллельных (одновременных) вычислений на примере решения системы линейных алгебраических уравнений с двухдиагональной матрицей

Рассмотрим систему линейных алгебраических уравнений с двухдиагональной матрицей:

p j + i x j + i ^- q j x j = j ^,0 ^ ^j ^ L ^{- i};

P 0 x 0 = f н

Решение (1) может быть вычислено последовательно по формулам:

x о = f ) / p o ;

x j + i = q j x j / P j + i ⁺ f j + i / P j + i ^,0 ^ ^j ^ L ^{- 1}.

Будем считать, что формулы (2) реализуют устойчивые вычисления, то есть относительно q _; / p _/+1,0 <  j <  L — i потребуем, чтобы абсолютные значения этих величин не превышали 1.

К системам линейных алгебраических уравнений с двухдиагональной матрицей сводятся конечноразностные аппроксимации смешанных задач для систем гиперболических дифференциальных уравнений, см. например, [5].

В соответствии с параметрическим способом организации параллельных (одновременных) вычислений, сеточный интервал изменения индекса (0, L ), разобьем M точками на M + 1 подинтервал вспомогательного разбиения, каждый из которых содержит N + 1 узлов. Значения искомых функций в точках M примем за параметры и будем обозначать далее через { а к } , к = 0 , „., M = 1 . Подразбиение исходного сеточного интервала порождает 2 способа индексации сеточных функций — с помощью глобального индекса и с помощью указания номера подинтервала вспомогательного разбиения, в котором находится соответствующий узел. Номер подинтервала вспомогательного разбиения при втором способе индексации узла и соответственно значения сеточной функции в этом узле будем указывать его добавлением в квадратных скобках:

gi = g,[k], j = 0,1, 2.....L i = 0, 1, 2, …, N + 1;

k = 1, …, M , где индексы j , i и к связаны естественным образом.

На каждом подинтервале [ к • ( N + 1), ( к + 1) • ( N + 1)] вспомогательного разбиения решим 2 вспомогательные задачи:

Р + 1 [ k ] u %[ k ] - q^ k ]и ⁽¹\ [ к ] = 0, i = 0,1...

Р 0 ^{[ к ] и (1)} 0 ^{[ к ]} = ¹

и

Р + 1 [ к ] u ⁽⁰⁾ i ₊ i [ к ] - q[ к ] и ⁽⁰\ [ к ] = f ₊ i [ к ], i = 0,1... Р о [ к ]и    к ] = 0.

Отметим, что решения вспомогательных задач могут быть выполнены независимо друг от друга и это предоставляет возможность параллельных (одновременных) вычисление этого этапа. По терминологии параметрической прогонки Яненко этот этап носит название - «предрешение исходной задачи». Заметим, что верхний индекс, указанный в круглых скобках, определяет граничное значение вспомогательного решения задач (3) и (4).

В силу линейности исходной задачи (1) на каждом подинтервале вспомогательного разбиения [ к • ( N + 1), ( к + 1) • ( N + 1)] имеет место представление:

Х -[ к ] = и ⁽¹⁾ [ к ] Х к ( n + 1) + и _ ⁽⁰⁾ [ к ], i = 0,1,...

в виде суммы решений однородного и неоднородного решений задач, связанных с (1). Подставляя (5) в оставшиеся разностные уравнения при к = 1, 2, 3, ... получаем двухточечные конечно-разностные уравнения относительно параметрических неизвестных:

Р кт ^Х [ ^к ] — q km - U^^к — ^1] x Г ^к — 1 1 = f m ⁺ Ч кт - 1 ^{u (} - 1 ^{[ к} — ^1], ^к = L" m ,

которые могут быть решены после этапа предрешения задачи.

Таким образом, решение исходной задачи может быть осуществлено в 3 этапа. После вспомогательного подразбиения расчетного интервала выполняем предрешение задачи с помощью (3), (4). Затем формируем модуль обмена (6) и вычисляем с его помощью параметрические переменные. Далее по представлению (5) вычисляем искомое решение во всех узлах расчетной области. На этапах предрешения задачи и восстановления искомого решения возможны параллельные (одновременные) вычисления. Возможно, и рекуррентная реализация способа - когда параметрические неизвестные вычисляются тоже с помощью параметрического подхода.

Отметим положительные стороны параметрического подхода.

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

При выделении параметрических неизвестных равенство числа узлов на подинтервалах вспомогательного разбиения не является обязательным требованием. Это может использоваться для балансировки загрузки вычислительных устройств. Назовем эту возможность «управляемым параллелизмом алгоритма».

Предрешение задачи с небольшим числом узлов

При реализации параметрической прогонки мы располагаем возможностью выбора числа узлов на подинтервалах вспомогательного разбиения исходной задачи. Мы можем ограничиться весьма малым числом узлов - 2, 3 и т. п. на этапе предрешения задачи. В этом случае предрешение задачи может быть выполнено аналитически без привлечения компьютера. Рассмотрим случай, когда подинтервал вспомогательного разбиения состоит из 3 узлов - двух граничных и одного внутреннего. Очевидно, в этом случае модуль обмена будет содержать параметры - искомые неизвестные только с четными индек- сами узлов сетки. Сформируем модуль обмена – разностные уравнения для параметрических переменных – для случая редукции неизвестных в узлах сетки с нечетными индексами. Для упрощения записи и вычислений исходную разностную задачу нам будет удобно записывать в следующей форме:

x0 задано, x1 - a0x0 = f0,                                                    (7)

Xt - (Ц-iXi-i = fi-1, так что индекс первого слагаемого – значения искомой сеточной функции будет совпадать с номером уравнения. В качестве параметрических неизвестных возьмем неизвестные с четными номерами индекса:

^x 0 ,^x 2 ,^x 4 X 2_k

Для описания этапа предрешения задачи введем для сеточной функции u t [ к ] локальные индексы i и k , где i = 0, 1, 2 – локальный индекс на подинтервале вспомогательного разбиения и k = 1, 2, … - номер подинтервала вспомогательного разбиения. На каждом подинтервале вспомогательного разбиения рассмотрим две вспомогательные задачи:

«00)[к] = 0,

U10)[/C] - Й2(к-1)и00)[^] = /2(k-1) и u01) [fc] = 1,

u11) [fc] - d2(k-i)u01)[k] = 0, где верхние индексы (0) и (1) определяют краевое условие в начальном узле подинтервала вспомогательного разбиения. Очевидно:

U10)[fc] = f2(k-1),

U 11^{) [}fc^] = d 2(k-1) .

Поэтому на подинтервале вспомогательного разбиения имеет место представление:

Xt+1 = U(0)[fc] + U(1)[fc] • Xi-1.

Подставляя это представление в разностное уравнение с индексом 2 к , получаем модуль обмена:

^X 2k - ^ 2k-1 ^ 2k-2 ^X 2k-2 = / 2k-1 ^{+ a}2k-1 f 2k-2 .                      ⁽⁸⁾

Полученный модуль обмена для нечетной редукции совпадает с видом уравнения циклической редукции, полученной элементарным приемом исключения неизвестных исходной задачи, сначала нечетных, потом с индексами, некратными 4, затем 8 и т.д. Для иллюстрации этого умножим разностное уравнение исходной задачи с нечетным индексом i на a i и сложим его с уравнением четного индекса i + 1. Получим уравнение, совпадающее с предыдущим. Однако для реализации алгоритма полной циклической редукции требуется, чтобы общее число неизвестных исходной задачи равнялось степени 2. Параметрический подход не ограничен таким условием и, более того, возможно различное число исключаемых неизвестных на этапе предрешения задачи. Действительно, возьмем в качестве параметрических неизвестные исходной задачи в узлах сетки с индексами, кратными 3. Это будет означать, что на этапе предрешения задачи исключаются неизвестные в узлах с индексами 1, 2, 4, 5, 7, 8 и т.п. В этом случае подинтервал вспомогательного разбиения будет содержать 2 внутренних узла, а весь подинтервал вспомогательного разбиения составлять 4 сеточных узла. Предрешение задачи в этом случае заключается в решении вспомогательных задач:

u(0)[fc] = 0;

г^И — l3(k-1)l°°)IAI - f3(k-1);

^l 2°^{) [}k] ^{— l} 3(k-1)+1 ^l 1°^{) [}fc] ^- f 3(k-i)+i ;

то есть l2°)[^] - f3(k-1)+1 + l3(k-1)f3(k-1)

и

u°1)[k] - 1;

u11)[/c] - a3(k-1)U°1)[k] - 0;

u21)[/c] - d3(k-1)+1u11)[k] - 0;

то есть:

^l 21' И ^{- l} 3(k-1)+1 ^l 3(k-1) -

Воспользовавшись представлением:

^3k - u^M + u21)[k] • %3k-3

получаем разностные уравнения, формирующие модуль обмена:

^% 3k ^{— l} 3k-1 ^l 3k-2 ^l 3k-3 ^% 3k-3 ^{- f} 3k-1 ^{+ l} 3k-1 ^f 3k-2 ^{+ l} 3k-1 ^l 3k-2 ^f 3k-3 .

Здесь в качестве индексов параметрических неизвестных взяты их индексы исходного разбиения. Аналогично может быть получены разностные уравнения, формирующие модуль обмена для случая параметрических переменных в узлах сетки с индексами, значения которых кратно 4:

% 4k ^{— l} 4k-1 ^l 4k-2 ^l 4k-3 ^l 4k-4 ^ 4k-4 ^- / 4k-1 ^{+ l} 4k-1 / 4k-2 ^{+ l}4k-1^l4k-2^f4k-3 ^{+ l} 4k-1 ^l 4k-2 ^l 4k-3 / 4k-4 ^.

Замечание. Разбиение исходного интервала на подинтервалы вспомогательного разбиения, содержащие одинаковое число узлов, в параметрическом подходе не является обязательным. К узлу, в котором записывается разностное уравнение относительно параметрического неизвестного (модуль обмена), могут примыкать подинтервалы вспомогательного разбиения с разным числом узлов сетки.

Реализация параметрического способа организации параллельных (одновременных) вычислений на примере решения системы линейных алгебраических уравнений с трехдиагональной матрицей

Рассмотрим в качестве исходной задачу отыскания решения системы линейных алгебраических уравнений с трехдиагональной матрицей. Например, к подобной задаче сводится отыскание решения конечно-разностной аппроксимации краевой задачи для уравнения второго порядка:

Л u = ( au- ) - bu = f ( x ), u ( 0 ) = ^ , u ( 1 ) = ^ .                   (9)

Разностные уравнения будем записывать в виде системы линейных алгебраических уравнений:

i° - ц°;

a2u2 — d1u1 + a1u° - fj,;

l3l3 — ^2l2 + l2l1 - f2;

^aL^UL - d L—l ^ L—l + ^aL-1^UL-1 = ^f L-1^;

U _L = H l , полученным стандартным образом из (9).

Обозначим через Ʌ[ k ] сужение разностного оператора на подинтервал вспомогательного разбиения под номером k . Это сужение по глобальному сеточному оператору Ʌ определяется только во внутренних точках подинтервала вспомогательного разбиения при определенных заданных краевых значениях сеточных функций. На каждом подинтервале вспомогательного разбиения выполним предрешение задачи, то есть найдем решения вспомогательных задач:

Л[к]и1°[к] = 0,u°°[k] = 1,u^°k] = 0;

Л[к]и°1[к] = 0,u°4k] = 0,u°°1k] = 1;

Л[к]н°°[к] = fi[k],u°°[k] = 0,u^°k] = 0, здесь узлы подинтервала вспомогательного разбиения с индексами, равными 0 и ^[к] соответствуют граничным узлам подинтервала вспомогательного разбиения - узлам основной сетки, где неизвестные приняты за параметрические. Для формирования модуля обмена - разностных уравнений относительно неизвестных, принятых за параметрические, воспользуемся представлением:

и[ k ] = au 10[ к ]+а и,'*[ к ]+и,l0[ к ];                            (10)

i = 0,1.....АШ

Подставляя представления (10) в разностные уравнения, соответствующие индексам параметрических переменных, получаем модуль обмена:

«кт+ХЧк + 1]«k+i — (dkm — akm+iW°°[k + 1] - akm^wM-JkDak + akmu1°[k]ak-i =  (11)

= fkm - akm+1u1°[k + 1] — «kmu^°k]-1[^].

Здесь в отличии от (8) у параметрических неизвестных взяты индексы, привязанные к номерам подинтервалов вспомогательного разбиения.

Вопрос об устойчивости прогонки для определения { а к } , к = 1 , ... M . был положительно решен в [6]. В работе была доказана устойчивость прогонки для определения { а к } , к = 1 , ... M, как следствие, устойчивости прогонки для исходной задачи.

Представляет определенный интерес рассмотреть случай, когда число узлов на подинтервале вспомогательного разбиения невелико, а именно - 1 или 2. В этом случае предрешение задачи может быть выполнено аналитически: вручную или с привлечением системы символьной алгебры. Такой крайний случай, на наш взгляд, повышает степень параллелизма алгоритма.

Для формирования модуля обмена в узле распараллеливания под номером к (этот номер не совпадает со значением индекса узла исходной сетки!) в подинтервалах вспомогательного разбиения, примыкающих к узлу необходимо выполнить предрешение задачи и воспользоваться представлением (10).

Для случая одного внутреннего узла на подинтервале вспомогательного разбиения имеем вспомогательные задачи:

u°°[k] = 1,м2°[^] = 0;

a 2k U 1^0[ kA - d 2k—i U 1° [k ^] + a 2k-i U °° = 0

или или

или

Далее

Далее

_u10 _ « ^2k-1 /

1          /a2k-1

Uq1^] = 0, и°Ч

« 2k U ⁰¹M - d 2k-1 U i^{1 [}k^] + « 2k-1 W ⁰¹ = 0

=    ₂

u°0[k] = 0,u00[k] = 0;

« 2k « 2 ^{0 [} k ^] - d 2k-1 U i⁰M + « 2k-1 ^ o⁰ = / 2^-1

<_

-

/ 2^-1 /

' ^2k-1'

Так как в случае одного внутреннего узла на подинтервале вспомогательного разбиения модуль обмена в узле под номером к имеет вид:

a2k+1U°1[k + 1]ak+1 - (d2k — «2k+1W°0[fc + 1] - «.kU0^1)»/,- + «2kU°0Mak-1

= /2k - «2k+1^°0[^ + 1] - «2ku°0M и с помощью найденных значений на этапе предрешения задачи получаем вид модуля обмена в узле под номером k:

« 2k + 1 « 2k+2 /^ 2k+1^ak + 1 ^{- (}^ 2k

-

« 2k+1     « 2k

-

^ 2k+1   ^ 2k-1

■) a k + « 2k « 2k-1 ^ k-1

= / 2k ^- « 2k + 1 /d 2k + 1 f 2k + 1

-

« 2k / d. 2k-1 f 2k-1 -

Приведем вид вспомогательных задач на этапе предрешения и вид модуля обмена для случая 2 внутренних узлов подинтервала вспомогательного разбиения. В этом случае под параметрами принимаются переменные в узлах исходной сетки со значениями кратными 3:

uo0[^] = 1,u10[^] = 0;

«3k-1U10[^] - d3k-2U°0M + «3k-2U°0M = 0;

«3kU10[fc] - d3k-1U°0[k] + «3k-1U°0[k] = 0

или

10 11,1 _ «3k-2 ^ 3k-1/

^- « ²k-1 ^);

^{U 1 [k]}              /(d 3k-2 d 3k-1

10[r,] _ « 3k-2 « 3k-1 /

²                  4d 3k-2 d 3k-1 ^- « 3k-1 )

Далее

uo1[^] _ 0, u01M _ 1;

«3k-1U21[^] - d3k-2U01M + «3k-2U°1[fc] _ 0;

«3kU01[fc] - d3k-1U°1[fc] + «3k-1U01M _ 0

или

„011/, I — ^азк ^ зк-2 /

^{№          /}(^ зк-2 ^ зк-1

^{— а 2}к-1 ^);

и 01 гц = « зк « зк-1 /

²              ' №к-2 ^ зк-1 ^— « зк-1 )

или

Далее

u00[fc] = 0, u00[fc] = 0;

азк-1и00[к] — ^зк—2^00[^] + азк-2и00[к] = /зк-2;

а зк и ⁰⁰И — ^ зк-т^М + й зк-Х⁰^ ] = / зк-1

00 гц — ^—(а 3к-1 / зк-1 + ^ зк-1 / зк-2⁾/

^U1 ^                            ^/(^ зк-2 ^ зк-1

—

« 2к

-1

);

00 гц — ^—(a3k-1 f 3k-2 + ^ зк-2. / зк-1 ⁾ /

'                                      ^/(^ зк-2 ^ зк-1

—

« зк-1 )

и с помощью найденных значений на этапе предрешения задачи получаем вид модуля обмена в узле под номером k :

^a 3fc+t^3k+2 _____^     _ ( ^__ ^а зк+1^зк+2 ___ ^а зк-1 ^а зк^зк-2    J ^  _|__ ^a 3fc ^a 3fc- 1______^     _

(^ зк+1 ^ зк+2 -а |к+2 )  ^к+1   \ ^зк   (^ зк+1 ^ зк+2 -д |к+2 )   (d 3k-i d 3k-2 -a |_k-i )/ ^к   (d 3_k-i d 3_k-2 -a 2_k_i )  ^к-1

г   ,         ^зк+2/зк+1^+азк+2/зк+2 ,          d 3k-2 / 3k_1 ^+a 3k_1 / 3k-2

7 зк ^{+ а} зк+1 r_a a „2    1 ^{+ ы} зк-1 (a a     „2    1.

⁽“ зк+1 “ зк+2 ^-Я 3к₊₂ ^{)              (}“ зк_2 “ зк_1 ^-Я 3к_₁ ⁾

Случай исключения каждого 4 неизвестного исходной задачи достаточно громоздкий, но тоже поддается «ручному» рассмотрению.

Заключение

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