CC BY
37
АНАЛИЗ ПОГРЕШНОСТЕЙ ЧИСЛЕННОГО РЕШЕНИЯ ЗАДАЧИ О РАСПРОСТРАНЕНИИ ЭЛЕКТРОМАГНИТНОГО ИЗЛУЧЕНИЯ В РАДИАЛЬНО-СИММЕТРИЧНОМ ВОЛНОВОДЕ
А.А. Белоусов, А.В. Гаврилов, А.А. Дегтярев Самарский государственный аэрокосмический университет
Аннотация
В работе решается задача о распространении монохроматической электромагнитной волны в радиально-симметричном волноводе, окруженном проводящей оболочкой. Для решения задачи построена разностная схема Кранка-Николсона. Исследование погрешности проведено в следующих двух частных случаях, допускающих аналитическое решение: а) в случае неучета эффекта самовоздействия среды [1] и постоянства коэффициента преломления; б) в случае параболической зависимости коэффициента преломления в радиальном сечении.
Введение
Рассмотрим процесс распространения монохроматической электромагнитной волны в среде, который может быть описан следующим уравнением [1]:
dU
-+-
k ■ П2П I ,2
AU--— |U| ■U = 0, (1)
дг 2 • к • п0 2 • п0
где и - напряженность электрического поля, к -волновое число, п0 - показатель преломления среды, пнл - изменение показателя преломления среды под действием поля распространяющейся волны, г - направление распространения волны.
Рассмотрим радиально-симметричный волновод, окруженный проводящей оболочкой. Для него уравнение (1) вместе с начальными и граничными условиями в цилиндрической системе координат примет вид:
dU
dz 2 ■ к ■ n0 2 ■ n0
0 < r < R, 0 < z < L,
UL=0=v(r) 0<r<R
U|r R = 0, 0 < z < L.
AU —^\U\2 ■U = 0,
(2)
Здесь R - радиус волновода, L - его длина, i//(r) - функция, описывающая напряженность электрического поля на входе в волновод, а напряженность является функцией координат U = U (z, r).
Далее рассматривается частный случай задачи, не учитывающий изменение показателя преломления под действием поля волны, то есть n нл = 0. Также будем считать, что n0 = n = const.
Подставив в (2) волновое число к = 2п/Л , где Л - длина волны, и записав оператор Лапласа Д r в полярной системе координат, получим окончательную формулировку задачи:
dU XI
-+ — I 2
dz 4nn ( dr
d 2U 1dU +--
r dr
= 0,
0 < r < R, 0 < z < L,
UL=0=v(r) 0<r<R,
U|r=R = 0, 0 < z < L, 5 2U
(3)
dU + _X__
dz 2nn dr2
= 0, r = 0.
Построение разностной схемы
Построим разностную схему для решения (3). Запишем для этого простейшие явную и неявную схемы для уравнения задачи, используя вспомогательный узел на сетке [3].
Пусть к2 = Ь/К - шаг сетки по направлению оси волновода, К - число узлов сетки по этому направлению, их координаты при этом гк = кк2,
к = 0,К.
Пусть также кг = К^ - шаг сетки по направлению, перпендикулярному оси волновода, 3 -число узлов сетки по этому направлению, их координаты при этом г^ = ]НГ, ] = 0,3 .
Тогда явная и неявная схемы примут, соответственно, следующий вид:
к+— к U 2 — ик i \
1 1 + C ( +ßk )= 0,
0,5 ■ h
.к+1 .. 2
1-+ C (ак+ +ßк+ )= 0,
0,5 ■ hz v 1 1 '
где C =-!-X--комплексный коэффициент,
4nn
к = ик+1 — 2ик + ик—— а 1 = 2
(4)
(5)
hi
(6)
в к =
и]+1 и1 -1 2Г]кг
(7)
+ и
к+1 СИ2_ 7+12И2
1+-
2 7
= -и
7-1
СИ 2И
а ик - сеточный аналог искомой функции напряженности.
Из (4) и (5), добавив начальные и граничные условия, получим разностную схему Кранка-Николсона для исходной задачи (3):
ик+1 - ик
С +—I
2 1
4+1
а" + аТ1 + вк +
вк+1)
= 0,
7 = 1, J -1, к = 0, К -1
и0 = у, 7 = 0, J
к+1 к [ к+1 7/ . - 1! .
- + 2С
А,
к+1
к Л
(8)
= 0,
к = 0,К -1 и^ = 0, к = 0К
где У] ).
Представим исходную задачу (3) и разностную схему (8) в операторной форме:
Ьи = / , (9)
ЬИиИ = /и . Очевидно, что
/ь /.
(10)
(11)
Нетрудно показать, что погрешность аппроксимации такой схемой исходной задачи будет:
||5/и|| = тах{0(А2,И2), 0, 0(И2,И2), 0} . (12)
Очевидно, что
5/и-
0.
(13)
Из (11) и (13) следует по определению, что разностная схема (8) аппроксимирует исходную систему (3) с порядком 2 относительно И, и Иг.
Также можно показать устойчивость разностной схемы (8) [4].
Так как сходимость следует из аппроксимации и устойчивости, то построенная схема Кранка-Николсона (8) сходится, причем, как видно из (12),
погрешность 0(И2г,И2г) .
характеризуется
выражением
+ и7
, КС 1
-и
]+1
СИ 2И
+27
(14)
При к = 0 получим для ] = 1, J -1:
7-1 2Й?
1 -1[1 -4т 1 + и1
2 7
1 - ИС
1
, СИ, ( 1 1 0 СИ2 ( 1 1 + и7-Н 1 +- 1 = -и 0-1-Н 1--1 +
7+1 2И2 У 27) 71 2ИГ2 У 27 )
+и
, И,С
1 + -£гг-
- и 0+1С1''+-
Значения функции и(г, г)
27
г=0
(15)
известны, т. е.
выражения в правых частях системы полученных уравнений можно считать величинами, зависящими только от индекса 7 :
ё0 =-и0-! ^ 7 1 2И?
]
1 1 + и 0
2 7 1 7
, И,С 1 +
0 СИ, [ 1 - и 0+!-Н 1 +-
7 2И2 У 27
7 = 1, J -1.
(16)
Матрица коэффициентов правых частей системы уравнений относительно функции напряженности электрического поля представляет собой трех-диагональную матрицу, т.е. для решения данной системы уравнений можно применить метод прогонки.
Дополним систему уравнением для 7 = 0 (преобразуем третье уравнение схемы (8)):
(
1 -
2ИС
Л
+ и
1 2СИг
= и0
1 +
2И,С
- и.
2СИ,
ё0 = и 0
00
1-
2И,С
и2
2СИ,
и2
(17)
(18)
Реализация разностной схемы Рассмотрим алгоритм решения системы (8). Выделим в левой части первого уравнения системы напряженность электрического поля на слое к +1 , а в правой части - на слое к; получим следующее выражение для 7 = 1, J -1:
к+1
СИг 17-1 2И2
1-
1
27
+ и
к+1
1 - ИС
Используя краевое условие разностной схемы
и° = 0 для 7 = J -1, получим:
СИ 2И
2,! - — 21 2^ -1)
+ и,
1 - ИС
17 2 2ИГ У 2( J -1)
+ и,
1 И,С
1 + -ТТ-
(19)
И
+
И
И
г
+
И
И
И
И
и
И
И
И
И
- и
и
! - 2
И
И
+
И
0 0 ChÄ_ uj-1 _ uJ-2 -, 2 2h2
1 --
1
2( J -1)
+ u,
i , ^Л i hzC 1 + ^ h2
(20)
Учитывая изложенные выше преобразования, получим следующую систему из 3 -1 уравнений с 3 -1 неизвестными, которая решается методом прогонки:
(
1-
2hC
\
+u
1 2Chz
= d0
1
Ch,
uj-1—2" I 1 - — 1 + u1 j12h2 l 2 j1 j
1 - hC
1
Ch,
+ uj+1—2r\ 1 + — 1 = d 0, j+12h2 l 2 j) j
j = 1, J -1
(21)
Ch
V - 2
+ U ,
z11 -
1
! h,C 1 ^
2( J -1)
Значение функции при ] = 3 получим из краевого условия:
= 0.
(22)
Нетрудно заметить, что матрица данной системы обладает диагональным преобладанием, что означает, что метод прогонки будет устойчивым.
Запишем общий вид системы на слое
к = 0, К -1:
U0
1 --
2h C
v
k+1 Chz Uj-1 2hr2
+ u,
+u
j+1
Chz 2h2
1
2 j
1
+ — 2j
1 2 Chz
h2
r
f
k+1 1 - l
j
-d)
= <,
hC
j = 1, J -1,
(23)
Ch
2\1 -"
1
2ht l 2(J -1)
Kf hr
+uJ-;\1 | = dJ-„
uk+1 = 0.
Решая методом прогонки систему (23) для к = 0, К -1, получим значения функции напряженности электрического поля и на слое к = 0, к = 1 и так далее, каждый раз используя решение на предыдущем слое. Действуя таким образом, получим значения искомой функции напряженности электрического поля для всех узлов сетки.
Аналитическое решение задачи Для проверки корректности работы построенной разностной схемы найдем аналитическое решение системы (3). В данном случае это возможно сделать. Решим для этого следующее уравнение:
Г д 2и 1 ди 1 - + —
dU + Äi dz 4rni
dr2
r dr
= 0.
(24)
Решение будем проводить методом разделения переменных. Сделаем подстановку (25) в (24), получим (26).
и (г, г ) = )• м>(г ) (25)
= -г{^(г) + ^(г)\ = М (26)
) /Я ^(г д г \
Решением уравнения
dv(z) = l/uÄ dz • v(z) 4nn является функция:
v(z)= С^хр^^П,j .
Уравнение
w"(r ) + 1 w'(r )+ jUw(r )= 0 r
(27)
(28)
(29)
является уравнением Бесселя [5]. В данном случае его решение имеет вид:
w(r) = J0 (4м • r)• С2, (30)
где J 0 (х) - функция Бесселя нулевого порядка [6]. Воспользуемся краевым условием:
Ar=R = С2 J0 (МR)= 0
(31)
Пусть у k - нули функции Бесселя J0. Тогда
R
Таким образом
Uk (z )= Ck • exP Тогда
ад
U(r,z)=£Ck • exp
f■ Äyl ^
l-/jLrr z
4nnR 2
J0 \ RR Yk
k=1
( 0 2 А ■ ÄYk
1—z 4rniR 2
•J0 \ R Yk
Воспользуемся начальным условием: UL=0 =^(r)= ZCk •J01RYk] .
(32)
(33)
(34)
(35)
То есть Ск - коэффициенты разложения функции у/(г) в ряд Фурье по базису функций Бес-
г
селя J0\ —Yk |:
R
Ck =
R 2J1
Ч^)r •dr . (36)
+
u
h
h
+
h
+
u
u
j -2
Реализация аналитического решения
Как видно из выражений (34) и (36), численная реализация данных выражений порождает ряд сложностей:
1) численное нахождение значений у к, J0 и J1;
2) численное интегрирование в (36);
3) приближенное вычисление значения бесконечного ряда в (34).
Значения у к являются табличными и могут быть взяты из справочников, либо рассчитаны с необходимой точностью [5, 6]. Для вычисления J 0 (х) и (х) могут быть использованы приближенные формулы необходимой точности [6].
Для общности применения реализации аналитического решения будем считать функцию у(г) заранее неизвестной и заданной не аналитически, а сеточным аналогом у = у(7). Значения всех остальных подынтегральных функций в выражении (36) могут быть вычислены с известной точностью в любой точке, в том числе и в узлах сетки. Таким образом, это задача численного интегрирования, в которой невозможно увеличение точности за счет измельчения шага интегрирования (т.к. размер сетки задан заранее). Поэтому необходимо использовать метод интегрирования, дающий достаточно точный результат при заданном числе узлов. В данной работе использовались квадратурные формулы интегрирования Ньютона-Котеса порядка 11 [7].
Вычислить численно бесконечный ряд вида (34) при неизвестной заранее функции у(г) не представляется возможным, поэтому целесообразно заменить его конечной суммой первых Ш его членов:
Ш { 1..2 \
UW ( z )=Z Ck exP
k=1
4nnR
2
r
Jo\ RYk
(37)
Число W выбирается из соображений минимизации вычислительной погрешности, возникающей из-за увеличения числа операций, также не обладающих абсолютной точностью, при увеличении числа членов ряда, что приводит к накоплению ошибки. Заметим, что указанный выше способ вычисления значений функций Бесселя дает не совсем корректные результаты при неасимптотических значениях аргумента, что приводит к сильному ухудшению результата при увеличении количества используемых членов ряда.
В качестве критерия для выбора W могут
служить вычисляемые значения U w (ri ,0), для которых известно точное значение y/(ri). Оценка точности проводится по формуле Д ^ = maxL(r-) - Uw ,0) . Тогда W = arg min Д w.
i=0, P 1 w
Поиск минимума, вообще говоря, можно проводить любым способом, т.к. время, затрачиваемое на эту операцию при численной реализации, даже при неоптимальном способе поиска на порядки
меньше времени, затрачиваемом на вычисление значений функции напряженности.
Заметим также, что полученная оценка погрешности для иШ (г ,0) может служить оценкой по-
грешности для UW (r, z)
zi, т.к.
exp
l-/Jkrr z
4mR 2
= 1.
Волноводы с неоднородным распределением показателя преломления Рассмотрим волновод, распределение показателя преломления которого неоднородно в направлении, перпендикулярном оси волновода. Заметим, что в уравнении (1) по-прежнему пнл = 0, но теперь п0 = n(rconst. Построенная нами разностная схема (8) применима и в таких случаях, потребуются лишь дополнительные вычисления при определении коэффициентов в системах в процессе решения.
Рассмотрим наиболее простой вид симметричного изменения показателя преломления n(r), а именно по параболическому закону:
l(r ) = V
(
1 - 2S
2
(38)
где п1 - значение показателя преломления на оси волновода, а 5 и а являются параметрами изменения показателя преломления [2], причем 5 > 0 и 0 < а < Я .
Рассмотрим случай, когда
¿1-1 << 1.
а,
(39)
Тогда (38) можно представить в следующем виде:
( / \ 2 Л
n(r) = п1
1 -\\
(40)
Подставив (40) в (1), решим его методом раз-
деления переменных: U (z, r ) = A(r )B(z ),
2ki • B-& = -1- ArA(r) = —u
B(z) n(r) A(r) ' Решением уравнения 2ki • B (z) = — u,
B(z) ^
очевидно, является функция
B(z )=c 4 uz ).
(41)
(42)
(43)
(44)
Преобразовав второе уравнение из (42), получим:
d2 A 1 dA
+--+ n(r)iuA = 0 .
dr2 r dr
Будем искать решение в следующем виде:
( ' \ 2 Л
A(r ) = X (r )exp
(45)
(46)
r
n
a
2
r
z
a
В результате ряда замен переменных и преобразований, потребовав выполнения равенства
1 = пхц8
, (47)
Ь4
4а2
(это возможно в силу отсутствия ограничений на ц и Ь ), получим уравнение Лагерра порядка
т =
8
а решением его будет функция
X (г )= Ьт
Г 2г 2 А
Ь2
(48)
(49)
где Ьт (х) - многочлен Лагерра порядка т [6].
Таким образом, решением уравнения (1) при условиях (38, 39) будет следующая функция:
Г / N2А
I г
ехр
и = С ехр|ф |Ь
Г 2г 2 А
Ь2
(50)
Экспериментальное исследование разностной схемы
Для проверки работы разностной схемы на вход волновода подавались различные распределения напряженности ц/(г), а потом анализировались различия между аналитическим решением и решением, полученным с помощью разностной схемы.
ЦуСГЪ ¥(г ) = 301 -К Ут
причем порядок корня
т зафиксирован. Очевидно, что решение примет вид:
Г л..2 А
и (г, г) = ехр
Яу т
4ппЯ
2
3 01 К У т
(51)
Заметим, что распределение интенсивности |и|2 не будет зависеть от г и будет в точности совпадать с начальным.
Далее в таблице 1 приведены результаты тестирования работы разностной схемы при различных параметрах и т = 2 . В колонке «%» указано максимальное по всем узлам сетки отклонение рассчитанной интенсивности от начальной по отношению к максимуму начальной интенсивности, в процентах.
Проведем теперь исследование для реальной функции распределения:
Г / \ 2 А
¥
(г ) = I(г )^_Р_ехр
(52)
где Р - мощность излучения, а а - эффективный радиус. Расчеты проводились для Р = 0,5 и а = 0,1- К. Результаты представлены в таблице 2. В колонке «%» указано максимальное по узлам сетки
отклонение сеточного решения от аналитического по модулю в квадрате в процентах относительно интенсивности максимального значения аналитического решения в этом слое (г = гк ).
Таблица 1. Результаты тестирования для решения (51)
п Я, мкм К , мкм Ь , м к, мкм К, мкм %
1 0,63 50 0,1 1,5625 1,5259 0,14
1,5 0,63 50 0,1 1,5625 1,5259 0,14
1,5 0,63 50 1 2,5 15,259 0,38
1,5 0,63 150 10 6 152,59 0,23
1,5 0,63 1000 30 10 1000 0,01
Таблица 2. Результаты тестирования для функции распределения (52)
п Я, мкм К , мкм Ь , м кг, мкм К, мкм %
1,5 0,63 50 0,01 0,167 10 6,66
1,5 0,63 50 0,01 0,25 2 4,21
1,5 0,63 1000 1 5 200 0,57
1,5 0,63 1000 1 1,25 1000 0,01
Заметим, что до определенной степени эффект несовпадения решений удается нивелировать измельчением сетки по радиусу волновода: уменьшение шага сетки по радиусу кг в 4 раза даже при увеличении шага сетки по длине в 5 раз снизило погрешность вычисления на 2 порядка. Но подобный подход существенно увеличивает время расчетов и имеет предел, связанный с усилением вклада вычислительной погрешности при слишком мелком разбиении.
На рис. 1 показано распределение интенсивности напряженности электрического поля на выходе волновода при разных способах решения для следующих параметров: длина волновода Ь = 0,01 м, радиус К = 50 мкм, показатель преломления п = 1,5 , длина волны Я = 0,63 мкм.
----- К=0,166 мкм, к=10мкм
------- К=0,333 мкм, К=3,3 мкм
----К=0,1 мкм, 1%=5 мкм
- аналитическое решение
Рис. 1. Интенсивность напряженности электрического поля на выходе волновода
г
в) 0,002
г) ,=0,01
Рис. 2. Интенсивность напряженности электрического поля в волноводе
На рис. 2 для этих же параметров показано изменение распределения интенсивности в процессе прохождения электромагнитного излучения по волноводу.
Для проверки работы схемы в случае волновода с неоднородным распределением показателя преломления естественно определить постоянную интегрирования в выражении (50) из начальных условий. Будем использовать в качестве начальных условий следующую функцию:
[ / \ 21
У(г ) = Ьт
[ 2г 2 1
ь2
ехр
Очевидно, что при этом решение примет вид
[ / Л21
и ( г ) =
ехр
2к
[ 2г 21
ь2
ехр
и( г) = и(0, г)
(53)
(54)
(55)
то есть можно провести проверку разностной схемы по интенсивности электромагнитного излучения, как это уже проделывалось ранее для случая однородного показателя преломления в волноводе. Но для этого определим сначала ряд параметров.
Зафиксируем т, тогда из (47) и (48) опреде-
Например, если взять т = 10 , то
лим ц и Ь2 215
Ь =-
(42)25 К б
и ц = -—'-—-. Кроме того, для обеспече-
та
ния корректности краевых условий потребуем, чтобы выполнялось равенство:
2Я2
(56)
где ^ 1 - один из корней полинома Лагерра Ьт (х). Возьмем, например, 1 = 8, тогда
^ 8 «16,2792578313781. Далее в таблице 3 приведены результаты тестирования схемы при параметрах п1 = 1,5; Х = 0,63 мкм; а = 0,000025; Ь = 0,01 м; И, = 10 мкм.
Таблица 3. Результаты тестирования для модели волновода с неоднородным распределением показателя преломления
5 Я , мкм Иг, мкм %
0,1 49,219 0,164 49,056
0,01 155,644 0,519 19,123
0,001 492,19 1,64 0,121
Заметим, что в действительности не совсем корректно утверждать, что при заданных в таблице 3 параметрах выполняется равенство (55), т.к. оно было получено в предположении (39). Но условие (39) для первых строк таблицы 3 явно несправедливо, поэтому утверждать, что схема действительно дает погрешность по интенсивности такую, как указано в таблице, нельзя. Действительно, анализ графиков интенсивности электромагнитного излучения в волноводе показывает, что интенсивность в процессе распространения излучения по волноводу осциллирует около начального значения, при этом положение минимумов и максимумов напряженности примерно сохраняется.
Но уже для третьей строки таблицы 3 условие (39) выполняется, равенство (55) тоже, что и объясняет столь высокую точность вычислений.
Проверка порядка сходимости
Пусть иИ - сеточное решение при конкретных значениях Иг и И, , а и - аналитическое с точностью выше, чем у сеточного. Погрешностью будем называть величину
е = \\ии -[иII , (57)
где норма имеет вид:
т = тах I
'и 1=0,! 1 к=0, К
(58)
Для рассматриваемой задачи должно быть справедливо выражение:
е=аИ2г + РН1, + о{к2г, И,2 ). (59)
И, вследствие (59), должно наблюдаться:
11т е = е, = РИ2
Нг
11т е = ег = ссИ Т
К ^0
(60)
Результаты исследования данных пределов приведены в таблице 4. Эксперимент проводился для параметров п = 1,5 ; X = 0,63 мкм; Я = 50 мкм; Ь = 0,01 м; И, = 33,3 мкм; Иг = 0,25 мкм и начально-
2
а
2
Ь
го распределения напряженности вида (52) с прежними параметрами. Отметим, что столь большие абсолютные значения погрешностей имеют причиной еще большие значения самой функции.
Таблица 4. Исследование пределов (60) для начального распределения напряженности вида (52)
Результаты исследование величин а и в приведены в таблице 5. Эксперимент проводился для параметров п = 1,5 ; X = 0,63 мкм; Я = 50 мкм; Ь = 0,01 м; к2 = 1,66 мкм; кг = 0,0083 мкм.
Опираясь на данные, приведенные в таблицах 4 и 5, можно определить, что сходимость является квадратичной.
Заключение
В работе построена разностная схема для решения задачи о распространении монохроматической
электромагнитной волны в радиально-симметричном волноводе, окруженном проводящей оболочкой, без учета влияния волны на преломляющие свойства вещества волновода. Также задача решена аналитически для двух случаев (в том числе для неравномерного распределения показателя преломления среды волновода), что позволило проверить корректность работы построенной разностной схемы. Результаты проверки показали, что данная схема вполне может быть использована для изучения процесса прохождения излучения через волновод. На расстояниях порядка 1 метра погрешность численного метода характеризуется величинами порядка 5%.
Благодарность
Работа поддержана грантом Президента РФ № НШ-1007.2003.1 и российско-американской программой «Фундаментальные исследования и высшее образование» («ВИНЕ»).
Литература
1. Виноградова М.Б., Руденко О.В., Сухоруков А.П. Теория волн // М.: Наука, 1979.
2. Адамс М. Введение в теорию оптических волноводов // М.: Мир, 1984.
3. Тихонов А.Н., Самарский А.А. Уравнения математической физики //Учебное пособие для университетов, М.: Наука, 1972.
4. Годунов С.К., Рябенький В.С. Разностные схемы (введение в теорию) // Учебное пособие М.: Наука, 1977.
5. Корн Г., Корн Т. Справочник по математике // М.: Наука, 1970.
6. Справочник по специальным функциям с формулами, графиками и математическими таблицами // Под ред. Абрамовица М., Стиган И., М.: Наука, 1979.
7. Сметанникова Е.Н. Численные методы математического анализа // Учебное пособие Самара: СГАУ, 2002.
К, мкм (к), 108 К, мкм ег (кг), 109
25 8,96 0,1 1,059
8,34 3,52 0,05 1,046
4,17 2,68 0,033 1,044
2,08 2,46 0,025 1,043
1,78 2,44 0,02 1,042
1,31 2,42 0,016 1,042
1 2,40 0,014 1,042
0,66 2,39 0,0125 1,042
Таблица 5. Исследование параметров а и в в формулах (60)
К, мкм в(к), 1018 К, мкм а(кг), 1021
25 1,215 0,1 3,86
20 1,465 0,05 4,05
16,7 1,617 0,033 4,2
