ВУЗ: Не указан
Категория: Не указан
Дисциплина: Не указана
Добавлен: 06.04.2021
Просмотров: 2130
Скачиваний: 4
УДК 517.968
ПРИМЕНЕНИЕ МЕТОДА НЕПОЛНОЙ
КРЕСТОВОЙ АППРОКСИМАЦИИ ДЛЯ
ЧИСЛЕННОГО РЕШЕНИЯ ЗАДАЧИ
ДИРИХЛЕ ДЛЯ УРАВНЕНИЯ
ГЕЛЬМГОЛЬЦА
М.Ю. Талтыкина
Вычислительный центр ДВО РАН
Россия, 680000, Хабаровск, Ким-Ю-Чена 65
E-mail:
taltykina@yandex.ru
А.А. Каширин
Вычислительный центр ДВО РАН
Россия, 680000, Хабаровск, Ким-Ю-Чена 65
E-mail:
elomer@mail.ru
Ключевые слова:
задачи Дирихле, уравнение Гельмгольца, быстрые ме-
тоды, мозаично-скелетонный алгоритм, неполная крестовая аппроксимация
Рассматриваются быстрые методы численного решения систем линейных
алгебраических уравнений с плотно заполненными матрицами, аппрокси-
мирующих интегральные уравнения задачи Дирихле. С целью снижения
времени работы программы, предназначенной для решения задачи Ди-
рихле для уравнения Гельмгольца, предлагается использовать мозаично-
скелетонный алгоритм на основе неполной крестовой аппроксимации. При-
ведены результаты экспериментов, демонстрирующие эффективность пред-
лагаемого метода.
Введение
Задачи Дирихле для уравнения Гельмгольца являются классическими задача-
ми математической физики. Аналитические решения таких задач могут быть най-
дены лишь в случае простейших областей. Вследствие этого, большое развитие по-
лучили методы их численного решения. Численное решение предполагает предва-
рительное построение дискретного аналога исходной задачи. При его построении
необходимо учитывать, что искомые решения зависят от трех пространственных
переменных и при больших действительных волновых числах являются быстро ос-
циллирующими функциями. Кроме этого, решения внешних краевых задач отыски-
ваются в неограниченных областях, где должны удовлетворять условиям излучения
на бесконечности. Отмеченные свойства приводят к тому, что дискретные аналоги
Сборник материалов XXXVII Дальневосточной Математической Школы-Семинара
имени академика Е.В. Золотова, Владивосток, 8 – 14 сентября 2013 г.
221
задач Дирихле, построенные методами конечных элементов или конечных разно-
стей, предъявляют весьма высокие требования к ресурсам компьютера. С вычисли-
тельной точки зрения более эффективной представляется дискретизация исходных
задач, сформулированных в виде эквивалентных им граничных интегральных урав-
нений. В этом случае неизвестные функции отыскиваются на компактных границах
областей, что особенно важно при численном решении внешних краевых задач, так
как при этом существенно понижается их вычислительная сложность. Численное
решение интегральных уравнений приводит к решению систем линейных алгебра-
ических уравнений (СЛАУ). Вычислительная сложность их решения прямыми ме-
тодами имеет оценку
O
(
n
3
)
, где
n
– порядок СЛАУ. Однако свойства получаемых
матриц таковы, что использование для поиска приближенных решений обобщенно-
го метода минимальных невязок (GMRES) [1] позволяет понизить эту сложность до
O
(
n
2
)
. Но даже в этом случае требования к ресурсам компьютера остаются высоки-
ми. Эта проблема стимулировала развитие быстрых методов. Называя метод «быст-
рым», обычно имеют в виду, что его сложность составляет
o
(
n
2
)
при
n
→ ∞
. После
нахождения приближенных решений интегральных уравнений решения исходных
задач при помощи интегральных представлений достаточно просто восстанавлива-
ются в любой точке пространства.
1.
Постановка задачи
Рассмотрим трехмерное евклидово пространство
R
3
с ортогональной системой
координат
ox
1
x
2
x
3
. Пусть в этом пространстве имеется произвольная замкнутая
липшицева поверхность
Γ
, разделяющая его на внутреннюю область
Ω
i
и внешнюю
область
Ω
e
=
R
3
\
Ω
i
Проблема 1.
(внешняя задача Дирихле для уравнения Гельмгольца). Найти функ-
цию
u
e
(
x
)
∈
H
1
(Ω
e
)
, удовлетворяющую уравнению Гельмгольца
∆
u
e
+
k
2
e
u
e
= 0
,
x
∈
Ω
e
,
(1)
граничному условию
γu
e
(
x
) =
f
(
x
)
,
(2)
и условию излучения на бесконечности
∂u
e
/
∂
|
x
| −
ik
e
u
e
=
o
|
x
|
−
1
,
|
x
| → ∞
.
(3)
Здесь
k
— волновое число,
f
∈
H
1
/
2
(Γ)
— известная функция,
H
1
(Ω
e
) =
{
u
:
||
u
||
2
H
1
(Ω
e
)
=
||
u
||
2
H
(Ω
e
)
+
||∇
u
||
2
H
(Ω
e
)
<
∞}
,
H
1
/
2
(Γ) =
{
u
:
||
u
||
H
1
/
2
(Γ)
= inf
γν
=
u
||
ν
||
H
1
(Ω
e
)
<
∞}
.
Теорема 1.
Для любой функции
f
∈
H
1
/
2
(Γ)
существует единственное решение
внешней задачи 1 из пространства
H
1
(Ω
e
)
.
Решение задачи 1 будем искать в виде потенциала простого слоя
u
e
(
x
) = (
S
e
q
e
) (
x
)
≡ h
G
e
(
x,
·
)
, q
e
i
Γ
,
x
∈
Ω
e
,
(4)
Ядром интегрального оператора (4) является фундаментальное решение уравнения
Гельмгольца, поэтому функция
u
e
удовлетворяет уравнению (1) в соответствующей
Сборник материалов XXXVII Дальневосточной Математической Школы-Семинара
имени академика Е.В. Золотова, Владивосток, 8 – 14 сентября 2013 г.
222
области и условию излучения на бесконечности (3). Эта функция будет решением
задачи 1, если подобрать плотность
q
e
так, чтобы
u
e
удовлетворяла граничному
условию (2). Таким образом, задача 1 сводится к граничному тождеству
h
S
e
q
e
, µ
i
Γ
=
h
f, µ
i
Γ
∀
µ
∈
H
−
1/2
(Γ)
.
(5)
Теорема 2.
Пусть
Im(
k
e
)
>
0
или
k
2
e
не является собственным значением задачи
∆
u
+
k
2
e
u
= 0
,
x
∈
Ω
i
,
γu
= 0
.
(6)
Тогда уравнение (5) корректно разрешимо в пространстве
H
−
1
/
2
(Γ)
и формула (4)
дает решение задачи 1.
2.
Описание метода
Численное решение интегральных уравнений приводит к системам линейных
алгебраических уравнений с плотно заполненными матрицами. Существуют мето-
ды приближенного быстрого матрично-векторного умножения для таким матриц.
Метод считается быстрым, если его сложность составляет
o
(
n
2
)
, где
n
— размер-
ность матрицы. Возможность построения быстрого метода может быть основана на
обнаружении и использовании некоторой особенности возникающих плотных мат-
риц. В качестве особенности можно выделять в матрице блоки, для которых может
быть построена малоранговая аппроксимация. Рассмотрим следующие известные
методы построения быстрых алгоритмов:
•
мультипольные разложения;
•
кластеризация граничных элементов;
•
интерполяция на регулярную (иерархическую) сетку;
•
мозаично-скелетонный метод.
Первые три метода являются операторными или безматричными, то есть матрица
в явном виде не участвует в операции умножения. Поскольку операторные методы
построения аппроксимаций основаны на специальном представлении исходного опе-
ратора, для их практического применения в расчетах требуется переработка всех
этапов алгоритма, начиная от процедуры дискретизации и заканчивая обработкой
результата решения линейной системы. С вычислительной точки зрения более при-
влекательными являются методы, использующие непосредственно блоки аппрокси-
мируемой матрицы. К таким методам относится, например, мозаично-скелетонный
алгоритм. Для его применения надо изменить лишь процедуру построения матрицы
системы и операцию матрично-векторного умножения. Рассмотрим основные этапы
мозаично-скелетонного алгоритма.
1. Построение дерева кластеров
2. Построение списка блоков
3. Аппроксимация блоков.
Сборник материалов XXXVII Дальневосточной Математической Школы-Семинара
имени академика Е.В. Золотова, Владивосток, 8 – 14 сентября 2013 г.
223
На этапе построения дерева кластеров погрузим область
Ω
i
в куб со стороной
α
и, поделив, каждое ребро куба пополам, разделим его на 8 подкубов. Проделывая
такую же операцию с подкубами, получим
иерархическое разбиение
исходного куба.
По дереву кластеров, матрица разбивается на иерархический набор блоков разного
размера, каждый из которых отвечает взаимодействию кластеров. Блоки, отвечаю-
щие взаимодействию достаточно удаленных друг от друга кластеров точек («даль-
няя» зона), могут быть приближены матрицей малого ранга. Блоки, полученные
взаимодействием рядом лежащих точек, определяют «ближнюю» зону. Для каждо-
го блока
A
∈
R
m
×
n
решается задача поиска малоранговой аппроксимации
||
A
−
U V
T
||
6
ε,
(7)
причем ранг матрицы
U V
T
должен быть существенно меньшим, чем ее размеры.
Задача малоранговой аппроксимации может быть решена разными способами, на-
пример, на основе метода Гаусса с выбором ведущего элемента по некоторому шаб-
лону. В работе [2] изложен метод, где в качестве такого шаблона выбирается строка
и столбец. Задача состоит в том, чтобы получить приближение матрицы
A
другой
матрицей
˜
A
r
, являющейся суммой
r
одноранговых матриц
u
p
v
T
p
/max A
(называе-
мых также скелетонами) по формуле (8)
˜
A
k
=
r
X
α
=1
u
α
v
T
α
max
A
(
i
k
, j
k
)
,
(8)
где
u
— строка матрицы
A
,
v
— столбец матрицы
A
,
max A
— элемент, стоящий
на пересечении
u
и
v
. Аппроксимация, построенная по формуле (8), называется
неполной крестовой аппроксимацией. Аппроксимация является достаточно точной,
если выполняется условие
(
n
−
k
)
||
u
k
v
T
k
||
F
6
ε
||
˜
A
k
||
F
.
(9)
Для проверки критерия 9 не требуется знать всех элементов матрицы
A
, а доста-
точно только вычисленных
k
столбцов и строк.
3.
Численные эксперименты
В работе [5] дано теоретическое обоснование существования малоранговых ап-
проксимации для осцилляционных ядер интегральных уравнений, в частности, для
фундаментального решения уравнения Гельмгольца. Программа, использующая мо-
заично-скелетонный алгоритм, написана на языке Fortran 95 на основе готовой
программы, предназначенной для решения краевых задач для уравнений Лапла-
са и Гельмгольца. Блоки «дальней» зоны аппроксимируются неполной крестовой
аппроксимацией с точностью
ε
= 10
−
7
. Блоки «ближней» зоны считаются точно. В
качестве итерационного алгоритма решения СЛАУ используется GMRES. Аппрокси-
мация «дальней» зоны и расчет «ближней» происходит один раз до начала итерации
GMRES. Умножение на каждой итерации делится на два этапа. Блоки «ближней»
зоны умножаются на вектор обычным матрично-векторным умножением. Блоки
«дальней зоны» для умножения на вектор используют специальную функцию, учи-
тывающую специфику хранения таких блоков. Данная функция требует
2
n
ариф-
метических операции для умножения блока размерностью
n
×
n
на вектор размера
Сборник материалов XXXVII Дальневосточной Математической Школы-Семинара
имени академика Е.В. Золотова, Владивосток, 8 – 14 сентября 2013 г.
224
n
. Мозаично-скелетонный алгоритм сравнивается с исходной программой, где мат-
рица считается построчно на каждой итерации GMRES и умножается на вектор
обычным матрично-векторным умножением. Задача 1 сформулирована в виде ин-
тегрального уравнения на границе области
Γ
, где
Γ
– единичная сфера с центром в
начале координат. Граничное условие:
f
(
x
) = exp(
ikx
3
)
,
x
∈
Γ
.
(10)
Волновое число равно 1. На рис. 1 приведен график зависимости времени выполне-
ния программ от размерности матриц.
0.1
1
10
100
1000
10000
510
998
2042
3998
8150
15974
32598
Время работы программы, с
Размерность матриц
Исходная программа
Мозаично−скелетонный алгоритм
Рис. 1.
Время выполнения программы
Красным цветом показаны результаты работы мозаично-скелетонного алгорит-
ма, синим — исходной программы. По результатам численных экспериментов для
данной задачи получено ускорение (в среднем в 4 раза) в работе программы. Оцен-
кой правильности работы алгоритма является сравнение погрешностей вычисления
плотностей интегральных уравнений. В таблице 1 приведены значения погрешно-
стей вычисления плотностей. Оба метода удовлетворяют необходимой точности рас-
чета плотностей.
Таблица 1.
Таблица погрешностей вычисления плотностей
Размерность матриц
Исходная программа
Мозаично-скелетонный
алгоритм
510
9
.
905
·
10
−
3
9
.
905
·
10
−
3
998
4
.
978
·
10
−
3
4
.
98
·
10
−
3
2042
2
.
407
·
10
−
3
2
.
41
·
10
−
3
3998
1
.
221
·
10
−
3
1
.
222
·
10
−
3
8150
5
.
959
·
10
−
4
5
.
959
·
10
−
4
15974
3
.
031
·
10
−
4
3
.
032
·
10
−
4
32598
1
.
486
·
10
−
4
1
.
49
·
10
−
4
На рисунке 2 показано время умножения «ближней» и «дальней» зон на век-
тор. Время умножения «ближней» зоны имеет сложность
O
(
n
2
)
, время умножения
«дальней» зоны растет почти линейно. Этот результат существенно экономит время
на каждой итерации GMRES. Итак, можно сделать вывод, что применение моза-
ично-скелетонного алгоритма на основе неполной крестовой аппроксимации эффек-
тивно для решения задачи Дирихле для уравнения Гельмгольца.
Сборник материалов XXXVII Дальневосточной Математической Школы-Семинара
имени академика Е.В. Золотова, Владивосток, 8 – 14 сентября 2013 г.