CC BY
28
Вычислительные технологии
Том 10, № 5, 2005
СХОДИМОСТЬ ПАРАЛЛЕЛЬНОГО АЛГОРИТМА ДЛЯ КОМПАКТНЫХ СХЕМ В НЕОДНОРОДНЫХ ОБЛАСТЯХ*
В. И. ПААСОНЕН
Институт вычислительных технологий СО РАН, Новосибирск, Россия
e-mail: paas@ict.nsc.ru
In the article, we studies the convergence of the parallel algorithm with high-order boundary conditions for boundary value problems in inhomogeneous areas. We studies the general problem and symmetric problem without convective terms.
Введение
Известно, что многие дифференциальные уравнения в частных производных второго порядка, в том числе с переменными коэффициентами или нелинейные, в простых и сложных координатных системах могут быть аппроксимированы с четвертым (иногда с шестым) порядком относительно шагов пространственной сетки на компактном шаблоне, не выступающем за пределы 3 х 3 х ... х 3 — точечного шаблона (см. работы [1, 2] и библиографию в них). Для краевых задач, отличных от простой задачи Дирихле, необходимо дополнительно решить проблему аппроксимации с адекватной точностью краевых условий, в том числе внутренних, на границах раздела сред. Это можно сделать различными способами, например, аппроксимируя закон сохранения в балансной ячейке в окрестности границы [3, 4] или приближая граничное условие без использования самого дифференциального уравнения [5, 6].
Автор придерживается мнения, что с различных точек зрения (независимости от характера уравнений, универсальности формул для любого порядка точности, пригодности для любых типов граничных условий) для достижения адекватной точности в целом наиболее удобно использовать односторонние разностные аналоги первых производных [7, 8] в граничных условиях. Необходимая точность при этом может быть достигнута увеличением числа точек в односторонней аппроксимации потока. Очевидно, этот подход годится и для внешних граничных условий второго-третьего рода, и для условий равенства потоков на границах раздела различных сред. Такой симбиоз компактных схем и многоточечных аппроксимаций граничных условий предложен автором в работах [7, 8]. В первой из них основное внимание уделено традиционной последовательной реализации, а во второй разработана параллельная технология, представляющая собой обобщение алгоритма
* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 05-01-00146-а).
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2005.
распараллеливания прогонки [9]. Данная работа является продолжением этого цикла и посвящена теоретическому обоснованию параллельного алгоритма.
1. Система уравнений
Прежде всего представим форму односторонних многоточечных разностных аппроксимаций первых производных, которые будут использоваться в дальнейшем при записи разностных граничных (внешних и внутренних) условий. Аппроксимируем первую производную на равномерной сетке с шагом к примененая одностороннюю разделенную разность на произвольном множестве равноотстоящих точек
А* 1 V*
Т = к ^ ак Т* •
к=0
Здесь Т, — оператор сдвига (Т^и(х) = и(х + к)), а ак — коэффициенты в + 1-точечного разностного оператора. Если они удовлетворяют линейной алгебраической системе уравнений
У^уакк1 = 511, I = 0, •••, в, к=0
(5ц — символ Кронекера), то погрешность аппроксимации первой производной на достаточно гладких функциях составляет величину О (к8), причем ввиду линейной независимости строк матрицы (Вандермонда) системы — это есть максимально возможный порядок аппроксимации при данном в. Коэффициенты многоточечной разности "вперед" могут быть вычислены по правилу Крамера в явном виде
(_1) к+1 «
ак = -—к-С, к = 0, а0 = ак,
к=1
который понадобится в дальнейшем при доказательстве корректности параллельного алгоритма. Здесь С8: — число сочетаний из в по к.
Разность "назад" порядка в имеет противоположную ориентацию и те же коэффициенты, но с обратным знаком
А-8 = -к £ акТ-. = - к £ акТ-.
к=0 к=0
Пусть х0, х1, •••, хг — возрастающая последовательность значений пространственной координаты х, представляющих собой координаты границ однородных слоев Шк = {хк-1 < х < хк} (к = 1, •••, г), внутри которых пространственный шаг кк сохраняет постоянное значение, вообще говоря, различное для различных слоев.
Предположим, что в таком многослойном пакете решается система уравнений вида
а,гЩ-1 + ЬгМг + СгЩ+1 = 1Сг, 0 < % < N % = ¡к, к = 1, •••, Г - 1, с внешними граничными условиями при % = ¡0 = 0 и % = ¡г = N
^о«о - к- азиз = Ф0, - кТ^^ азиN-з = Ф1
1 з=0 г з=0
8
и с внутренними граничными условиями равенства потоков на границах раздела слоев пакета
Л 5 Л 5
а+ ^а= ^ = 1к к = 1>->г-1-
к ]=0 к+1 ]=0
Внешние граничные условия записаны так, что их тип регулируется параметрами V и ф на левой и правой границах.
Итак, разностное уравнение является неоднородным: во внутренних узлах слоев имеем трехточечные соотношения (это могут быть, например, компактные аппроксимации повышенной точности или приближения обычной точности), а на внешних и внутренних границах использованы описанные выше многоточечные односторонние аппроксимации потоков. Мы намеренно не конкретизируем порядок аппроксимации схемы внутри слоев. Ясно, что при любом порядке выбором "длины" граничного условия в можно добиться необходимой гармонии в порядках. Исследование алгоритма будем проводить при произвольном в. При этом будем предполагать, что "длинные" аппроксимации потоков не выходят за пределы одного слоя. Этого всегда можно добиться уменьшением шага сетки внутри каждого слоя так, чтобы произведение вк не превосходило ширины слоя.
2. Параллельный алгоритм
Предположим, что значения решения на границах слоев уже известны. Обозначим их 'к, (к = 0,...,г). Введем в каждом слое Шк = {хк-1 < х < Хк} локальную нумерацию узлов (т = 0,...,М, М = тк), обозначим фрагмент вектора решения и в слое шк через ук и представим его в виде специальной линейной комбинации
Ук = Рк 'Шк-1 + Як 'к + %к,
где %к — решение задачи Дирихле для неоднородного уравнения с нулевыми граничными условиями
ат%т-1 + Ьт + Fm, %М 0
а рк и як — решения задач Дирихле для соответствующего однородного уравнения со специальными граничными условиями (1, 0) и (0,1),
атРт-1 + ЬтРт + Стрт+1 = 0 р0 = 1, Рм = 0
атЯт-1 + ЬтЯт + СтЯт+1 = 0, Яю = 0 Ям = 1-
Здесь Ет = ^к+т. Аналогично выражаются коэффициенты трехточечных уравнений при локальной нумерации в слое, например ат = щк +т.
Приведенные выше три вспомогательные краевые задачи решаются, например, методом прогонки, причем соответствующие алгоритмы реализуются взаимно независимо для каждого из г слоев пакета, что позволяет проводить вычисления параллельно. Заметим, что линейная комбинация ук решений вспомогательных задач рк, як и %к специально составлена так, что она удовлетворяет исходному неоднородному уравнению с правильными
к к граничными условиями у0 = 'к-1 и Ум = 'к.
Подставляя после решения вспомогательных задач комбинации ук во все "длинные"
соотношения исходной системы уравнений, т. е. во все внутренние и внешние граничные
условия, получим замкнутую систему уравнений для вычисления значений ' решения
на границах слоев. Система оказывается трехдиагональной, так как по предположению потоки аппроксимируются по узлам в пределах слоя. В самом деле, подставляя выражение решения у1, соответствующее первому слою, в соотношение на левой границе, получим
- — ^ а(р]'0 + <]Ш1 + г1) = фо.
Т1 • п
3=0
Очевидно, что после приведения подобных равенство превращается в линейное уравнение для двух неизвестных ш0 и ш1. Из правого граничного условия аналогично устанавливается связь для двух неизвестных — шг-1 и шг, а из соотношения на к-й границе раздела сред вытекает трехточечное соотношение, связывающее значения искомого решения Шк-1, 'к и на трех последовательно расположенных границах слоев. Это следует из того, что все значения решения, участвующие в аппроксимации потока слева, относятся к левому слою и выражаются через Шк-1 и Шк, а все слагаемые в аппроксимации потока справа выражаются через Шк и Шк+1.
Таким образом, значения Шк на границах слоев определяются из системы уравнений с трехдиагональной матрицей
Актк-1 + Бк-шк + ОкШк+1 = Ск, к = 1,..., г - 1, BоWо + С0^1 = ^0, Аг Шг-1 + Бг 'г = ^1,
где
, Ак к ?-, Ак к , Ак+1 Ак = азРы-з, Бк = т У. аз <ы-з + 7- У. азР
Тк „ Тк п Т к+1 Г-Г
з=0
з=0
к+1 з
з=0
Ск
Ак+1 к к+1
Е
з=0
аз 3+1-
°к = - 3
3=0
к А
аз гЫ-3 - к
к+1
к+1
Е
з=0
аз 3+1
Б0 = »0 - аз Р1, С
з=0
Т0 аз з & = ф0 + аз г1
з=0
з=0
Аг = - аз'РГх-з, Бг = »1 - аз-з, = ф1 + ТГ^
з=0
кг
з=0
кг
аз гм-з •
з=0
Размерность полученной системы уравнений не растет при стремлении шага сетки к нулю, оставаясь постоянной, а именно равной общему числу границ (г + 1) в слоистом пакете. Казалось бы, это обстоятельство освобождает от необходимости ее исследования. Между тем известны матрицы, даже при малой размерности плохо обусловленные. В частности, это известная матрица Гильберта или примеры, сконструированные Годуновым [10]. Уже при размерности менее десяти решения систем с такими матрицами не могут считаться надежными при обычном машинном представлении чисел. Поэтому считаем целесообразным проанализировать устойчивость последней задачи, хотя ее размерность конечна.
8
3. Устойчивость параллельного алгоритма
Будем исследовать систему с точки зрения диагонального преобладания при постоянных в пределах каждого слоя коэффициентах.
Заметим, что любой трехточечный разностный оператор (в нашем случае — действующий на верхнем слое внутри области однородности) с точностью до множителя можно представить в виде
В Е + ^ - Т- Ть - 2Е +
В = Е + аит-—--^^т
2к к2
где а и 7 суть положительные веса в аппроксимации конвективных и диссипативных членов. Поэтому коэффициенты разностного уравнения выражаются в виде
к к
а =---а, Ь =1 + 2а, с =--а,
2 ' 2
где
ит цт
К = аТ> а =
безразмерные взвешенные конвективный и диссипативный параметры. Дискриминант характеристического уравнения
положителен, поэтому его корни
Р±
cp2 + bp + а = 0
1 + 2а ± V1 + 4а + к2 2а — к
вещественны и различны. Ясно, что при соотношении шагов т ~ h конвективный параметр ограничен, тогда как диссипативный в пределе стремится к бесконечности, а при т ^^
h2
конвективный стремится к нулю при ограниченном диссипативном. Иначе говоря, в обоих случаях масштаб конвективного параметра ниже, чем диссипативного, поэтому будем считать, что | к | < 2а. Это требование вполне аналогично ограничению на разностное число Рейнольдса. Для нас важно, что оно обеспечивает положительность обоих корней p+ и p- характеристического уравнения. Заметим, что из двух корней p+ — больший корень. Кроме того, простой анализ показывает, что
p+ > 1, p- < 1.
По теореме Виета произведение корней p+p- = а/c. В частности, в симметричном случае а = c =1.
Наконец потребуем, чтобы трехточечные соотношения системы в каждом слое подчинялись условию диагонального преобладания b > |c| + |а|. Это требование гарантирует устойчивость любой задачи Дирихле для любого отдельно взятого слоя. В частности, гарантируется устойчивость всех приведенных выше вспомогательных задач Дирихле, составляющих содержание первого (параллельного) этапа вычислений с потенциально бесконечномерными матрицами.
Таким образом, при сделанных предположениях проблема исследования устойчивости параллельного алгоритма в целом сводится к изучению устойчивости сформулированной
выше конечномерной задачи, связывающей значения искомого решения на границах слоев. Коэффициенты системы , В к, С к, соответствующие условию на к-й границе, запишем (опуская далее индекс к) в виде
Л = - В = Ь- + ®+ь+, С = £) +с+.
Здесь
а- = ^ армч , ь- = ^ а дм_ч, ъ+ = ^ азр+ , с+ = X] а я+ >
где суммирование ведется от 0 до в, а знаки "—" и " +" позволяют отнести величины к слою слева и справа от данной границы. Для получения достаточного критерия устойчивости потребуем выполнения условия диагонального преобладания, которое ввиду взаимной независимости отношения \/к слева и справа сводится к совпадению знаков Ь+ и Ь_ и преобладанию модулей Ь+ и Ь_ над с+ и а_ соответственно. Заметим, что такие неравенства должны выполняться на каждой границе, поэтому можно сгруппировать их более удобным для решения способом — не по границам, а по слоям, т. е. для левой и правой границ данного слоя. При этом для внешних границ не требуется записывать особых неравенств. Для условия первого рода соответствующее неравенство снимается, для второго остается неизменным, и для условия третьего рода только усиливается. Система неравенств запишется в виде
Ь+ > с+, Ь- > а-,
где
а- = X азРм-3' Ь- = X ат-3' Ь+ = X азРз, с+ = X а3 Я3, а векторы р и д суть решения наших вспомогательных задач для текущего слоя.
Воспользуемся тем, что решения р и д вспомогательных задач Дирихле для однородного уравнения в каждом слое могут быть найдены точно в виде линейных комбинаций степеней корней р± характеристического уравнения. С учетом граничных условий имеем
= рМ р- — рМ р+ = р+ — р-
Рз пм пм , Яз пм пм' р+ р- р+ р-
где М — число шагов в слое, а 3 — локальный номер узла нем. Введем в рассмотрение полином специального вида
ЗД = X
3=0
аз г3,
где — коэффициенты многоточечной односторонней аппроксимации производной. Заметим, что данный полином может быть выражен также в компактной интегральной форме. Подставим в него явные выражения для коэффициентов . В результате получим простую зависимость
= X (—1)3+1 С1 (г3 — 1) = / ¿(—1)з+1 С33х3-Чх = / 1 — (1 — ^ ¿х. з=1 3 1 3=1 1 Х
8
Данная функция играет важную роль — через ее значения выражаются коэффициенты Ь±, с+ и а- и, следовательно, все выражения в неравенствах системы.
В самом деле, подставляя значения Рз, <з в выражения для коэффициентов, получим
= Ф3(р+) - Ф3(р-) = рЫФв(1/р+) - РыФ.(1/р-)
С+ „Ы Ы , Ь-
рЫ - рЫ ' рЫ - рЫ
= рЫФ.(р-) - рЫФ8(р+ ) = Ф,(1/р_) - Ф,(1/р+) ( )Ы
Ь+ пЫ _ пЫ , а- пЫ _ пЫ (р-р+) .
р+ - р- р+ - р-
Очевидно, что Ф8(р-) и Ф8(1/р+) безусловно отрицательны, так как подынтегральное выражение положительно, а верхний предел интегрирования меньше нижнего. Аналогично устанавливается положительность Ф8(1/р-) и Ф8(р+) при всех нечетных в. Отсюда при нечетных порядках вытекает отрицательность Ь+ и Ь-, а вместе с тем и совпадение знаков слагаемых центрального коэффициента Б разностного уравнения.
При четных порядках в числителях выражений для Ь+ и Ь- уменьшаемые отрицательны, а вычитаемые могут быть любого знака. В этом случае отрицательность Ь+ и Ь- имеет место при одновременном выполнении условий
рЫФs(1/р+) < рЫФв(1/р-), рЫФв(р-) < рыФ8(р+).
Будем предполагать, что либо в нечетно, либо для выбранного соотношения конвективного и диссипативного безразмерных параметров выполняются сформулированные условия. Множество удовлетворяющих им решений непусто, что можно определенно утверждать, учитывая отрицательность левых частей обоих неравенств и положительность Ф8(г), по меньшей мере, на промежутке (1, 2]. Справедливо следующее утверждение, дающее критерий корректности параллельного алгоритма.
Теорема 1. При сделанных предположениях выполнение неравенств
рыФ3(р+) - рЫФз(р-) > |Фв(р+) - Фв(р-)|, рыФв(1/р-) - рЫф.(1/р+) > (р-р+)Ы|Ф-(1/р-) - Ф-(1/р+)|
является критерием диагонального преобладания матрицы системы уравнений для определения значений решения на границах слоев.
Доказательство. Непосредственная подстановка выражений коэффициентов а-, Ь±, С+ в условия диагонального преобладания после упрощений и сокращений положительных множителей приводит к сформулированным неравенствам.
Рассмотрим особо симметричный случай а = С. Как отмечалось выше, по теореме Виета р+ = р, р- = 1/р, поэтому
Ь = Ь = рЫФа(1/р) - 1/рЫФ,(р) Ь+ рЫ - 1/рЫ ,
= = Фа(р) - Фа(1/р)
С+ а- рЫ - 1/рЫ .
Поскольку коэффициенты попарно совпадают, неравенства системы оказываются идентичны, каждое из них приводится к условию
р2Ы [Фз(1/р)]2 > [Ф^(р)]2.
Предположим, шаг сетки по пространству h уменьшается согласованно с временным шагом (или итерационным параметром) т так, что диссипативный безразмерный параметр a остается постоянным. Тогда из последнего неравенства можно сделать следующие выводы.
Теорема 2. Если при некотором шаге h неравенство выполняется, то оно выполняется при любом более мелком шаге.
Доказательство. Ввиду предположения о постоянстве диссипативного параметра при h ^ 0 коэффициенты системы, а следовательно, и корни характеристического уравнения, не меняются. Единственной переменной величиной в неравенстве является число шагов в слое M, растущее с уменьшением шага h. Очевидно, множитель в левой части неравенства возрастает, являясь натуральной степенью числа р = р+ > 1. Если неравенство выполняется при некотором M, то при всех больших значениях оно также выполняется.
Теорема 3. Если неравенство выполняется в частном случае M = s, то оно выполняется при любом шаге.
Доказательство. Справедливость утверждения следует из теоремы 2, если заметить, что случай M = s соответствует максимально возможному шагу в слое, когда "длинное" граничное условие записывается по всем узлам слоя вплоть до соседней границы, т. е. занимает всю сеточную область слоя.
Теорема 3 дает простой критерий для анализа корректности системы при любом шаге
|рвФв(1/р)| > |Фв(р)|.
Легко видеть, что в левой части неравенства под знаком модуля находится полином с теми же коэффициентами, как полином в правой части, но взятыми в обратной последовательности. Это обстоятельство делает неравенство полиномиальным, облегчая его анализ, в том числе с применением операционных вычислительных сред, таких как Matlab или Maple.
В заключение заметим, что в общем случае при изменении шага сетки корни характеристического уравнения не остаются постоянными, так как подбором соотношения шагов невозможно одновременно стабилизировать a и к, поэтому в критерии теоремы 1 все величины зависят от шага. Однако в пределе при h ^ 0 в силу приведенного выше рассуждения о соотношении порядков диссипативного и конвективного безразмерных параметров второй из них пренебрежимо мал, поэтому неравенства, сформулированные в теореме 1, приближаются к критерию теоремы 2.
Список литературы
[1] Самарский А.А. Схемы повышенного порядка точности для многомерного уравнения теплопроводности // Журн. вычисл. математики и мат. физики. 1963. Т. 3, № 5. С. 812-840.
[2] Паасонен В.И. Компактные схемы для систем уравнений второго порядка с конвективными членами // Вычисл. технологии. 1998. Т. 3, № 1. С. 55-66.
[3] Коскин П.И. Схема повышенной точности для уравнения теплопроводности с разрывными коэффициентами // Численные методы механики сплошной среды: Сб. науч. тр. / АН СССР. Сиб. отд-ние. ВЦ; ИТПМ. 1979. Т. 10, № 2. С. 85-96.
[4] Ильин В.П. Балансные аппроксимации повышенной точности для уравнения Пуассона // Сиб. мат. журнал. 1996. Т. 37, № 1. С. 151-169.
[5] Вллиуллин А.Н., Сафин Р.И., Паасонен В.И. О схеме расщепления с повышенным порядком аппроксимации краевых задач для уравнения Пуассона // Численные методы механики сплошной среды: Сб. науч. тр. / АН СССР. Сиб. отд-ние. ВЦ; ИТПМ. 3, № 1, 1972, 17-25.
[6] Паасонен В.И. Разностные схемы высокого порядка точности для краевых задач в неоднородных областях // Тр. Междунар. конф. по вычисл. математике. Ч. II. Новосибирск, 2004. C. 574-579.
[7] V.I. Paasonen. Compact difference schemes for inhomogeneous boundary value problems // Russ. J. Numer. Analys. and Math. Modell. 2004. Vol. 19, N 1. P. 65-81.
[8] Паасонен В.И. Параллельный алгоритм для компактных схем в неоднородных областях // Вычисл. технологии. 2003. Т 8, № 3. С. 98-106.
[9] Яненко Н.Н., Коновалов А.Н., Бугров А.Н., Шустов Г.В. Об организации параллельных вычислений и распараллеливании прогонки // Численные методы механики сплошной среды: Сб. науч. тр. / АН СССР. Сиб. отд-ние. ВЦ; ИТПМ. 1978. Т. 9, № 7. С. 136-139.
[10] Godünov C.K. Лекции по современным аспектам линейной алгебры. Унив. сер. Т. 12. Новосибирск: Науч. книга, 2002. 202 с.
Поступила в редакцию 30 августа 2005 г.
