Многосеточный метод решения двумерной задачи дифракции Текст научной статьи по специальности «Математика»

■ъ

ФИЗИКО-МАТЕМАТИЧЕСКИЕ НАУКИ

УДК 317.9:51*6 н.А. ЛУГОВАЯ,

С.А. ТЕРЕНТЬЕВ

Омский государственный университет им. Ф.М. Достоевского

МНОГОСЕТОЧНЫЙ МЕТОД РЕШЕНИЯ ДВУМЕРНОЙ ЗАДАЧИ ДИФРАКЦИИ

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

1. Постановка задачи

В работе рассматривается задача дифракции электромагнитного поля на неоднородном локальном теле. Предполагается, что среда и первичное поле имеют плоскую симметрию, и вдоль оси Ог характеристики среды не меняются. В этом случае задача распадается на две независимые скалярные задачи: для Е. в случае ^-поляризации поля и ддя Н. в случае/-/-поляризации поля [1].

Допустим Ц — двумерная область, содержащая неоднородность и ограниченная кусочно-гладкой кривой Ляпунова Г, — внешняя к Г область. Проводимость ег' и магнитная проницаемость в области П, являются функциями координат (х,у), а проводимость а', и магнитная проницаемость /,!с в области П. предполагаются постоянными. Временная зависимость выбрана в виде е"°' с круговой частотой со.

Вводя обозначения: ;/ = £., у = (л в случае Е-поля-ризации поля и и = Н., у - с' в случае //-поляризации поля, исходную формулировку задачи записываем в виде:

■ Л"„

к-1. ч„

Ди„

к: и»

-V-V», -Г, У,

1 ди„

О

и,. = и,,

К дп

1 ди, У, дп

в а, в а,,

на Г.

(1)

Здесь ие — поле в области , и1 - поле в области П,, и„ - известное первичное поле. Аномальное поле иа =ас-и0 удовлетворяет условиям излучения на бесконечности. Коэффициенты к^ и к] определяются формулами: к] = \(оц{а\, к~ = коцр] ■

Воспользовавшись леммой Лоренца для вспомогательных потенциалов [3] либо формулой Стрэт-тона-Чу [1], приходим к интегро-дифференциаль-ной формулировке исходной задачи:

- V

V -

V

= О

Г,

(2)

К 8и, 2У, &»«

дт 5л„

Г, д"

Здесь пит- нормальное и касательное направления в точке М.аЛцИ г„ - нормальное и касательное направления в точке Ма гладкости контура Г, ,МВ) = ~ Кп(г'Кг)12п ~ Функция Грина для оператора Гельмгольца первого уравнения (1) в Л2, К (■) - функция Макдональда, г — расстояние междуточками М и М0 в евклидовой норме. При совпадении точек М и Ма функция Грина имеет логарифмическую особенность. Интегрирование в (3) следует понимать в смысле главного значения Коши, поскольку ядро (7С,г0) имеет сильную особенность.

Аналогично можно получить несколько иную формулировку задачи, которой удобней воспользоваться при численном решении:

- V —Vu, Г,

— «i = О

Y,

íf

Y¡ дп дип

а- дп„

(4)

(5)

ди, Тт

Здесь Г' — кривая, являющаяся границей замкнутой подобласти области О,, а п„ и г0 - нормальное и касательное направления к ней в точке М0. Несмотря на то, что интегро-дифференциальное граничное условие (5) является уравнением первого рода относительно нормальной производной, задача легко ре-гуляризируется выбором расстояний между контурами Г и Г', а схема вычисления интегралов упрощается. Поле в области Г2е может быть вычислено согласно формуле Грина:

и,(М0) = ив(Ч)"

дп

ке контура. Коэффициенты разложения локально-сглаживающих сплайнов во внутренних узлах были выбраны в виде:

В =--!-у -(-¿V _ 1 у С)

В точках, выходящих за границы отрезков, значения функций и и V экстраполировались многочленами третьей и второй степени соответственно.

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

Точки коллокации выбираются на контуре Г' по направлению внутренней нормали отточек Г*^. В качестве контура Г' использовался прямоугольный контур, расположенный в области Г2, достаточно близко к контуру Г. Выписьшая (5) во всех точках коллокации М,, приходим к системе уравнений:

УВ(у) ч 0А(и) = (8)

где (т^-) - известная сеточная функция, определенная в точках коллокации, а (у) и (и) - искомые сеточные функции, определенные на и Г* соответственно. Матрицы Л и В построены с учетом вышеприведенных формул (7), а элементы матриц V и И/имеют вид:

V= \у,(VмС{М,М,)ММ,))ии]П)<ии, (9)

(10)

у\М) дп '

Интегро-дифференцйальные граничные условия обобщаются и на случай слоистой вмещающей среды [1].

2. Дискретизация задачи

Численный метод решения задачи опишем для формулировки (4)-(5) в предположении, что О, является прямоугольной областью, что не изменит общности постановки. Для простоты изложения переобозначим функции и = м, на и V = на Г. Построим треугольное разбиение области, покрыв ^ прямоугольной сеткой шага И = (А,, Л,), и разбив каждый прямоугольник на два треугольника параллельными диагоналями. Через ф* обозначим множество узлов сетки принадлежащих , через Г11 - множество узлов сетки принадлежащих Г, через Г*: - множество центральных точек отрезков разбиения контура Г .

Для дискретизации интегро-дифференциального условия (5) воспользуемся методом коллокации. На границе области полагаем:

и1г58 £ ~ £ Ачч Ч-ч* (6(

ЫГ" '

где систему координатных функций {тк} строим в виде В-сплайнов третьей степени, экстремальные точки которых принадлежат Гл, а систему {Цы/2} - в виде В-сплайнов второй степени, экстремальные точки которых принадлежат Г*2. Поскольку, в данном случае, контур Г является лишь КС'-гладким контуром, то такие сплайны строятся на каждом С'-отрез-

дт

Для вычисления интегралов по отрезку, ближайшему к точке коллокации, применяются квадратурные формулы Гаусса с достаточно большим четным количеством узлов. Симметричность суммирования относительно точки коллокации обеспечивает достаточно точное вычисление интеграла, близкого к сингулярному. Для вычисления интефалов по остальным отрезкам разбиения применяются квадратурные формулы Гаусса с небольшим количеством узлов.

Для разностной аппроксимации эллиптического уравненияД4) воспользуемся методом конечных элементов [2]. Решение представляем в виде линейного сплайна на треуг ольном разбиении области:

Ф,У) ~ <РкАх'У)- (11)

Здесь <ри - функции, линейные на каждом треугольнике разбиения, равные единице в узле (А,/) и нулю в остальных узлах. Используя вариационную постановку для уравнения (4), приходим к системе линейных алгебраических уравнений:

X "и \ —Ч<РиЧ<Рт,Л<1у -

(t.l)en

У,

~ Е_ ии\—<Ри<Р„Л<*У = )'v^dC. . Y¡ г

{к.1)еО

(12)

—И

(т,п) е П

Поскольку подынтегральные функции в левой части равенства (12) в общем случае могут быть разрывные, то для численного подсчета данных интегралов применяется итерационный метод интегрирования. Идея данного метода состоит в том, что треугольники исходного разбиения делятся на четыре более мелких, тё, в свою очередь, делятся анало-

Погрешность решения задачи дифракции

Таблица 1

1 1 1 Е-поляризация Н-поляризаиия

< т 1 Ли - % на О, Дн\ ■ и ^ % на Г Д ди дп ■ ди дп % на Г / ^ . Ди;, ' % на О. " 2 - и , % на Г Дди дп и;, 2 % на Г . ди дп у^

а Ь а Ь а ь а Ь А Ь а Ь

0 0.073 0.018 0.046 0.011 0.268 0.097 0.073 0.018 0.046 0.011 0.268 0.097

1 0.061 0.015 0.043 0.010 0.235 0.085 0.061 0.015 0.043 0.010 0.235 0.085

2 0.062 0.015 0.052 0.012 0.235 0.085 0.062 0.015 0.052 0.012 0.235 0.085

4 4 3 0.092 0.023 0.079 0.018 0.460 0.168 0.092 0.023 0.079 0.018 0.460 0.168

4 0.122 0.030 0.121 0.029 0.386 0.139 0.122 0.030 0.121 0.029 0.386 0.139

5 0.181 0.043 0.197 0.046 0.754 0.273 0.181 0.043 0.197 0.046 0.754 0.273

6 0.270 0.064 0.295 0.069 0.800 0.289 0.270 0.064 0.295 0.069 0.800 0.289

1 0 0.063 0.016 0.040 0.010 0.095 0.031 2.296 1.170 0.339 0.171 2.324 1.464

1 0.042 0.011 0.026 0.007 0.054 0.018 2.003 0.990 0.190 0.093 0.721 0.483

2 0.015 0.004 0.016 0.004 0.058 0.020 1.651 0.852 0.124 0.059 0.323 0.211

4 16 3 0.054 0.013 0.060 0.014 0.292 0.104 1.234 0.666 0.105 0.041 0.384 0.185

4 0.122 0.030 0.129 0.031 0.381 0.137 0.972 0.535 0.169 0.048 0.402 0.156

5 0.220 0.056 0.219 0.053 0.703 0.254 0.567 0.298 0.240 0.061 0.711 0.259

6 0.321 0.081 0.324 0.079 0.761 0.274 0.493 0.211 0.340 0.084 0.763 0.275

0 0.002 6е-4 0.003 7е-4 0.085 0.027 0.198 0.101 0.006 0.003 1.911 1.362

1 0.002 0.001 0.001 6е-4 0.004 0.002 1.377 0.677 0.121 0.055 0.498 0.337

2 0.013 0.005 0.016 0.005 0.053 0.019 1.541 0.796 0.112 0.054 0.304 0.199

16 64 3 0.051 0.012 0.061 0.015 0.291 0.104 1.213 0.655 0.103 0.039 0.381 0.183

4 0.481 0.132 0.282 0.074 0.381 0.137 1.069 0.547 0.318 0.087 0.402 0.156

5 0.316 0.084 0.251 0.063 0.703 0.254 0.619 0.306 0.280 0.073 0.712 0.259

6 0.725 0.196 0.449 0.114 0.761 0.274 0.874 0.289 0.485 0.124 0.763 0.275

0 0.002 4е-4 0.002 6е-4 0.117 0.041 0.184 0.093 0.071 0.034 2.606 1.714

1 0.002 0.001 0.003 7е-4 0.021 0.007 4.110 2.151 3.408 1.802 5.938 4.043

2 0.010 0.002 0.015 0.003 0.055 0.020 2.555 1.302 1.859 0.942 4.015 2.790

64 16 3 0.052 0.012 0.061 0.015 0.292 0.105 1.876 0.955 1.073 0.541 2.930 1.963

4 0.199 0.052 0.156 0.039 0.382 0.137 1.406 0.683 0.828 0.394 1.649 0.984

5 0.238 0.061 0.225 0.055 0.704 0.254 0.753 0.363 0.474 0.208 1.312 0.684

6 0.379 0.098 0.337 0.083 0.761 0.274 0.583 0.253 0.423 0.151 0.942 0.425

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

При аппроксимации интегралов в правой части (12) применялись квадратурные формулы второго порядка точности с линейной интерполяцией V по узлам из Г';/2. Для внутренних точек отрезков границы эти квадратуры принимают вид:

Г ^

Для угловых точек границы:

+ 7 V><

(13)

(14)

ближенное решение задачи ищется за счет сглаживания невязки и рекурсивного списка на более редкую сетку. В связи с этим через Ф обозначим пространство функций и*, определенных на сетке П . Для простоты изложения предположим, что Л, =2-С, и Л„ = 2~" Су.

В данной работе при квазистационарной постановки задачи (Яе А,2 = 0) в качестве оператора, эффективно гасящего негладкую компоненту невязки, использовался метод Зейделя. Для восстановления значений сеточных функций на более частую сетку испол^зовадся оператор линейной интерполяции Р1' :Ф -> Ф такой, что

м2Л, если / и _/ четны

(ц. ,.1 )

если I нечетно, _/ четно если / четно, 1 нечетно

, если ; и _/ нечетны

Где, в зависимости от (рг, в качестве Л, и используются как /?г, так и А,.. Выражая из системы (8) значения функции V на Г*, через значения и на Гл и подставляя полученные зависимости в правую часть системы (12), приходим к системе линейных алгебраических уравнений с неизвестными значениями функции и на О :>/;< + 2и = /. Здесь 3 - оператор левой части равенства (12), а оператор 2 действует только на граничные точки, оставляя на месте внутренние точки области. Построенная система решается циклически многосеточным методом [2].

3, Многосеточный метод

Многосеточный метод представляет собой итерационный процесс, на каждой итерации которого при-

а при переходе^ бол^е частую сетку использовался оператор ()'' :Ф Ф , сопряженный к Рь, определяемый равенством:

2к /лЛ Л I / Л , V Л ч

=0 ",,=ТК + 1 Щ,) (15)

Систему, построенную на сетке , запишем в виде

I"и"=/"■ (16)

Одна итерация многосеточного метода состоит в следующем.

Сделав к «3-4 итераций методом Зейделя, имеем приближение йк. Из равенства = - /А определяется невязка /•''на сетке <3 . Решив задачу вида ¿Аи>'' = г'1, мы найдем поправку к приближенному решению и' и определим точное решение задачи

Рис. 1. Модель неоднородности для тестовой задачи

1.5

-0.5

0.5

1.5

|о.„Л>| |ДЛг,;с>|

Рис. 5. Поведение и в зависимости

от .V при у = V,, =1.1

Ли \ п,

пе

Рис. 2. Модель неоднородности, на которой демонстрируется сходимость итерационного метода

1 Ь-03

Е-04 :

----

1.Е-05

-1.5 -0.5 0.5

х

Рис. 6. Поведение и -р^- в зависимости

1.5

от х при у = уи = 5

1.Е-02

К1 К.1

Рис. 3. Поведение щ и -^т" в зависимости от х при у = у„ = 1.1

| |л й|„/Дг|

Рис. 7. Поведение и | Д| в зависимости

от X при у = у1} = 5

.Е-02

1.Е-03

1 .Е-04

-1.5

-0.5

03

1.5

\п1„!сх\ ¡Д пг/гУ

Рис. 4. Поведение |л ! и

Сии!с\\ IП1ц1гу\

в зависимости от х при у = у{] — 1

\<'1<Я1<У\ п^/сг]

Рис. 8. Поведение ^ и в зависимости от х при у = _у0 = 5

(16) как и1' =11* -иЛ Поскольку невязка в данном случае является гладкой функцией, то с небольшой погрешностью поправку можно найти, решая аналогичную задачу на более редкой^сетке. Для этого строится оператор ¿2А на сетке . При этом, значения интегралов, составляющих оператор ./, пере-считываются одной итерацией итерационного метода интегрирования, а оператор 2 перестраивается согласно формуле (15). Преобразовав систему к виду ¡}\у2л=0У, (17)

и решив ее приближенно, определяем приближенное значение м'2Л. Взяв новое приближенное решение и1 = и* - /}А)г2л, сглаживаем невязку к итерациями метода Зейделя, с противоположным порядком перебора узлов.

Для решения задачи (17) на сетке О используется та же самая процедура с нулевым нача^ным приближением и вспомогательной сеткой О . Таким образом, рекурсивный спуск продолжается до тех пор, пока число узлов на очередной вспомогательной сетке не станет совсем незначительным, на которой система легко решается прямым методом.

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

К сожалению, итерационный процесс, эффективно гасящий негладкую компоненту невязки системы (8), (12) в целом, в данной работе построен не был. Поэтому, при многократном решении задачи на сетке с п узлами в предположении, что число граничных узлов равно 0(я|/2), трудоемкость алгоритма составляет 0{п + ~), где р - число вариантов возбуждающего поля, а д - число вариантов электромагнитных параметров неоднородности.

4. Вычислительный эксперимент

Вычислительный эксперимент проводился при постоянной магнитной проницаемости среды (//, = /ис) ■ В качестве тестовой задачи рассматривалась задача дифракции на круге 0 = {(х,у)\ (л:-0.!)2+(.у-0.2)2<0.92 Ь размещенном в прямоугольнике П,- = [-1,11 х [-1.5, 1.5] (рис. 1). Проводимость среды была выбрана в виде кусочно-постоянной функции, при этом , при (.г, у) е О

при {х,у)<£В-В качестве первичного поля рассматривались цилиндрические гармоники 1т{-1к/)со$,т{в-ви), где /„(■) — модифицированная функция Бесселя, а полярная ось была направлена под углом #„=-0.2 радиан к оси Ох, и. Точное решение такой задачи строится методом разделения переменных. Численно задача решалась, при расстоянии между соответствующими сторонами контуров Г и Г' равным 0.25 длины отрезков разбиения Г.

В таблице 1 приведены относительные погрешности решения задачи в норме на сетках размернос-тьюа) (25+ 1)х(2'Ч 1)иЬ) (2"+ 1)х(27+ 1) в зависимос-

ти от номера т цилиндрической гармоники и электромагнитных параметров среды Д и Лс.

Кроме того, экспериментально была продемонстрирована сходимость метода в случае неоднородности достаточно произвольной формы. В прямоугольнике С1, = [-0.5, 0.5]х [-1, 1] выделялся сегмент 0 = {(х,у) | х > -0.4, (* + 0.5)2+/<1 } (рис. 2). Проводимость в области была постоянной функцией и к] = I , а в области Г2, проводимость являлась кусочно-постоянной функцией и

/с2, при (,г,)/)гО

при (х,у)ей, х < у при (х,у)е Д х>у В качестве первичного поля рассматривалась плос-

кая волна падающего поля и0(х,у) = е

— 0-'кЛУ-Уа'>

. На ри-

сунках с 3 по 8 приведены результаты для Н-поляри-зации поля и электромагнитных параметрах среды Л,, = 4096, А,, = 16 и Ле = 256. Задача решалась на сетках с количеством узлов (2Л Ч- 1) х (2Ь -Ь 1), (2" + 1)х(27 + 1) и (27 + 1)х(2е + 1). Линия, обозначенная символом V, показывает разницу при вычислении на сетках размерностью (25+ 1) х (2'* -Ь 1) и (2*'+ 1)х(27+ 1), линия, обозначенная символом ' •', -разницу при вычислении на сетках размерностью (2" + 1)х(27+ 1)и(27+ 1)х(2в+ 1). Линия, обозначенная символом ' +', показывает модуль аномального поля, вычисленного на сетке размерностью (27+ 1)х(2в+ 1). Результаты на графиках изображены в логарифмическом масштабе.

Библиографический список

1. В.И. Дмитриев, Е В. Захаров. Интегральные уравнения в краевых задачах электродинамики. - М.: МГУ, 1987. - 167 с.

2. Е.Г. Дьяконов. Минимизация вычислительной работы. -М.: Наука, 1989, - 272 с.

3. С. А. Терентьев. Построение интегральных уравнений для задач дифракции методом вспомогательных потенциалов. ■ В кн.: Численные методы и задачи оптимизации. Томск: ТГУ, 1983, с. 49-58.

4. Н. А. Луговая, С.А. Терентьев. Численное решение двумерной задачи дифракции. // Методы оптимизации и их приложения: Труды XIII Байкальской международной школы-семинара. Том 3. Обратные и некорректные задачи прикладной математики. — Иркутск. 2005, с. 153-158.

5. Н А. Луговая Многосеточный метод решения двумерной задачи дифракции. // Труды II Международной конференции студентов и молодых ученых «Перспективы развития фундаментальных наук». - Томск: ТПУ, 2005, с. 239-241.

ТЕРЕНТЬЕВ Сергей Александрович, доцент, кандидат физ.-мат. наук, кафедра математического моделирования, математический факультет. ЛУГОВАЯ Наталья Александровна, аспирант, кафедра математического моделирования, математический факультет.

Дата поступления статьи в редакцию: 04.08.2006 г. © Терентьев С.А., Луговая H.A.