CC BY
10УДК 517.518.87
Решение интегральных уравнений на многопроцессорных вычислительных системах
М.Д. Рамазанов*, Д.Я. Рахматуллин, Л.З. Валеева, Е.Л. Банникова
Институт математики с ВЦ УНЦ РАН 450077 Россия, Уфа, Чернышевского, 1121
Получена 24.06.2008, окончательный вариант 28.07.2008, принята к печати 10.03.2009 Рассматриваются уравнения Фредгольма по областям различных форм. Приводится описание алгоритма приближенного вычисления интегралов с использованием многопроцессорной техники.
Ключевые слова: многопроцессорные системы, приближенное вычисление интегралов, многопроцессорная техника.
1. Введение
Мы рассматриваем уравнения Фредгольма второго рода
u(x) - J dyK(x,y)u(y) = f (x), (1)
Q
записанные для функций двух переменных в ограниченной области произвольной формы,
О С R2. Граница области Г = дО, функции K(x,y), f (x) и решение предполагаются глад-
кими, М раз непрерывно дифференцируемыми по своим аргументам. Норма интегрального оператора в пространстве C(О) считается малой, чтобы можно было применить метод последовательных приближений, max J dy\K (x,y)\ = в < 1. Тогда существует единственное
x£Q Q
решение - непрерывная функция, u G C(О), и
u(x) = lim us(x), Us(x) = f(x)+ dyK(x,y)us-i(y), u0 = f(x). (2)
s^<x, J
Из самого вида интегрального уравнения (1) следует гладкость этого решения, u G CM(О.). А итерационное равенство (2) позволяет обосновать равномерную ограниченность CM(О) —норм последовательности {us(x)\s = 1, 2...}.
KIIcm < ||us-i Ус • max f dy V D K (x,y)\ + \\f ||cm < K.
x£Q J i^i
Q \a<m\
Это дает возможность в вычислениях интегрального члена итерационного равенства (2) использовать кубатурные формулы с равномерными оценками в CM(О) — пространствах. Мы применим решетчатые кубатурные формулы с ограниченным пограничным слоем порядка М (ОПС-формулы). При наложенных ограничениях они являются асимптотически
* Corresponding author E-mail address: ramazanov@yandex.ru
1 © Siberian Federal University. All rights reserved
оптимальными в нормах пространств Шрт и Ст с любым т < М (в конкретной, одной из эквивалентных, нормировке и, конечно, с т > п/р). Порядок приближения интегрального слагаемого этими формулами - 0(Нт) с соответствующей примененному пространству наилучшей константой при этом порядке.
Отметим также, что эти формулы условно ненасыщаемы для всех т < М, то есть формулы остаются асимптотически оптимальными при любой промежуточной гладкости т € (0, М), и мы можем работать с ними для К(х,у) и ](х) из классов Ст, получая соответственно решения этой же гладкости.
В общем случае область произвольной формы задается равенством П = {х|Ф(х) > 0}. Мы предполагаем, что Ф € См(М2), множество {х|Ф(х) ^ 0} является ограниченным, а ВФ(х) = 0 при Ф(х) = 0, что обеспечивает гладкость границы. Для удобства программной реализации мы накладываем дополнительные условия |ДФ(х)| = 1 при Ф(х) = 0, достигая его заменой первоначального Ф(х) на Ф(х)/^Ф(х)У
Решетку узлов кубатурной формулы возьмем кубической {НЩк € Ъп, Н € М+, Н ^ 0}. Интегральное слагаемое итерационного процесса (2) заменяется кубатурной формулой
Нп ек (Н)К (Н],Нк)ив(Нк)
р(Кк,0.)КЬК
с некоторым Ь ^ 4М.
ОПС-свойство этой формулы означает, что для р(Нк, М2 \ П) ^ ЬН коэффициенты равны единице, ек(Н) = 1. Заметим, что если / € См(П), К € См(П х П), то все коэффициенты формулы можно взять равными единице, при этом пограничный слой будет отсутствовать. Тогда и решение получится из класса См(П). Мы применяем этот вариант в вычислительных экспериментах при оценке объема работы по вспомогательным расчетам.
Выбранный нами алгоритм последовательных приближений позволяет полностью распараллелить расчеты на каждой отдельной итерации. Поэтому мы можем применить многопроцессорные вычислительные системы. Конкретно расчеты проводились на МВС-15000ВМ и МВС-50К Межведомственного Суперкомпьютерного Центра РАН.
Идея применения решетчатых кубатурных формул для приближения интегралов по областям произвольных форм заключается в локализации вычислений при помощи разбиения области на более простые и преобразовании границ этих локальных областей в плоские - координатные гиперплоскости. Алгоритм локализации (предложенный А.Н. Игнатьевым [1]) и реализующая алгоритм программа автоматического разбиения области описаны в п. 3. Эта программа составлена Л.З. Валеевой [2]. А программа для вычисления интегрального слагаемого составлена Е.Л. Банниковой [3], ее описание и анализ работы приведены в п. 4. Алгоритм приближенного вычисления интегралов дан Д.Я. Рахматуллиным [4], описание приведено в п. 2. Следует отметить, что в описании алгоритма Д.Я. Рахматуллина наложено условие выпуклости области интегрирования. Однако оно было использовано им только для упрощения разбиения области на подобласти более простых форм. Эта проблема локализации решается в общем виде в п. 3 без требования выпуклости области П = {х|Ф(х) > 0}, что и применено в п. 4 в алгоритме и программе численного решения интегральных уравнений. В основе решетчатых кубатурных формул лежит алгоритм формул С.Л. Соболева
- кубатурных формул с регулярным пограничным слоем (РПС-формул) [5] и алгоритм построения ОПС-формул для локальных областей М.Д. Рамазанова [6].
Метод последовательных приближений разбивает вычисление решения интегрального уравнения на отдельные временные такты. Это мешает полному распараллеливанию алгоритма (естественному для вычисления интеграла по ОПС-формулам). Впоследствии мы
намереваемся преодолеть эту трудность записью приближенных решений рядами Неймана с одновременными расчетами членами этих рядов, обеспечивающими заданную точность приближения решения уравнения (1).
В п. 5 мы подводим итоги исследований по результатам вычислительных экспериментов и даем общие рекомендации по применению разработанных программ.
2. Алгоритм и программа приближенного вычисления интегралов в п-мерном случае
2.1. Постановка задачи
Пусть задана ограниченная область П С Мп с гладкой границей Г € Ст.
Рассмотрим пространство Соболева ^), р > 1, т > п, периодических функций, интегрируемых с р-ой степенью вместе с производными до т-го порядка включительно, в одной из эквивалентных норм:
Як ■= д(х)е 2пгхк(1х,
где Q — единичный гиперкуб: Q := {х : х^ € [0; 1), г = 1,п}. Взяв это пространство за основу, построим новое пространство Ш^(П) с нормой
Требуется вычислить интеграл произвольной функции ] € Ш ^(П) по области П. Мы решаем эту задачу путем приближения интеграла решетчатыми кубатурными формулами с ограниченным пограничным слоем (ОПС-формулами). Поясним эти термины.
Решетчатыми кубатурными формулами называются кубатурные формулы, узлы которых совпадают с узлами некоторой решетки. Мы будем рассматривать последовательности формул, соответствующие последовательностям кубических решеток со сгущающимися узлами {Нк}ке2п, Н ^ 0 :
КкЕ4
Здесь за знак суммы вынесен масштабный множитель Нп, равный объему элементарной (наименьшей) ячейки, образуемой узлами решетки.
Кубатурная формула обладает ограниченным пограничным слоем, если ее коэффициенты равномерно ограничены по к и Н :
К£(1 ) = Нп^ ек(Н)1 (Нк), Н ^ 0.
3 Ь1 : эир Ук (Н) < Ь1,
к, К
а в глубине области все коэффициенты равны единице:
3 Ь2 : УН, к : р(Нк, Мп\П)^Ь2Н ^ ек(Н) = 1.
Для решения данной задачи создан алгоритм, позволяющий вычислить коэффициенты кубатурной формулы, приближающей точное значение интеграла с погрешностью порядка Нт.
2.2. Алгоритм
Зададимся произвольным целым числом М, таким что М > п. Ниже мы приведем алгоритм (см. [4-9]) вычисления кубатурной формулы, оптимальной по порядку на каждом из пространств Ш“(П) с т = , М^ , т.е. универсально оптимальной по порядку на этом
классе пространств.
Учитывая оценки
3С1 : >С1 \\и\{^Ш* >С1Ьт
для оптимальности по порядку достаточно иметь оценки сверху:
3 С2 : |Ю|(^(П))* <С2Ьт, Н ^ 0.
По следующей теореме универсально оптимальные по порядку кубатурные формулы с ОПС являются универсально асимптотически оптимальными.
Теорема (М. Д. Рамазанов). Кубатурная формула с ОПС асимптотически оптимальна на каждом пространстве из множества
{«Г™} т € (т.т),
р € ( р1 , р2 )
тогда и только тогда, когда она оптим,альна на каждом из них по порядку, при п < т1 < т2 < М, 1 < р1 <р2 < ж.
В технических целях наложим дополнительные ограничения на условия задачи. Пусть П — выпуклая область, лежащая внутри единичного гиперкуба вместе с замыканием: П С С Ш Q, а параметр Н стремится к нулю по одной из последовательностей, каждый член которой можно представить в виде N-1, где N € .
Алгоритм основан на сведении задачи интегрирования по области с гладкой границей к нескольким однотипным задачам интегрирования по единичному гиперкубу. Для нашей области всегда можно подобрать специальное разбиение единицы — конечный набор функций (будем называть их срезывающими), удовлетворяющий условиям:
1) они являются периодическими с основным периодом Q ив сумме составляют функцию, тождественно равную единице в фиксированной окрестности П области П, целиком лежащей в гиперкубе Q (периодически продолженной с основным периодом Q);
2) они принадлежат пространству См^);
3) носитель одной из срезывающих функций лежит полностью внутри области П, не пересекаясь с ее границей; пересечение носителя каждой из остальных срезывающих функций с границей д П непусто и проектируется взаимно однозначно и гладко (класса См) на одну из координатных гиперплоскостей;
Примем количество срезывающих функций за Т + 1. Обозначим их как
т
{¥т (х)}Т=0, $>т (х) 1
т=0
Функцию ^>о (х) подберем так, чтобы множество точек Ер0, где она равна единице, находилось не ближе, чем на некотором расстоянии £о от границы области:
р(Е^о, Мп\П)> £о, Е^о := {х : <ро(х) = 1}.
Конкретный способ построения функций разбиения единицы приведем ниже. Введем обозначения Тт := виррут П П и Гт := виррут П Г.
Интеграл dxf (х) можно представить в виде суммы интегралов:
т
/ dxf (х) = У dx Ут(x)f (х). иа _п
т=0 "т
Для удобства и универсальности программной реализации область П будем задавать неявно — как множество точек, где неотрицательна гладкая функция Ф(х) :
П:= {х : Ф(х)^0}, Ф(х) € См^).
При этом граница задается как
Г := {х : Ф(х) = 0}.
Пусть градиент функции Ф(х) на границе не обращается в нуль:
УФ(х)|г = 0.
Наложенных на Ф(х) условий достаточно для того, чтобы уравнение для каждого из участков границы в достаточно малой окрестности каждой точки границы можно было разрешить относительно координаты, задающей направление, ортогональное гиперплоскости, на которую мы гладко проектируем Гт :
Ут 3 2т , 3 ^т (xjт ) : Гт = {х € Тт : xjт Пт (xjт )},
где Xjт := (хл,...,х^т-1,х^т+1,...,хп).
Если каждой области Тт сопоставить функционал погрешности
К (х):= ХГт (х) - 53 ек (Н)6(х - hk),
Нк € Тт, к € Zn
то общий функционал погрешности можно собрать из таких функционалов так:
т
1а(х) --=^2 Ут (х)К(х).
т =0
При этом коэффициенты Ск(Н) функционала погрешности (х) получаются следующим образом:
ек (н) = ^2 Ут (х)ек (н).
т=0
Так как в определении функционала погрешности 1^£(х) каждый функционал 1^т (х) умножен на срезывающую функцию ут (х), мы можем в дальнейшем произвольно доопределять К (х) вне вирр ут (х), не изменяя итоговой суммы.
Для того чтобы функционал погрешности 1£(х) был оптимальным по порядку в (П), достаточно оптимальности по порядку каждого из функционалов Г^т (х). Необходимо подобрать для каждого т коэффициенты ек(Н) так, чтобы выполнялись оценки
\\К Н(^(а))*-°(НШ).
Приведем результат, полученный в предположении, что уравнение для границы Гт задается функцией, выражающей последнюю координату вектора х через остальные:
Гт := {х е Тт : Хп = ч(х')}, X := (х1,..., х„_1),
и выполняется хп^7(х■'), Ух е Тт.
Искомые коэффициенты:
ЕМ +1 1 М +1 ъ—1^^.т[п{кп, — а—1,М +1} ^Т"'т'т{кп — а—г,М +1} г -^г) ,
р=1 Р 1^г=1 п 2^г=1 "Угг 2-^т = 1 'итр, кп^2 + а
-к := '
0, кп^1 + а.
где — элементы матрицы, обратной к матрице, составленной из элементов определителя Вандермонда.
Эту формулу для ускорения вычислений лучше переписать в виде:
Г ЕТ=п{гп,М} п^—г ЕМ=0 фдг1, и>0
ст (Ь) = I ,
[ о, гГ1 < о
м 1 г
£ : к , 1п : кп а 2; Г&1у : ^1+1, З+1 , ^, ^ ^ ^ ^ ^ . 1 ^ЗР , ^ ^ .
Р=0 Р + З = 0
Для практической реализации алгоритма необходимо построить Т + 1 срезывающих функций {ут(х)}т=о, 5^т=о Ут(х)\хвй = 1, удовлетворяющих наложенным выше условиям.
Так как в алгоритме фактически используются лишь внутренние узлы области О, можно не следить за разбиением единицы вне области и не соблюдать условие периодичности срезывающих функций. То же касается и подынтегральной функции ](х) — она не обязана быть периодичной и определенной вне области О.
Разбиение единицы произведем в два этапа. Вначале построим центральную срезывающую функцию уо(х), затем построим разбиение единицы для гиперкуба Q: {фт(х)}т=1, 5^т=1 Фт(х)\хея = 1. И, наконец, получим окончательный вид срезываю-
щих функций: ут(х) := фт(х)(1 — уо(х)), т = 1..Т.
При этом ХГ=0 Уг(х) = Уо(х) + Е,Т=1 Ут(х) = Уо(х) + (1 — Уо(х)) ^Т=1 Фт(х) = 1, х е Q.
Нам понадобится вспомогательная функция С(г) :
0, г < 0
С(г) := { 1(0 ^(г)вг/ /о х(г)в;ь, х(г) = (г(1 — г))М, 0^ь < 1.
1, 1^Ь.
Центральную срезывающую функцию уо(х) зададим так:
Ф(х) £ 1
У0(х) := £
£2 — £1 £2 — £1
где Ф(х) — функция, задающая область О.
0, Ф(х)^£1
Тогда уо(х) = ^ р е (0,1), £ 1 < Ф(х) < £2.
1, £2<Ф(х).
Таким образом, конкретный вид функции уо(х) зависит как от способа задания области Q, так и от параметров е\ и е2.
Построим функции {фт (ж)}т=1. Потребуем, чтобы центр гиперкуба Q принадлежал области Q. Тогда, учитывая выпуклость области, можно построить ровно 2п срезывающих функций (по числу граней гиперкуба), удовлетворяющих условиям, наложенным на них выше. Далее будем считать T = 2п.
Для удобства будем нумеровать координатные оси от 0 до п — 1, а систему функций {фт (x)}2=i заменим на систему
{pdim,dir (x)}dim=0..n-1, dir=±1.
Первый индекс dim (от dimension) означает размерность, второй индекс dir (от direction)
— направление (отрицательное или положительное). Будем строить систему срезывающих функций так, чтобы функция 4[!dim,dir (x) «опиралась» на грань гиперкуба Q, ортогональную оси Xdim, и такой, чтобы Xdim = &е±Х для любой точки х, принадлежащей этой грани. Под выражением <^dim.,dir (х) опирается на грань» имеется в виду то, что соответствующая срезывающая функция (pdim,dir (х), определяемая из
<Pdim,dir(х) ■= Фdim,dir(х)(1 — фо(х)), dim = 0..п — 1, dir = ±1,
должна гладко и взаимнооднозначно проектировать вырезаемый ей участок границы области Q на координатную плоскость, содержащую эту грань.
Зададим вначале разбиение единицы для двумерного случая. Введем вспомогательные функции:
Ф(г,х) ■= £(хг)£((1 — t)x),
A_ 1(t, b, c) ■= —t + b, A1(t, b, c) ■= -t — (--—.
c c c
Рассмотрим плоскость {хо,х1}. Зададим разбиение единицы:
фо,-1(хо, хл) ■= ф(х 1, А-1(хо, -, c)), фо,1(хо,х1) ■= ф(х1, А1(хо, —, c))фl,—l(xо, х1) + ф 1д(хо,х1) ■=
■ =1 — (фо,-1(хо,х1) + фо,1(хо,х1)).
Основой для построения разбиения единицы в многомерном случае будут двумерные разбиения во всевозможных координатных плоскостях. В n-мерном пространстве таких
п(п— 1) « п
плоскостей ровно —^"2—- — число сочетаний из п по 2.
Каждую из плоскостей {xi, xj} будем разбивать двумя функциями:
$i,j ■ = rpi,-1(xi,xj) + 'pi,1(xi,xj ),
Фл ■ = 'pj-1(xi,xj) + 'pj,1(xi,xj).
Имеем:
Фо,1 + О, =1
Фо,п-1 + Фп-1,о =1
$1,2 + Ф2,1 =1
$1,п-1 + Фп-1,1 =1
Фп-2,п-1 + Фп-1,п-2 =1
Перемножив левые части этих уравнений, получим выражение с 2) слагаемыми. Каждое из этих слагаемых представляет собой функцию п переменных и может быть использовано в качестве одной из многомерных функций разбиения единицы. Однако количество слагаемых может быть очень велико — в десятимерном случае их 245 « 35 • 1012 штук. Сгруппируем их в п слагаемых следующим образом:
' |о : = УП-1,
Ф1 : = У П-и „ „ „ „
ф2 : = 1 - У0 - УФ° : = Фо,1 Фо,2 ... Фо,1,
< ф3 : = У 2 + У1 - 'У0 - Уз, , где У1 : = Ф1’0 Ф 1,2 ■■■ $1*’
г = 0..п - 1.
$п-1: = У0п-2 + УП-2 - УП-1 - УП-1,
Каждая из полученных функций ф^ легко разбивается на два слагаемых в зависимости от знака выражения (х^ - 0.5). Таким образом, получаем необходимые 2п функций фт разбиения единицы для гиперкуба Q, дающие искомую систему срезывающих функций. Для наглядности приведем их на одном графике для п = 2, М = 2, Ь = 6, с = 0.5 (рис. 1).
Рис. 1. Функции {ут}2п=1
2.3. Программа
Программа “СиЬаШ;” [10], предназначенная для вычисления многомерных интегралов по ограниченным выпуклым областям с гладкими границами, тестировалась при следующих параметрах:
1. п от 2 до 10.
2 I(х) = ЕГ=1 аХ ■
3. М от 2 до 6.
4. N от 10 до 105, где N - количество точек решетки на ребре единичного куба Q■ При этом шаг Н = N-1 ■
5. П = {х : Ф(х) = 0, Ф(х) = 1 — ^П=1 Сг(Хг — 0.5)^ }■
6. Р от 1 до 1000.
Для примера возьмем а = (2,1, 2,1,2,1), Ь = (2,4, 2, 42, 4),с =
(6.25, 39.0625, ■■■, 6.25, 39.0625), 3 = (2, 4, 2,4..., 2,4),
Точность вычислений рассчитывалась по устойчивости десятичных знаков в ответе при уменьшении параметра Н■ Независимыми параметрами являются п, М, N и Р. Приведем некоторые результаты экспериментов.
Таблица 1. п = 2, порядки теоретических и экспериментальных погрешностей
пса ■ -ш - ш -т - _шя ■ - ш я - ш - я -
50 | 3 2 2 1 ■ 50 нн 7 9 11
100 4 4 3 2 2 100 4 5 5 10 12
200 7 5 4 5 3 200 5 7 10 12 14
400 5 9 11 7 7 400 5 5 11 14 15
БОО 5 10 12 13 14 БОО 5 9 12 15 15
1500 9 11 13 15 15 1500 7 10 13 17 20
3200 10 12 15 15 17 3200 5 11 15 15 22
■5400 11 14 15 1В 1В ■5400 5 12 15 20 23
12500 12 15 17 15 17л 12500 9 13 17 21 25а
Таблица 2. п = 3, порядки теоретических и экспериментальных погрешностей
шив а ла -ш -т - _шя ш - я -в - ■ - я -
50 нн 3 2 2 1 1 »■ нн 7 9 11
100 4 4 3 3 3 100 4 5 5 10 12
200 5 7 5 4 4 200 5 7 10 12 14
400 5 5 5 5 7 400 5 5 11 И 15
БОО 9 9 10 9 9 БОО 5 9 12 15 15
1500 9 10 11 11 10 1500 7 10 13 17 20
Таблица 3. п = 5, порядки теоретических и экспериментальных погрешностей
----- - - -- -- -
23 4 3 3 3 2 25 3 5 5 7?
25 ■з ■з 1 3 2 | 25 5 7 9
50 4 4 4 3 3 50 4 5 7 9 11
75 5 4 4 3 2 75 4 5 Б 10 12
100 5 4 4 4 3 100 4 5 5 10 12
125 5 5 4 4 4 125 5 7 9 11 13
150 7 5 5 5 4 150 5 7 9 :: 14
175 7 7 5 5 4 175 5 7 9 12 14
200 7 7 5 5 5 200 5 7 10 12 14
В табл. 1-3 показаны порядки абсолютных погрешностей, полученных как при непосредственных вычислениях, так и при теоретической оценке. Как видно, результаты вычислений примерно соответствуют ожидаемым теоретически, уступая им, в основном, по двум причинам. Во-первых, при малом количестве точек N и большом M фактическая толщина пограничного слоя не обеспечивает нужных для обоснования теоретической точности свойств (вмещать 2M узлов). Во-вторых, при использовании типа long double с 18 значащими цифрами нет возможности увеличить точность выше 18-го порядка. Также не следует забывать, что при расчете теоретической погрешности мы не учитывали норму подынтегральной функции, которая может ухудшить результат.
Перейдем теперь к анализу быстроты счета и качества распараллеливания программы. Для этого введем понятия ускорения и эффективности программы:
Sp = TP ’
где Тр — время, за которое задача выполняется на Р процессорах.
На рис. 2 показаны отклонения экспериментальных ускорений Бр (темные ломаные) от идеальных ускорений (светлые прямые).
п=5, N=100, М=2
Sp
УШ 1
400 -300 -
3—f Ы
100 -
1 1 1
100 200 300 400 J00 600 700 800 900
• Ideal------Experimental |
Рис. 2. Ускорение вычислений
Как видно, эффективность распараллеливания постепенно падает с увеличением числа процессоров. Это связано, в первую очередь, с особенностями распределения вычислительной работы по процессорам. Каждый процессор получает примерно одинаковое количество узлов, однако вычислительная сложность в разных узлах различна. В результате процессоры не всегда загружаются равномерно, особенно при сложной форме области.
3. Построение разбиения области интегрирования
3.1. Постановка задачи
Проблема вычисления интегралов по многомерным областям решается путем приближения интеграла решетчатыми кубатурными формулами с ограниченным пограничным слоем
(ОПС-формулы [4,5]). Этот алгоритм предусматривает вычисление интеграла как суммы локальных интегралов по областям более простых форм, которые являются подобластями
О.
Пусть семейство {От }т=0 образует покрытые области О. Это взаимно пересекающиеся открытые области с кусочно-гладкими границами, удовлетворяющие главному свойству: часть границы Г, попавшая в От может быть задана явной формулой, выражающей одну из координат однозначной функцией класса См от остальных координат, Х8. = ^3 (XI, ...,Х5._г ,Х5.+1 ,...Хп).
Такую систему подобластей можно выбрать с помощью разбиения единицы. Это конечный набор финитных функций {ут}т=о, сумма которых тождественно равна единице: т
ут(х) = 1, для х € Q. У нас они принадлежат пространству См(^).
т=0
Для каждой функции ут строится своя кубатурная формула.
Таким образом, итоговая кубатурная формула для области О выглядит следующим образом: К^(Н) = Е Е С1 (Н)ут(Нк)/(Нк), От = О П вирр^у.
т=0 кН£^
То есть приближение интеграла кубатурной формулой сводится к приближению интеграла суммированием кубатурных формул более простого вида.
Наша цель - создать программу, которая автоматически разбивает области произвольных форм на нужные подобласти. Эта задача решена для произвольных двумерных ограниченных областей с гладкими границами.
3.2. Описание алгоритма
Дадим описание конструктивного алгоритма разбиения единицы для такой области О. В общем случае расположенная в единичном кубе Q = [0,1)” область О задается следующим образом: О = {х|Ф(х) ^ 0} С Q, Ф(х) € См([0,1]”). Ее граница Г={х|Ф(х) = 0}. Достаточным условием гладкости границы является условие ДФ(х)| = 0.
Для дальнейших построений будем использовать вспомогательную функцию £(£), определенную в п. 2.
Разбиение единицы проведем в два этапа. Сначала построим центральную срезывающую
функцию у>о(х), затем приграничные срезывающие функции. Получим окончательный вид
____ т
срезывающих функций: ут(х) = Фт(х)(1 — уо(х)), т = 1..Т и ^ ут = 1 в некоторой окрест-
т=0
ности О.
Центральную срезывающую функцию зададим как в п. 2: ^(х) = £( ------£1—). Па-
£2 £ 1 £2 £ 1
раметры е\ и £2 подбираются исходя из условия, чтобы функция у>0(х) была сосредоточена
в О и равна 0 в некоторой окрестности границы Г.
Мы решаем задачу для двумерных областей О.
Возьмем простейшую область - круг, лежащий в единичном квадрате, со следующими
значениями параметров е1=0.2, £2=0.5, а =0.3, в=0.05, М=3.
4
Наше разбиение единицы имеет следующий вид: ^(у) + Е Ут (у) = ^(у) + (1 —
т =1
—у0(у))(фг(у2) + Ф2(у2) + 'Фз(у) + Ф4(у)) = 1, где функции Фт(у) можно задать следующим образом: ф^) = £(У2 — ), Ф2Ы = £(—У2 — ), Фз(у) = £(У1 — ) • (1 —
—Ф1(У2) — Ф2(У2)), Ф4(у) = (1 — £(У1 — 1-в-в)) • (1 — Ф1 (У2) — Ф2(У2)), где а и в подбирают-
4
ся в зависимости от ширины пограничного слоя. Причем ^ Фт = 1 на носителе функции
т=1
Рис. 3. Область {у|Фі(у) > 0}
(1 — У0(Е)), тогда и £ ут(у) = 1, У € [0,1]2.
т=0
Таким образом, функции {ут}4=0 - искомые срезывающие функции для круга.
Если бы это был п-мерный шар, то количество срезывающих функций было бы равно 2п+1, где п - размерность пространства. Для круга их было 4 (вместе с центральной - 5).
Перейдем к построению разбиения единицы для произвольной двумерной области с гладкой границей.
Сделаем замену переменных х = |^ф(^)|. Она является гладким отображением приграничной полоски, в которой 1 — у0 (х) = 0, на единичную сферу, У = 1. А в переменных у применим функции Фт (у), £ Фт (у)1|у| = 1 = 1.
Приведем примеры разбиения единицы для некоторых невыпуклых областей. На рис. 4 изображено разбиение для области {х|Ф2(х) > 0} с границей Ф2(х) = 1 —
— /8(х1 — 0.5)2 + 8(х2 — 0.5)2 + 5)2+(Х—о.5)2)2 =0 с параметрами £1=0.2, £2=0.5, а
=0.3, в=0.05, М=3.
Рис. 4. Область Ф2(х) ^ 0
Следует отметить, что предложенный алгоритм получения разбиения единицы универсален. Он не зависит от формы и расположения области.
Например, повернем область, заданную на рис. 4, на 1 радиан. Тогда уравнение границы запишется в следующем виде:
+
Фз(х) = 1 -у/8(хі — 0.5)2 + 8(х2 - 0.5)2 + 2(х1 — 0.5)(х2 — 0.5) яіп(1) + ео8(1)((х1 — 0.5)2 — (х2 — 0.5)2)2
0.
((х 1 — 0.5)2 + (х2 — 0.5)2 )2
Заметим, что разбиение единицы автоматически подстроилось под новую форму области (относительно декартовой системы координат) (рис. 5).
Рис. 5. Область Ф3(х) ^ 0
Проиллюстрируем результат еще несколькими примерами разбиения невыпуклых областей более сложных форм. На рис. 6 изображена область с границей Ф4(х) = ^ — 625((хі —
0.5)2 + (х2 —0.7)2 )((хі —0.2)2 + (х2 —0.2)2)-((хі —0.4)4+(х2 —0.9)2)++(хі —0.4)2 + (х2 —0.6)2 =0
Рис. 6. Двусвязная область {х|Ф4(х) ^ 0}
Ф5(х) = уб — 150((хі—0.3) +(х2—0.7) )((хі—0.2) + (х2—0.2) ) = 0 На рис. 8 и 9 изображены соответственно области
Фб(х) — 16^/ (хі—0.5) + (х2 —0.5) — 5+
Рис. 7. Область {ж|ф5(ж) ^ 0}
Рис. 8. Область {ж|ф6(ж) ^ 0}
+ (х1-0.5)((х1-0.5)4 + -10(х1-0.5)2(х2-0.5)2 + 5(х2-0.5)4) = 0 ((х1-0.5)2 + (х2-0.5)2) 2
и
ф7(х) = ((х1-0.4)2 + (х2-0.4)2) - 136((х1-0.4)2 + (х2-0.4)2) - ((4х1-0.9)2+
+(4х2-0.9)2 - -1) • 9((4х1-0.9)2 + (4х2-0.9)2) = 0
3.3. Описание программы
По изложенному алгоритму была создана программа на языке программирования С+—Н, подробно закомментированная. Перед ее использованием необходимо провести вспомогательные действия: подсчитать производные границы области интегрирования П, которые необходимы для вычисления градиента.
Входными данными в программе являются гладкость функции М, ширина пограничного слоя а и в, параметры £1 и £2 и количество точек N на стороне единичного квадрата. Главный цикл программы пробегает по всем узлам решетки, и для каждой срезывающей функции рисуется свой график по точкам. Координаты точек записываются в соответствующий текстовый файл, а рисунок выводится на экран.
Рис. 9. Область {х|ф7(х) ^ 0}
4. Решение интегральных уравнений
4.1. Алгоритм
Напомним постановку задачи.
Пусть в области П С Я2 задана функция /(х) и функция К(х, у). Наша задача - найти функцию и(х), удовлетворяющую уравнению
и(х) - J К(х,у)и(у)ву = ](х),х € П С Я2. п
Здесь П - ограниченная замкнутая область с гладкой границей, лежащая внутри единичного квадрата (0, I)2, К € См(П х П) и ] € См(П).
В связи с тем, что изначально планировалось применение многопроцессорной вычислительной техники, необходимо было использовать такие алгоритмы численного решения интегральных уравнений, при распараллеливании которых не возникнут трудности. Для решения данного уравнения был выбран итерационный метод, позволяющий получить относительно простые численные алгоритмы решения интегральных уравнений. Предположим, что ЦК||с(пхй) = @ < 1. Можно применить метод последовательных приближений. За начальное приближение берется ио(х) = ](х). Последующие приближения строятся по
рекуррентной формуле ия+1(х) = f (х) + § К(х,у)и8(у)с1у, в = 1, 2,.... Последовательность
п
функций {ия(х)} является приближением к искомому решению данного уравнения в норме пространства С(17). Было показано, что решение интегрального уравнения и(х) обладает той же гладкостью, что ядро интегрального оператора и правая часть уравнения.
При численной реализации этого метода для вычисления интеграла на каждой итерации используется решетчатая кубатурная формула с ограниченным пограничным слоем (ОПС-формула) (см. п. 2). Для приближенного вычисления интеграла создан алгоритм (см. [5,7,8]), который позволяет вычислить коэффициенты кубатурной формулы, приближающей точное значение интеграла с погрешностью порядка (Н/а)т, т - гладкость интегрируемой функции, а - параметр локализации области интегрирования (см. п. 3).
Отметим, что выбор решетчатых кубатурных формул позволяет достичь необходимой универсальности и единообразия в интегрировании по областям произвольных форм вследствие того, что последовательности решеток не зависят от форм областей интегрирования.
А использование ОПС-формул позволяет использовать алгоритмы вычисления коэффициентов формул с хорошими и теоретически обоснованными аппроксимационными свойствами.
4.2. Описание программы
В этой части рассмотрена возможность применения решетчатых кубатурных формул к решению интегральных уравнений с произвольной формой области интегрирования. Исходя из изложенных алгоритмов была написана программа численного решения интегрального уравнения вида (1). Программа написана на языке С+—Н с использованием библиотеки параллельных функций MPI. Методом последовательных приближений вычисляется последовательность {us}, сходящаяся к точному решению данного интегрального уравнения. Процесс вычисления заканчивается при выполнении условия ||us — ws-i|| < е, где \\п\\ = max |u(x)|, е - заранее заданная точность.
X
Входные данные программы:
1. Область интегрирования П, которая задается в неявном виде: П = {ж|Ф(ж) ^ 0},
Ф(х) е CM(Q), П С (0,1)2.
2. Функция K(xi,yi,x2,y2), max K(x,y)| < 1.
x,y
3. Функция f (xi,x2).
4. Шаг кубической решетки h < 0.01, он задается как h = 1/N, где N - натуральное число.
5. Параметры срезывающих функций.
6. Параметр гладкости M.
7. Количество процессоров P.
Приближенным решением интегрального уравнения является результат последней итерации, который выводится в виде табличного значения функции us(x) в узлах x = hk е
е П С R2.
Программа тестировалась на суперкомпьютерах МВС-15000ВМ и МВС-50К Межведомственного Суперкомпьютерного Центра РАН.
Приведем результаты вычислительного эксперимента решения интегрального уравнения с выпуклой областью интегрирования [3]. Для того чтобы иметь возможность сравнения с точным решением уравнения, мы задавались функцией, которая должна быть точным решением, и по ней вычислялась правая часть уравнения. Вычислительный эксперимент проводился с этой правой частью.
Расчеты проводились со следующими данными. В качестве точного решения была взята функция u(x) = (xi — x2)5.
1. K(x, y) = (xiyi + x2y2)3.
2. f (x) = (xi — x2)5 — (0.00099x3 + 0.0007x2x2 — 0.00059x2xi — 0.00093x2).
3. Область интегрирования П = {x^i(x) ^ 0}, Ф1(x) = 1 — ((xi — 0.5)/0.4)2 —
— ((x2 — 0.5)/0.4)4, (рис. 10).
4. Шаг кубической решетки h = 1/N, N = 200, 300.
5. Гладкость M = 2, 3.
6. Параметры срезывающих функций , е2 = 0.5, а = 0.1, в = 0.05.
7. Количество процессоров 100-1000.
Теоретическая ожидаемая точность - O(hM), без учета влияния ширины пограничного слоя, в данном примере M =2, h = 1/200, ожидаемая точность - 10~5.
Точность вычислений рассчитывалась двумя способами. Сравнением с известным точным решением (погрешность обозначена е') и по правилу Рунге - по устойчивости десятич-
Рис. 10. Область Ф^х) ^ 0
ных знаков в ответе при уменьшении Н (погрешность обозначена е). В табл. 4 приведены результаты вычислительного эксперимента.
Таблица 4. Достигнутая точность при учтенной гладкости М = 2
к - номер итерации Н = 1/200 Н = 1/300
е е
1 0.0987 0.0988
2 0.000776 0.000778
3 0.0000268 0.0000269
4 0.00000310 0.00000310
В табл. 4 точность рассчитывалась по устойчивости десятичных знаков. При сравнении результата последней итерации с точным решением и(х) = (х1 — ж2)5 получили совпадение 5-6 знаков после запятой.
Следует отметить, что программа предназначена для интегральных уравнений вида (1) с произвольной формой области интегрирования. В связи с чем был проведен вычислительный эксперимент для невыпуклой области. Приведем результаты тестирования для области Ф2(х) ^ 0, Ф2(х) = 1 — -у/8(х1 — 0.5)2 + 8(х2 — 0.5)2+ +(2(х1 — 0.5)(х2 — 0.5) вт(1) + ео8(1)((х1 — 0.5)2 — (х2 — 0.5)2)2)/((х1 — 0.5)2 + (х2 — 0.5)2)2, изображенной на рис. 5, с параметрами срезывающих функций: е1 = 0.1, е2 = 0.5, а = 0.1, в = 0.05.
Точность рассчитывалась по устойчивости десятичных знаков. Получили, что теоретически ожидаемая точность (Н/а)м порядка 10~5 достигается на пятой итерации при N = 250, М = 3, а = 0.1.
Далее приведем анализ эффективности распараллеливания. В табл. 5 указано время работы программы на разном количестве процессоров и подсчитаны ускорение и эффективность.
Таблица 5. Время счета (с) при разном числе Р процессоров, Н = 1/200, М = 2.
Достигнутая точность 10~5
Р 10 20 30 40 50 100
Тр 230 127 93 75 63 36
Sp 10 18 24.7 36 42.6 63.8
Ер 1 0.9 0.8 0.8 0.72 0.63
На графике (рис. 11) изображена зависимость ускорения от количества процессоров.
p
Ideal
Experiment
Рис. 11. Зависимость ускорения от количества процессоров
5. Заключение
Суммируем выводы по описанным выше программам.
Создана программа приближенного интегрирования “СиЬаШ;”. Анализ скорости работы программы показал, что эффективность распараллеливания снижается с увеличением числа процессоров, и это связано с, вообще говоря, различной сложностью вычисления коэффициентов кубатурной формулы в разных узлах. Отметим, что даже в случае идеальной эффективности программы нет смысла увеличивать размерность пространства выше 10 и количество узлов выше 1012.
Разработан алгоритм автоматического разбиения двумерных областей произвольных форм с гладкими границами на подобласти и создана программа, реализующая этот алгоритм. Этот алгоритм позволяет построить разбиение и в п-мерном случае (даже для несвязных и не односвязных областей). Однако метод не подходит для областей с кусочногладкими границами.
Для наглядности показаны картинки только для двумерного случая, для размерностей больше 3 осложняется иллюстрация результатов, но алгоритм также работает. Результаты опубликованы в [2].
И наконец, оправдано применение ОПС-формул для численного решения интегральных уравнений. Применение итерационного метода в сочетании с решетчатой кубатурной ОПС-формулой позволяет достичь точности порядка 10~5 за 5-6 итераций. Отметим, что этот алгоритм решения интегральных уравнений позволяет хорошо распараллелить вычислительную программу для применения многопроцессорных вычислительных систем. На данный момент эффективность распараллеливания около 63 %. Это связано с тем, что количество последовательных операций становится сравнимым с количеством параллельных операций и значительная часть процессоров простаивает. В дальнейшем планируется повысить эту эффективность путем внесения изменений в алгоритм.
Работа выполнена при поддержке Программы № 14 Президиума РАН, Программы целевых расходов Президиума РАН «Поддержка молодых ученых» и Российского фонда фун-
даментальных исследований (№ 06-01-00597).
Список литературы
[1] Игнатьев А.Н. Универсальный алгоритм вычисления интегралов по ограниченным областям с гладкими границами. Кубатурные формулы и их приложения. Уфа: ИМВЦ УНЦ РАН, 1996. С. 21-31.
[2] Валеева Л.З. Применение локализации алгоритма кубатурных формул для решения интегральных уравнений. Сб. статей II Межд. конференции "Математическое и компьютерное моделирование естественнонаучных и социальных проблем". Пенза, 2008. С. 14-17.
[3] Банникова Е.Л. Программа численного решения интегральных уравнений Фредгольма второго рода “IntUr”. Свидетельство №10418 от 15.04.2008.
[4] Рахматуллин Д.Я. Интегрирование функций по выпуклым областям на многопроцессорных вычислительных системах: Автореф. дис... канд. физ.-мат. наук. - Уфа: ИМВЦ УНЦ РАН, 2006. - 22 с.
[5] Рамазанов М.Д. Теория решетчатых кубатурных формул с ограниченным пограничным слоем. Препринт 2007-1 - Уфа: ИМВЦ УНЦ РАН, 2007. - 105 с.
[6] Рамазанов М.Д. Лекции по теории кубатурных формул. Уфа: Изд-во БашГУ, 1973. 177 с.
[7] Соболев С.Л., Васкевич В.Л. Кубатурные формулы. Новосибирск: Изд-во ИМ СО РАН, 1996. 484 с. С. 254-263.
[8] Рахматуллин Д.Я. Вычисление интегралов по многомерным областям на многопроцессорных вычислительных системах // Вычислительные технологии. Т. 11, 3, 2006. С. 118-125.
[9] M.D.Ramazanov To the Lp-Theory of Sobolev Formulas //Siberian advances in mathematics. 1999. Vol. 9, N 1. P. 99-125.
[10] Рахматуллин Д.Я. Программа интегрирования по многомерным областям “CubaInt”. Свидетельство № 2007614331 от 10.10.2007.
Solving of Integral Equations on Multiprocessing Computing Systems
Marat D. Ramazanov, Djangir Y. Rakhmatullin, Leysan Z. Valeeva and Ekaterina L. Bannikova
Institute of Mathematics with Computing Centre 112 Chernyshevsky str., Ufa, 450077 Russia
It is described Fredholm equations for the various forms of domains. It is reduced the description of the algorithm of approximate calculation of integrals using multiprocessing technology.
Key words: multiprocessing systems, approximate calculation of integrals, multiprocessing technology.
