Содержание к диссертации
Введение
1 Дифференциальная и интегральная формулировки за дач дифракции. 18
1.1 Интегро-дифференциальное уравнение 18
1.2 Тензорная функция Грина прямоугольного резонатора 23
1.3 Эквивалентность дифференциальной и интегральной формулировок 28
1.4 Векторное сингулярное интегральное уравнение 37
2 Исследование векторного сингулярного интегрального уравнения задачи дифракции 39
2.1 Основные положения теории сингулярных интегральных уравнений 39
2.2 Анализ интегральных уравнений 47
2.3 Обобщенные решения в случае задачи с потерями 57
2.4 Обобщенные решения в случае задачи без потерь 58
3 Численный метод решения задачи дифракции 66
3.1 Метод Галеркина 66
3.2 Сходимость метода Галеркина 67
3.3 Построение схемы Галеркина 69
3.4 Программная реализация численного метода 73
3.4.1 Об оптимизации вычисления матрицы СЛАУ 74
3.4.2 О решении серии задач дифракции на различных телах 75
3.4.3 Решение задачи дифракции на магнитоэлектрическом теле 77
3.5 О выборе кубатурной формулы 78
3.5.1 Кубатурная формула прямоугольников 78
3.6 Численное тестирование метода выделения особенности функции Грина 82
3.7 Метод решения задачи на многопроцессорных системах 84
3.7.1 Численные расчеты по методу Галеркина 86
3.7.2 Статистика расчетов на многопроцессорных кластерах 87
Приложение 92
- Тензорная функция Грина прямоугольного резонатора
- Анализ интегральных уравнений
- Построение схемы Галеркина
- Метод решения задачи на многопроцессорных системах
Введение к работе
Задачи дифракции электромагнитных волн на рассеивателях с конечной проводимостью имеют широкий диапазон применений в физике плазмы, микроволновой технике, медицине и в других областях. В качестве примера можно привести задачи проектирования микроволновых печей (как бытовых, так и больших промышленных), обладающих заданными параметрами. Задача дифракции на неоднородном теле произвольной формы имеет тесную связь с вопросами о влиянии электромагнитного излучения на биологические объекты. Решения задач дифракции на анизотропных телах могут быть использованы при исследовании характеристик плазменных образований. В связи с этим важное значение приобретает разработка математических моделей задач дифракции и эффективных численных методов их решения.
Точное решение задачи, представляющее рассеянное поле во всем пространстве в аналитическом виде, удается получить только для весьма узкого класса объектов и при отсутствии экранов. Один из первых результатов представлен в работе [1], где рассматривается двумерная задача возбуждения диэлектрического цилиндра круглого сечения плоской волной. Решение находится методом разделения переменных и представляет собой ряд по тригонометрическим и цилиндрическим функциям (ряд Рэлея), который сходится при всех значениях koR, где ко — 2тг/Х~ волновое число, А— длина падающей волны, Я— радиус цилиндра. Если k$R « 1, то ряд сходится достаточно быстро. При koR» 1 сходимость ряда может быть улучшена при помощи преобразования Ватсона [2]. В резонансном случае (т.е.
для цилиндров с диаметром порядка длины волны) ряд Рэлея сходит-
ся плохо, и для его суммирования необходимо применять специальные методы [3]. Методом разделения переменных получено решение для диэлектрического шара [4] и для некоторых других тел простой формы. В отдельных случаях этим методом можно получить точные
'" решения и для более сложных конфигураций рассеивателя.
Для большинства задач дифракции получить решение в аналитическом виде не удается, и требуется применение численных методов. Диапазон частот, в котором рассматривается задача, определяет
Ф подбор численного метода для ее решения. В задачах, где характер-
ные размеры много больше длины падающей волны (коротковолновый диапазон), находят успешное применение асимптотические методы [5]. Для резонансного диапазона (когда длина волны и характер-
Ф ные размеры суть величины одного порядка) разработан ряд мето-
дов, учитывающих специфику решаемых задач (поведение парамет
ров среды, геометрия тела) и допускающих эффективную численную
реализацию. Наиболее естественными для решения задач дифракции
щ. рассматриваемого диапазона длин волн являются проекционные ме-
тоды и методы интегральных уравнений.
Проекционные методы типа метода конечных элементов нахо
дят применение для численного решения трехмерных задач дифрак-
* ции. Однако прямое применение метода конечных элементов встре-
чает ряд трудностей. Во-первых, краевая задача для системы уравнений Максвелла не является эллиптической, поэтому "не работают" стандартные схемы доказательства сходимости проекционных методов [15]. Во-вторых, для получения приемлемой точности расчета поля в теле с диэлектрической проницаемостью є и 10-^20єо (например,
если тело, в основном, состоит из воды) необходимо выбирать доста-точно мелкую сетку, что влечет также выбор мелкой сетки и в объеме вне тела (выбор же сетки разного масштаба внутри и вне тела ведет к неверным результатам). А это, в свою очередь, учитывая трехмерный векторный характер задачи, приводит к разреженным матри-
* цам очень больших порядков в методе конечных элементов. В неко-
торых случаях эффективность проекционных методов удается повысить посредством выбора специальной системы аппроксимирующих функций. Так, в некоторых задачах, лучше использовать не тригоно-
4 метрическую систему, а, например, систему локальных сплайнов.
В последнее время метод конечных элементов, а также некоторые другие реализуются в комплексах программ, предназначенных для моделирования микроволновых процессов внутри абсолютно про-
Ф водящих структур и, в частности, в прямоугольных резонаторах —
это учебные ( MAGNA/TDM, EMA3D, JMAG, CONCEPT Я и другие) и промышленные [Agilent HFSS Designer 5.5, ANSYS/EMAG, Microwave Studio 2, EMPIRE 2.2, QuickWave 1.9) пакеты программ.
ве тел внутри прямоугольного резонатора с абсолютно проводящей поверхностью. Резонатор возбуждается монохроматической волной СВЧ-диапазона, источник поля - волновод, соединенный с резонато-
л ром через отверстие. Тело характеризуется большой относительной
диэлектрической проницаемостью (например, рассматривается случай є = 65 + 20г). Список решаемых задач следующий: отыскание собственных частот в присутствие тел с потерями, распределение поля от заданного источника, определение рассеянной энергии поля и коэффициента удельного поглощения. Наиболее распространенными
7 математическими методами для решения этих задач являются метод
конечных элементов, метод конечно-разностных схем, метод моментов. Однако, несмотря на достаточно высокую цену предлагаемых программных продуктов, результаты, анонсируемые их производителями, сильно разнятся и не могут считаться надежными [6]. По-следнее объясняется, прежде всего, отсутствием до последнего времени интереса у разработчиков ПО к данному кругу задач. Кроме того, указанные пакеты являются лишь модификациями программ для решения задач, связанных с техникой связи, и не имеют хороше-
ф го теоретического обоснования.
До недавнего времени для задач дифракции резонансного диапазона не существовало универсального метода решения всех перечисленных выше задач, приемлемого с точки зрения вычислительных
Ф ресурсов, и для решения более общих задач приходилось вносить в
их постановку различные упрощения, например, считать рассеива-тель телом вращения или заменять непрерывное распределение параметров среды кусочно-постоянным. Впервые объемные интеграль-
гф ные уравнения были выведены в [8]. В работах [9, 10] были приве-
дены объемные сингулярные интегральные уравнения относительно вектора электрического поля, описывающие трехмерные задачи дифракции на неоднородных диэлектрических рассеивателях. К досто-
щ инствам этих уравнений следует отнести простоту и универсальность
учета неоднородности любого типа, т.е. анизотропии.
В работах А.Б. Самохина [11, 12] эти уравнения получили наиболее подробное исследование. Был построен численный метод [13, 14], позволяющий решать за приемлемое время трехмерные за-дачи дифракции на неоднородных диэлектрических телах с харак-
8 терными размерами до нескольких длин волн с использованием персональных компьютеров. Опишем кратко результаты, полученные в этих работах.
Рассмотрен следующий класс задач. Пусть область Q С R3, ограниченная поверхностью S, характеризуется произвольно распределенной тензор-функцией (#), обращающейся вне Q в скалярную константу Єо свободного пространства. Магнитная проницаемость во всем пространстве предполагается постоянной. Требуется определить электромагнитное поле, возбуждаемое в данной среде внешним монохроматическим полем, представляющим собой либо плоскую волну, либо поле заданного стороннего тока. Вводя электрический ток поляризации, применяя известные формулы векторного потенциала и теорему о дифференцировании интегралов со слабой особенностью [22], для рассматриваемой задачи получено [11] следующее векторное интегральное уравнение относительно электрического поля:
-И>.р./ ( I ^ - /) E(y),grad х) grad xG{r)dy,
где G(r) = exp(ikor)/(4nr) - фундаментальное решение трехмерного уравнения Гельмгольца в свободном пространстве, соответствующее временной зависимости ехр(—iwt). При достаточно общих предположениях о компонентах тензора проницаемости, оператор уравнения является ограниченным в пространстве jLa(Q), и в [12] получена оценка для его нормы.
Уравнение (1) относится к классу многомерных сингулярных интегральных уравнений, общая теория которых построена в [22]. Ос-
новные результаты этой теории получены для уравнений на многообразиях без края, поэтому для того, чтобы применить их к (1), необходимо предварительно определить продолжение оператора уравнения в пространство ^(R3). В [11] такое продолжение построено, исследована его взаимосвязь с исходным оператором, и в явном виде выписан символ продолженного оператора. Установлен следующий критерий нетеровости оператора (1).
Теорема В1. Пусть компоненты тензора диэлектрической проницаемости непрерывны во всем пространстве. Тогда для того, чтобы оператор уравнения (1) был нетеров в L-2,{Q), необходимо и достаточно выполнения условия
х,в,ф
> 0, (2)
он(8,ф)ац(Є,ф)єу{х)
где <Хі{6,ф)— декартовы координаты точек единичной сферы, у— компоненты тензора є в декартовой системе координат.
Отметим, что требование непрерывности компонент тензора диэлектрической проницаемости во всем пространстве является принципиальным при использовании теории многомерных сингулярных интегральных уравнений, поскольку эти функции входят в выражения для символа оператора (1). С использованием теоремы В1 и факта эквивалентности исходной дифференциальной задачи и системы (1), в [11] получена следующая теорема существования и единственности решения задачи дифракции на неоднородном диэлектрике.
Теорема В2. Пусть декартовые составляющие исходного поля и компоненты тензора диэлектрической проницаемости непрерывны по Гельдеру во всем пространстве, выполняется условие (2), а тензор ^((#) — е*{х]) положительно определен почти всюду
10 в Q. Тогда решение системы (1) существует, единственно в L2{Q)
и является классическим решением задачи дифракции, т.е. удовлетворяет поточечно уравнениям Максвелла и условиям излучения.
Для изотропной среды требование непрерывности диэлектри
ческой проницаемости во всем пространстве удается снять. В [12] для
* системы (1) в явном виде построен регуляризующий оператор, и те-
орема существования и единственности решения задачи дифракции доказана при выполнении неравенства
Ітє(х) > Ітє0. (3)
В [11] также проведено исследование системы (1) при минимальном требовании ограниченности компонент тензора диэлектрической проницаемости. При условии
(є(х) - е*(х) - 2%1гп во I) /(2%) > 0, х Є Q,
обобщающем (3) на анизотропный случай, доказана теорема существования и единственности решения (1) в Li2{Q).
Подробное исследование задачи дифракции в свободном про-
щ странстве на сложных диэлектриках с применением метода интег-
ральных уравнений проведено А.С. Ильинским, А.Б. Самохиным и Ю.Ю. Капустиным в работах [30]-[32].
Метод интегральных уравнений находит применение и при ре-
^ шении других классов задач, например - в волноведущих структурах.
В работе Ю.А. Еремина, В.И. Ивахненко [33] рассматривается проблема моделирования царапины на стенке волновода. Исходная дифференциальная задача сводится с использованием тензорной функции Грина к интегральному уравнению, для приближенного ре-шения которого построен и обоснован численный метод.
В работе Карчевского Е.М. [34] метод сингулярных интегральных уравнений успешно применен при исследовании проблемы о собственных модах волновода с размытой границей.
До недавнего времени главной областью применения реальных резонаторов была техника СВЧ. К основным вопросам, рассматривавшимся исследователями, можно отнести следующие: определение собственных частот и волн для резонаторов разной формы; исследование характера волновых процессов в резонаторах; изучение влияния на поля потерь в реальных резонаторах и определение добротности; исследование диэлектрических, оптических и квазиоптических резонаторов. Эти и ряд других вопросов исследованы в монографиях [36], [37], [38].
Наиболее полно математический аппарат для решения задач электродинамики применительно к резонаторам исследован в работах В.В. Никольского [39] - [41].
В [42] рассматриваются следующие задачи.
Исследуется вопрос о свободных колебаниях в полом резонаторе, задача сводится к трехмерному уравнению Гельмгольца. Для прямоугольного и цилиндрического резонаторов выписаны собственные частоты и представлены полные поля собственных колебаний (с.303-310).
Далее автором изучается вопрос о вынужденных колебаниях полого резонатора. Ставится задача найти распределение электромагнитного поля в объеме резонатора при заданных электрических и магнитных источниках (токи и падающая волна). Математическая постановка задачи формулируется в терминах уравнений Максвелла и граничных условий для полей (с. 388). Решение представляется в
12 виде разложения в бесконечный ряд по системам собственных Е- и
Н -полей.
В настоящей работе для исследования задачи дифракции в прямоугольном резонаторе используется метод сингулярных интеграль-
,^, ных уравнений, который используется и для теоретического исследо-
вания задач дифракции, и при построении численных методов. Применение данного метода исследования вызвано спецификой проблем изучаемой задачи. Переход от дифференциальных уравнений к интегральным осуществляется введением токов поляризации с использова-ниєм функции Грина или при помощи методов теории потенциала. Граничные условия для решений обеспечивается свойствами специальным способом выбираемых функций Грина, и их проверка не вызывает затруднений. Таким образом, исследование задачи сводится к изучению свойств ядер и символов уравнении с привлечением существующей теории, что является несомненным преимуществом метода интегральных уравнений. Решение задач дифракции численными методами сводится к решению систем линейных алгебраических уравнений. Для методов, основанных на интегральных уравнениях, размерность систем получается меньше: 104 — 105 (для метода конечных элементов - 106 — Ю8), поэтому для решения систем можно использовать даже прямые методы. Кроме того, интегралы вычисляются по области неоднородности, а не по всему объему резонатора, что также можно отнести к преимуществам данного метода. С другой стороны, применение ИУ приводит к заполненным матрицам с более сложным характером их элементов, что существенно сказывается на
ф- времени заполнения матрицы СЛАУ.
13 Данная работа содержит следующие основные результаты
Для задач дифракции на неоднородном анизотропном диэлектрическом теле получено объемное векторное сингулярное интегральное уравнение. Исследована функция Грина прямоугольного ре-зонатора. Установлена эквивалентность исходной дифференциальной и интегральной формулировок.
Для векторного сингулярного уравнения установлен критерий нормальной разрешимости в пространстве интегрируемых с квадратом функций. Доказаны теоремы существования и единствен-
^ ности решения задачи в случаях, когда диэлектрическая проницае-
мость является комплексной или вещественной тензор-функцией.
3. Предложен численный метод решения задачи, пригодный
для решения задач дифракции в СВЧ-диапазоне волн. Доказана схо-
димость метода. Исследована возможность оптимизации алгоритма;
предложен метод решения задачи на многопроцессорных кластерах.
Работа состоит из введения и трех глав. В первой главе рассмотрены дифференциальная и интегральная формулировки задачи дифракции. Во второй главе проводится исследование объемного векторного сингулярного интегрального уравнения, рассмотрен вопрос о существовании и единственности решения задачи. Третья глава работы посвящена численному методу решения задачи, его точной
'- формулировке и обоснованию, вопросам оптимизации вычислений и
реализации метода для многопроцессорных вычислительных систем. В работе принята сквозная нумерация формул.
В первой главе задача дифракции сформулирована в диф-
л ференциальной и интегральной форме, установлена эквивалентность
14 этих формулировок.
В первом пункте проведен формальный вывод интегро-дифференциального уравнения для векторов электрического и магнитного полей с применением теории векторного потенциала и с применением тензорной функции Грина. Сформулировано определение классического решения и получены достаточные условия, при которых любое решение интегро-дифференциального уравнения является классическим решением задачи в дифференциальной формулировке.
В пункте 2 выписана тензорная функция Грина полого прямоугольного резонатора, исследованы ее свойства. Предложен конструктивный метод выделения ее особенности. Доказано представление функции Грина в виде суммы сингулярного и гладкого во всем объеме резонатора слагаемого.
В третьем пункте для решений интегро-дифференциальных уравнений проверяется выполнение условий на стенках резонатора и на границе разрыва сред. Устанавливается эквивалентность интегро-дифференциальной и исходной формулировок задачи дифракции.
В четвертом пункте из интегро-дифференциального уравнения получена система объемных сингулярных интегральных уравнений.
Вторая глава посвящена исследованию объемного интегрального уравнения для электрического поля в пространстве L\{Q), основанному на теории многомерных сингулярных интегральных уравнений [22]. Построено продолжение оператора в пространство Z|(R3). Доказано, что исходное уравнение и уравнение с новым оператором одновременно являются нормально разрешимыми.
В первом пункте приведены краткие сведения из теории многомерных сингулярных интегральных уравнений: даны определения
15 сингулярного интегрального оператора и его символа. Сформулиро-
вано определение фредгольмовости интегрального уравнения и сис
темы уравнений. Приведены несколько достаточных признаков фред
гольмовости сингулярного оператора, основанных на свойствах сим
вола.
* В пунктах 2 и 3 рассматриваются классические и обобщенные
решения задачи дифракции на теле с прозрачной границей - тензор
диэлектрической проницаемости является всюду непрерывным. На
основе результатов исследования задачи дифракции в свободном про-
*. странстве 1 доказываются теоремы о фредгольмовости сингулярного
интегрального оператора.
Пункт 4 второй главы посвящен исследованию задачи дифрак
ции на диэлектрике без потерь и с разрывной проницаемостью. Ис-
ф> ходная краевая задача рассматривается в обобщенной постановке [17].
На основе результатов исследования краевой задачи в пространстве
L\[Q)i 2 используя эквивалентность исходной и интегральной фор
мулировок, доказывается теорема о существовании и единственности
ж решения сингулярного интегрального уравнения.
Третья глава посвящена построению и обоснованию числен
ного метода для отыскания приближенного решения задачи дифрак
ции. Доказывается теорема о сходимости метода Галеркина для ре-
- шения интегрального уравнения. Выведено достаточное условия схо-
димости метода, накладывающее ограничение не тензор диэлектрической проницаемости.
'А.Б, Самохик. Интегральные уравнения и итерационные методы в электромагнитном рассеянии. - М.: Радио и Связь, 1998.
2М.Ш. Бирман, М.З. Соломяк. .2-теория оператора Максвелла в произвольных областях. //
ф_ Успехи математических наук. - 1987. - Т.42. - Вып.6. - С.61-75.
В первом пункте приведены общие сведения из теории проек-ционных методов: сформулированы определение сходимости метода Галеркина для операторного уравнения и теорема о сходимости метода.
Пункт 2 посвящен исследованию метода Галеркина для объем-
* ного сингулярного интегрального уравнения задачи дифракции. До-
казывается теорема о сходимости метода Галеркина для решения ин
тегрального уравнения. Выведено достаточное условия сходимости
метода, накладывающее ограничение не тензор диэлектрической про-
л ницаемости.
В пункте 3 рассматривается конкретная реализация метода Га
леркина в случае, когда тело представляет собой анизотропный неод
нородный параллелепипед. Для базисных функций выписаны точные
1^1 выражения, исследованы их свойства. Для отыскания приближенных
решений метод сводится к решению системы линейных алгебраичес
ких уравнений; даны точные формулы для элементов расширенной
матрицы этой системы.
ш Четвертый пункт главы 3 посвящен некоторым аспектам про-
граммной реализации численного метода: рассмотрен вопрос опти
мизации вычислений за счет симметрии тензорной функции Грина;
исследована возможность быстрого решения серии задач дифракции
ш, на телах сложной формы и различными значениями диэлектрической
проницаемости, а также задач дифракции на магнитных телах.
В пункте 5 исследуется проблема приближенного вычисления шестимерных интегралов в методе Галеркина. Предложено использование квадратурной формулы прямоугольников. Исследуется вопрос о точности вычисления интегралов, а также о влиянии числа узлов
квадратуры на сходимость численного метода.
В шестом пункте третьей главы приведены результаты численного тестировании конструктивного метода выделения особенности функции Грина. Исследована область применимости исходных формул для компонент функции Грина. Показано, что предложенный конструктивный метод эффективен при вычислении значений функции Грина, особенно вблизи особых точек, где применение исходных представлений вообще невозможно.
Пункт 7 третьей главы посвящен параллельному методу решения задачи дифракции. Предложена эффективная реализация численного метода на многопроцессорных системах. Получены новые численные результаты решения задачи дифракции в резонаторе. Приведены статистические данные о работе кластеров НИВЦ МГУ, участвовавших в вычислительном эксперименте, подтверждающие высокую эффективность предложенного параллельного алгоритма.
Тензорная функция Грина прямоугольного резонатора
Одним из ключевых моментов в исследовании уравнений (14) явля-ется изучение свойств тензора GE. ДЛЯ ограниченных пространств и, в частности, для прямоугольных резонаторов, компоненты функции Грина являются решением уравнения Гельмгольца и удовлетворяют следующим граничным условиям, обеспечивающим исчезновение тангенциальных составляющих напряженности электрического поля на стенках резонатора: д В [26] показано, что при произвольном распределении источников тензор GE - диагональный, а его компоненты имеют вид В этих выражениях 7шп = у () + () — Аг, при этом квадратного корня выбирается так, чтобы 1тупт 0. (см. [21], стр.55). Ряды (16) сходятся экспоненциально, если тз ф уз- Если х% = уз» то ряды (16) расходятся даже при вычислении в точке регуляр ности функции Грина. Можно поступить следующим образом. Пусть, например, х\ ф у\. Перейдем к представлению функции Грина в ви де тройных рядов по тригонометрическим функциям, а затем свер нем ряд по индексу, соответствующему переменным xi,ух. Новые Ф ряды будут содержать множители вида сіг/шл і ch7nm(c — yi) или sh7nmxl зп7пт(с — 1)) и будут экспоненциально сходится. В случае, когда х — у\ -э- 0 ряды требуют дополнительного исследования. Теорема 1.1. Тензор Грина GE допускает представление где матрица-функция (тензор) g Є C(Q X Р) и g Є С(Р X Q). Доказательство. Запишем GlE с выделенной особенностью в при х = у: где функция 0і Є C(Q X Р) (см. [21], стр. 132). Аналогичные пред- ставлення справедливы и для функций G%, GE. Отсюда, и в силу симметрии функций Грина G (x y) — G1g(y,x), (m = 1,2,3) получим утверждение теоремы. Такое представление функции Грина удобно для теоретического исследования математической модели задачи дифракции, но не-пригодно для численных расчетов, так как оно не дает явного пред- ставлення гладкого остатка д(х,у). Ниже на примере компоненты G1 будет изложен конструктивный метод выделения особенности, позволяющий корректно вычислять значения функции Грина вблизи особых точек. Пусть для функции / по узлам /ц(п + I), h2{m -f j) (l.j — —1,0,1), rnm := f — Lnm. Выражения для Сц имеют следующий вид: где при больших п,т, поскольку rnm() выражается через третьи производные функции /, имеющие асимптотику 0( 4) при — оо. Оценки равномерны по х. Просуммировав интегралы 1пт, получим где сумма ряда понимается как предел "сферических" частичных сумм f, Бл - круг радиуса i с центром в нуле.
Таким образом имеем Из последнего соотношения находим Выражая в (21) экспоненциальные функции через sin и cos, придем к следющему представлению ряда Н Вернемся теперь к выражениям (16). Проведя ряд преобразований, получим: и разобьем в соответствии с этим представлением ряд (4) на два слагаемых и X)" При больших п,тп функция /1 имеет асимптотику e 5lnm, где S $а 0. Это свойство, обеспечивающее равномерную сходимость ряда , нарушается лишь в случаях: хз + уз = 2с (для первого слагаемого в квадратных скобках) и жз + Уз = 0 (для второго слагаемого). Однако эти случаи можно исключить, так как мы предположили, что тело в резонаторе не касается его стенок. Рассмотрим теперь сумму Воспользовавшись формулами для произведения косинусов и синусов, получим Учитывая полученные соотношения, получим выражение для ?, с выделенной особенностью: где функция дг непрерывно дифференцируема в Р. Аналогично используем описанный выше метод и для выделения особенностей компонент G\ и G\ тензорной функции Грина. Заметим, что приведенные выше рассуждения не могут служить обоснованием того, что любое решение исходной дифференциальной задачи (4) будет также решением системы (14). Действительно, при выводе представлений (14) использовалось предположение о том, что потенциалы удовлетворяют условию (12). Условие /\ф = iu Q ф, обычно накладываемое на произвольную функцию ф из (10), не гарантирует существования требуемой пары потенциалов, а лишь утверждает, что если найдется хотя бы одна такая пара, то и все другие потенциалы, получаемые при помощи формулы (10), будут удовлетворять условию (12). Покажем обратное, т.е. что любое решение интегро-дифференциального уравнения будет также и решением задачи (4). Предварительно уточним, какими свойствами гладкости обладают рассматриваемые решения. Для операции ротора, входящей в уравнения Максвелла, можно дать несколько определений. Пусть А(х) - некоторое векторное поле. Определение rot A = [Vi- ] предполагает существование в рассматриваемой точке первых производных. Определение ротора, основанное на формуле Стокса не содержит в своей записи производных, однако можно указать непрерывную по Гельдеру функцию, для которой ротор, определяемый по формуле (26), не существует. Например, рассмотрим поле А = Є2А2, где А2(хі,Х2,хз) = я? (0 р 1) при х\ 0 и А2 = 0 при х\ 0. В силу неравенства №—а? {Ь—а)р, справедливого для любых Ь а 0, находим, что \А(х") - где G - произвольная ограниченная область. Взяв в качестве V в (26) шар с центром в начале координат, получим, что Тем не менее, формула (26) позволяет распространить определение ротора на несколько более широкий класс, чем дифференцируемые в заданной точке вектор-функции.
Пусть векторное поле представляет собой градиент функции, непрерывно дифференцируемой в окрестности рассматриваемой точки хо : А(х) = grad у(х). Возьмем в формуле (26) область V в виде шара с центром в хо и поместим в эту точку начало прямоугольной декартовой и сферической систем координат. Запишем выражение для ротора по формуле (26): Для вычисления последнего интеграла используем представление градиента в сферической системе координат и выражение декартовых координат вектора через сферические координаты и выполнив необходимые преобразования, находим, что декартовы «% компоненты ротора равны нулю, что согласуется с известным тож- деством rot grad 7 — 0 для дважды дифференцируемых скалярных функций. Заметим, что для дифференцируемых векторных полей все перечисленные выше определения ротора эквивалентны. Понимая под ротором векторного поля оператор, определяемый формулой (26), вве-дем следующее определение классического решения. Определение 1.1. Пусть компоненты тензора диэлектрической проницаемости - непрерывные по Гельдеру функции в Q, ток JE непрерывная по Гельдеру в Р финитная вектор-функция. Будем называть классическим решением дифференциальной задачи (4) поле {Е,#}, представимое в Q и P\dQ в виде суммы дифференцируемого поля и градиента непрерывно дифференцируемой функции, удовлетворяющие уравнениям Максвелла (4) при х Є Р \ 8Q, условиям сопряжения и условиям на поверхности резонатора (5). В качестве множества допустимых функций, входящих в опре деление классического решения дифференциальных задач, обычно выбирают класс Ср, где р - максимальный порядок дифференциаль- ных операторов, используемых в постановке задачи. Приведенное вы ше определение классического решения допускает функции, которые не являются, вообще говоря, дифференцируемыми в каждой точке рассматриваемых множеств. Такое расширение множества допусти- jt мых решений обусловлено двумя причинами. Как будет показано ни- же, в задаче о возбуждении электромагнитного поля заданными токами решение не является дифференцируемым во всем объеме, если компоненты стороннего тока непрерывны по Гельдеру (но не дифференцируемы) в области своей локализации. Кроме того, при исследовании задачи с использованием сингулярных интегральных уравнений, нельзя гарантировать дифференцируемости получающихся решений даже в случае гладкого распределения параметров неоднородности. Можно сделать предположение, что условие представимости вектор-функции в виде суммы дифференцируемого поля и непрерывного градиента вьщеляет максимально широкий класс функций, для которых уравнения Максвелла выполняются поточечно.
Анализ интегральных уравнений
Гильбертово пространство трехмерных векторов будем обозначать R. В трехмерной области Q, которая может совпадать с пространством R3, зададим вектор-функции U(x), х Є Q и определим функциональное пространство вектор-функций Щ{0) со скалярным произведением Построенное пространство является гильбертовым пространством интегрируемых с квадратом трехмерных вектор-функций. Это про- " странство будет основным при исследовании уравнений (39). Перепишем уравнение (39) в виде Здесь г/0) т), Ді, І?2, Яз - матричные операторы в R , определяемые формулами Здесь Г (ж, J/), Гі(х, у), Г2(х, у) - операторы, определенные в (40), удовлетворяющие соотношениям Теорема 2.1. Интегральный оператор, отвечающий оператору Г2(ж,у), компактен в L\{Q) Доказательство. Для компактности интегрального операто-ра с ядром К(х, у) в пространстве Щ (Q) достаточно выполнение уело- Для исследуемого оператора это условие выполняется в силу непрерывности тензора диэлектрической проницаемости и гладкости оператора І2(а;, у) (см. теорему 1.1). Утверждение доказано. Далее будем изучать уравнение (59), опираясь на исследования оператора интегрального уравнения без последнего компактного слагаемого (будем называть его основным оператором задачи): Q Основной оператор удобен для исследования - как будет показано ниже, матричные операторы Г(х,у) и Гі(ж, у) симметричны, и для них удобно выписываются сопряженные операторы. Пусть /3 - матричный оператор в пространстве т—мерных векторов Rm. Через (г обозначим оператор в Rm, определенный соотношением где (3 пі, pni - компоненты матричных операторов /г и /3 в декартовой системе координат. Имеет место равенство Здесь / - компоненты матричного оператора /? , т.е. оператора, со-пряженного к (3 в Rm; черта сверху в (65) обозначает комплексное сопряжение. где U - произвольный элемент R3. Для исследования основного оператора векторного сингулярного уравнения будем использовать теорию разрешимости уравнений на многообразиях без края. Рассмотрим следующее уравнение, определенное в Щ{8?) Здесь p(x) - дифференцируемая финитная функция в L(R3), такая, что р(х) = 1, х Є Q.
Эту функцию необходимо ввести для того, чтобы первый из интегральных операторов в (68) был компактным в Из структуры уравнений (63) и (68) видно, что, если т){х) = 0 при х Є R3 \ Q, то (при условии существования) решения этих уравнений совпадают в области Q. Следовательно, уравнение (68) - одно из возможных способов продолжения (63) из пространства L\(Q) в 3(и»). Обозначим через Во и В операторы уравнений (63) и (68) со- Лемма 2.1. Если f}(x) = 0 при х Є R3 \ Q, то из кетеровосттш оператора уравнения (68) в Z(R3) следует нетеровостъ оператора уравнения (63) в L\{Q) и обратно. Доказательство. Учитывая, что г] (х) = 0 при х Є R3 \ Q, из структуры операторов В и В получаем, что множество нулей Кег(В ) оператора В состоит из тех элементов ф, связанных с Из (73) видно, что размерность подпространств нулей Кег(В ) и Кег(В ) одинакова. Из (69) следует связь между нулями операторов BQ И В, откуда выводим, что размерности подпространств Кег(В) и Кег(Во) также совпадают- Таким образом IndB — ITMIBQ. Пусть оператор В нормально разрешим в (R3). Рассмотрим уравнение BQV = VQ, причем VQ J- Кег(В ). Определим Щ Є 3(R3) следующим образом: Тогда, учитывая (73), получаем, что Uo JL Кег(В ) и, следовательно, уравнение BU = UQ разрешимо. Представим решение V в виде U = V + W, где W(x) = 0 при х Є R3 \ Q, a V{x) - 0 при х Є Q. Тогда из структуры операторов В, Во, учитывая, что 7}(х) = О при х Є R3 \ Q, получаем, что W решение уравнения BQW VWQ, Значит оператор BQ нормально разрешим Пусть теперь оператор BQ нормально разрешим в Щ{Я). Рас- Л смотрим уравнение BU = UQ, причем Щ X. Кег(В ). Представим UQ в виде щ = VQ + Wo, где Wo(:c) = 0 при х є R \ Q, a Vft(x) = 0 при х Є Q. Тогда, учитывая (73), получаем, что Wo _1_ Кег(Во) и значит уравнение BQW = VWQ разрешимо. Из структуры операторов BQ и В ясно, что решение BU = UQ определяется формулой Теорема 2.2. Пусть заданы непрерывная в R3 тензор-функция (х), отличная от постоянного тензора є І только в области Q. Для того, чтобы оператор уравнения (63) был нетеров в Щ{Я)у необходимо и достаточно выполнение условия Ё an{e,tp)am(6, ір)єпт(х) где ап(в, р) - декартовы координаты точек единичной сферы, определенные в (29). Доказательство. Найдем матричный символ оператора В, действующего в Щ{кг). Для этого определим сначала матричный символ оператора В . Последний представим в виде суммы Вс + i?s, где Вс - компактный оператор в Щ(К3), & Bs - сингулярная часть В , которая определяется первым и последним членами в правой части (70) и имеет следующую структуру: Из (78), (79) видно, что характеристики операторов Bfn выра жаются только через трехмерные сферические функции второй сте- пени. Поэтому для определения символа нужно, согласно теории син гулярных операторов, соответствующую характеристику умножить на коэффициент 72,3 = —47г/3. Таким образом, символы операторов Bfn определяются формулами Введем следующие ограничения. Пусть тензорная функция ё{х)) отличная от постоянного тензора 6QI только в области Q, всю- 4t ДУ непрерывна, в том числе при переходе через границу 0Q. Тогда из (80) и теоремы С5 находим, что матричный символ оператора В представляется в виде где элементы матрицы определяются формулами Детерминант det Ф можно записать в виде Здесь суммирование проводится по всем повторяющимся индексам. Пусть ж — любая перестановка на множестве целых чисел (1,2,3). Тогда из соображений симметрии следует, что если коэффи-циенты CijUt Qjkiml равны нулю при некоторых значениях индексов, то и коэффициенты Cx(ijki)i Ci:{ijkimi) также равны нулю.
Принимая это во внимание и проводя соответствующие выкладки с учетом (79), можно показать, что в формуле (83) коэффициенты Cijkh Cijkiml равны нулю при всевозможных значениях индексов. Окончательно получим Теперь, используя теорему (С6) для оператора В} а также (84) и лемму 2.1, получим утверждение теоремы. Если тело Q изотропно, то есть є(х) = є(х)І, то условие (76) принимает более простой вид Это условие можно записать в векторной форме: для любого ненулевого вектора ft с действительными компонентами и любого х Є Q. Выражения (86) удобны, если тензор диэлектрической проницаемости задан в криволинейной системе координат. Из теоремы 2.1 следует, что при нарушении условий (76)-(86) рассматриваемый оператор не является нетеровым, а значит не имеет замкнутую область значений и, следовательно, не является корректно разрешимым. Для задач дифракции на изотропном локально-неоднородном теле можно получить условия фредгольмовости оператора уравнения (59). А именно, верно утверждение: Теорема 2.3. Пусть задана всюду непрерывная скалярная +. функция є(х), отличная от постоянной є$ только в области Q. Тог- да оператор уравнения (59) будет фредгольмовым в Ь\{0] если Доказательство. Согласно замечанию к теореме С6, оператор В, определяемый (69), будет фредгольмовым, если нижние грани модулей миноров матрицы Ф строго положительны: Последнее из условий (88) выполняется, так как имеет место (87), а значит и (85). Учитывая (82) и то, что тп(х) = є(х)бтп, получим, что первые два условия (88) приводится к следующему выражению где —1 а 1. Ясно, что, если имеет место (87), то (89) выполняется для всех х Є Q. Далее из леммы 2.1 следует утверждение теоремы. Полученные условия фредгольмовости более жесткие, чем условия нетеровости (85). В частности, (87) не выполняется для плазмы с отрицательными значениями диэлектрической проницаемости. Отметим, что приведенные условия являются только достаточными. Теорема 2.4. Пусть компоненты тензора диэлектрической проницаемости всюду непрерывны, є = 6QI вне Q и выполнено условие (76). Кроме того, Іт Єо 0, Irn Q 0, и тензор (є — є )/(2і) положительно определен почти всюду в Q. Тогда оператор интегрального уравнения (59) фредголъмов в Щ(Я), и его решение является классическим решением исходной задачи дифракции. Доказательство. Согласно [25], для решений однородных уравнений Максвелла в свободном пространстве, удовлетворяющих условиям излучения, справедливо следующее интегральное тождество где Ф 0 имеет физический смысл потока энергии на бесконечности и обращается в ноль только при наличии поглощения в свободном пространстве, т.е. если Imepsa 0 или Іт[і$ 0.
Построение схемы Галеркина
Вернемся теперь к вопросу о конкретной реализации схемы Галеркина для рассматриваемой задачи дифракции. Будем формулировать метод не для сингулярного интегрального уравнения (39), а для интегро-дифференциального уравнения (14). Этот подход оказывается эффективным в силу более удобного представления интегралов. Будем предполагать, что матрица 1 " — -Л обратима в Q / В1 _ /J Loo(Q)j Ї - единичная матрица. Введя обозначения: перейдем от (14) к следующему уравнению Представим это уравнение в виде системы трех скалярных уравне- нии: Определим компоненты приближенного решения J следующим образом: где fl - базисные функции-" крышки", существенно зависящие лишь ОТ переменной Х{. Ниже проводится построение функций /jj. Будем считать, что Q - параллелепипед: Q — {х : а\ Функции /ыт /fm, зависящие от переменных Х2 и х% соответственно, определяются аналогичными соотношениями. Так как (125) то каждая компонента вектора приближенного решения обращается в нуль на двух гранях параллелепипеда Q (такой выбор базисных функций позволяет избежать возникновения поверхностных интегралов). Отметим, что построенное множество базисных функций удовлетворяет требуемому условию аппроксимации в L% [24]. Введем сквозную нумерацию для функций Д : Расширенную матрицу для нахождения неизвестных коэффициентов ajt, bk, Ck удобно представить в блочной форме: элементы колонок Bk и матриц Ац определяются из соотношений: Преобразуем второе скалярное произведение в (129). Применяя к внешнему интегралу формулу интегрирования по частям, получим: Поверхностные интегралы отсутствуют в силу условия (125). Введем следующую вспомогательную функцию Производные функции G связаны с производными компонент тензорной функции Грина ( формулами: Используя последние соотношения и применяя ко внутреннему интегралу формулу интегрирования по частям, получим: или, записывая явно скалярные произведения по носителям базисных функций: Примечание - предположение о том, что Q — параллелепипед, сделано лишь для удобства построения схемы метода Галеркина.
Предложенный численный метод может быть реализован и для более сложных структур: можно рассматривать тела с кусочно-гладкой границей (см. определение 2.6 главы 2), многосвязные тела . В этих случаях наибольшую сложность представляет задача перенумерации базисных функций. Кроме того, необходимо следить за тем, чтобы их носители лежали в области неоднородности - в противном случае нарушается условие существования тензора . С целью проверки предложенных численных методов решения задачи дифракции в прямоугольном резонаторе предложенные алгоритмы были запрограммированы на языке программирования Fortran 90. Система программирования - Microsoft Fortran Powerstation 4.0 для операционной системы Windows 95. Программа состоит из следующих основных подпрограмм (перечислены в порядке выполнения): 1. подпрограмма вычисления элементов матрицы СЛАУ и записи их в файлы; 2. подпрограмма заполнения матрицы СЛАУ; 3. подпрограмма решения СЛАУ методом Гаусса; 4. подпрограмма табулирования приближенного решения; 5. подпрограмма построения графиков приближенного решения. Так как заполнение основной матрицы СЛАУ для рассматриваемой задачи является трудоемким, возникает необходимость в оптимизации счета. Довольно существенно вычисления можно сократить, учитывая характер функции Грина, присутствующей в выражениях для интегралов. А именно, для функции Грина верны следующие соотношения: Рассмотрим второе скалярное произведение в (133): Введя обратно (для удобства записи) тройную нумерацию индексов для базисных функций с учетом свойств (134) получим: Подобные же равенства верны и для измененных вторых и третьих индексов. Таким образом достигается уменьшение времени вычислений в 8 раз. Для второго скалярного произведения в (133) рассуждения аналогично повторяются. Решение задачи дифракции методом интегральных уравнений с применением численного метода Галеркина является трудоемкой задачи, однако предложенный алгоритм можно эффективно применять даже для решения целой серии задач на телах различной формы и харак-теризующейся разными значениями диэлектрической проницаемости, что можно отнести к существенным достоинствам метода интегральных уравнений.. Элементы расширенной матрицы представляют собой сумму нескольких интегралов: трехмерного и шестимерных, причем вычис- ление последних - наиболее трудоемкая задача. В первой из перечисленных выше подпрограмм вычисляются второе и третье скалярные произведения формулы (133), не зависящие от тензора (то есть - от электрических параметров области неоднородности). Первое же скалярное произведение представляет собой тройной интеграл, и его вычисление проводится на этапе заполнения матрицы СЛАУ. Важно отметить, что вычисление выражений (137) на несколько порядков быстрее счета записываемых в файл интегралов. Это позволяет быстро решать несколько задач для случаев областей неоднородности с различными электри- ческими характеристиками (при неизменной геометрии этих областей). Описанный принцип можно применять и для решения задач дифракции на телах с различной геометрией. Пусть геометрические параметры резонатора Р выбраны. Вве- дем прямоугольную область Q С Р, граница которой расположена очень близко к стенкам резонатора. Выберем очень мелкую сетку, введем на Q базисные функции и вычислим и запишем в файл все шестимерные интегралы в соответствии со схемой Галеркина.
Метод решения задачи на многопроцессорных системах
Для получения численных результатов решения задачи дифракции эффективным является применение методов параллельных вычисле- щ ний на многопроцессорных кластерах и распределенных вычисли- тельных системах. Неизбежность использования подобных алгоритмов в электродинамике вызвана большим объемом вычислительной работы. , Например, при решении краевой задачи методом конечных эле- ментов проблема сводится к решению систем линейных алгебраических уравнений высокого порядка (106 — 108). Матрицы в данном случае получаются разреженными, а их элементы имеют простой вид и щ во многих случаях могут быть вычислены аналитически. В данном случае очень трудоемким является решение СЛАУ. При решении задачи методом Галеркина наибольшую трудность представляет не решение системы, а ее заполнение. Это обусловлено невысоким (в сравнении с методом конечных элементов) порядком системы уравнений, но более сложным характером элементов матрицы СЛАУ, которые представляются через шестимерные интегралы; в выражениях для этих интегралов присутствуют компоненты функции Грина, которые представляются в виде двойных рядов (в случае свободного пространства, когда функция Грина имеет простой вид, интегралы вычисляются значительно быстрее). Время заполнения основной матрицы СЛАУ составляет Эбиспользования методов параллельных вычислений. Решение же СЛАУ может быть проведено на обычной однопроцессорной системе. Поясним суть алгоритма параллельных вычислений. Предположим, что расширенная матрица СЛАУ состоит из п строк и в нашем распоряжении р процессоров. Распределим задание на заполнение матрицы следующим образом: 1-ый процессор заполняет строки с номерами 1 ... \Щ 2-ой процессор заполняет строки с номерами ИМ + 1 ...2-:4 р— 1-ый процессор заполняет строки (р—2) +1 ... (р— 1) После того как заполнение матрицы завершено, вызывается процедура SCALAPAK. Если процедура SCLAPAK недоступна, поступаем следующим образом: по завершении заполнения матрицы всеми процессорами расширенная матрица полностью "собирается" и передается одному из процессоров, на котором СЛАУ решается одним из распространенных методов. Предложенный алгоритм был запрограммирован на языке Fortran 90, для распараллеливания алгоритма использовалась систе- ма параллельного программирования Message Passing Interface [27], [28] в реализации MPI Chameleon для операционной системы Linux.
Отладка и запуск программы проводились на двух платформах: 1. 4-х процессорный кластер в Региональном суперкомпьютерном центре Пензенского государственного университета. Узлы на ос-нове процессора Alpha21264A, 667MHz, сеть Myrinet 2 Gbit/s. 2. Кластер рабочих станций НИВЦ МГУ ("sky"). Узлы на основе Pentium-Ill, 850 MHz, коммутатор Fast Ethernet. В качестве примера использования метода интегральных уравнений и отыскания его приближенных решений методом Галеркина рассмотрена следующая задача. В кубическом резонаторе со стороной 1 расположен однородный диэлектрический куб [0.2; 0.б]3 с комплексной проницаемостью є — (6,0.1). Волновое число ко = 1. Компоненты решения однородной задачи Е определены соотношениями: На рисунках A. 1 - A. 10 приложения А приведены графические представления модуля приближенного решения задачи в сечениях области неоднородности: на каждом из рисунков переменные х\ и х% пробегают значения от 0.2 до 0.6, а значения переменной гпз фиксированы (эти значения приведены в примечаниях к рисункам). Число базисных функций в методе Галеркина равно 336. Описанный алгоритм параллельных вычислений был програм- мно реализован, с его помощью была решена описанная задача дифракции. При этом было увеличено число базисных функций (общее их число составило 1911) и была повышена точность вычисления интегралов расширенной матрицы системы линейных уравнений. В приложении А приведены графики, представляющие модули приближенных решений в сечениях области неоднородности. Темные участки на рисунках соответствуют малым значениям модуля напряженности \Е\, светлые - большим значениям. На каждом из рисунков переменные х\ и а?2 пробегают значения от 0,2 до 0.6, а значения пе- ременной 3 фиксированы (эти значения приведены в примечаниях к рисункам). Приведенные в рисунках графические данные показывают, что ф полученные приближенные решения в достаточной степени симмет- ричны.
Кроме того, амплитуда напряженности электрического поля уменьшается внутри области неоднородности. Такое поведение решения имеет физический смысл, так как решается задача дифракции на , теле с поглощением (диэлектрическая проницаемость имеет ненуле- вую мнимую часть), и в этом случае имеет место небольшое затухание поля внутри области неоднородности. 3.7.2 Статистика расчетов на многопроцессорных кластерах В численном эксперименте участвовали 4 вычислительных кластера НИВЦ МГУ AQUA, LEO, SCI и SKY, объединенных в единую сеть. "Нагрузка" процессоров заданиями осуществлялась следую щим образом. Если процессор не был занят выполнением "посторон- них" (например, системных) заданий, то ему передавалась команда на вычисление элемента расширенной матрицы с номером /, J По окончании вычисления элемента процессор производил запись элемента в формате / J A(I, J) в общий файл и получал следующее задание. Ниже приведена общая статистика, описывающая работу кластеров при решении задачи дифракции. Общие данные: Время эксперимента: 643 часа 26 дней 19 часов Обработано заданий: 3 653 832 Суммарное время работы узлов: 34 582 часов 4 года Суммарное время расчета заданий: 34 115 часов 4 года Эффективность передачи данных: 98.64% Примечание. Эффективность передачи данных есть отношение времени на вычисление элементов матрицы к общему времени решения задачи на кластере. Приведенное значение говорит о том, что время на пересылку служебной информации и результатов вычислений очень мало в сравнении со временем счета элемента (менее 2-х процентов) . Столь высокую эффективность работы системы можно объяснить тем, что вычисления элементов производятся независимо, общие данные для процессоров отсутствуют, и разным процессорам не приходится обмениваться никакой информацией. Этот момент является очень существенным, так как часто при разработке алгоритмов параллельных вычислений приходится сталкиваться с проблемой обмена данными между процессорами. В рассматриваемой же задаче эта проблема отсутствует. Рачет заданий Минимальное время: 0:00:01 Максимальное время: 0:02:44 Среднее время: 0:00:42 Примечание. Время расчета заданий - время вычисления одного элемента расширенной матрицы на одном процессоре. Сильное различие имеет две причины. Во-первых, элементы правой части матрицы вычисляются быстрее (так как представляются трехмерными интегралами), чем элементы основной матрицы - шестимерные интегралы, зависящие от функции Грина. Во-вторых, кластеры собраны на процессорах разных классов (Pentium III и Pentium IV), что также сказывается на времени вычисления элементов.
