Алгоритмическое и программное обеспечение для численного решения задач электромагнитного рассеяния на диэлектрике тема диссертации и автореферата по ВАК РФ 05.13.18, кандидат наук Сотникова, Наталья Юрьевна

• Сотникова, Наталья Юрьевна
• кандидат наук
• 2013

Оглавление диссертации кандидат наук Сотникова, Наталья Юрьевна

ОГЛАВЛЕНИЕ

Введение

1. Задача электромагнитного рассеяния на диэлектрике

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

1.2. Расчетная вычислительная сетка

2. Итерационные методы для решения задач электромагнитного рассеяния

2.1. Решение задач методом простой итерации

2.2. Обобщение метода простой итерации

2.3. Оптимальный итерационный параметр для сходимости МОПИ

2.4. Метод релаксации

3. Разработка алгоритмов и программ для решения интегральных уравнений задач рассеяния волн с помощью МОПИ

3.1. Алгоритмы определения оптимального итерационного параметра для различных практически важных случаев области локализации спектра матрицы перехода

3.2. Описание подпрограмм определения оптимального параметра

3.3. Описание подпрограммы МОПИ для решения систем линейных уравнений

4. Численные исследования

4.1. Двумерные скалярные задачи рассеяния

4.2. Двумерные векторные задачи рассеяния

4.3. Трехмерные векторные задачи рассеяния

4.4. Исследование задач рассеяния на диэлектрике с помощью борновского приближения

4.5. Двумерные краевые задачи электростатики

Заключение

Литература

Приложение 1

Приложение 2

Введение диссертации (часть автореферата) на тему «Алгоритмическое и программное обеспечение для численного решения задач электромагнитного рассеяния на диэлектрике»

Введение

Диссертационная работа посвящена разработке алгоритмического и программного обеспечения для численного решения задач двумерного и трехмерного рассеяния и поглощения электромагнитных волн рэлеевского и резонансного диапазонов на неоднородном диэлектрическом теле. Указанные задачи являются объектом исследования многих специалистов. Это обусловлено тем, что подобные задачи имеют большое значение для практики, например, при проектировании диэлектрических антенн, при взаимодействии электромагнитного поля с биологическими объектами, при дистанционном обнаружении дислокаций в кристаллах и так далее. Математическое моделирование взаимодействия электромагнитного поля с диэлектрическими структурами значительно облегчает разработку устройств, а в ряде случаев является единственно возможным, например, при взаимодействии поля с биологическими объектами.

Указанные задачи описываются интегральными уравнениями относительно электрического поля в диэлектрическом теле. Для численного решения интегральные уравнения необходимо дискретизировать и свести к системе линейных алгебраических уравнений (СЛАУ) относительно неизвестных значений электрического поля в узлах расчетной сетки. Поскольку размерность СЛАУ N очень велика, особенно для трехмерных задач, то использование прямых методов типа метода Гаусса практически невозможно, так как это требует порядка N3 арифметических операций, а значит возможно использование только итерационных методов.

Указанные задачи изучались в работах отечественных ученых - А.Б. Самохина, С.П. Куликова, Ю.Г. Смирнова, В.И. Дмитриева, A.C. Ильинского, Е.В. Захарова [1-7] и других отечественных и зарубежных специалистов. В данной работе мы основываемся на полученных ими результатах.

При использовании итерационных методов решения СЛАУ наблюдается значительный выигрыш по числу арифметических операций и объему оперативной памяти, требуемых для реализации алгоритмов. По сравнению с методом Гаусса время решения задачи будет в N / т раз меньше, где т — требуемое число итераций. Таким образом, решение СЛАУ размерностью И» 1000 становится реально осуществимой задачей для персональных компьютеров.

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

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

В соответствии с целью работы поставлены следующие научные задачи:

• Исследование спектра рассматриваемых интегральных уравнений для различных классов задач.

• Алгоритмы нахождения оптимальных итерационных параметров для МОПИ.

• Разработка алгоритмов и подпрограмм, основанных на МОПИ, для решения СЛАУ, возникающих после дискретизации интегральных уравнений.

• Численные решения двумерных скалярных, а также двумерных и трехмерных векторных задач.

Научная новизна работы определяется следующими положениями:

• Численно исследованы спектральные картины интегрального оператора и их закономерности использованы для построения оптимального итерационного алгоритма.

• Определены явные формулы и алгоритмы для нахождения оптимального итерационного параметра для сходимости МОПИ.

• Исследована точность и область применения первого приближения МОПИ (оптимизированное борновское приближение), которое существенно использует информацию о спектре интегрального оператора. Показана эффективность этого метода, как в рэлеевском, так и в ближнем резонансном диапазонах длин волн.

В первой главе диссертации представлена подробная постановка задачи рассеяния электромагнитных волн на неоднородном диэлектрическом теле, а также представлены алгоритмы и подпрограммы расчетных моделей.

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

В третьей главе диссертации описаны алгоритмы нахождения оптимального итерационного параметра МОПИ для различных областей локализации спектра матрицы перехода на комплексной плоскости. Представлено описание подпрограмм определения оптимального параметра и подпрограммы, реализующей МОПИ для решения СЛАУ.

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

биологических структурах и тканях, двумерных и трехмерных векторных задач рассеяния. Представлены численные исследования задач рассеяния на диэлектрике с помощью борновского приближения. Также приводятся результаты решения двумерных краевых задач электростатики.

Основные положения, выносимые на защиту:

• Комплекс программ, разработанный в среде МаШсас!, позволяющий получать решения задач электромагнитного рассеяния на диэлектрических структурах.

• Алгоритмы нахождения оптимальных итерационных параметров для МОПИ для различных практически важных случаев области локализации спектра матрицы перехода.

• Комплекс алгоритмов и программ на основе МОПИ для решения интегральных уравнений.

• Численные результаты, полученные с помощью разработанного пакета программ по решению двумерных и трехмерных задач рассеяния электромагнитных волн.

Практическая ценность работы заключается в следующих положениях:

• Исследования показали практическую возможность моделирования электромагнитного рассеяния на двумерных и трехмерных локально неоднородных объектах.

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

• Оптимизированное борновское приближение может служить предобу-славливателем для нестационарных итерационных методов.

• Результаты исследований вошли в работы по проектам: «Исследование путей создания экспериментального образца аппаратно-программного комплекса для решения больших вычислительных задач (с приложениями в области электромагнетизма и гидроакустики)» (2013 г.), который проводился по Государственному контракту с Министерством образования и науки Российской Федерации; «Математическое моделирование взаимодействия электромагнитного поля со сложными трехмерными структурами на основе сингулярных объемных интегральных уравнений» ведомственной научной программы «Развитие научного потенциала высшей школы» (2009 г.). Практическая ценность выполненных разработок подтверждена актами о внедрении, представленными в приложениях 1 и 2.

Достоверность полученных результатов подтверждается математическими выкладками, а также проведенными численными исследованиями.

Апробация работы. Основные положения и результаты работы докладывались на:

• Международном симпозиуме «Progress In Electromagnetics Research Symposium» (Москва 2009 г.) - одном из крупнейших ежегодно проводящихся симпозиумов по применению и прогрессивному развитию теоретической и прикладной электродинамики и разнообразных ее приложений во всех частотных диапазонах; симпозиумы PIERS проводятся при поддержке Академии электромагнетизма с 1989 г.;

• Международной научно-практической конференции «Актуальные проблемы и перспективы развития радиотехнических и инфокоммуникацион-ных систем» (Москва 2013 г.);

• Международной научно-практической конференции «Современные информационные технологии и ИТ-образование» (Москва 2011 г.);

• Конкурсе прикладных разработок и исследований в области компьютерных технологий «Компьютерный континуум: от идеи до воплощения» (Москва 2011 г.) - работа отмечена дипломом;

• Международных конференциях «Информационно-вычислительные технологии в науке» (Москва 2007, 2008, 2009 гг.);

• Научно-практических конференциях «Современные информационные технологии в управлении и образовании» (Москва 2010, 2011 гг.);

• Научно-технических конференциях МГТУ МИРЭА (Москва 2007, 2008, 2009, 2011 гг.);

• Конкурсе «Лучшая научная работа студентов и молодых ученых» (Москва 2008 г.) - работа отмечена дипломом 3-й степени.

1. Задача электромагнитного рассеяния на диэлектрике

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

Пусть в свободном пространстве существует некая область Q, внутри которой магнитная проницаемость ¡л равна магнитной проницаемости свободного пространства /л0. Диэлектрическая проницаемость е внутри области задана достаточно произвольным распределением. Вне области диэлектрическая проницаемость постоянна и равняется е0.

Задача рассеяния электромагнитных волн на неоднородном диэлектрическом теле описывается интегральными уравнениями в области относительно полного электрического поля Е.

Интегральные уравнения электромагнитного рассеяния имеют вид:

• в двумерном скалярном случае -

Ё(х) = Ё\х) + кЩ [е(у) -1 }Ё{у)С{г) аУ,

е

• в двумерном и трехмерном векторных случаях -

Ё(х) = Ё\х) + и[е(х) -1]£(» + кЦ[Б{у) -1 ]Ё(у)С(г)с1у +

е

+ ър. I (ИдО -1 Щу), ра^роЛ 0{г) йу, (1.1.1)

е

где Ё{х) - неизвестное поле;

Ё° (х) — заданное первичное поле от сторонних источников в отсутствии диэлектрика;

б - комплекснозначное распределение относительной диэлектрической проницаемости по локально неоднородному телу <2;

(/(г) = Н^2\к0г) - функция Грина в двумерном пространстве, или

(7(г) =--функция Грина в трехмерном пространстве, где г=|х-у| -

4 • я • г

расстояние между точкой наблюдения х и точкой истока (интегрирования) у; и = -1/2 - для двумерного векторного случая, или V — -1/3 - для трехмерного векторного случая;

обозначает сингулярные интегралы, то есть интегралы по площади

(для двумерных задач) или объему (для трехмерных задач) за вычетом бесконечно малого круга (для двумерных задач) или шара (для трехмерных задач) в окрестности точки у=х.

Приведем уравнение (1.1.1) в трехмерном случае к виду, удобному для дальнейших вычислений, с помощью введения тензорной функции Грина:

т = -\{е{р)-№{рУ\-о{к5))+ / к\е{Ч)-\)ТШ¥1] , (1.1.2)

где Г = (1 + Г2УУ)С = в((\-1(кЯу1 ~(кК)-2)1 + (3(кЯу2 +31(кКу] -1)М). (1.1.3) Здесь С-С{К) - функция Грина в трехмерном пространстве, I -

единичныи тензор,

R =

rp-rq

, R = R/R, RR - тензор второго ранга, элементы

которого в декартовой системе координат имеют вид:

RR = -V R

if \2 > (Лг) АхАу Ах Az

AyAx (Ay) AyAz AzAx AzAy (Az)

(1.1.4)

у

Здесь Ах = xp-xq, Ay = yp-yq, Az = zp-zq\ для цилиндрической сетки

Ах = xm = rm cos ((pPiq), Ау = ум = fp q Sin ((pp<q).

Интеграл в смысле главного значения заменяется интегралом по области VjQ.p§, где Q.p S — шар исчезающе малого радиуса 8 вокруг точки наблюдения р. Известно, что

lim f FEdva = 0 ,

*-*°а /а

где Q.pr - шар любого радиуса г >5 вокруг точки р. Поэтому в (1.1.2) при

интегрировании по ячейке сетки при q = р остается взять численно интеграл по области между объемом ячейки и объемом шара, вписанного в нее.

Аналогичное (1.1.2) выражение получается и в двумерном случае Н-поляризации. Отличие касается первого слагаемого, коэффициент при котором и = -1/2. Тензорная функция Грина в этом случае имеет вид

где G = G(R) - функция Грина в двумерном пространстве; Hq2\JcR), H[2\kR)

А А

- функции Ханкеля 2-го рода нулевого и первого порядка; RR - тензор второго ранга с четырьмя элементами, аналогичными (1.1.4).

Отметим, что дискретным аналогом воздействия оператора на вектор Е является в трехмерном случае воздействие матрицы Т порядка 3N на столбец значения поля в узлах расчетной сетки. Основной частью диагональных

элементов ац матрицы Т является -^-(£„-1), и это значение не изменяется с

измельчением сетки.

Рассеянное поле в дальней зоне определяется на основе (1.1.2), причем в (1.1.2) первое слагаемое отсутствует ввиду отсутствия особенности в

подынтегральном выражении при р eV :

E^=\k2(s(q)-l)TEdVq (1.1.5)

v

Тензор Г = (I + r2VV)G представлен в (1.1.3), (1.1.4).

Поэтому, окончательно в трехмерном случае для поля в дальней зоне получим выражение

Е

scat

1

ATTR

рх\k2(e(q)-W(q)exp(ik(v,rq)dVq

Здесь Е - поле уже в основной системе координат, связанной с телом, и имеет, вообще говоря, все три отличные от нуля компоненты, а поле Е5Са' имеет угловые компоненты в сферической системе координат.

Аналогично, в двумерном случае Н-поляризации в дальней зоне:

Е.

scat

= ^ yj-^jj PxJ*2 (e(q) ~ 1)Е (q) exp (ikrq cos (q> - tpq )dSq

Диаграмма направленности (рассеяния) в двумерном и трехмерном случаях определяется соответственно как

ДН(^) = lim IkR

Я-> 00

Е

scat

2 »

ДЩв,<р)=\т AkR2

R-* оо

Е

. scat

Е"

Если источник поля - падающая плоская волна, то в векторном двумерном случае получаем

1

ДВД =

4 к

I х J к2 (e(q) - l)E(^) exp(ikrq cos(<p - (pq )dSq

(1.1.6)

а в трехмерном случае

ДН (в,<р) =

Алк

р х J кг (s(q) - l)E(g) ехр01(р, г, )dVq

(1.1.7)

1.2. Расчетная вычислительная сетка

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

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

Рис. 1.2.1. Цилиндрическая и декартовая системы координат.

Остановимся более подробно на цилиндрической расчетной модели. В этом случае для расчетов была принята равномерная цилиндрическая расчетная сетка, состоящая из КЯ круговых слоев в сечении плоскостью 2=сош1 и из К2.

слоев вдоль переменной Z с одинаковым шагом h = R/KR между слоями и между узлами расчетной сетки, где R - радиус цилиндра.

Общее число узлов расчетной сетки в сечении Z=const для равномерной сетки равно NR=3*KR , где KR - количество круговых слоев. Так, на рис. 1.2.2,а представлена сеточная картина в сечении Z=const при KR=6 и NR=108, на рис. 1.2.2,6 -KR=20 и NR=1200.

а) б)

Рис. 1.2.2. Цилиндрическая расчетная сетка в сечении Z=const.

В трехмерной области общее количество узлов сетки составляет где КУ, - количество слоев вдоль переменной Z.

Для формирования равномерной сетки в цилиндрических координатах задаем массивы радиусов (г) и углов (ф):

гфк(11,1,кос1) :=

к <- О

Гог 1е 0..1 - 1

гг. (2Л + 1)— 1 2

71

бф, 2—

3-(24+ 1)

Гог je0..3-(2.i+l)-l

Г, <Г- гг. к 1

<--71 + — +

кч-к+ 1 г 1Г код = 1 ф К код = 2 к кой = 3

к — количество ячеек; I - номер максимального слоя; щ - радиус ¡-го слоя; Ь - шаг по радиусу;

5ф1 — разница между двумя соседними ячейками по углу в одном слое; I - номер слоя; j - номер ячейки в слое; гк - массив радиусов;

фк - угол между осью х и радиус-вектором из центра окружности в центр ячейки.

Формирование равномерной сетки в объемных цилиндрических координатах на всем цилиндре происходит в 2 этапа:

1) задаем координаты в сечении (хх^ уэд) и координаты по ъ {ът^ъ);

2) переходим к декартовым координатам каждого узла цилиндрической

сетки.

} := 0,1.. Иху- 1 кЗ :=0..КЗ - 1

хх := гусоБ^) Ь 1 1 22кЗ 2 + ~ КЗ-— 2

N := Иху-КЗ ЛЛ/Л 7 I := 0.. N - 1 /л

'"тобО, Иху) У1:= ^0(1(1, Ыху) ^ := гг С \ } йтшс -

Т\[ху- количество узлов сетки на окружности; Г) - радиус ячейки; ((^ — угол ячейки;

КЗ - количество узлов на длине цилиндра; И - шаг сетки;

N — общее количество узлов цилиндра; хь у], ъ\ — координаты центра ячейки.

В случае декартовой расчетной модели формирование декартовых координат каждого узла происходит следующим образом:

N :=К1-К2-К3

к1 :=0..К1 - 1 т1 :=0.. N - 1 Ь

"ЪГ'

:=- + к1-Ь - К1- —

к2:=0..К2 - 1 т2:= 0.. N - 1

уу^ :Л + к2.Ь -

кЗ :=0..КЗ - 1

тЗ :=0.. N - 1

Ь Ъ

22^-.= -+ кЗ-Ь-КЗ—

Хт1 ' ^тоёСпй.К!)

1гипс

то(!(т2,К1К2) К1

2 о :=22: Г Л

тЗ ( тЗ 4

Ьгипс -

К1-К2

N - общее количество узлов сетки; к1, к2, кЗ - счетчики; К1 - количество точек разбиения по х; К2 - количество точек разбиения по у; КЗ - количество точек разбиения по ъ\ Ь - шаг сетки;

ххтЬ УУт2? 2гтз - вспомогательные координаты (в сечении и по длине); Хшь Уш2? 2т3 - координаты центра ячейки.

Таким образом, далее задача состоит в нахождении поля в узлах расчетной сетки.

Воздействие интегрального оператора на поле Е заменяется произведением матрицы перехода на вектор неизвестных значений Е в узлах сетки, что приводит к СЛАУ из 2И (для двумерных задач) или ЗИ (для трехмерных задач) уравнений, где N - количество узлов расчетной сетки в области тела. Матрица этой СЛАУ - комплексная, симметричная, неэрмитова.

В двумерном случае матрица перехода представляет собой:

В =

Вп вл

>12

\В2\ Вг2; Блоки матрицы перехода в двумерном случае:

В11

т, п

у(Ет - 1) + (ет " 1)^-~я-а-Н2(1,а) - 1 (еп - 1)-А11 (т,п)-1-а о&ег^Бе

Л

у

1:Г т = п

В21

т ,п

О ¡Г т = п

(сп - 1)-А21 (т,п)-71 -а"2 оЛет^е

В12

т, п

В22

т, п

О т = п

(еп - 1^-А12 (ш,п)-тг-а2 оШегтоБе

~ 0 + (£ш " 0( "7^-а-Н2(1,а) - 1

V Х

(бп - 1)-А22 (т,п)-7г-а оШепмзе

Л

И ш-

Составляющие матрицы перехода в двумерном случае:

All(m, п) := —

ч

Н2(0!Рт5П)-

H2(l,pm,n)

+ -4V

H2(0,pmjn)-2

А12(т,п):=-- Н2(0,рш>п)-2

Pm,n у Н2(15Рш,п)1 (Хш-ХпНУт~Уп)

Н2(1,рт,п)Л (хт"Л)

(Рт,п)"

'т.п

Ч

'т,п

А21(т,п) := А12(т,п)

А22(т,п) := — 4

Г

H2(0,pmjn)-

V

'm,n J

Н2(0,рт?п)-2

Н2(1,рт5П)Л (уш-Уп)

V

Рт,п

(рт3пУ

В трехмерном случае матрица перехода представляет собой:

^ Ви В\2 В\з^

В

Вг\ Вгг Вгъ чВъ\ Въг Взз;

Блоки матрицы перехода в трехмерном случае:

В11

m,n

В12

m,n

В22

m,n

В13

m,n

B23

m,n

у-(ет l) if n

(en - l)-Al l(m,n) h otherwise 0 if m = n

(en - l)-A12(m,n)-h3 otherwise - l) if m= n

(cn - l)-A22(m,n) h3 otherwise 0 if m = n

(sn - l)-A13(m,n)-h otherwise 0 if m = n

(en- l)-A23(m,n) h3 otherwise

B21 := B12 m,n m,n

B31 B13 m,n m,n

B32 := B23 m,n m,n

B33

m,n

-1

(em ~ l) if m~ n

(en - l)-A33(m,n)-h3 otherwise

В := Bll m,n m,n

В ЛТ:=В12 m,n+N m, n

Bm 7 xj := B13

m,n+2-JN m, n

В XT :=B21

В \т XT := B22 m+N,n+N m, n

В -vT _ := B23 m+N,n+2-N m,n

В _XT :=B31

В

'm+2-N,n+N

:=B32

m,n

^m+2-N,n+2-N '

Составляющие матрицы перехода в трехмерном случае:

2

+

А11(ш,п) := О

т,п

А12(п,т)

т,п

А22(п, ш) :=С

т,п

А13(п,т) := в

ш,п

А23(п,т) :=в

т,п

А33(п,т) := С

т,п

( Р ш, п) 3

2 Рт,п За

+

- 1

(рщ,п) 3

2 Рт,п

ЗЛ

+--1

(Рш,„)2 Рт'П 3-1

+

- 1

(рщ,п) 3

2 Рт,п ЗЛ

(х -х V

V т п)

(Рт,п)2

(х -х). 1т п/

(Уш"УпГ

(Рт,п)2

(к -к)-

V т п)

1 -1—— + ——

У ~ У т ' п

(Рт,п)2.

Рт'П (Рш,пГ| А21(п,т) := А12(п,т)

+

1 -1

1 - 1--+ -

ъ — ъ т и

(Рт,п)"

+

- 1

1(Рт,п)

3

2 Рт,п ЗЛ

(Рш,п )*

+

- 1

Рт,п

(Ут ~ Уп)

(г -г)2 V т п)

(Рт,п)2

г -г т п

(Рт,п)'

Рт'° (Рш,пГ А31(п ,т) := А13(п,т)

А32(п,т) :=А23(п,т)

1 -1 1 - I--+-

Рт,п

(Рт,п)2_[

Для численного решения СЛАУ используем итерационные методы ввиду их явного преимущества перед прямыми методами. В результате решения

получаем значения поля Е в узловых точках. Зная их, можно рассчитать поле в любой точке пространства, а также получить диаграмму направленности в дальней зоне.

2. Итерационные методы для решения задач электромагнитного рассеяния

Интегральное уравнение (1.1.1) даже до дискретизации представлено в виде, удобном для применения итерационного процесса, а именно Е = Еперв + ХАЕ (есть падающая волна и интегральный оператор).

Итерационные алгоритмы, в особенности быстро сходящиеся, имеют ряд преимуществ над прямыми методами, например:

• количество арифметических операций ~т-И , вместо

— число

итераций);

• отсутствие накопления ошибок в процессе итераций со сжимающим оператором;

• пониженные требования к оперативной памяти ЭВМ;

• если даже приходится детально исследовать спектр задачи для построения быстро сходящегося итерационного процесса то, однажды его построив, можно затем многократно использовать для расчетов с различными источниками - правыми частями уравнений.

Особенно эти преимущества заметны для задач с большими матрицами N»1000. Решение комплексной СЛАУ с N=1000 стандартным методом МаШсас! на ЭВМ Р-4 2ГТц занимает около 3 минут машинного времени, в то время как решение той же системы быстро сходящимся итерационным методом с числом итераций т ~ 10..20 требует всего около 1..2 секунды.

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

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

Одним из наиболее популярных стационарных итерационных методов является метод простой итерации, а также его модификации.

2.1. Решение задач методом простой итерации Запишем исходное уравнение электромагнитного рассеяния (1.1.1) в виде, более удобном для вычисления итераций

х = Тх + ф. (2.1.1)

Здесь X - неизвестный вектор, ф - заданный вектор правой части, Т - заданная матрица коэффициентов, так называемая матрица перехода, используя которую вычисляются итерации по следующим формулам

х(А:+1) = Тх{к) +ф, к-0,1,2,... (2.1.2)

Для случая Т = Е - А выражения (2.1.2) представляют собой запись метода простой итерации.

Если матрица Т постоянна (не зависит от номера итерации к), то такой итерационный процесс называется стационарным. Если матрица Т легко вычисляема и может быть явно представлена в законченном виде до начала вычислений, то такой итерационный процесс называется явным.

Процесс простой итерации может быть эквивалентно записан также в виде ряда по степеням оператора Т, то есть в виде, так называемого, ряда Неймана

х = / + Т/ + Т2/ + ...= ^Т1/.

/ = 0

00

(2.1.3)

Пусть X*- точное решение, строго удовлетворяющее исходному уравнению х =Тх + /, а АхК -хК -х - ошибка на к -м шаге итерации. Тогда из (2.1.2) получаем выражение для соотношения ошибок на к+\ п к-м шаге _ т£х(кК Значит для нормы ошибок имеем соотношения:

Дх(*+1) < \Т\\ Ах^

< Т

к+1

Д3с(0)

Отсюда следует достаточное условие сходимости процесса простой итерации: ||г|| = 8 < 1.

Действительно, тогда

<71

к+1

А3с(°)

= 5

к+1

Дх<°> « Ах<°>

(2.1.4)

Матрица перехода с = <1 называется сжимающей, а процесс (2.1.2)

для него сходящимся, т.к. ошибка убывает с каждым шагом независимо от её начальной величины.

Спектральным радиусом матрицы (конечномерного оператора) Т называется р(Т) = тах|#/|, где в. - собственные числа оператора Т [27]. /

Для любой нормы справедливо соотношение р(Т) < ||г|.

Доказывается [11, 24], что необходимым и достаточным условием сходимости процесса простой итерации (2.1.2) является условие

Р(Т)<\, (2.1.5)

при этом итерации сходятся не хуже геометрической прогрессии со знаменателем q = р(Т).

Условие (2.1.5) является, как правило, сильным ограничением для непосредственного применения метода простой итерации (2.1.2) для решения заданной СЛАУ. Выбор иного оператора Т с другим спектром при эквивалентности исходной системе уравнений может существенно расширить область сходимости процесса простой итерации:

= Тх{к) к = 0,1,2,... (2.1.6)

В качестве условия выхода из вычислительного процесса при достижении заданной точности решения £, как показано ниже, можно принять:

<

1-д

б , где д - спектральный радиус Т или какая-либо оценка дру-

гой нормы Т.

Покажем практический способ выхода из процесса итераций, гарантирующий достижение заданной точности вычислений в общем случае процесса простой итерации со знаменателем д. Будем полагать, что ошибка на п-ой ите-

рации вычислена с точностью е, то есть

Ах!

<е. Контролю же в процессе

вычислений поддаётся величина |Лхя|, где Ахи =х(,7) -х(п разность векторов

решения на двух соседних итерациях. Установив связь между этими величинами, мы получим возможность проводить вычисления с заданной точностью. Заметим, что х{п+к) - х(и) -» Ах* при к —» оо. Из неравенства треугольника по-

лучим

< \Ахп+к + Ахп+к_х +... + Ах„+11 < дк|Ахи| + дк~х| ДхЛ_,| +... + д\Ахп| = Я(1-Як)\

= д(1 + д + ... + дк-')\Ахп Тогда при к —» со имеем

1-д

Ах"

1 -9

К1-

Таким образом, требование

Ч

(2.1.7)

обеспечивает заданную точность вычислений.

В неявных методах простой итерации матрица перехода не может быть представлена в явном, легко вычисляемом виде. Исходное уравнение тем или иным способом представляется в виде

Вх = Сх + (р. (2.1.8)

Здесь В - матрица, для которой вычисление обратной матрицы В~1 представляет вычислительные трудности и не может быть представлено в аналитическом виде.

Процесс итераций строится аналогично вышеприведенному

Вх^=Сх^т)+(р. (2.1.9)

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

Т{ЭК6)=В~1С (2.1.10)

Результат процесса (2.1.9) такой же, как и процесса (2.1.2) при использовании в качестве матрицы перехода Т эквивалентной матрицы однако процессы (2.1.9) и (2.1.2) существенно различные.

В виде (2.1.9) могут быть записаны как неявные, так и явные итерационные процессы, но для явных процессов матрица В~1 легко вычисляема. Так, если В = Е- единичная матрица, то С -Т = Е- А и процесс известен, как явный метод простой итерации Ричардсона.

Если В = И, то С -Б- А = -{Ь + и), и матрица перехода

Т = -0~Х{Ь + и) (2.1.11)

легко аналитически вычисляется (так как для обратной матрицы элементы

известны и равны с1..=\ /а-, при i = j и нулю при / Ф _/). Процесс (2.1.2) с

1,1 1,1

матрицей (2.1.11) называется явным методом простой итерации Якоби. Здесь исходная матрица А представляется в виде суммы нижней треугольной матри-

цы Ь, диагональной матрицы И и верхней треугольной матрицы и. Таким образом, матрица перехода Якоби имеет нулевую главную диагональ.

В случае В = Ь + И и С = -II эквивалентная матрица перехода есть

Т(э кв) =_(£ + £>)"!£/. (2.1.12)

Обратная матрица здесь не является легко вычисляемой и в вычислениях не присутствует. Процесс (2.1.9) в этом случае называется неявным методом Зейделя.

Список литературы диссертационного исследования кандидат наук Сотникова, Наталья Юрьевна, 2013 год

Литература

1) Самохин А.Б. Интегральные уравнения и итерационные методы в электромагнитном рассеянии. —М: Радио и связь, 1998. —160 с.

2) Самохин А.Б. Объемные интегральные уравнения: методы и алгоритмы. М.: МИРЭА, 2011 - 96 с.

3) Смирнов Ю.Г. Математические методы исследования задач электродинамики. Пенза, Изд. Пензенского университета, 2009.

4) Математические модели прикладной электродинамики. Под. Ред. В.И. Дмитриева, А.С. Ильинского. М.: МГУ, 1984 - 216 с.

5) Куликов С.П. Итерационные методы и алгоритмы линейной алгебры. Учебное пособие — М.: МИРЭА, 2011 —168 с.

6) Куликов С.П., Самохин А.Б. Численное решение интегрального уравнения электромагнитного рассеяния: от 1D скалярного до 3D векторного случая // Электромагнитные волны и электронные системы. 2010. Т.15. №10. С.32-48.

7) A.B.Samokhin. Integral equations and iteration methods in electromagnetic scattering, VSP, Utrecht, The Netherlands, 2001, - 101 p.

8) Никольский B.B. Электродинамика и распространение радиоволн. М.: Наука, 1978.

9) Воеводин В.В., Тыртышников Е.Е. Вычислительные процессы с тепли-цевыми матрицами. - М.: Наука, 1987. - 234с.

10)Хейгеман J1., Янг Д. Прикладные итерационные методы. Пер. с англ.-М.: Мир, 1986.-448 с.

11) Самарский А.А., Гулин А.В. Численные методы. Учебное пособие для вузов. М.: Наука. 1989. - 432 с.

12)Волошанюк А.В., Воронина Н.Ю., Куликов С.П., Сафиулина А.А. Алгоритмы поиска оптимального значения комплексного параметра в модифицированном методе простой итерации. 56 Научно-техническая конференция, посвященная 60-летию МИРЭА. Сборник трудов. 4.2. Физико-математические науки. / МИРЭА. - М: 2007. С. 11-17.

13)Волошанюк А.В., Воронина Н.Ю., Куликов С.П., Сафиулина А.А. Использование алгоритмов поиска оптимального значения комплексного параметра модифицированного метода простой итерации для численного решения интегрального уравнения рассеяния. 57 Научно-техническая конференция. Сборник трудов. Часть 2. Физико-математические науки. / МИРЭА. - М: 2008. С.60-66.

14) Воронина Н.Ю., Куликов С.П. Алгоритмическое и программное обеспечение для численного решения интегрального уравнения рассеяния в трехмерном случае. 58 Научно-техническая конференция. Сборник трудов. Часть 2. Физико-математические науки./МИРЭА. - М: 2009. С.35-41.

15)S.P. Kulikov, N.Y.Voronina. Numerical Analysis of Scattering and Absorption Problems of Electromagnetic Waves of a Mobile Communication Range on Non-uniform Biological Structures. Progress In Electromagnetics Research Symposium Proceedings, 2009, pp. 1905-1909.

16) Воронина Н.Ю., Куликов С.П. Численное исследование задачи поглощения электромагнитных волн диапазона мобильной связи на неоднородных биологических тканях. Современные информационные технологии в управлении и образовании. Сборник научных трудов. М: ФГУП НИИ «Восход», 2010. 4.2. С.47-55.

17) Куликов С.П., Сотникова Н.Ю. Численное решение краевой задачи для уравнения Пуассона с использованием метода сопряженных градиентов и других итерационных методов. 60 Научно-техническая конференция. Сборник

Сотникова Наталья Юрьевна. Диссертация на соиск. уч. ст. канд. техн. наук. Москва 2013

трудов. Часть 2. Физико-математические науки./ МИРЭА. - М: 2011. С.54-58.

18) Куликов С.П., Сотникова Н.Ю. Численное решение краевой задачи для уравнения Пуассона с использованием метода релаксации и других итерационных методов. Современные информационные технологии в управлении и образовании: Сборник научных трудов. - М: ФГУП НИИ «Восход», 2011. 4.2. С.103-108.

19) Куликов С.П., Сотникова Н.Ю. Оптимизированное борновское приближение для оценки электромагнитного поля дифракции на диэлектрике. Современные информационные технологии и ИТ-образование. Сборник избранных докладов научно-практической конференции. Под ред. проф. В.А. Сухом-лина.-М.: ИНТУИТ.РУ, 2011. С.811-823.

20) Куликов С.П., Сотникова Н.Ю. Оптимизированное борновское приближение для векторного электромагнитного рассеяния на диэлектрике. — Электромагнитные волны и электронные системы. 2012. Т.17. №3. С. 11-17.

21) Сотникова Н.Ю. Численное решение задач электромагнитного рассеяния методом оптимальной простой итерации. — Актуальные проблемы и перспективы развития радиотехнических и инфокоммуникационных систем. Сборник научных трудов международной научно-практической конференции. -М.: МГТУ МИРЭА, 2013. 4.2. С.285-288.

22)Самохин А.Б. Численные методы решения многомерных интегральных уравнений математической физики с ядрами зависящими от разности аргументов. — Радиотехника и электроника, №2, 2005.

23) Самарский А. А. Введение в численные методы. Учебное пособие для вузов. 3-е изд., стер. - СПб: Издательство «Лань», 2005. - 288с.

24) Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы.— М.:Бином. Лаборатория знаний, 2006. — 636 с.

25) Вержбицкий В.М. Основы численных методов.—М.гВысшая школа, 2002. — 840 с.

26) Райе Дж. Матричные вычисления и математическое обеспечение.— М.: Мир, 1984.—264 с.

27)Воеводин В.В., Кузнецов Ю.А. Матрицы и вычисления.— М.:Наука. — 320 с.

28) Куликов С.П. Разработка и математическое обоснование алгоритмов и программ библиотеки для решения систем линейных алгебраических уравнений большой размерности стационарными итерационными методами // Теоретические вопросы вычислительной техники программного обеспечения. Межвузовский сборник научных трудов. — М.: МИРЭА, 2008. — с.78-94.

29) Куликов С.П. Строгая теория сходимости стационарных итерационных процессов на примере решения СЛАУ с трехдиагональной несимметричной теплицевой матрицей. // Сб. трудов 3-ей междунар. конф. «Параллельные вычисления и задачи управления». — М.:ИПУ им. В.А.Трапезникова, 2006, с.1208-1233.

30)Самохин А.Б. Метод сингулярных объемных интегральных уравнений для задач электромагнитного рассеяния на сложных трехмерных диэлектрических структурах. — Электромагнитные волны и электронные системы, 2007, 12, №8, с.7-14.

31) Dielectric Properties of Tissues and Biological Materials. A critical review. Crit.Rev.Biomed. Eng., Vol.17, pp. 25-104, 1989.

32) M. Okoniewski, M.A.Stuchly. A Study of the Handset Antenna and Human Body Interaction. IEEE Trans. On Microwave Theory and Techniques, Vol. 44, No. 10, 1996.

33) American National Standart safety levels with respect to human exposure to radio frequency electromagnetic fields, 300 kHz to 100 GHz, NY:IEEE, 1991.

Сотникова Наталья Юрьевна. Диссертация на соиск. уч. ст. канд. техн. наук. Москва 2013

34) Saad Y. Iterative Methods for Spars Linear Systems. SIAM. 2003. - 439 p.