Программирование, 2023, № 2, стр. 46-53

СИМВОЛЬНО-ЧИСЛЕННАЯ РЕАЛИЗАЦИИ МЕТОДА ГАЛЕРКИНА ДЛЯ ПРИБЛИЖЕННОГО РЕШЕНИЯ ЗАДАЧИ ВОЛНОВОДНОЙ ДИФРАКЦИИ

Д. В. Диваков^a, b, *, А. А. Тютюнник^a, b, **

^a Российский университет дружбы народов
117198 Москва, ул. Миклухо-Маклая, д. 6, Россия

^b Объединенный институт ядерных исследований
141980 Московская обл., г. Дубна, ул. Жолио-Кюри, д. 6, Россия

^* E-mail: divakov-dv@rudn.ru
^** E-mail: tyutyunnik-aa@rudn.ru

Поступила в редакцию 31.08.2022
После доработки 16.10.2022
Принята к публикации 30.10.2022

EDN: MFUGDU
DOI: 10.31857/S0132347423020097

Полный текст (PDF)

Аннотация

В работе построена символьно-численная реализация метода Галеркина для приближенного решения задачи волноводной дифракции на стыке двух открытых планарных трехслойных волноводов. Метод Гелеркина реализован в системе компьютерной алгебры Maple с использованием символьных манипуляций, основа программной реализации – символьно-численная процедура scprod, реализующая численный расчет скалярных произведений метода Галеркина на основе символьных выражений. Использование символьных манипуляций позволяет ускорить расчет интегралов в методе Галеркина благодаря однократному символьному расчету типовых для задачи интегралов вместо многократного численного интегрирования.

1. ВВЕДЕНИЕ

В работе рассматривается построение символьно-численной реализации метода Галеркина в системе компьютерной алгебры для приближенного решения задачи волноводной дифракции на стыке двух открытых планарных трехслойных волноводов. Самая ресурсозатратная процедура в методе Галеркина – это вычисление большого числа скалярных произведений элементов функционального пространства. Каждое скалярное произведение представляет собой интеграл от произведения двух функций, причем для планарных открытых волноводов эти функции являются кусочно-гладкими собственными функциями задачи Штурма–Лиувилля, заданными в символьно-численном виде: в символьном виде задана фундаментальная система решений (ФСР) дифференциального уравнения с постоянным коэффициентом в каждом слое волновода, а в численном виде заданы только коэффициенты разложения по ФСР. Различные коэффициенты при ФСР определяют разные собственные функции. Если численные константы заменить на символьные, то интегралы от произведения собственных функций можно вычислить однократно в каждом слое волновода в виде символьного выражения, а после подставлять в полученные символьные выражения численные значения констант. В таком подходе вычисление интегралов сведено к минимуму для ускорения численных расчетов.

Описанный подход реализован в системе компьютерной алгебры Maple для решения задачи волноводной дифракции на стыке двух открытых планарных трехслойных волноводов. Собственные функции вычисляются в символьно-численном виде методом волнового сопряжения, который реализован в виде процедуры wave_conj в Maple. Используя собственные функции применяем процедуру scprod для численного расчета матрицы скалярных произведений метода Галеркина. Матрица скалярных произведений используется для формирования системы уравнений, которая решается численно в Maple методом LinearSolve пакета LinearAlgebra.

С помощью разработанного символьно-численного метода исследована сходимость метода Галеркина для задачи волноводной дифракции.

2. МОТИВАЦИЯ

Корректная постановка задачи волноводной дифракции на стыке двух планарных открытых волноводов может быть сформулирована следующим образом. Имеется стык двух планарных открытых волноводов, который для корректной постановки задачи помещаем в объемлющий закрытый волновод аналогично подходу [1].

В таком случае компонента ${{E}_{y}}\left( {x,z} \right) = u\left( {x,z} \right)$ удовлетворяет уравнению Гельмгольца [2, 3]

(1)

$\left( {\frac{{{{\partial }^{2}}}}{{\partial {{x}^{2}}}} + \frac{{{{\partial }^{2}}}}{{\partial {{x}^{2}}}} + k_{0}^{2}{{n}^{2}}\left( {x,z} \right)} \right)u \equiv (\Delta + k_{0}^{2}{{n}^{2}})u = 0,$

условиям на границе объемлющего закрытого волновода [2, 3]

(2)

${{\left. u \right|}_{{x = \pm L}}} = 0,$

и условиям сопряжения на стыке волноводов [4]

(3)

${{\left. {\left[ u \right]{\kern 1pt} } \right|}_{{z = 0}}} = 0,\quad {{\left. {\left[ {\frac{{\partial u}}{{\partial z}}} \right]{\kern 1pt} } \right|}_{{z = 0}}} = 0,$

где ${{\left. {\left[ f \right]} \right|}_{{s = a}}} = {{\left. f \right|}_{{s = a - 0}}} - {{\left. f \right|}_{{s = a + 0}}}$. Кусочно-постоянный коэффициент n²(x, z) в уравнении Гельмгольца определяет показатель преломления в рассматриваемой области: ${{\left. {{{n}^{2}}(x,z)} \right|}_{{z < 0}}} = n_{1}^{2}(x),$ ${{\left. {{{n}^{2}}(x,z)} \right|}_{{z > 0}}} = n_{2}^{2}(x)$, где

(4)

$n_{l}^{2}\left( x \right) = \left\{ \begin{gathered} n_{c}^{2},\quad x > {{h}_{l}} \hfill \\ n_{f}^{2},\quad 0 < x < {{h}_{l}} \hfill \\ n_{s}^{2},\quad x < 0 \hfill \\ \end{gathered} \right..$

Применяя метод разделения переменных к уравнению Гельмгольца в области $z < 0$ и в области $z > 0$ получаем общий вид решения в каждой из областей в виде разложения по нормальным модам

(5)

$\begin{gathered} {{\left. u \right|}_{{z < 0}}} = \sum\limits_{j = 0}^\infty {{A}_{j}}{{\varphi }_{j}}\left( x \right){{e}^{{i{{k}_{0}}{{\beta }_{j}}z}}} + \sum\limits_{j = 0}^\infty {{R}_{j}}{{\varphi }_{j}}\left( x \right){{e}^{{ - i{{k}_{0}}{{\beta }_{j}}z}}}, \\ {{\left. u \right|}_{{z > 0}}} = \sum\limits_{j = 0}^\infty {{T}_{j}}{{\psi }_{j}}\left( x \right){{e}^{{i{{k}_{0}}{{\gamma }_{j}}z}}} + \sum\limits_{j = 0}^\infty {{S}_{j}}{{\psi }_{j}}\left( x \right){{e}^{{ - i{{k}_{0}}{{\gamma }_{j}}z}}}, \\ \end{gathered} $

где ${{\varphi }_{j}}\left( x \right)$ и ${{\beta }_{j}}$ ортонормированные собственные функции и собственные значения задачи Шурма–Лиувилля при l = 1

(6)

$\left\{ \begin{gathered} {v}{\kern 1pt} '{\kern 1pt} '\; + k_{0}^{2}n_{l}^{2}\left( x \right){v} = k_{0}^{2}{{\eta }^{2}}{v}, \hfill \\ {{\left[ {v} \right]}_{{x = 0,{{h}_{l}}}}} = 0,\quad {{\left[ {{v}{\kern 1pt} '} \right]}_{{x = 0,{{h}_{l}}}}} = 0, \hfill \\ {v}\left( { \pm L} \right) = 0, \hfill \\ \end{gathered} \right.$

а ${{\psi }_{j}}\left( x \right)$ и ${{\gamma }_{j}}$ – ортонормированные собственные функции и собственные значения приведенной выше задачи Шурма–Лиувилля при l = 2.

Постановка задачи дифракции формулируется следующим образом: в области $z < 0$ имеется одна падающая m-я мода с заданным амплитудным коэффициентом $A$ и счетное количество отраженных волноводных мод, характеризующихся амплитудными коэффициентами отражения ${{R}_{j}}$ [2]

(7)

${{\left. u \right|}_{{z < 0}}} = A{{\varphi }_{m}}\left( x \right){{e}^{{i{{k}_{0}}{{\beta }_{m}}z}}} + \sum\limits_{j = 0}^\infty {{R}_{j}}{{\varphi }_{j}}\left( x \right){{e}^{{ - i{{k}_{0}}{{\beta }_{j}}z}}},$

а в области z > 0 имеется счетное количество прошедших волноводных мод, характеризующихся амплитудными коэффициентами прохождения ${{T}_{j}}$ [2]

(8)

${{\left. u \right|}_{{z > 0}}} = \sum\limits_{j = 0}^\infty {{T}_{j}}{{\psi }_{j}}\left( x \right){{e}^{{i{{k}_{0}}{{\gamma }_{j}}z}}}.$

Условие ${{\left. {\left[ u \right]} \right|}_{{z = 0}}} = 0$, учитывая вид решения (7) и (8), принимает следующий вид

(9)

$\sum\limits_{j = 0}^\infty {{R}_{j}}{{\varphi }_{j}}\left( x \right) - \sum\limits_{j = 0}^\infty {{T}_{j}}{{\psi }_{j}}\left( x \right) = - A{{\varphi }_{m}}\left( x \right),$

а условие ${{\left. {\left[ {\frac{{\partial u}}{{\partial z}}} \right]} \right|}_{{z = 0}}} = 0$ с учетом (7) и (8) принимает следующий вид

(10)

$\sum\limits_{j = 0}^\infty {{\beta }_{j}}{{R}_{j}}{{\varphi }_{j}}\left( x \right) + \sum\limits_{j = 0}^\infty {{\gamma }_{j}}{{T}_{j}}{{\psi }_{j}}\left( x \right) = {{\beta }_{m}}A{{\varphi }_{m}}\left( x \right).$

Задача волноводной дифракции на стыке двух планарных открытых волноводов, помещенных в объемлющий закрытый волновод, состоит в определении амплитудных коэффициентов отражения ${{R}_{j}}$ и прохождения ${{T}_{j}}$ по заданным характеристикам волноводного перехода, а также по заданной амплитуде A падающей m-й моды.

3. МЕТОДЫ

3.1. Приближенный метод решения задачи волноводной дифракции

Систему из двух уравнений (9), (10) решаем приближенно, для этого в каждой из бесконечных сумм оставляем конечное число слагаемых: в каждой сумме оставляем только те слагаемые, для которых выполняется ${{\beta }_{j}} > 0$ в первой сумме и ${{\gamma }_{j}} > 0$ во второй сумме соответственно, то есть неэванесцентные моды. К усеченным суммам применяем метод Галеркина [5] и получаем систему линейных алгебраических уравнений следующего вида

(11)

$\left\{ \begin{gathered} \vec {r} - \Phi \vec {t} = - {{{\vec {e}}}_{m}}, \hfill \\ {{D}_{1}}\vec {r} + \Phi {{D}_{2}}\vec {t} = {{D}_{1}}{{{\vec {e}}}_{m}}, \hfill \\ \end{gathered} \right.$

где $\vec {r} = {{\left( {{{R}_{0}},{{R}_{1}}, \cdots ,{{R}_{N}}} \right)}^{T}}$ и $\vec {t} = {{\left( {{{T}_{0}},{{T}_{1}}, \cdots ,{{T}_{M}}} \right)}^{T}}$ – векторы искомых величин, ${{D}_{1}} = {\text{diag}}\left\{ {{{\beta }_{j}}} \right\}_{{j = 0}}^{N}$, D₂ = diag$\left\{ {{{\gamma }_{j}}} \right\}_{{j = 0}}^{M}$, ${{\vec {e}}_{m}}$ – вектор, m-я компонента которого равна A, а остальные компоненты нулевые. Матрица $\Phi $ определена следующим образом

(12)

$\Phi = \left( {\begin{array}{*{20}{c}} {\left( {{{\varphi }_{0}},{{\psi }_{0}}} \right)}&{\left( {{{\varphi }_{0}},{{\psi }_{1}}} \right)}& \cdots &{\left( {{{\varphi }_{0}},{{\psi }_{M}}} \right)} \\ {\left( {{{\varphi }_{1}},{{\psi }_{0}}} \right)}&{\left( {{{\varphi }_{1}},{{\psi }_{1}}} \right)}& \cdots &{\left( {{{\varphi }_{1}},{{\psi }_{M}}} \right)} \\ \vdots & \vdots &{}& \vdots \\ {\left( {{{\varphi }_{N}},{{\psi }_{0}}} \right)}&{\left( {{{\varphi }_{N}},{{\psi }_{1}}} \right)}& \cdots &{\left( {{{\varphi }_{N}},{{\psi }_{M}}} \right)} \end{array}} \right),$

где $\left( {{{\varphi }_{p}},{{\psi }_{q}}} \right) = \int_{ - L}^L {{{\varphi }_{p}}} \left( x \right){{\psi }_{q}}\left( x \right)dx$.

Численно решаем систему линейных алгебраических уравнений (11) и определяем приближенные значения амплитудных коэффициентов отражения и прохождения. Численное решение системы не представляет трудностей, более ресурсоемкая операция в данном случае – это численный расчет матрицы Φ, каждый элемент которой представляет собой интеграл от произведения двух собственных функций. Используем возможности системы компьютерной алгебры Maple по символьному расчету интегралов для символьно-численной реализации метода Галеркина чтобы сократить общее время расчетов. Для этого в первую очередь отыскиваем собственные функции задачи (6) в символьно-численном виде методом волнового сопряжения.

3.2. Символьно-численная реализация метода волнового сопряжения

Определяем собственные функции и собственные значения задачи Штурма–Лиувилля (6) методом волнового сопряжения: в каждой подобласти с постоянным показателем преломления записываем общее решение соответствующего обыкновенного дифференциального уравнения с постоянными коэффициентами

(13)

${v}_{\alpha }^{{''}} + k_{0}^{2}(n_{\alpha }^{2} - {{\eta }^{2}}){{{v}}_{\alpha }} = 0,\quad \alpha = c,f,s,$

в виде разложения по фундаментальной системе решений ${{{v}}_{\alpha }}\left( x \right) = A_{\alpha }^{1}{v}_{\alpha }^{1}\left( x \right) + A_{\alpha }^{2}{v}_{\alpha }^{2}\left( x \right)$ с неизвестными коэффициентами $A_{\alpha }^{{1,2}}$. Общее решение ${{{v}}_{\alpha }}\left( x \right)$ для $\alpha = c,f,s$ подставляем в граничные условия и условия сопряжения из (6). Получается однородная система линейных алгебраических уравнений относительно неизвестных величин $A_{\alpha }^{{1,2}}$, $\alpha \, = \,c,f$, s, коэффициенты которой зависят от спектрального параметра ${{\eta }^{2}}$

(14)

$M({{\eta }^{2}})\vec {A} = \vec {0}.$

Система (14) имеет нетривиальное решение только тогда, когда $\det M({{\eta }^{2}}) = 0$. Каждому корню $\eta _{j}^{2}$ этого уравнения соответствует однородная система ${{M}_{j}}{{\vec {A}}_{j}} = \vec {0}$, где ${{M}_{j}} = M(\eta _{j}^{2})$. Решив систему, то есть отыскав вектор ${{\vec {A}}_{j}}$, мы определяем искомые коэффициенты j-й моды $A_{{\alpha ,j}}^{{1,2}}$, $\alpha = c,f,s$ тем самым полностью определяя и собственную функцию задачи, соответствующую собственному значению $\eta _{j}^{2}$.

Численное отыскание корней нелинейного уравнения $\det M({{\eta }^{2}}) = 0$ осуществляется методом деления отрезка пополам. Приближенное решение системы ${{M}_{j}}{{\vec {A}}_{j}} = \vec {0}$ с численно вырожденной матрицей M_j осуществляется с помощью численного отыскания собственного вектора матрицы M_j, отвечающего почти нулевому по модулю собственному значению. Алгоритм реализован в численном виде в системе компьютерной алгебры Maple с использованием метода Eigenvectors библиотеки LinearAlgebra [6], псевдокод алгоритма приведен ниже.

Algorithm 1: Algorithm wave_conj for numeric calculation ${{\eta }_{j}}$, ${{\vec {A}}_{j}}$

Input:${v}_{c}^{1}(x)$, ${v}_{c}^{2}(x)$, ${v}_{s}^{1}(x)$, ${v}_{s}^{2}(x)$, ${v}_{f}^{2}(x)$ – symbolic expressions for each layer defined in eqs (15), (16), $h$ – numeric value for thickness of the waveguide layer, $L$ – numeric value for thickness of the “Dirichlet box”, ${{n}_{c}}$, ${{n}_{s}}$, ${{n}_{f}}$ – numeric values for refractive indices
Output: array of ${{\eta }_{j}}$, array of ${{\vec {A}}_{j}}$
generate matrix for $\eta \in ({{n}_{s}};{{n}_{f}})$
`GenerateMatrix`${{M}_{1}}(\eta )$
using ${v}_{c}^{1}(x)$, ${v}_{f}^{2}(x)$, ${v}_{s}^{1}(x)$ and conditions (6)
generate matrix for $\eta \in ({{n}_{c}};{{n}_{s}})$
`GenerateMatrix`${{M}_{2}}(\eta )$
using ${v}_{c}^{1}(x)$, ${v}_{f}^{2}(x)$, ${v}_{s}^{2}(x)$ and conditions (6)
generate matrix for $\eta \in (0;{{n}_{c}})$
`GenerateMatrix`${{M}_{3}}(\eta )$
using ${v}_{c}^{2}(x)$, ${v}_{f}^{2}(x)$, ${v}_{s}^{2}(x)$ and conditions (6)
${{\eta }_{j}}$$ \leftarrow $`RootOf` (`Determinant`$({{M}_{1}}(\eta ))$),
$\eta \in ({{n}_{s}};{{n}_{f}})$, $j = 1,..,{{N}_{1}}$
${{\eta }_{j}}$$ \leftarrow $`RootOf` (`Determinant`$({{M}_{2}}(\eta ))$),
$\eta \in ({{n}_{c}};{{n}_{s}})$, $j = {{N}_{1}} + 1,..,{{N}_{2}}$
${{\eta }_{j}}$$ \leftarrow $`RootOf` (`Determinant`$({{M}_{3}}(\eta ))$),
$\eta \in (0;{{n}_{c}})$, $j = {{N}_{2}} + 1,..,{{N}_{3}}$
for$j = 1,..,{{N}_{1}}$do
		`e,v $ \leftarrow $ Eigenvectors $({{M}_{1}}({{\eta }_{j}}))$;`
		`indx_min $ \leftarrow $ min[index](abs$ \sim $(e));`
		`A_j${\kern 1pt} \leftarrow $ Column(v, indx_min);`
for$j = {{N}_{1}} + 1,..,{{N}_{2}}$do
		`e,v $ \leftarrow $ Eigenvectors $({{M}_{2}}({{\eta }_{j}}))$;`
		`indx_min $ \leftarrow $ min[index](abs$ \sim $(e));`
		`A_j${\kern 1pt} \leftarrow $ Column(v, indx_min);`
for$j = {{N}_{2}} + 1,..,{{N}_{3}}$do
		`e,v $ \leftarrow $ Eigenvectors $({{M}_{3}}({{\eta }_{j}}))$;`
		`indx_min`$ \leftarrow $`min[index](abs$ \sim $(e));`
		`A_j`${\kern 1pt} \leftarrow $`Column(v, indx_min);`
return$\left[ {{{A}_{j}}} \right],\left[ {{{\eta }_{j}}} \right],{\text{ }}j = 1,...,{{N}_{3}}$

3.3. Символьно-численная реализация метода Галеркина

Общее решение уравнения (13) при $\eta > {{n}_{\alpha }}$ имеет следующий вид

(15)

${v}_{\alpha }^{1}\left( x \right) = A_{\alpha }^{1}{{e}^{{ - {{k}_{0}}\sqrt {{{\eta }^{2}} - n_{\alpha }^{2}} x}}} + A_{\alpha }^{2}{{e}^{{{{k}_{0}}\sqrt {{{\eta }^{2}} - n_{\alpha }^{2}} x}}},$

а при $\eta < {{n}_{\alpha }}$

(16)

${v}_{\alpha }^{2}(x)\, = \,A_{\alpha }^{1}{\text{sin}}({{k}_{0}}\sqrt {n_{\alpha }^{2} - {{\eta }^{2}}} x)\, + \,A_{\alpha }^{2}{\text{cos}}({{k}_{0}}\sqrt {n_{\alpha }^{2} - {{\eta }^{2}}} x).$

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

(17)

$\begin{gathered} \int\limits_{ - L}^L {{\varphi }_{p}}{{\psi }_{q}}dx = \int\limits_{ - L}^0 {{\varphi }_{{s,p}}}{{\psi }_{{s,q}}}dx + \int\limits_0^{{{h}_{1}}} {{\varphi }_{{f,p}}}{{\psi }_{{f,q}}}dx + \\ \, + \int\limits_{{{h}_{1}}}^{{{h}_{2}}} {{\varphi }_{{c,p}}}{{\psi }_{{f,q}}}dx + \int\limits_{{{h}_{2}}}^L {{\varphi }_{{c,p}}}{{\psi }_{{c,q}}}dx, \\ \end{gathered} $

где зависимость от переменной $x$ опущена для краткости записи. Отметим, что ${{\varphi }_{{\alpha ,p}}}\left( x \right)$ и ${{\psi }_{{\alpha ,q}}}\left( x \right)$, $\alpha = c,f,s$ принимают либо вид (15), либо вид (16) в зависимости от положительности/отрицательности выражений ${{\beta }_{p}} - {{n}_{\alpha }}$ и ${{\gamma }_{q}} - {{n}_{\alpha }}$. Поэтому для численного расчета любого интеграла вида (17) достаточно иметь только 3 символьных выражения для определенных интегралов:

(18)

$\begin{gathered} {{I}_{{1,1}}}\left( {{{A}_{1}},{{A}_{2}},{{B}_{1}},{{B}_{2}},a,b,\beta ,\gamma } \right) = \\ \, = \int\limits_a^b ({{A}_{1}}{{e}^{{ - \beta x}}} + {{A}_{2}}{{e}^{{\beta x}}})({{B}_{1}}{{e}^{{ - \gamma x}}} + {{B}_{2}}{{e}^{{\gamma x}}})dx, \\ \end{gathered} $

(19)

$\begin{gathered} {{I}_{{1,2}}}\left( {{{A}_{1}},{{A}_{2}},{{B}_{1}},{{B}_{2}},a,b,\beta ,\gamma } \right) = \\ \, = \int\limits_a^b ({{A}_{1}}{{e}^{{ - \beta x}}} + {{A}_{2}}{{e}^{{\beta x}}})\left( {{{B}_{1}}\sin \left( {\gamma x} \right) + {{B}_{2}}\cos \left( {\gamma x} \right)} \right)dx, \\ \end{gathered} $

(20)

$\begin{gathered} {{I}_{{2,2}}}\left( {{{A}_{1}},{{A}_{2}},{{B}_{1}},{{B}_{2}},a,b,\beta ,\gamma } \right) = \\ \, = \int\limits_a^b \left( {{{A}_{1}}\sin \left( {\beta x} \right) + {{A}_{2}}\cos \left( {\beta x} \right)} \right) \times \\ \, \times \left( {{{B}_{1}}\sin \left( {\gamma x} \right) + {{B}_{2}}\cos \left( {\gamma x} \right)} \right)dx. \\ \end{gathered} $

Используем возможности системы компьютерной алгебры Maple, а именно функцию int для символьного интегрирования, а также функции factor, expand для упрощения символьных выражений [6]. Полученные в системе компьютерной алгебры символьные выражения ${{I}_{{1,1}}},\;{{I}_{{1,2}}},\;{{I}_{{2,2}}}$ используются для расчета интегралов по формуле (17).

Описанные символьно-численные методы реализованы в системе компьютерной алгебры Maple.

Algorithm 2: Algorithm ipq for numeric calculation of integral (18), (19) or (20)

Input:$\beta $, $\gamma $ – numeric values for eigenvalues, ${{n}_{1}}$, ${{n}_{2}}$ – numeric values for refractive indices, a, b – numeric values for integration limits, ${{A}_{{1,2}}}$, ${{B}_{{1,2}}}$ – numeric values for amplitudes
Output: numeric value of sub-integral
if$\beta > {{n}_{1}}$$\& $$\gamma > {{n}_{2}}$then
		res $ \leftarrow $ subs numeric values of
		${{A}_{{1,2}}},{{B}_{{1,2}}},a,b,\beta ,\gamma $ in ${{I}_{{1,1}}}$
if$\beta < {{n}_{1}}$$\& $$\gamma < {{n}_{2}}$then
		res $ \leftarrow $ subs numeric values of
		${{A}_{{1,2}}},{{B}_{{1,2}}},a,b,\beta ,\gamma $ in ${{I}_{{2,2}}}$
if$\beta < {{n}_{1}}$$\& $$\gamma > {{n}_{2}}$then
		res $ \leftarrow $ subs numeric values of
		${{A}_{{1,2}}},{{B}_{{1,2}}},a,b,\beta ,\gamma $ in ${{I}_{{1,2}}}$
if $\beta > {{n}_{1}}$$\& $$\gamma < {{n}_{2}}$then
		res $ \leftarrow $ subs numeric values of
		${{B}_{{1,2}}},{{A}_{{1,2}}},a,b,\beta ,\gamma $ in ${{I}_{{1,2}}}$
return res

3.4. Реализация символьно-численных алгоритмов в Maple

Реализация символьно-численного метода решения волноводной задачи состоит из 3-х символьно-численных алгоритмов, реализованных в системе компьютерной алгебры Maple. Метод волнового сопряжения, описанный в разделе 3.2, реализован в системе компьютерной алгебры Maple в виде процедуры wave_conj, псевдокод которой приведен выше, см. алг. 1. Символьная часть метода состоит в формировании матрицы ${{M}_{{1,2,3}}}(\eta )$ для каждого из трех диапазонов $\eta \, \in \,({{n}_{s}},{{n}_{f}})$, $\eta \, \in \,({{n}_{c}},{{n}_{s}})$ и $\eta \, \in \,(0,{{n}_{c}})$. Численно отыскиваем значения ${{\beta }_{j}}$, решая уравнение ${\text{det}}{{M}_{{1,2,3}}}(\beta )$ = 0 и значения ${{\vec {A}}_{j}}$, решая систему (14) для левого волновода и, аналогично, ${{\gamma }_{j}}$ и ${{\vec {B}}_{j}}$ для правого волновода.

Далее численные значения ${{\beta }_{j}}$, ${{\vec {A}}_{j}}$ и ${{\gamma }_{j}}$, ${{\vec {B}}_{j}}$ в качестве аргументов передаются в процедуру scprod (см. алг. 3), которая вычисляет матрицу скалярных произведений (12), опираясь на заранее полученные символьные выражения для интегралов типа (18), (19) и (20) с помощью процедуры ipq (см. алг. 2).

После численного расчета матрицы (12) решаем систему (11) с помощью метода LinearSolve библиотеки LinearAlgebra.

Algorithm 3: Algorithm scprod for numeric calculation matrix defined in (12)

Input:${{\beta }_{j}}$, ${{A}_{j}}$, $j = 1,...,N$ and ${{\gamma }_{j}}$, ${{B}_{j}}$, $j = 1,...,M$ – numeric values for eigenpairs, calculated by `wave_conj` algorithm, ${{h}_{1}}$, ${{h}_{2}}$, $L$ – numeric values for thicknesses, ${{n}_{c}}$, ${{n}_{s}}$, ${{n}_{f}}$ – numeric values for refractive indices
Output: Matrix $F$ of scalar product defined in (12)
forp = 1, …,Ndo
	forq = 1, …, Mdo
		$I_{{pq}}^{{(cc)}} \leftarrow $`ipq`(${{A}_{1}}$, ${{A}_{2}}$, ${{B}_{1}}$, ${{B}_{2}}$, ${{h}_{2}}$, $L$, ${{n}_{c}}$, ${{n}_{c}}$, ${{\beta }_{p}}$, ${{\gamma }_{q}}$)
		$I_{{pq}}^{{(cf)}} \leftarrow $`ipq`(${{A}_{1}}$, ${{A}_{2}}$, ${{B}_{3}}$, ${{B}_{4}}$, ${{h}_{1}}$, ${{h}_{2}}$, ${{n}_{c}}$, ${{n}_{f}}$, ${{\beta }_{p}}$, ${{\gamma }_{q}}$)
		$I_{{pq}}^{{(ff)}} \leftarrow $`ipq`(${{A}_{3}}$, ${{A}_{4}}$, ${{B}_{3}}$, ${{B}_{4}}$, $0$, ${{h}_{1}}$, ${{n}_{f}}$, ${{n}_{f}}$, ${{\beta }_{p}}$, ${{\gamma }_{q}}$)
		$I_{{pq}}^{{(ss)}} \leftarrow $`ipq`(${{A}_{5}}$, ${{A}_{6}}$, ${{B}_{5}}$, ${{B}_{6}}$, $ - L$, $0$, ${{n}_{s}}$, ${{n}_{s}}$, ${{\beta }_{p}}$, ${{\gamma }_{q}}$)
		${{F}_{{pq}}} \leftarrow I_{{pq}}^{{(cc)}} + $$I_{{pq}}^{{(cf)}}$ + $I_{{pq}}^{{(ff)}} + I_{{pq}}^{{(ss)}}$
returnF

4. РЕЗУЛЬТАТЫ

Используем разработанный символьно-численный алгоритм для приближенного решения задачи дифракции на стыке двух волноводов методом Галеркина. Исследуем для начала точность метода Галеркина в зависимости от параметра L, для этого оцениваем невязку $\delta $, определяющую баланс энергии

(21)

$\delta = \sqrt {\sum\limits_{j = 0}^N {{{\left| {{{R}_{j}}} \right|}}^{2}} + \sum\limits_{j = 0}^M {{{\left| {{{T}_{j}}} \right|}}^{2}}} - 1.$

В зависимости от величины $L$ изменяется число базисных функций, а именно $N = 52$, $M = 52$ при $L = 5\lambda $, $N = 152$, $M = 151$ при $L = 15\lambda $ и $N = 248$, $M = 249$ при $L = 25\lambda $. Для дальнейших расчетов выбираем значение $L = 25\lambda $, при котором величина $\delta $ достигаем минимума. Приведем далее численные результаты решения задачи дифракции на стыке двух волноводов, где толщина ${{h}_{1}} = 1.5\lambda $ фиксирована, а толщина второго волновода h₂ изменяется от $1.5\lambda $ до $3.7\lambda $ при $\lambda = 0.55$ μm, показатели преломления слоев определны следующим образом ${{n}_{c}} = 1.0$, ${{n}_{s}} = 1.47$, ${{n}_{f}}$ = 1.565. Оцениваем численно изменение абсолютных значений амплитуд прошедших волноводных мод $\left| {{{T}_{0}}} \right|$ и $\left| {{{T}_{1}}} \right|$.

Величины абсолютных значений амплитуд прошедших волноводных мод $\left| {{{T}_{0}}} \right|$ и $\left| {{{T}_{1}}} \right|$ характеризуют перераспределение энергии между первыми двумя волноводными модами. Численно оценим также суммарную интенсивность ${{E}_{1}}$, которая приходится на моды с ${{n}_{s}} < {{\gamma }_{q}} < {{n}_{f}}$ и суммарную интенсивность ${{E}_{2}}$, приходящуюся на моды с $0 < {{\gamma }_{q}} < {{n}_{s}}$. Величины ${{E}_{1}}$ и ${{E}_{2}}$ определены следующим образом

(22)

5. ОБСУЖДЕНИЕ

В работе построена символьно-численная реализация метода Галеркина для решения задачи волноводной дифракции на стыке двух планарных трехслойных волноводов. Багодаря тому, что собственные функции задачи (6) вычисляются в символьно-численном виде можно применить метод волнового сопряжения, который реализован в виде процедуры wave_conj в Maple. В системе Maple также реализована процедура scprod для численного расчета матрицы скалярных произведений метода Галеркина на основе вычисленных собственных функций. Процедура scprod использует рассчитанные и упрощенные формулы для типовых интегралов (18)–(20), возникающих в задаче дифракции. Другими словами процедура интегрирования осуществляется три раза для расчета интегралов (18)–(20), далее полученные символьные выражения упрощаются и используются в численных расчетах.

Важно отметить, что несмотря на использование символьных манипуляций, матрица скалярных произведений $\Phi $ из (12) вычисляется приближенно, из-за того, что собственные функции ${{\varphi }_{p}}$ и ${{\psi }_{q}}$ определяются приближенно. Каждый элемент диагональных матриц собственных значений D_{1, 2} системы (11) известен с точностью 10^–15 – так как собственные значения являются корнями нелинейного уравнения, которые без труда отыскиваются методом деления отрезка пополам с точностью 10^–15. При такой точности вычисления собственных значений собственные функции отыскиваются с точностью $a \times {{10}^{{ - 13}}}$, где a < 10 – см. например [8]. Для расчета матричных элементов $\Phi $ из (12) необходимо вычислить интеграл от произведения ортонормированных функций ${{\varphi }_{p}}$ и ${{\psi }_{q}}$, ошибка при расчете такого интеграла будет порядка $2a \times {{10}^{{ - 13}}}$.

Символьно-численная реализация метода Галеркина позволяет однократно вычислить однотипные интегралы, ускоряя тем самым последующие численные расчеты, что особенно актуально для векторных задач – см. [9]. В настоящей работе рассматривается скалярная задача, которую, однако, необходимо решать многократно, как для исследования сходимости метода Галеркина (см. рис. 2), так и для исследования самой задачи дифракции на различных стыках. Стоит отметить, что использование систем компьютерной алгебры характерно не только для методов Галеркина [10] и Ритца [11], но и для других приближенных методов решения волноводных задач [12].

Рис. 1.

Стык двух планарных открытых волноводов, помещенных в объемлющий закрытый волновод с границами $x = \pm L$.

Рис. 2.

График сходимости величины $\delta $ для стыка с параметрами ${{h}_{1}} = 1.5w,$ ${{h}_{2}} = 1.8w,$ ${{n}_{c}} = {\text{1}}{\text{.0,}}$ ${{n}_{s}} = {\text{1}}{\text{.47,}}$ ${{n}_{f}} = {\text{1}}{\text{.565}}$, $A = 1$, m = 0, где длина волны $\lambda = 0.55$ μm.

Рис. 3.

График изменения величин $\left| {{{T}_{0}}} \right|$ и $\left| {{{T}_{1}}} \right|$ в зависимости от толщины ${{h}_{2}}$ для стыка с параметрами ${{h}_{1}} = 1.5\lambda $, ${{n}_{c}} = {\text{1}}{\text{.0}}$, ${{n}_{s}} = {\text{1}}{\text{.47}}$, ${{n}_{f}} = {\text{1}}{\text{.565}}$, $A = 1$, m = 0, где длина волны $\lambda = 0.55$ μm.

Корректная постановка задачи дифракции основана на использовании подхода, предложенного А.Г. Свешниковым для задач волноводной дифракции [7], схожему с подходом О.И. Толстихина для задач квантовой механики [1]. Подход состоит в помещении рассматриваемого волноводного стыка открытых волноводов в объемлющий закрытый волновод (“ящик Дирихле”), границы которого $x = \pm L$ удалены от рассматриваемой структуры на “достаточно большое” расстояние. В работе величина этого расстояния $L$ определена из численного эксперимента: исследовано поведение невязки $\delta $, описывающей численное “нарушение” закона сохранения энергии (21). Полученные результаты (см. рис. 2) демонстрируют немонотонное поведение невязки, которое характерно для метода Галеркина в волноводных задачах – см. [9]. Получено значение $L = 25\lambda $, обеспечивающее минимальное значение невязки $\delta $, где $\lambda $ – длина волны.

Использование “ящика Дирихле” позволяет аппроксимировать непрерывный спектр исходной задачи для открытых планарных волноводов [13–15] дискретным спектром благодаря ввдению искуственных границ $x = \pm L$. Этот подход позволяет сформулировать корректную задачу дифракции и применить для ее исследования метод Галеркина, а также приближенно оценить какое количество энергии приходится на моды непрерывного спектра исходной задачи (см. рис. 4, 5 ). На рис. 4 представлен график величины ${{E}_{1}}$, которая описывает суммарную интенсивность, приходящуюся на моды, аппроксимирующие дискретный спектр исходной задачи в зависимости от толщиты ${{h}_{2}}$. Когда ${{h}_{2}} = {{h}_{1}}$ стык отсутствует, поэтому $\left| {{{T}_{0}}} \right| = \left| A \right| = 1$ (см. рис. 3) – единственная волноводная мода, налетающая на стык из области z < 0, продолжает распространяться вдоль волновода с неизменной амплитудой. В данном случае перераспределение энергии отсутствует, о чем свидетельствует значение величин ${{E}_{1}} = 1$ на рисунке 4. Величина ${{E}_{2}}$ на рисунке 5 описывает суммарную интенсивность, приходящуюся на моды, аппроксимирующие непрерывный спектр исходной задачи в зависимости от толщиты ${{h}_{2}}$. Когда ${{h}_{2}} = {{h}_{1}}$ стык отсутствует, поэтому ${{E}_{2}} = 0$.

Рис. 4.

График изменения величины ${{E}_{1}}$ в зависимости от толщины ${{h}_{2}}$ для стыка с параметрами ${{h}_{1}} = 1.5\lambda $, ${{n}_{c}} = {\text{1}}{\text{.0}}$, ${{n}_{s}} = {\text{1}}{\text{.47}}$, ${{n}_{f}} = {\text{1}}{\text{.565}}$, $A = 1$, m = 0, где длина волны $\lambda = 0.55$ μm.

Рис. 5.

График изменения величины ${{E}_{2}}$ в зависимости от толщины ${{h}_{2}}$ для стыка с параметрами ${{h}_{1}} = 1.5\lambda $, ${{n}_{c}} = {\text{1}}{\text{.0}}$, ${{n}_{s}} = {\text{1}}{\text{.47}}$, ${{n}_{f}} = {\text{1}}{\text{.565}}$, $A = 1$, m = 0, где длина волны $\lambda = 0.55$ μm.

Первые две моды, которым соответствуют амплитуды ${{T}_{0}}$ и ${{T}_{1}}$, аппроксимируют направляемые моды исходной задачи – моды дискретного спектра исходной задачи. При увеличении толщины ${{h}_{2}} > {{h}_{1}}$ происходит перераспределение энергии между модами, о чем свидельствует рисунок 3: амплитуда ${{T}_{0}}$ уменьшается по абсолютному значению, а амплитуда ${{T}_{1}}$ – увеличивается. Перераспределение энергии происходит не только между модами дискретного спектра исходной задачи, о чем свидетельствуют графики на рис. 4 и 5 : при увеличении толщины ${{h}_{2}}$ суммарная интенсивность мод, аппроксимирующих дискретный спектр исходной задачи, уменьшается от 1.0 до примерно 0.984 при увеличении ${{h}_{2}}$ от $1.5\lambda $ до $1.9\lambda $ (см. рис. 4), а суммарная интенсивность мод, аппроксимирующих непрерывный спектр исходной задачи, увеличивается от 0.0 до примерно 0.016 при увеличении ${{h}_{2}}$ от $1.5\lambda $ до $1.9\lambda $ (см. рис. 5). При дальнейшем увеличении толщины ${{h}_{2}}$ до значения $2.3\lambda $ суммарная интенсивность мод, аппроксимирующих дискретный спектр исходной задачи, увеличивается до примерно 0.995 (см. рис. 4), а суммарная интенсивность мод, аппроксимирующих непрерывный спектр исходной задачи, уменьшается до примерно 0.005 (см. рис. 5). Стоит также отметить, что суммарная интенсивность всех отраженных мод не более 0.001, поэтому графики этих значений не приводятся.

Резюмируя вышесказанное, перераспределение энергии между модами, аппроксимирующими дискретный спектр исходной задачи, существенно – при увеличении толщины ${{h}_{2}}$ абсолютное значение амплитуды ${{T}_{0}}$ уменьшается практически в 2 раза, а абсолютное значение амплитуды ${{T}_{1}}$ увеличивается от 0.0 до 0.7. Доля энергии, которая приходится на моды, аппроксимирующие непрерывный спектр исходной задачи, составляет не более 1.6% для рассматриваемой структуры (см. рис. 5).

6. ЗАКЛЮЧЕНИЕ

В работе построена символьно-численная реализация метода Галеркина, которая позволяет исследовать задачу волноводной дифракции на стыке двух планарных волноводов в “ящике Дирихле”. Использование символьных манипуляций позволяет ускорить расчет интегралов в методе Галеркина благодаря однократному символьному расчету типовых для задачи интегралов вместо многократного численного интегрирования, что позволяет проводить ресурсоемкие расчеты на персональном компьютере для задачи дифракции с количеством собственных функций, исчисляемым сотнями.

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

Список литературы

Tolstikhin O.I., Ostrovsky V.N., Nakamura H. Siegert Pseudo-States as a Universal Tool: Resonances, S Matrix, Green Function // Physical Review Letters. 1997. V. 79. № 11. P. 2026–2029.
Sveshnikov A.G. The basis for a method of calculating irregular waveguides // USSR Computational Mathematics and Mathematical Physics. 1963. V. 3. № 1. P. 219–232.
Eremin Y.A., Sveshnikov A.G. Study of scalar diffraction at a locally inhomogeneous body by a projection method // USSR Computational Mathematics and Mathematical Physics. 1976. V. 16. № 3. P. 255–260.
Delitsyn A.L. On the completeness of the system of eigenvectors of electromagnetic waveguies // Computational Mathematics and Mathematical Physics. 2011. V. 51. № 10. P. 1771–1776.
Sveshnikov A.G. A substantiation of a method for computing the propagation of electromagnetic oscillations in irregular waveguides // USSR Computational Mathematics and Mathematical Physics. 1963. V. 3. № 2. P. 413–429.
Mathematics-based software and services for education, engineering, and research https://www.maplesoft.com/
Свешников А.Г. Неполный метод Галеркина // ДАН СССР. 1977. Т. 236. № 5. С. 1076–1079.
Диваков Д.В., Тютюнник А.А. Символьное исследование спектральных характеристик направляемых мод плавно-нерегулярных волноводов // Программирование. 2022. № 2. С. 23–32.
Tiutiunnik A.A., Divakov D.V., Malykh M.D., Sevastianov L.A. Symbolic-Numeric Implementation of the Four Potential Method for Calculating Normal Modes: An Example of Square Electromagnetic Waveguide with Rectangular Insert // Lecture Notes in Computer Science. 2019. V. 11661. P. 412–429.
Виницкий C.И., Гердт В.П., Гусев А.А., Касчиев М.С., Ростовцев B.А., Самойлов В.Н., Тюпикова Т.В., Чулуунбаатар О. Символьночисленный алгоритм вычисления матричных элементов параметрической задачи на собственные значения // Программирование. 2007. Т. 33. № 2. С. 63–76.
Зорин A.В., Севастьянов Л.А., Третьяков Н.П. Компьютерное моделирование водородоподобных атомов в квантовой механике с неотрицательной функцией распределения // Программирование. 2007. Т. 33. № 2. С. 50–62.
Диваков Д.В., Тютюнник А.А. Символьное исследование собственных векторов для построения общего решения системы ОДУ с символьной матрицей коэффициентов // Программирование. 2021. 1. С. 11–24.
Shevchenko V.V. Spectral decomposition in eigen- and associated functions of a nonselfadjoint problem of Sturm–Liouville type on the entire axis // Differ. Uravn. 1979. V. 15. № 11. P. 2004–2020.
Gevorkyan M.N., Kulyabov D.S., Lovetskiy K.P., Sevastyanov A.L., Sevastyanov L.A. Waveguide modes of a planar optical waveguide // Mathematical Modelling and Geometry. 2015. V. 3. № 1. P. 43–63.
Sevastianov L.A., Egorov A.A., Sevastyanov A.L. Method of adiabatic modes in studying problems of smoothly irregular open waveguide structures // Physics of Atomic Nuclei. 2013. V. 76. № 2. P. 224–239.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Программирование

Архивы выпусков Информация о журнале Отправить рукопись в журнал