CC BY
64
ИССЛЕДОВАНИЕ ПОГРЕШНОСТИ РАЗНОСТНОГО РЕШЕНИЯ ОДНОНАПРАВЛЕННОГО УРАВНЕНИЯ ГЕЛЬМГОЛЬЦА МЕТОДОМ ВЫЧИСЛИТЕЛЬНОГО ЭКСПЕРИМЕНТА
Дегтярёв А.А.1, Козлова Е.С.2 1Самарский государственный аэрокосмический университет имени академика С.П. Королёва (национальный исследовательский университет),
2Институт систем обработки изображений РАН
Аннотация
Рассматривается вопрос об исследовании сходимости разностной схемы для однонаправленного уравнения Г ельмгольца. Основное внимание уделено экспериментальному исследованию скорости сходимости сеточного решения в равномерной и среднеквадратической нормах для линейного, а также нелинейного варианта уравнения Гельмгольца, учитывающего эффект «самовоздействия». Расчёты проводились на суперкомпьютере «Сергей Королёв» с использованием специально разработанных параллельных алгоритмов.
Ключевые слова: коэффициент преломления, эффект «самовоздействия», фокусировка, конечно-разностная схема, погрешность аппроксимации, исследование сходимости, метод прогонки, параллельный алгоритм, суперкомпьютер.
Введение
Высокая универсальность метода конечных разностей и связанные с ним широкие возможности для разработки параллельных вычислительных алгоритмов делает его эффективным средством решения многомерных краевых задач для дифференциальных уравнений, описывающих процессы распространения света в волноводах, что объясняет стабильный интерес исследователей к этой теме. При практическом применении метода конечных разностей важным является вопрос о фактической сходимости сеточных решений, получаемых с помощью разностной схемы, к решению исходной краевой задачи.
Численному решению краевых задач рассматриваемого класса, а также экспериментальной проверке сходимости разностных схем посвящён ряд работ [1 - 5 и др.]. Например, в работе [1] отмечено существенное влияние разрывности поля коэффициента преломления на скорость сходимости разностной схемы для уравнения Г ельмгольца. В работах [2 - 4] приведены результаты экспериментального исследования сходимости разностных схем для одно- и двумерных уравнений Гельмгольца с разрывным коэффициентом преломления, которые показали хорошее согласие экспериментальной скорости сходимости с теоретической. Однако следует отметить, что приведённые в этих работах результаты экспериментального исследования выполнены для весьма частных тестовых примеров, которые далеко не всегда могут быть формально «перенесены» на задачи с модификацией параметров, правых частей, начальных условий и т.п.
Настоящая работа посвящена экспериментальному исследованию методической погрешности конечно-разностного решения трёхмерного однонаправленного уравнения Гельмгольца, описывающего распространение электромагнитного излучения в волноводе прямоугольного сечения. Основное внимание уделено экспериментальному исследованию влияния на скорость сходимости сеточного решения таких параметров, как величина разрыва Аи = п1 - п2
поля показателя преломления среды, а также величина параметра нелинейности среды гн / е0.
Проведение такого исследования потребовало применения весьма мелких сеток, что в случае трёх пространственных переменных сопряжено с большими объёмами вычислительных затрат и необходимой оперативной памяти. Для проведения расчётов были разработаны параллельные вычислительные алгоритмы, основанные на методе расщепления с использованием для решения разностных СЛАУ методов встречной и правой прогонок с циклическим разбиением области. Для решения нелинейной задачи использовалась процедура итерационного уточнения. Экспериментальные исследования сходимости производились с помощью суперкомпьютера «Сергей Королёв» в Самарском государственном аэрокосмическом университете.
1. Теоретические основы
Математическая модель
Рассмотрим обобщённое волновое уравнение для векторного поля электрической напряжённости [6]:
-У[1п ц]х[Ух Е ]-У[у[1 п е]-Е ]-
д2Е (1)
дt2
-ДЕ + єц^^ = 0,
где ц - магнитная проницаемость, Генри/м; Е - напряжённость электрического поля, В/м; є - диэлектрическая проницаемость, Фарад/м.
В случае учёта эффекта «самовоздействия» диэлектрическая проницаемость может быть представлена в следующем виде [7]:
е = е„ +Є,- Е
(2)
где гср - среднее значение диэлектрической проницаемости, Фарад/м; гн - коэффициент изменения диэлектрической проницаемости под воздействием
поля распространяющейся волны; сивность излучения, мкДж.
Е = I - интен-
Будем использовать относительные величины магнитной проницаемости Д и диэлектрической проницаемости е:
Ц = ^ = 1, ^0
є
є=
є
є
є
є« +є,'1 , ~t т - -(I,-,
----------=-------+—I = є+—-|E
є0 є0 єі
(3)
(4)
J0
Также учтём, что при описании явлений для оптического (видимого, ультрафиолетового, инфракрасного) диапазона частот вместо относительной диэлектрической проницаемости рассматривается
коэффициент преломления я=л , который будем полагать функцией пространственной точки. Учитывая (2) - (4), запишем волновое уравнение в виде:
где с =
д2 E
n +-0- El |є0^0 -д~Г -
-v[v[ln є]- E ] - ДЕЕ = 0, 1
(5)
- скорость света.
■\/є0^0
Для получения однонаправленного уравнения Гельмгольца используем замену вида:
Е ( x, y, z, t ) = E( x, y, z)exp|i(rat + ф + kz)] , (б)
где ю = 2пс IX0 - угловая частота; ф - начальная фаза; k = 2nn IX0 - волновое число в заданной точке; k = 2nn IX0 - усреднённое значение волнового числа; Х0 - длина волны в вакууме; n - среднее по пространственной области значение коэффициента преломления.
Далее будем предполагать:
n (x ) = n + n (x), (l)
где n (x) - малое уклонение от значения n . Подставив (б) в уравнение (5), получим:
Л
n2 +-М E є„ 1 1
E д 2E є0^0Е —-Sz
„ r дЕ -2- f д2 д
—2ik---+ k E — I —— +-------------
Sz ^дг 3y
.2 ^2 +—— |e—
(8)
—vIv
[v[ln n2 ]-E] = 0.
Сделаем предположение о «медленности изменения огибающей» вдоль оси г [6], т.е.:
д 2Е
Sz 2
дЕ
Sz
(9)
В результате получим однонаправленное уравнение Г ельмгольца, учитывающее эффект «самовоз-действия»:
2ik -дЕ = ( k2 — k2 )Е —
Sz
д2
2п
є 1-І2 --ЧЕ Е —
д2
(10)
, +—^ |Е—v[v[inn2]-Е].
дт2 Sy2 1 [ [ ] ]
В случае ТЕ-поляризации (Е(x,y,z) = (0,Ey,0)) последнее уравнение представится в виде:
SEy=^_ Sz 2k
—у+^T |e +(k2 — k2 )e y +
Sx2 Sy2 1 y v > y
г 2
2п є 1 12
— Е Е
[Xo J є \ y y Ь0
(11)
Далее при записи формул будем опускать индекс в записи Еу, подразумевая под Е соответствующую компоненту вектора напряжённости электрического поля.
Математическая модель распространения электромагнитного излучения в волноводе прямоугольного сечения, ограниченном идеально проводящей оболочкой, будет иметь вид:
SE = j_ Sz 2k m є.
(k" —k~2 )Е+(І7+|HE
(12)
—F (Е),
\п е0
х 6 [0,К], У £ [0,1у], г £ [0,4];
Е|г=0 = у(х,у), X £ [0,1х], у £ [0,1у];
Е1=0 = ЕЦ =0, уе(0, 1у), ге(0,4];
Еу=0 =Еу=1у = 0, х £ (0,4 ), г е(0, 4 ]
где 4, 1У, 4 - ширина, высота и длина волновода соответственно, мкм; F(Е) = |Е|2 Е - функция «само-
воздействия»; у (х,у) - напряжённость электрического поля на входе в волновод, В/м.
Разностная схема Заменим область
О = {0 < х < 1х } х {0 < у < 1у } X {0 < г < 1г }
непрерывного изменения аргументов для поля Е(х,у, г) равномерной сеткой:
ц, ={x
*h ={xl = ihx , І = 0 I }x
}x{ yj = jhy , j = 0 J }x
x{z, = sh ,s = 0.
S>,
(1З)
и 1* и 1у и 1*
где ^ = —, hy = —, hг = — - шаги дискретизации.
I J S
Для численного решения нелинейной задачи (12) используем следующую конечно-разностную схему с итерационным уточнением:
2
X
є
0
0
2
+
+
2
—ю
Л,
.к * г
[и+" Г- и
0,5-2 2k
+Л-и;+(к2 - її2 )[и;+* ]
^ №Г), [Ц+* Т = и
л ' '
+ 1П —
х0п є0
т — 0, М; і — 1, I -1; ] — 1, J -1;
и ,+1 - и ,+/г г г
у 4 —-=1 Л и"
2к[ -х ]
0,5-
+л и;;1 +
(14)
+(Ч - *: К *]+ХП £-? (и *)
Х0П -0
] — V -1; і — 1, I -1; , = 0, £ -1;
и;1=0=^(х,у.), і—01, ]=017;
и^ = и‘\ — 0, , = 1, £, ] — 1,7 -1;
]11=0 ]І і—І •’
и!.\ — и , = 0, , = 1,£, /= 1,I -1,
]1]=0 ]1]=7
где и; - сеточная функция, взятая в узле (і,], 5),
и* . - 2и, + и 1,.
Л и" — ._______._____і-1.
- и — h2 ’
и5+. - 2и; + и ; .
Л и" — і+1___________________-_______і
к и
к
- разностные операторы
Лапласа по переменным х и у соответственно; т -номер итерации.
В случае линейного варианта задачи (12), т.е. если £к = 0, необходимость в итерационном уточнении не возникает и схема (14) фактически совпадёт со схемой переменных направлений Писмена-Рекфорда [8].
Исследуем свойство аппроксимации схемы (14). Прежде всего, запишем выражение для вспомогательной сеточной функции и^, которое нетрудно получить из первого и второго уравнений системы (14):
^+Х=\ и - ит ]+1 И+иг). (15)
Подставив в (15) вместо сеточной функции и точное решение Е(х, у, г) задачи (12) и проведя разложение с помощью формулы Тейлора в окрестности точки (х,у],г^), получим:
к к2 к
и— 1Е + Е"^ + Е" ^-Е'" +
. 1 2 2 22 4 ^ 8к
+о (-2-2)}
(16)
^■;-УП\ = {Е + ДЕ}х у _
,у.,2, хі,у., 2
где и5+/г - результат подстановки в (15) функции
Е (х, у, 2) вместо и ; ДЕІ — Е" — + О (к2).
У У ’ ’ Х ,У] ,2, 2 2 У2>
Запишем невязку 5/- для первого уравнения схемы (14) на точном решении краевой задачи (12):
5/—
2
х,, У., 2; к
і 2к
— -(и;+* -е( х,, у.,2))-
- І [л Л+х +Л у е( х,, у.,2;)+
(17)
(Ї2 - к2 и']-^2^F и).
2є„
Следует заметить, что функция F(E) не является дифференцируемой по комплексному аргументу Е. Поэтому воспользуемся следующим представлением:
F и+‘О—і? (е+де)!| х,, Уі .2, =
—{1е+дЕ 2 х(е+де)}| =
— {[ІЕ2 +|д^|2 + 2 (Re Е- Re ДЕ +
+1тЕ- 1т ДЕ)](Е +ДЕ)}| —
, у.,2,
— {|Е2 -е+|де| 2 -де+ІЕ 2 - ДЕ+|ДЕ|2 -Е+
+2 (Re Е - Re ДЕ+Іт Е - Іт ДЕ) х
(18)
((Е+ДЕ)1,, у., 2, = ? (Е)+?;(Е) к2Т+О (—2) •
X, У], г^
Далее, используя разложение функций, входящих в (17), кроме F (и^+'^2), по действительным переменным х, у, г и учитывая (18), получим:
5/—
х,, у., 2,
к К
— Е' +Е" -Е'" ^ -
4ї
-Е" -4-Е"' ^-Е" —L(ї2 -ї2)х
2ї 2хх 4ї уу 2ї 2ї ' . '
х| Е + Е„— +ЕІ-^|-^їП0£н? (Е)-
2
4
(19)
2є,
ІїП0Є^ (Е) |+О (-3, -2-2, -2, -2, -2-2, -2).
2є,
Продифференцируем первое уравнение системы (12) по переменной 2:
д 2Е і
с22 2ї
(к2 - 12 )дЕ + ^_^^ +^_
' ' д2 1 дх2д2 ду2д2
Е
(20)
И £ , , ч
+—-^'(Е).
Х0П -0
Если в выражении (19) учесть первое уравнение системы (12) и формулу (20), то будем иметь:
5/-
‘<х,,. у., 2, — О (-’,-2,-'2 ), і — ",1 - >•
(21)
. — 1,7 -1, к — 0, К -1.
Невязки 5/"/ = 0, S = 2,6, т.е. начальное и граничные условия разностной схемы, точно аппроксимируют исходные дифференциальные уравнения.
+
+
2
+
Таким образом, в случае линейной задачи в предположении достаточной гладкости решения погрешность аппроксимации схемы (14) будет квадратичной по всем шагам дискретизации. Для нелинейной задачи квадратичность погрешности аппроксимации будет иметь место в случае сходимости итерационного процесса при т ^ да.
2. Экспериментальное исследование сходимости
Линейное уравнение Гельмгольца с постоянными коэффициентами
Пусть и
- численное решение схемы (14),
полученное на сетке с шагами дискретизации Кх, Ку, hz. Если решение Е (х, у, 7) задачи (12) достаточно гладкое и имеет место (21), а также если схема (14) устойчива, то для погрешности метода сеток можно записать:
| = акХ + рКу2 + А2 + О (К, hA2, КК ) , (22)
где £ - погрешность численного решения и|к К К ■
Если взять более мелкую сетку, например с шагами дискретизации Кх = Кх / 2, Ку = Ку / 2, К7 = К7 / 2,
и рассчитать решение ^, то погрешность £
2 ’ 2 ’ 2
будет иметь вид:
Г К 1 2 Г Ку 1 2 Г к 1 2
а х + Р у + У 7 +
2 2 2
+О
( г к 1 3 К Г к 1 2 К 1 Г 2 Л
V _ 2 _ ’ 2 _ 2 _ 7 2 2 У
(23)
Очевидно, что погрешности £ и £ будут связаны соотношением:
(24)
Таким образом, следует ожидать уменьшения погрешности численного решения приблизительно в четыре раза.
Рассмотрим следующее начальное условие:
у(х, у) = sin | — х | sin
( Л —
Ту
V ‘у У
(25)
Оно является гладким и представляет собой собственные функции оператора Лапласа. В этом случае краевая задача (12) при постоянном коэффициенте преломления и при отсутствии эффекта «самовоз-действия» (т.е. £н = 0) имеет следующее аналитическое решение:
Е (х, у, 7) = ехр
I
2кп
(( л2 ( Л2^
~ —
V
V ‘у У
(26)
х sin
На рис. 1-2 представлены графики интенсивности комплекснозначной амплитуды, полученные с помощью разработанной схемы (14), а также рассчитанные по формуле (26). Расчёты проводились при следующих значениях параметров: Х=0,63 мкм,
1х = 10 мкм, 1у -J= 32, £ = 300.
10 мкм, 17 = 6 мкм,
= 3,6, I = 32,
Рис. 1. Распределение интенсивности вдоль оси Ох при у = 5 мкм, 7 = 6 мкм, полученное аналитически (кривая 1) и с помощью схемы (кривая 2)
I, мкДж
Рис. 2. Распределение интенсивности вдоль оси 07 при х = 5 мкм, у = 5 мкм, полученное аналитически (кривая 1) и с помощью схемы (кривая2)
Как видно из графиков, численное и аналитическое решения совпадают с высокой степенью точности, при этом мода сохраняет свой вид на рассмотренном участке волновода, что хорошо согласуется с теорией [9].
Проведём исследование скорости сходимости конечно-разностной схемы (14), меняя шаги согласно вышеописанному методу, и проследим за величиной погрешности в равномерной:
тах
1=0,/; j=0^; 1=0,Я
/ \|2 | Г
Е(х, у1,7.) - И
£6 I ( М"
, .03 .таЫЕ( х • у‘ ,'■■ I
и среднеквадратической нормах:
(27)
I / \|2 | |'
Е(х^у1 ^71 - И
(„ Л ( Л £ = -=0,/
— х) sin V ‘х У — Ту V 1у У '"зпё /=
/ • J • я
тах
(28)
п
2
7
х
На рис. 3 представлен график зависимости погрешности численного решения от номера сетки D (количества интервалов разбиения). Соответствующие значения погрешностей приведены в табл. 1. Расчёт проводился при следующих значениях параметров:
Тогда
X -п --
- 0,63 мкм, 1Х --
- 3,6, £ = 300. Ъу-10г10,%
10 мкм, 1У = 10 мкм,
4 = 2000 мкм,
Рис. 3. Графики погрешностей в среднеквадратической (кривая 1) и равномерной (кривая 2) нормах Таблица 1. Значения погрешностей
^ш-10, % £вх10"10, % D
0,5060 2,3800 1 3х3 х3
0,2390 0,9280 2 6х6 х 6
0,0376 0,1440 3 12х12 х12
0,0102 0,0377 4 24х24 х24
0,0028 0,0111 5 48х48 х48
Из таблицы 1 видно, что как в среднеквадратической, так и в равномерной норме скорость сходимости численного решения к точному близка к квадратичной. Стоит отметить, что погрешность сеточного решения близка к машинному нулю, что свидетельствует о высокой точности этого решения.
Линейное уравнение Гельмгольца с пространственным распределением коэффициента преломления Рассмотрим случай, когда коэффициент преломления среды зависит от пространственных переменных и представляет собой ступенчатую функцию:
п ( х ) =
1* и 1Х и
п, — + Ь < х < — - Ь;
1 2 2
п., — - Ь < х < — + Ь,
2 2 2
(29)
где Ь - ширина, а Ап = п1 - п2 - высота ступеньки.
В данном случае проблематично получить аналитическое решение задачи (12). Следовательно, вычислить погрешность сеточного решения так же, как это было и в предыдущем случае, не представляется возможным. Поэтому для получения информации о скорости сходимости численного метода поступим следующим образом. Запишем соотношение, связывающее точное и сеточное решения, предполагая сходимость последнего квадратичной:
= Е(х,-, У1, Ъ) + айХ + РйУ + т^2 +
/ ч — — — (30)
+0 (й3 йД2 йй2), i = 0,1,1 = 0, J, 5 = 0, £.
иц
ч\Ъх А а
и
у А ^ ^ 2’ 2 ’ 2
= Е( Х, Уу, ) + с
Г А 1 2 Г А, 1 2
Х + Р +
1_ 2 ] 2
+У
Г А 1 2 Г Г А 1 3 А Г А 1 2 А Г 1 2 ^
+о
2 V 2 2 2 2 2 У
(31)
і = 0, I, у = 0, J, 5 = 0, £.
В качестве характеристики погрешности численного решения возьмём следующую величину:
П =
ик а А -
х у 2* 2’ 2
и
ах А •'
4 И
+рАу2 +УК + о ( а3, ад2, аА).
-(“АХ2
(32)
Очевидно, что для сетки, в 2 раза более мелкой, будем иметь:
и|А: А К -[ і\а к к
2 ' 2 ' 2 4* 4*4
3
и\к К К 16 и|к К к
2* 2*2 Р 2* 2*2
-(аАх2
(33)
+РА2 +УА2 + о (А3, НАІАА).
Таким образом, следует ожидать выполнения приближённого соотношения:
П * (1/4)п . (34)
Вычисление погрешностей п и п проводится по значениям сеточного решения на двух сетках, одна из которых вдвое мельче другой. При этом разность значений сеточной функции вычисляется в узлах крупной сетки.
Проведём исследование скорости сходимости схемы (14) на тестовых примерах для двух случаев, различающихся высотой ступени Аи. Подадим на вход волновода собственные функции оператора Лапласа (27) и рассчитаем решение на сетках различной мелкости. На рис. 4-5 приведены графики зависимости погрешности численного решения от номера сетки, рассчитанные при следующих параметрах: X = 0,63 мкм, п1 = 3,6, Ь = 2 мкм, 1Х = 10 мкм, 1у = 10 мкм, 4 = 6 мкм.
Соответствующие значения погрешностей приведены в табл. 2 и 3.
Таблица 2. Значения погрешностей, п2=3,58
Пс,х10-3, % ги % D
8,6437 5,1646 1 60х8 х6
4,4839 6,3085 2 120х16 х12
1,7268 9,5504 3 240х32х24
0,2662 3,2471 4 480х64 х48
0,0437 1,6857 5 960х 128 х96
0,0079 0,8926 6 1920х256 х192
0,0015 0,4460 7 3840х512 х484
0,0003 0,2247 8 7680х1024 х768
3
+
+
Таблица 3. Значения погрешностей, п2=3,5
Пс,х10-3, % ги % D 8хШ
6,4882 9,7011 4 480х64 х48
4,1272 7,1334 5 960х 128 х96
2,3142 5,7083 6 1920х256 х192
1,1735 3,7718 7 3840х512 х484
0,5758 2,5396 8 7680х1024 х768
1 2 3 4 5 6 7 0
Рис. 4. Графики погрешностей в среднеквадратической норме при п2=3,58 (кривая 1) и п2=3,5 (кривая 2) г\,%
Рис. 5. Графики погрешностей в равномерной норме при п2=3,58 (кривая 1) и п2=3,5 (кривая 2)
Анализируя полученные результаты, можно сделать вывод, что сеточное решение конечно-разностной схемы (14) сходится.
Однако скорость сходимости сильно зависит от высоты ступени коэффициента преломления. Так, в случае п2 = 3,58 скорость сходимости близка к квадратичной в среднеквадратической норме, а в случае п2 = 3,5 она на порядок ниже. В равномерной норме в обоих случаях сходимость близка к линейной, хотя для случая с большим разрывом в коэффициенте преломления схема сходится значительно медленнее.
Результаты моделирования, представленные на рис. 6 - 7, были рассчитаны на следующих параметрах: X = 0,63 мкм, 1Х = 10 мкм, 1У = 10 мкм, 12 = 6 мкм, п1 = 3,6.
Из графиков видно, что чем сильнее разрыв в коэффициенте преломления, тем сильнее осцилляции сеточного решения в окрестности точки разрыва, которые медленно стабилизируются с измельчением
сетки. Стоит отметить, что в случае п2 = 3,58 использование сеток более высокой мелкости, чем сетка 4000 х 150 х 20, уже не даёт каких-либо визуальных изменений численного решения, в то время как в случае п2 = 3,5 визуальные изменения решения значительны, хотя неуклонно уменьшаются.
I, мкДж
Рис. 6. Численное решение при п2 = 3,5, у = 5 мкм, ъ = 6 мкм, полученное на сетке:
2000 х 150 х 20 отсчётов (кривая 1),
1000 х 100 х 10 отсчётов (кривая 2),
6000 х 200 х 30 отсчётов (кривая 3)
I, мкДж
Рис. 7. Численное решение при п2 = 3,58, у = 5 мкм, ъ = 6 мкм, полученное на сетке:
2000 х 150 х 20 отсчётов (кривая 1),
1000 х 100 х 10 отсчётов (кривая 2),
6000 х 200 х 30 отсчётов (кривая 3)
Нелинейное уравнение Гельмгольца с учётом эффекта самовоздействия Поскольку на данный момент решение краевой задачи для нелинейного однонаправленного уравнения Г ельмгольца не получено, для исследования погрешности и скорости сходимости схемы (14) применим ту же технологию, что и для линейной задачи со ступенчатым коэффициентом преломления. Подадим на вход волновода гауссов пучок:
р
2 2 х + У
(35)
Результаты численного моделирования представлены на рис. 8 - 9. Расчёты производились при следующих значениях параметров: п1 = 3,6, X = 0,63 мкм, 1Х = 10 мкм, 1у = 10 мкм, 4 = 6 мкм, I = 30, J = 30, £ = 4000, т = 5, £н / £0 = 0,5.
2
а
I, мкДж
£ мкм
X, мкм
Рис. 8. Распределение на входе в волновод
,, I, мкДж , __"-"Ч
10 "V х, мкм
Рис. 9. Эффект фокусировки в точке ъ = 6 мкм
Численное моделирование процесса распространения электромагнитной волны в нелинейной среде, проведённое с помощью разработанной разностной схемы, позволило наблюдать на вышеприведённом графике явление самофокусировки излучения (т.е. смещения энергии излучения к центру пучка), что вполне согласуется с теорией электромагнитного излучения [10]. Однако на рис. 9 наблюдаются ложные пики, которые исчезают с дальнейшим измельчением сетки. На рис. 10 приведён график установившегося поля решения, полученного в процессе измельчения сеток. Расчёты были проведены при следующих значениях параметров: п1 = 3,6, X = 0,63 мкм, 1Х = 10 мкм, 1у = 10 мкм, 4 = 6 мкм, I = 100, 3 = 100, £ = 6000, т = 5, ен / е0 = 0,5.
I, мкДж
На рис. 11 и в табл. 4 приведены результаты исследования погрешности численного решения и скорости сходимости схемы (14).
Л
Рис. 11. Графики погрешностей в среднеквадратической (кривая 1) и равномерной (кривая 2) нормах Таблица 4. Значения погрешностей
Пск, % Цр, % Б 8*1^
0,6474 1,5081 1 128х8 х8
0,2832 0,9631 2 256х16 х16
0,1236 0,6116 3 512х32х32
0,0544 0,4054 4 1024х64 х64
0,0255 0,2317 5 2048х128 х128
0,0118 0,1562 6 4096х256 х256
0,0056 0,0957 7 8192х512 х512
,х, мкм
Рис. 10. Эффект фокусировки в точке ъ = 6 мкм
Анализируя полученные результаты, легко заметить явление сходимости разностной схемы. Однако скорость сходимости не совпадает с теоретической, т. е. не является квадратичной и близка к линейной. Попытки увеличить количество итераций уточнения не привели к улучшению сходимости. Так, к примеру, различие сеточных решений для т = 4 и т = 5 характеризуется величиной 10-7 % и быстро убывает на следующих итерациях.
Приведённые результаты потребовали проведения расчётов на весьма мелких сетках. Эти вычисления были выполнены с использованием разработанных параллельных алгоритмов и привлечением высокопроизводительной техники.
3. Параллельные алгоритмы
Вследствие огромной трудоёмкости (десятки миллиардов и более вычислительных операций) для решения задачи были разработаны параллельные алгоритмы расчёта. Первый алгоритм основан на решении СЛАУ методом встречных прогонок с циклическим разбиением, второй - на решении СЛАУ методом правой прогонки с циклическим разбиением области [11, 12]. Для учёта зависимости правой части от решения на предыдущем слое, в которой задействован разностный оператор Лапласа, необходимо разбивать расчётную область на участки с перекрытием в три узла, что является особенностью разработанных алгоритмов.
Приведём результаты исследований для нелинейной задачи. При решении данной задачи применяется итерационное уточнение. В процессе расчёта на вход волновода подавались собственные функции оператора Лапласа (25) и использовались следующие значения параметров: п1 = 3,6, X = 0,63 мкм, 1Х = 10 мкм, 1у = 10 мкм, 1Ъ = 6 мкм, т = 5, £н / £0 = 0,5.
На рис. 12 представлено ускорение Талг, полученное за счёт применения параллельных алгоритмов. Расчёты производились на суперкомпьютере с использованием сеток, указанных в табл. 5.
Талг
Рис. 12. Ускорение параллельных алгоритмов: правой прогонки на двух процессорах (кривая 1), правой прогонки на четырёх процессорах (кривая 2), встречных прогонок на четырёх процессорах (кривая 3)
Таблица 5. Размеры сеток
D SxIxJ D SxIxJ
1 100x10 x10 9 1700x160 x160
2 250x20 x20 10 1900x180 x180
3 450x40x40 11 2100x200 x200
4 650x60 x60 12 2500x240x240
5 850x100x100 13 2900x280 x280
б 1100x120x120 14 3300x320 x320
7 1300x140x140 15 3700x360 x360
8 1500x160x160 16 4100x400 x400
Из графиков видно, что наилучший результат показал алгоритм правой прогонки для четырёх процессоров, достигнув ожидаемого ускорения в четыре раза. Несколько отстаёт от него алгоритм встречных прогонок.
На рис. 13 представлено ускорение Ттех, которое определяется как отношение времени работы программы на персональном компьютере ко времени работы программы на суперкомпьютере для соответствующих реализаций последовательного и параллельных алгоритмов.
Легко видеть, что суперкомпьютер в разы быстрее производит вычисления. В случае решения нелинейной задачи с итерационным уточнением применение последовательного алгоритма даёт ускорение примерно в 32 раза по сравнению с ПК. В случае применения параллельных алгоритмов наблюдается падение ускорения в направлении от крупных сеток к мелким, однако оно остаётся весь-
ма значительным (не менее, чем в 18 раз). Различия в ускорениях для алгоритмов распараллеливания на четыре и два процессора объясняется конфигурацией ПК (на основе Intel® Core(TM) 2 Duo CPU, 2 Гб оперативной памяти). Следует отметить, что решение такой задачи на ПК с использованием сеток мельче, чем 1100 х 120 х 120 отсчётов, является проблематичным вследствие нехватки оперативной памяти.
Ттех_______|_______(________|_______
65 С ' ' ' 1
ч\
55 \ \ -ф—4
45
35
25
15__________,_______.______,_________:___
1 2 3 4 5 0
Рис. 13. Ускорение, полученное за счёт использования суперкомпьютера «Сергей Королёв», для алгоритмов: правой прогонки на двух процессорах (кривая 1), правой прогонки на четырёх процессорах (кривая 2), встречных прогонок на четырёх процессорах (кривая 3), последовательного на одном процессоре (кривая 4)
Заключение
Получены следующие результаты:
1. Экспериментально подтверждена сходимость схемы (14) для рассмотренных тестовых примеров. При этом для случая линейной краевой задачи (12) с постоянными коэффициентами, т.е. при отсутствии эффекта «самовоздействия» (Ен = 0) и постоянстве коэффициента преломления, скорость сходимости оказалась близка к квадратичной как в равномерной, так и в среднеквадратической норме, что вполне согласуется с результатами теоретического исследования.
2. Установлено, что в случае пространственной зависимости коэффициента преломления, описываемого разрывной функцией (29), скорость сходимости существенно зависит от величины разрыва Ап = п1 - п2. Так, с ростом Ап эта скорость в среднеквадратической норме фактически падает на порядок: от квадратичной до линейной.
3. При решении нелинейной задачи, учитывающей эффект «самовоздействия» (ен / е0 = 0,5), скорость сходимости схемы (14) оказалась близка к линейной. Увеличение количества шагов т итерационного уточнения не позволило существенно повысить скорость сходимости схемы.
4. Использование специальных параллельных алгоритмов и высокопроизводительной техники позволило существенно (в 18 - 67 раз) уменьшить время расчётов для решения рассматри-
\ і
-*■-5 -Ф—4
NCv. \ ■
ваемой задачи, а также одновременно повысить на один-два порядка точность искомого решения за счёт измельчения сетки, задействовав при этом несущественную часть всех мощностей суперкомпьютера.
Благодарности
Работа выполнена при поддержке ФЦП «Научные и научно-педагогические кадры инновационной
России» (госконтракт № 14.740.11.0016), грантов
Президента РФ поддержки ведущих научных школ
(НШ-4128.2012.9) и молодого кандидата наук (МК-
3912.2012.2), и гранта РФФИ (12-07-00269).
Литература
1. Gustafsson, B. Time compact difference methods for wave propagation in discontinuous media / B. Gustafsson, P. Wahlund // J. Sci. Comput. - 2004. - V. 26 - P. 272-293.
2. Baruch, G. Fourth-order schemes for time-harmonic wave equations with discontinuous coefficients / G. Baruch, G. Fibich, S. Tsynkov, E. Turkel // Commun. Comput. Phys. - 2009. - V. 5 - P. 442-455.
3. Feng, X. High order compact finite difference schemes for the Helmholtz Equation with discontinuous coefficients / X. Feng, Z. Li, Z. Qiao // Journal of Computation Mathematics. - 2011. - V. 29(3) - P. 324-340.
4. Feng, X. A high-order compact scheme for the onedimensional Helmholtz equation with a discontinuous coefficient / X. Feng // International Journal of Computer Mathematics. - 2012. - V. 1 - P. 1-7.
5. Злотник, И.А. Семейство разностных схем с приближёнными прозрачными граничными условиями для обобщённого нестационарного уравнения Шрёдингера в полуполосе / И.А. Злотник // Журнал вычислительной математики и математической физики. - 2011. -Т. 51(3) - С. 384-406.
6. Неганов, В.А. Линейная макроскопическая электродинамика. Т. 1. / В.А. Неганов, С.Б. Раевский, Г.П. Яровой. - М.: Радио и связь, 2000. - 509 с.
7. Виноградова, М.Б. Теория волн / М.Б. Виноградова,
О.В. Руденко, А.П. Сухоруков. - М.: Наука, 1979. - 385 c.
8. Самарский, А.А. Численные методы математической физики / А.А. Самарский, А.В. Гулин. - М.: Научный мир, 2003. - 316 c.
9. Борн, М. Основы оптики / М. Борн, Э. Вольф. - М.: Наука, 1973. - 720 с.
10. Адамс, М. Введение в теорию оптических волноводов / М. Адамс. - М.: Мир, 1984. - 512 с.
11. Ортега Джеймс, М. Введение в параллельные и векторные методы решения линейных систем / М. Ортега Джеймс. - М.: Мир, 1991. - 364 с.
12. Головашкин, Д.Л. Параллельные алгоритмы решения сеточных уравнений трёхдиагонального вида, основанные на методе встречных прогонок / Д.Л. Голо-вашкин // Математическое моделирование. - 2005. -Т. 17, № 11. - С. 118-128.
References
1. Gustafsson, B. Time compact difference methods for wave propagation in discontinuous media / B. Gustafsson, P. Wahlund // J. Sci. Comput. - 2004. - V. 26 - P. 272-293.
2. Baruch, G. Fourth-order schemes for time-harmonic wave equations with discontinuous coefficients / G. Baruch, G. Fibich, S. Tsynkov, E. Turkel // Commun. Comput. Phys. - 2009. - V. 5 - P. 442-455.
3. Feng, X. High order compact finite difference schemes for the Helmholtz Equation with discontinuous coefficients / X. Feng, Z. Li, Z. Qiao // Journal of Computation Mathematics. - 2011. - V. 29(3) - P. 324-340.
4. Feng, X. A high-order compact scheme for the onedimensional Helmholtz equation with a discontinuous coefficient / X. Feng // International Journal of Computer Mathematics. - 2012. - V. 1 - P. 1-7.
5. Zlotnik, I. Family of finite-difference schemes with approximate transparent boundary conditions for the generalized nonstationary Schrodinger equation in a semi-infinite strip / I. Zlotnik // Computational Mathematics and Mathematical Physics. - 2011. - V. 51(3) - P. 384-406. - (In Russian).
6. Neganov, V.A. Linear macroscopic electro-dinamics / V.A. Neganov, S.B. Raevsky, G.P. Yarovoi. - Moscow: “Radio I svyas” Publisher, 2000. - V. 1. - 509 p. - (In Russian).
7. Vinogradova, M.B. The theory of waves / M.B. Vinogradova, O.B. Rudenko, A.P. Sukhorukov. - Moskow: “Nau-ka” Publisher, 1979. - 358 p. - (In Russian).
8. Samarsky, A.A. Numerical methods of mathematical physics / A.A. Samarsky, A.B. Gulin. - Moscow: “Nauch-niy mir” Publisher, 2003. - 316 p. - (In Russian).
9. Born, M Fundamentals of optics / M. Born, E. Wolf. -Moscow: “Nauka” Publisher, 1973. - 720 p. - (In Russian).
10. Adams, M. Introduction to the theory of optical waveguides / M. Adams. - Moscow: “Mir” Publisher, 1984. -512 p. - (In Russian).
11. Ortega James, M. Introduction to parallel and vector solution of linear systems / M. Ortega James. - Moscow: “Mir” Publisher, 1991. - 364 p. - (In Russian).
12. Golovashkin, D.L. Parallel algorithms for solving tridiagonal finite-difference equations based on method of counter-runs / D.L. Golovashkin // Mathematical modeling. - 2005. - V. 17, N 11. - P. 118-128. - (In Russian).
INVESTIGATION OF ACCURACY OF NUMERICAL SOLUTION OF THE ONE-WAY HELMHOLTZ EQUATION BY METHOD OF COMPUTATIONAL EXPERIMENT
A. A. Degtyarev1, E. S. Kozlova2
1 S.P. Korolyov Samara State Aerospace University,
2 Image Processing Systems Institute of the RAS
Abstract
The question of investigating the convergence of the finite-difference scheme for one-way Helmholtz equation is considered. The main attention is paid to the evaluation of the speed of convergence of the numerical solutions in a uniform and standard norms for linear as well as nonlinear variant of
the Helmholtz equation that takes into account the “self-influence” effect. The calculations were performed on a supercomputer “Sergey Korolev” using specially developed parallel algorithms.
Key words: guidelines refractive index, “self-influence” effect, focusing, finite-difference scheme, accuracy of approximation, research of convergence, sweep method, parallel algorithm, supercomputer.
Сведения об авторах
Дегтярёв Александр Александрович, 1953 года рождения. В 1977 году с отличием окончил Куйбышевский авиационный институт, ныне - Самарский государственный аэрокосмический университет имени академика С.П. Королёва (СГАУ), по специальности «Прикладная математика». В 1985 году защитил диссертацию на соискание степени кандидата технических наук. В настоящее время работает доцентом кафедры технической кибернетики СГАУ. Область научных интересов: численные методы и компьютерное моделирование в физике.
E-mail: aadest@,smail.com .
Aleksandr Aleksandrovich Degtyaryov (b. 1953) graduated with honour (1977) from the S.P. Korolyov Kuibyshev Aviation Institute (presently, Samara State Aerospace University, SSAU), majoring in Applied Mathematics. He received his PhD in Technical sciences (1985). At present time he is holding a position of Associate Professor at SSAU’s Technical Cybernetics sub-department. Research interests include numerical methods and computer simulation in physics.
Козлова Елена Сергеевна, магистр прикладной математики и информатики, аспирант кафедры технической кибернетики Самарского государственного аэрокосмического университета имени академика С.П. Королёва (СГАУ). Область научных интересов: дифракционная оптика, численные методы.
E-mail: [email protected].
Elena Sergeevna Kozlova, Master of Mathematics and Computer Science. Currently studies at Samara State Aerospace University. Research interests are diffractive optics and numerical methods.
Поступила в редакцию 26 октября 2011 г.
