Нахождение констант распространения методом Крылова при расчете мод фотонных волноводов Текст научной статьи по специальности «Математика»

НАХОЖДЕНИЕ КОНСТАНТ РАСПРОСТРАНЕНИЯ МЕТОДОМ КРЫЛОВА ПРИ РАСЧЕТЕ МОД ФОТОННЫХ ВОЛНОВОДОВ

Я. О. Шуюпова, В. В. Котляр Институт систем обработки изображений РАН Самарский государственный аэрокосмический университет им. С.П. Королева

Аннотация

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

Введение

Задача отыскания векторных мод волновода методом согласованных синусоидальных мод (ССМ-метод) [1-3, 5] состоит из нескольких этапов, среди которых наиболее вычислительно сложным является поиск констант распространения, связанный с решением нелинейной задачи на собственные значения матрицы большой размерности. В оригинальной работе по ССМ-методу [2] поиск констант распространения осуществлялся с помощью метода нахождения нулей (корней) некоторой функции, имеющей, в то же время, разрывы второго рода вблизи искомых нулей. Последнее обстоятельство приводит к пропуску нулей при их автоматическом обнаружении при разбиении области определения аргумента функции на конечное число отрезков.

В данной работе для поиска констант распространения мод волноводов применен итеративный метод Крылова решения нелинейной задачи на собственные значения [4]. Численно показана эффективность применения метода Крылова: найдено в два раза больше собственных значений при одинаковом разбиении на интервалы, чем обычным методом поиска нулей функции.

Основные положения ССМ-метода

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

ССМ-метод позволяет решать как скалярную, так и векторную задачи отыскания мод волновода, что делает его среди прочих известных методов [69] универсальным средством исследования, применимым и для слабонаправляющих, и для фотонных волноводов.

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

x(m) < x < x(m+1) n(x, y) = n("'m) = const, (n) " (n+n' (1)

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

Гт)(х, у) =

™ (2) = ХЖЧх^Су) +^)(х)^)(у)],'

р=е,И к=1

где Г представляет собой любую из электрических или магнитных компонент поля моды, а внешняя сумма соответствует суммированию по поляризациям: ТЕ- р = И, ТМ- р = е . Рассмотрим подробнее локальные моды, входящие в выражение (2). Локальная х-мода имеет вид:

u % Ч x) =

k (m) V V у

</ 'cos

[*£ (x - x(m))]-

,( m,l)

(3)

k ( m)

лxpk

■sin[k« (x-x(m>)],

Величины -Рр"-1 и Г^р"', входящие в выражение

(2), определяются через параметры волокна и локальные у-моды.

Локальные у-моды - есть непрерывные функции, заданные в каждой прямоугольной ячейке с номером п столбца разбиения с номером т выражениями:

?(т)

ф(y) = Ф^1 )cos [kyn) ((y - y(n)))]-

k( n.1)

kín)

■sin

[ kyn)(

y - y(n))],

y (y) =^sn'l)cos [k(n) (y -y(n))]+

№ sinl

k<n

in [kyn)(

y -y(n))],

(4)

(5)

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

Эти условия описываются нелинейной задачей на собственные значения относительно константы распространения кг:

Л(кг )и = 0, (6)

где матрица Л(к2), элементы которой нелинейно зависят от параметра кг, имеет достаточно сложную блочную структуру [1-3]. Размерность этой матрицы зависит от сложности системы однородных областей, описывающей сечение, и точности аппроксимации по формуле (2).

Метод Крылова решения нелинейных задач на собственные значения

На этапе решения нелинейной задачи на собственные значения (6) оригинальный метод согласованных синусоидальных мод [2] был модифицирован итеративным методом Крылова [4,5], позволяющим точно рассчитывать собственные значения даже для матриц большой размерности.

Интерполируя нелинейный матричный оператор Л(к2) между двумя произвольными точками ст и ц следующим образом:

к _ст ц_к

Л к) «Л(кг) = ^-Л(ц) + ^-^ Л(ст), (7)

ц _ ст ц _ ст

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

Л(кг )и = 0, которую можно представить в виде:

Цк+1 СТЛ(цк) +Цк Цк+' Л(ст)

(8)

Пк = 0, (9)

Цк _ст Цк _ст

где цк+2 _ приближенное значение искомого кг, получаемое на очередном шаге итерации. Введем обозначение:

0 =

Ц к+1 _Ц к

(10)

ц к+1 _ст

Тогда уравнение (9) можно представить как: [Л (ст)_1Л (ц к) _01]ик = 0. (11)

Пусть ц _ некоторое начальное приближение для собственного числа кг матрицы Л(кг), пусть ст _ некоторое фиксированное значение вблизи ц, тогда итерационная процедура уточнения собственного значения к2 состоит в следующем. На к _ м шаге решаем относительно 0 :

Л(ст)-1 Л(цк)и = 0и , (12)

рассчитываем новую оценку для собственного зна-

чения кг:

Итерации повторяются до тех пор, пока последовательность оценок {цк } не сойдется.

Таким образом, используя метод Крылова, можно отыскать собственные значения кг задачи (6), для каждого из которых можно сформировать по формуле (2) поля компонент векторной моды.

Комбинированный метод решения нелинейной задачи на собственные значения

В работах [1, 2] рассматривался другой метод решения нелинейной задачи на собственные значения (6). Назовем его методом нулей функции. Суть метода состоит в следующем. Выберем произвольный вектор V с ненулевыми компонентами, например, единичный. Будем решать уравнение

Л(кг )и' = V (14)

относительно и'. Для разных значений параметра к2 будем получать разные решения уравнения (14). Определим функцию

/(к2) = 1/и'р , (15)

где и' есть р _ ая компонента вектора и'. В окрестности искомых значений к2 функция /(кг) является непрерывной функцией скалярного аргумента, и ее нули есть искомые к2.

Поиск нулей /(кг) осуществляется стандартными методами математического пакета МайаЪ.

Однако разрывы функции /(к2), которые возникают в местах нулей функции и' р , существенно

затрудняют поиск констант распространения методом нулей. Поэтому даже при достаточно малом шаге дискретизации, значительно меньшем необходимого для отделения соседних нулей функции /(кг),

существует вероятность пропуска корней из-за разрывов расположенных вблизи нулей. Такая проблемная ситуация показана на рис. 1, где в районе корня кг = 6,9951 мкм-1 функция имеет разрыв.

ЯК)

140 120 100 80 60 40 20 0 -20

Г

0

цк+1 =цк +— (цк _ст). 1 _ 0

(13)

6,992 6,994

Рис. 1. Пример разрыва функции / (кг ) вблизи нуля

6,996 6,998 к, мкм1

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

Рассмотрим возможности совместного использования метода (14)-(15) с методом Крылова (7)-(13). Пусть метод Крылова используется для индикации наличия корня на некотором интервале, где дальнейшее уточнение значений корней ведется при помощи функции (15). При заданной точности разделения корней, выберем шаг дискретизации в методе Крылова таким образом, чтобы минимизировать время расчета

t(hk) = Mtki + Intsi ^ min. (16)

В (16) M = L/hk - количество интервалов, анализируемых по методу Крылова, где:

L - длина интервала, на котором осуществляется поиск собственных значений,

tki - время на выполнение одной итерации метода Крылова,

I - предполагаемое число корней на интервале длины L ,

n = hk/h - число интервалов, рассматриваемых на интервале длины hk при уточнении корней с шагом h ,

tsi - среднее время на выполнение одной итерации обнаружения нуля функции f (kz) на интервале длины h .

Оптимальным значением шага дискретизации hk является значение:

hopt = \Lhtki

Its

(17)

Расчет констант распространения собственных мод фотонного волновода различными методами

Для модели фотонного волновода, показанной на рис. 2 проводился расчет констант распространения нескольких собственных мод.

123456789

X, мкм

Рис. 2. Модель фотонного волновода светлым показаны области со значением показателя преломления

Для аппроксимации по формуле (2) использовалось 30 локальных мод, что, согласно исследованию, проведенному в работе [1], достаточно для получения оценок констант распространения мод с абсолютной ошибкой не более 10-4.

В соответствии с ССМ-методом была получена задача (6), где элементы матрицы Л(кг) нелинейно

зависят от параметра к2, искомые значения которого

обращают Л(кг) в вырожденную числовую матрицу.

Для заданной структуры сечения и количества локальных мод равного 30, размерность матрицы Л(кг) составила 480 х 480 .

На интервале единичной длины поиск констант распространения осуществлялся тремя способами: методом определения нулей функции (15), методом Крылова и комбинированным методом с оптимальным шагом дискретизации.

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

12

8

4 0

□ Комбинированный метод

□ Метод нулей функции Ш Метод Крылова

JZL

ГШ

ш

темным - с n

(m,n) _

= 1

2 5 10 20 50 100 200 500 1000

Рис. 3. Диаграмма зависимости количества найденных констант распространения нескольких первых мод фотонного волновода от величины обратно пропорциональной точности разделения корней 1/ h

Как следует из диаграммы, метод Крылова существенно лучше двух других обнаруживает искомые значения параметра kz, но при этом, как видно из графика на рис. 4, является самым затратным по времени при заданном значении h . Так, при расчете 13 корней методом Крылова с h = 10-3 на персональном компьютере Intel Pentium IV с частотой 3,2 Ггц и оперативной памятью 2 Гб было затрачено около 3 часов. Заметим (рис. 3), что метод Крылова находил те же 13 корней , но за время в 8 раз меньшее. В тоже время метод нулей функции и комбинированный метод имеют примерно одинаковую частоту обнаружения корней. Незначительное преимущество комбинированного метода в области малой точности разделения корней объясняется случайно лучшим позиционированием интервала длины h , на котором производятся вторичные итерации метода нулей функции, относительно корня. При этом комбинированный метод является более эффективным по времени, чем метод нулей функции, как следует из рис. 4.

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

Г" ___—- i—1 -i

у0"

w f

» Комбинированный метод -А- Метод нулей функции Метод Крылова

0

200

400

600

800 1000 1Лг, мкм'1

Рис. 4. График зависимости натурального логарифма времени в секундах, затраченного на расчет нулей, от величины обратно пропорциональной точности разделения корней 1/ И

Благодарности

Работа выполнена при поддержке российско-американской программы «Фундаментальные ис-

следования и высшее образование» (грант CRDF

RUX0-014-SA-06), а также грантов РФФИ 05-0850298 и 07-07-97601.

Литература

1. Котляр В.В., Шуюпова Я.О. Расчет пространственных мод оптических волноводов с неоднородным поперечным сечением методом согласованных синусоидальных мод // Компьютерная оптика, 2003. № 25 P. 41-48.

2. Sudbo A.S. Film mode matching: a versatile method for mode film calculations in dielectric waveguides // Pure Appl. Opt., 1993. V. 2. P. 211-233.

3. Sudbo A.S. Improved formulation of the film mode matching method for mode film calculations in dielectric waveguides // Pure Appl. Opt., 1994. V. 3. Рр. 381-388.

4. Ruhe A., Zapiski Nauchnih Seminarov, Steklov Mathematical Institute, 2000. Pp.176 - 180.

5. Котляр В.В., Шуюпова Я.О. Расчет векторных мод оптического волновода // Компьютерная оптика, 2005. № 27 P. 89-94

6. Hadley G.R., Smith R.E. Full-vector waveguide modeling using an iterative finite-difference method with transparent boundary conditions // J. Lighhtwave Techn., 1995. V. 13. N. 3. P. 465-469.

7. Rogge U., Pregla R. Method of lines for the analysis of dielectric waveguides // J. Lightwave Techn., 1993. V. 11. P. 2015-2020.

8. Sztefka G., Nogling H.P. Bidirectional eigenmode propagation for large refractive index steps // IEEE Photonics Techn. Left., 1993. V. 5. P. 554-557.

9. Rahman B.M.A., Davies J.B. Finite-Elements Solution of integrated optical waveguides // J. Lightwave Techn., 1984. V. 2. P. 682-687.