CC BY
29
УДК 517.6
С. Н. Куприянова
ЧИСЛЕННЫЙ МЕТОД РЕШЕНИЯ НЕЛИНЕЙНОЙ ЗАДАЧИ НА СОБСТВЕННЫЕ ЗНАЧЕНИЯ ДЛЯ НЕОДНОРОДНОГО ВОЛНОВОДА
Аннотация. Разработаны и обоснованы два численных метода нахождения собственных значений и соответствующих им собственных функций. Результаты обобщены на случай неоднородного диэлектрического волновода кругового сечения с заполнением нелинейной средой, выраженной законом Керра. Ключевые слова: неоднородный волновод, закон Керра, задача на собственные значения, дисперсионные соотношения.
Abstract. Two numerical methods for solving nonlinear eigenvalue problem are proposed and justified. The results are generalized for the case of non-homogeneous dielectric waveguide of circular cross section filled with a dielectric exhibiting a local Kerr-type nonlinearity.
Keywords: non-homogeneous waveguide, Kerr-type nonlinearity, eigenvalue problem, dispersion relations.
Рассмотрим задачу о собственных волнах цилиндрического диэлектрического волновода. Пусть все трехмерное пространство R заполнено изотропной средой без источников с диэлектрической проницаемостью £l = const. В эту среду помещен цилиндрический диэлектрический волновод неоднородного заполнения с образующей параллельно оси Oz и поперечным
2 2 21
сечением W :={х: Xi + х2 < R i. Будем предполагать гармоническую зависимость полей от времени в виде [1]
Е (x, y, z, t) = E+ (x, y, z)cos at + E - (x, y, z)sin at;
Н(x, y, z, t) = Н + (x, y, z) cos at + Н - (x, y, z) sin at,
где a - круговая частота; Е , Е +, Н , Н - вещественные искомые функции. Образуем комплексные амплитуды полей E(x, y, z), H (x, y, z):
Е = Е + + iE -; н = Н + + iH -.
Везде ниже множители cos at, sin at будем опускать.
Пусть диэлектрическая проницаемость е внутри цилиндра определяется по закону Керра [2].
Среда предполагается изотропной и немагнитной, ц = jiQ .
Требуется отыскать поверхностные волны, распространяющиеся вдоль образующей волновода, т.е. собственные волны структуры. Электромагнитное поле Е, Н удовлетворяет системе уравнений Максвелла
rot H =-iaeE; (1)
rot Е = /юцИ, (2)
условиям непрерывности касательных составляющих поля Ит и Ет при переходе через границу волновода и условиям экспоненциального затухания поля на бесконечности.
Перейдем к цилиндрической системе координат (р, ф, z). Тогда уравнения Максвелла примут вид
1 Эеz ЭЕф
р Эф Эz
ЭЕр ЭЕг
Эz Эр ф
1 э ( Е ) 1 ЭЕр
- ч-(РЕФ)-----------
р Эр Y Р Эф
1 ЭИ2 ЭИф
■ = /юцИр; (3)
= z'fflvH ф; (4)
“ ~(рЕф) - - ^ = /юцИ2 ; (5)
р Эф Эz р
ЭИр ЭИ2
= -/юєЕр; (6)
Эz Эр
= -іюєЕф; (7)
1 Э 1 ЭИр
- — (рИф)----------р = -/югЕ2 . (8)
р Эр Y р Эф
В случае ТЕ-поляризации предположим, что Е = {0; Еф; 0J,
И ={Ир; 0; Hz J. В результате уравнения (3)-(8) приведутся к виду
ЭЕф
—э^ = /юцИр; (9)
1 А(рЕф) = /юцИ2 ; (10)
р Эр
1 ^=0, (И)
р Эф
= -/юеЕф; (12)
ЭИр эи2
Эz Эр _ <Ш^ф
1 эир
------^ = 0. (13)
р Эф
Из (11) и (13) следует, что Hz = Hz (р, z) и Ир = Ир (р, z) не зависят от ф . Из уравнений (9) и (10) находим
1 ЭЕф 1 1 Э
Ир= - ---------; Hz = ------------— (рЕф). (14)
к /юц Эz /юц р Эр Y
Э 1 Э Э2Еф 2
—(- — (рЕф)) + —2^^^ = 0. (15)
Эр Р Эр ^ Эz2 ^
Решение задачи будем искать в форме осесимметричных волн Еф (р, у, z) = и (р, у)е1^2, где у - вещественная постоянная распространения волны. Будем предполагать, что и (р, у) - вещественная функция.
Таким образом, (15) можно быть переписано в виде
(1 Л 2 2
— (ри)' ) + (ю ец-у )и = 0, (16)
р )
где производная означает дифференцирование по р. Во внешней области, учитывая, что е = е1, получаем уравнение Бесселя:
и" + 1и—1-и + к2и = 0, р>Я , (17)
р р2
,22 2 где к = ю е1^-у .
Внутри волновода, где е = е20 + а ^ е21 (р) + |Е|21, получаем кубическое
нелинейное дифференциальное уравнение второго порядка:
11? 3
и Н— и----г—Н к и Н ^(и Н е21 (р)и) = 0, 0 <р < Я , (18)
р р2и
2 ,2 2 2 где а = ю ац , к =ю е2Ц - у .
Условия сопряжения на поверхности волновода преобразуются к виду
[ЕФ ]р=Я = 0 и [Hz ]р=Я = 0, что дает
[и ]р=Я = 0,
[и ]р=К = 0, (19)
где [и ]р_^ = и(Я - 0) - и(Я + 0) - скачок предельных значений функции в точке Я . Спектральным параметром задачи является у .
Сформулируем теперь краевую задачу на собственные значения Р, к которой свелась исходная задача о распространяющихся поверхностных волнах цилиндрического волновода. Требуется отыскать ненулевую, ограниченную и непрерывно-дифференцируемую на полубесконечном интервале р > 0 функцию и(р) и соответствующие собственные значения у, такие что и(р) удовлетворяет уравнениям (17) и (18), условиям сопряжения (19) и условиям экспоненциального убывания функции и(р) на бесконечности при р ^ ^.
Запишем решение уравнения Бесселя (17) в виде
и = С1Н111)(кр), р>Я , (20)
где С1 - произвольная действительная неисчезающая константа, а Н® -функция Ханкеля.
Принимая во внимание условия излучения, выберем решение уравнения в форме
u = CXKX(\k\)р, р>Я ,
(21)
где К1 - функция Макдональда.
Условия излучения выполняются, потому что К1 (|к| )р ^ 0 экспоненциально при р ^ го .
Перепишем нелинейное уравнение (18) в виде
( - 1 Л
I 2 1 I / 3 \
(ри')' + k р------и + ар|и + Є2і(р)риj = 0
v Р / и рассмотрим линейное уравнение Бесселя:
ри + и +
12 1
k р------и
р
= 0.
Перепишем последнее в операторной форме:
Т 0 т d2 d (12 11
Lu = 0, L = р—2 + — + 1 Р-------
d р2 d р V р J
(22)
(23)
(24)
Используя стандартный метод, построим функцию Грина для краевой задачи
10 = -8(р-р0),
0р_ = 0'|р=я = 0 (0 <р0 < «)
в виде (см., например, [3])
п
N1 (kR) - Jl(kр<)Ni(kР>)
0 < р, ро < R,
где
р( = min{р, Ро } ,
р) = max{р, ро } . (25)
Функция Грина существует при таких значениях параметров, что
JlkR Ф 0.
Ьи + аВ(и) = 0, В(и) = ри + р£21 (р)и . (26)
Используя вторую формулу Грина [4]
Я я
|(иЬы - иЬх)йр = |(ри')'-и(р"и')'))р = Я(и'(Я)и(Я)-"и'(Я)и(Я)) (27)
0 0
и полагая и = О, получаем
Я
| (ОЬи - иЬ0)йр =
0
= Я(и(Я - 0)0(Я, р0) - 0'(Я, р0)и(Я - 0) = Яи(Я - 0)0(Я, р0). (28)
Используя уравнение (26), имеем
Я Я
|(ОЬи -иЬ0)йр = -а|ОВ(ы)йр + и(р0) (29)
0 0
и получаем интегральное представление решения и(р0) уравнения (18):
Я
и (р0) = а 10(р, р0)р(и3(р) + Ё21(р)и (р)) й р + Яи(Я - 0)0( Я, р0), 0 <р0 < Я . (30)
0
Принимая во внимание условия сопряжения
и( Я - 0) = и( Я + 0),
перепишем уравнение (30) в виде
Я
и (р0) = а| О (р, р0)р(и3(р) + е21(р)и (р))й р + /(р0), 0 <р0 < Я, (31)
0
где
/ (р0) = Яи'( Я + 0)0( Я, р0) (32)
и
0(Я,р0) =— ^р0). (33)
0 кЯ J{(kЯ)
При этом существенно, что /(р) не зависит от и. Из (31) и условий сопряжения и (Я - 0) = и (Я + 0) следует дисперсионное соотношение
Я
и (Я + 0) = а | О (р, Я )р(и 3(р) + Ё21(р)и (р))й р + Яи '(Я + 0)0 (Я, Я). (34)
Положим N (р, Po) = aG(p, Ро)Р и рассмотрим интегральное уравнение в C[0, R] (см. также [5]):
R
u(Po) = jN(р,Po)u3(p) + е21 (p)u(p))dр + f (Po) . (35)
0
Предполагается, что f е C[0, R] и J[(kR) Ф 0. Нетрудно видеть, что ядро N(р,Po) является непрерывной функцией в квадрате 0 < р, Po < R . Таким образом, справедливы следующие утверждения.
Утверждение 1. Линейный интегральный оператор
R
Nro = j N (р, Po)m(p)J р (36)
o
ограничен, вполне непрерывен в C [0, R] и
R
И = max Г|N(р,Po)dр . (37)
' ' Poe[0,R]0'
Утверждение 2. Нелинейный оператор R
F (u) = j N (р, Po)(u 3(р) + е 2i(p)u (р)) d р + f (po) (38)
o
является вполне непрерывным на каждом ограниченном в C[0,R] множестве.
Действительно, это следует из утверждения 1 и того, что нелинейный
оператор Bo(u) = u3(p) + E2i(p)u(p) ограничен и непрерывен вC[0, R].
Для расчета собственных значений и соответствующих им собственных функций предлагается два численных алгоритма.
Реализация первого из них предполагает использование итерационного процесса на основе интегрального представления собственной функции. При
использовании этого метода спектральный параметр у2 изначально не фиксируется, а для каждого его значения решается уравнение краевой задачи, решение которого затем подставляется в дисперсионное соотношение.
В формуле итерационного процесса
R
un+1 = (Po, у2) = а j G (р, Po)p(un + E2i(p)u„ )d р + f, n = 0,1,...
полагаем
»o(Po, у2) = f (Po) = Jpor
k2 J1( k2 R)
2
где у - решение g(у) = 0 .
где
и
Дисперсионное соотношение запишем в виде g (Я, у 2) = аЕ (Я, у2; и3),
g (Я, у2) = к2 ЯВДЯ) Jo(k2 Я) + ^ЯЗДЯ) Jl(k2 Я)
Я
Е (Я, у2; и3) = | й ррJl(k2P)(u3(р, у2) + £21(р)иИ (р, у2)).
0
Рассмотрим последовательность дисперсионных соотношений: g (Я,(у (и))2 ) = оЕ (Я,(у (и))2;и3).
Подстановка нулевой итерации
и0(р0,у2)
в приближенное дисперсионное соотношение дает уравнение
к 3 , Я
0(Я,у, а) = g(Я,у2) - а-3—---------Г йррJl4 (к2р) = 0,
к2 •А (к2Я) 0
решениями которого является соответствующее нулевое приближение (у(0))2. Первое приближение Ы1(р0,(у(0))2) функции поля внутри волновода
2 (0) 2 ы(р0, у ) достигается путем подстановки нулевых приближений (у^ ') и
Ы0(р0,(у(0))2) в формулу итерационного процесса:
'><р0.<у(0))2) = ТТ-Я; +
к2 Jl(k2 Я)
+а—
2
я ( N (к Я) > Я
^Я) ] Jl(-2р0)Гйрр/1(-2р)[ы03 (р,(у(0))2 ) + Е21(р)Ы0 (р,(у(0))2
р0
-а^Л^1(-2р0) Г йррJl(k2р) ы0 (,(у(0))2 ) + £21(р)ы0 (р,(у(0))" 20
Я
а-2Jl(k2Po) Г йрр^1 (-2р) ы0 (р, (у(0) )2 ) + £21(р)Ы0 (р, (у(0) )2
р0
Итак, функции
и, ((у(0))2)(у(0))2)= §-|
являются первыми приближениями собственных функций соответственно внутри и вне волновода.
Второй численный метод решения подразумевает полное решение поставленной краевой задачи на собственные значения на сетке с заданной точностью е .
Сначала выбирается интервал поиска собственных значений, потом определяется шаг деления рассматриваемого интервала и вычисляются промежуточные собственные значения в каждом узле. Для каждого такого значения параметра решается итерационное уравнение, пока собственная функция не будет определена с требующейся точностью. После этого каждая пара собственного значения и соответствующей собственной функции подставляется в дисперсионное уравнение на предмет проверки знака последнего в узле. Заключительным процессом является нахождение интервала разбиения, на концах которого дисперсионное уравнение меняет знак. Полагаем корнем уравнения среднюю точку этого отрезка.
Пусть собственные значения у ищутся на отрезке [А,, А2] . Введем на
этом отрезке сетку с узлами у(у) = А, + у (А2 - А,)/ N, у = 0,..., N, где N удовлетворяет условию А2 - А, < N8 , если собственное значение у требуется найти с точностью 8 . Вычисляем значения функции Ф(у) = g(А) -аЕ(А)(А = у )
в узлах у(у), причем при каждом у(у) решаем интегральное уравнение (3) с помощью итерационного алгоритма (1) с требуемой точностью. Далее определяем перемену знака в последовательности чисел Ф(у(у)). Если Ф(у(у)) Ф(у(у +1)) < 0, то приближенно полагаем у = (у(у) + у(у+1))/2 .
Алгоритм метода состоит из следующих шагов:
1. Выбираем значения параметров: радиуса волновода, коэффициента нелинейности, диэлектрической проницаемости среды внутри и вне волновода. В соответствии с выбранными параметрами определяются границы интервала поиска собственных значений спектрального параметра:
[А1т;А2т ], таких что А1; < у2 < А2;, 1 = 1,..., т .
2. Определяем шаг поиска решений уг- на каждом интервале:
И = А21^И,
N
где N - число разбиений интервала.
3. Вычисляем значение параметра у;- в каждом узле:
л А2; -А, .
у 1- = А1 +к,
где к = 1,..., N.
4. Для каждого значения параметра у;- вычисляем соответствующее значение функции и— с заданной точностью по формуле итерационного процесса:
R
Un+1 (Ро, Y2) = a j G (Р, Ро)р(«П + е21 (р)«и )d Р + f, n = 0,1,...
о
5. Для каждого вычисленного значения Yik и Uk в k-м узле вычисляем значение дисперсионного соотношения
(g -aF) (Yik)(uik),
которое затем сравниваем с нулем. Затем фиксируем два соседних узла, где происходит перемена знака дисперсионного отношения.
6. Полагаем корнем y° среднюю точку отрезка, на концах которого была зафиксирована перемена знака.
7. Вычисляем соответствующее значение собственной функции
u0k = u0k (Y°c) итерациями с требующейся точностью.
Список литературы
1. Никольский, В. В. Электродинамика и распространение радиоволн /
B. В. Никольский. - М. : Наука, 1978.
2. Смирнов, Ю. Г. Распространение электромагнитных волн в цилиндрических диэлектрических волноводах, заполненных нелинейной средой / Ю. Г. Смирнов,
C. Н. Куприянова // Журнал вычислительной математики и математической физики. - 2004. - Т. 44. - № 10. - С. 1850-1860.
3. Stakgold, I. Green's Functions and Boundary Value Problems / I. Stakgold. - Wiley ; N.Y., 1998.
4. Корн, Г. Справочник по математике для научных работников и инженеров / Г. Корн, Т. Корн. - М. : Наука, 1968.
5. Треногин, В. А. Функциональный анализ / В. А. Треногин. - М. : Наука, 1993.
Куприянова Светлана Николаевна
кандидат физико-математических наук, доцент, кафедра математики и суперкомпьютерного моделирования, Пензенский государственный университет
E-mail: [email protected]
Kupriyanova Svetlana Nikolaevna Candidate of physical and mathematical sciences, associate professor, sub-department of mathematics and supercomputer modeling,
Penza state university
УДК 517.6 Куприянова, С. Н.
Численный метод решения нелинейной задачи на собственные значения для неоднородного волновода / С. Н. Куприянова // Известия высших учебных заведений. Поволжский регион. Физико-математические науки. -2010. - № 2 (14). - С. 76-84.
