Программирование, 2021, № 1, стр. 11-24

СИМВОЛЬНОЕ ИССЛЕДОВАНИЕ СОБСТВЕННЫХ ВЕКТОРОВ ДЛЯ ПОСТРОЕНИЯ ОБЩЕГО РЕШЕНИЯ СИСТЕМЫ ОДУ С СИМВОЛЬНОЙ МАТРИЦЕЙ КОЭФФИЦИЕНТОВ

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

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

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

Поступила в редакцию 14.07.2020
После доработки 05.08.2020
Принята к публикации 13.09.2020

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

Аннотация

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

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

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

1. ВВЕДЕНИЕ

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

1.1. Символьный подход

Символьный подход состоит в ряде аналитических упрощений уравнений, выделяющих важные частные случаи, допускающие при этом исследование на символьном уровне [1, 2]. При таком подходе стараются получать символьные либо символьно-численные выражения, позволяющие качественно исследовать упрощенную модель.

1.2. Численный подход

В рамках численного подхода за постановкой задачи следует построение приближенного метода ее решения – либо разностного [35], либо какого-то иного, вроде методов Галеркина [68], Канторовича [9, 10] или метода конечных элементов [11, 12]. Численный подход требует от исследователя внимательного построения математического метода решения задачи и почти всегда исключает аналитическое упрощение представленной модели и угадывания вида решения. Однако в численном подходе можно строить методы, подходящие для исследования широкого класса задач с большим произволом параметров модели.

1.3. Цель исследования

Авторский подход относится к группе символьных методов решения уравнений Максвелла, отвечающих волноводному распространению в плавно-нерегулярных протяженных волноводных структурах [1, 2].

Основа подхода – представление электромагнитного поля в виде асимптотического ряда [13], в котором явно разделены амплитудная и фазовая части поля. Такое представление электромагнитного поля позволяет в нулевом приближении редуцировать уравнения Максвелла к двум алгебраическим уравнениям и системе четырех обыкновенных дифференциальных уравнений относительно амплитудных частей электромагнитного поля, причем матрица коэффициентов этой системы задана в символьном виде. Общее решение такой системы обыкновенных дифференциальных уравнений можно записать с помощью собственных векторов и собственных значений матрицы коэффициентов. Так как матрица коэффициентов задана в символьном виде, собственные значения и собственные векторы нужно по этой причине отыскивать в символьном виде.

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

2. МОТИВАЦИЯ

В рамках предлагаемого авторами подхода получается система обыкновенных дифференциальных уравнений следующего вида [2]

(2.1)
$\frac{{\partial{ \vec {u}}}}{{\partial x}} + A\left( {x,y,z} \right)\vec {u} = \vec {0},$
в которую входит производная искомой вектор-функции $\vec {u}$ только по x, однако коэффициенты системы есть кусочно-постоянные функции переменной x, а также функции от y, z, поэтому искомая величина $\vec {u}$ зависит от x, y, z.

Если рассмотреть систему в какой-нибудь области постоянства коэффициентов по x, то она будет иметь вид

(2.2)
$\frac{{\partial {{{\vec {u}}}_{j}}}}{{\partial x}} + {{A}_{j}}\left( {y,z} \right){{\vec {u}}_{j}} = \vec {0},$
то есть это система обыкновенных дифференциальных уравнений с постоянными по x коэффициентами, так как коэффициенты системы зависят только от y, z. Также важно отметить, что среди коэффициентов системы встречаются неизвестные функции от y, z, которые определяются из условия разрешимости системы граничных уравнений [2]. Система граничных уравнений формируется из записи граничных условий [2], следующих из условий сопряжения электромагнитных полей на всех поверхностях разрыва коэффициентов системы (2.1).

Для каждой системы вида (2.2), отвечающей каждой области постоянства коэффициентов по x, можно записать ее общее решение с помощью собственных векторов и собственных значений матрицы коэффициентов, в которые будут входить в виде символьных выражений неизвестные функции от y, z. Поэтому для записи общего решения необходимо в символьном виде решать задачу на собственные значения и собственные векторы квадратной матрицы, элементы которой заданы в символьном виде. Подобные возможности заложены в различных системах компьютерной алгебры.

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

В символьное представление собственных векторов могут входить неизвестные функции от y, z, поэтому необходимо, чтобы:

1) символьное представление собственного вектора имело смысл для всей области значений неизвестных функций: например, если функция находится в знаменателе выражения и может обращаться в ноль, то это выражение может устремиться к бесконечности и использовать его в общем случае некорректно;

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

Учитывая описанные особенности рассматриваемой системы дифференциальных уравнений вида (2.1), для построения эффективного численного алгоритма ее решения необходимо исследовать символьное представление собственных векторов ее матрицы коэффициентов с учетом требований 1) и 2) выше. Исследованию этой задачи и посвящена настоящая работа.

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

3. МЕТОДЫ И АЛГОРИТМЫ

3.1. Асимптотическое разложение решения

3.1.1. Вид решения. Электромагнитное поле представляем в виде асимптотического ряда [13], в котором явно разделены амплитудная и фазовая части поля [2]:

(3.1)
$\begin{gathered} \vec {E}(x,y,z,t) = \\ \, = \sum\limits_{s = 0}^\infty \,\frac{{{{{\vec {E}}}^{s}}(x;y,z)}}{{\mathop {\left( { - i\omega } \right)}\nolimits^{\gamma + s} }}exp\left\{ {i\omega t - i{{k}_{0}}\varphi (y,z)} \right\}, \\ \end{gathered} $
(3.2)
$\begin{gathered} \vec {H}(x,y,z,t) = \\ \, = \sum\limits_{s = 0}^\infty \,\frac{{{{{\vec {H}}}^{s}}(x;y,z)}}{{\mathop {\left( { - i\omega } \right)}\nolimits^{\gamma + s} }}exp\left\{ {i\omega t - i{{k}_{0}}\varphi \left( {y,z} \right)} \right\}, \\ \end{gathered} $
где k0 – волновое число, $\varphi \left( {y,z} \right)$ определяет фазу, а ${{\vec {E}}^{s}}(x;y,z),$ ${{\vec {H}}^{s}}\left( {x;y,z} \right)$ определяют амплитуду s-го порядка. В записи ${{\vec {E}}^{s}}\left( {x;y,z} \right),$ ${{\vec {H}}^{s}}\left( {x;y,z} \right)$ отделение x точкой с запятой означает следующее предположение: $\partial {{\vec {E}}^{s}}{\text{/}}\partial y,$ $\partial {{\vec {E}}^{s}}{\text{/}}\partial z,$ $\partial {{\vec {H}}^{s}}{\text{/}}\partial y,$ $\partial {{\vec {H}}^{s}}{\text{/}}\partial z$ являются малыми величинами. Другими словами справедливы следующие выражения для производных:
$\partial{ \vec {E}}{\text{/}}\partial y = - i{{k}_{0}}{{\varphi }_{y}}\vec {E},\quad \partial{ \vec {E}}{\text{/}}\partial z = - i{{k}_{0}}{{\varphi }_{z}}\vec {E},$
и аналогичные выражения:
$\partial cH{\text{/}}\partial y = - i{{k}_{0}}{{\varphi }_{y}}\vec {H},\quad \partial{ \vec {H}}{\text{/}}\partial z = - i{{k}_{0}}{{\varphi }_{z}}\vec {H},$
в которых φy и φz – частные производные $\varphi (y,z)$ по y и z соответственно.

3.1.2. Редукция уравнений Максвелла – общий случай. Уравнения Максвелла в нулевом приближении (s = 0) асимптотического разложения редуцируются к системе обыкновенных дифференциальных уравнений первого порядка [2]:

(3.3)
$\frac{{\partial{ \vec {u}}}}{{\partial x}} + A\left( {x,y,z} \right)\vec {u} = \vec {0},$
и двум дополнительным соотношениям
(3.4)
$\begin{gathered} E_{x}^{0} = {{\varepsilon }^{{ - 1}}}({{\varphi }_{z}}H_{y}^{0} - {{\varphi }_{y}}H_{z}^{0}), \\ H_{x}^{0} = - {{\mu }^{{ - 1}}}({{\varphi }_{z}}E_{y}^{0} - {{\varphi }_{y}}E_{z}^{0}), \\ \end{gathered} $
где искомая вектор-функция $\vec {u}$ содержит величины $\vec {u}(x;y,z) = \mathop {(\begin{array}{*{20}{c}} {E_{y}^{0}}&{H_{z}^{0}}&{H_{y}^{0}}&{E_{z}^{0}} \end{array})}\nolimits^T $, которые описывают распределение по x соответствующей компоненты поля в каждой точке (y, z). Матрица A определена следующим образом [2]
(3.5)
$\begin{gathered} A\left( {x,y,z} \right) = \\ \, = i{{k}_{0}}\left( {\begin{array}{*{20}{c}} 0&{ - \frac{{\varphi _{y}^{2}}}{\varepsilon } + \mu }&{\frac{{{{\varphi }_{y}}{{\varphi }_{z}}}}{\varepsilon }}&0 \\ { - \frac{{\varphi _{z}^{2}}}{\mu } + \varepsilon }&0&0&{\frac{{{{\varphi }_{y}}{{\varphi }_{z}}}}{\mu }} \\ { - \frac{{{{\varphi }_{y}}{{\varphi }_{z}}}}{\mu }}&0&0&{\frac{{\varphi _{y}^{2}}}{\mu } - \varepsilon } \\ 0&{ - \frac{{{{\varphi }_{y}}{{\varphi }_{z}}}}{\varepsilon }}&{\frac{{\varphi _{z}^{2}}}{\varepsilon } - \mu }&0 \end{array}} \right), \\ \end{gathered} $
где $\varepsilon = \varepsilon \left( {x,y,z} \right)$ и $\mu = \mu \left( {x,y,z} \right)$ есть кусочно-постоянные функции – диэлектрическая и магнитная проницаемости соответственно. В каждом слое волновода с постоянными диэлектрической и магнитной проницаемостями $\varepsilon ,\mu $ решение системы дифференциальных уравнений (3.3) можно представить в следующем виде [2], который получен с помощью функции Eigenvectors пакета LinearAlgebra в Maple [14]:
(3.6)
$\begin{gathered} \vec {u}\left( {x;y,z} \right) = A\left( {\begin{array}{*{20}{c}} q \\ { - i\varepsilon \eta } \\ 0 \\ {{{\varphi }_{y}}{{\varphi }_{z}}} \end{array}} \right){{e}^{{\gamma x}}} + B\left( {\begin{array}{*{20}{c}} { - i\mu \eta } \\ p \\ {{{\varphi }_{y}}{{\varphi }_{z}}} \\ 0 \end{array}} \right){{e}^{{\gamma x}}} + \\ \, + C\left( {\begin{array}{*{20}{c}} q \\ {i\varepsilon \eta } \\ 0 \\ {{{\varphi }_{y}}{{\varphi }_{z}}} \end{array}} \right){{e}^{{ - \gamma x}}} + D\left( {\begin{array}{*{20}{c}} {i\mu \eta } \\ p \\ {{{\varphi }_{y}}{{\varphi }_{z}}} \\ 0 \end{array}} \right){{e}^{{ - \gamma x}}}, \\ \end{gathered} $
где $q = \varphi _{y}^{2} - \varepsilon \mu $, $p = \varphi _{z}^{2} - \varepsilon \mu $, $\eta = \sqrt {\varphi _{y}^{2} + \varphi _{z}^{2} - \varepsilon \mu } $, γ = = k0η, а $A,B,C,D$ – неопределенные константы в каждой точке (y, z).

Замечание. Собственному значению γ = = ${{k}_{0}}\sqrt {\varphi _{y}^{2}\, + \,\varphi _{z}^{2}\, - \,\varepsilon \mu } $ отвечают два собственных вектора

(3.7)
${v}_{1}^{ + } = \left( {\begin{array}{*{20}{c}} q \\ { - i\varepsilon \eta } \\ 0 \\ {{{\varphi }_{y}}{{\varphi }_{z}}} \end{array}} \right),\quad {v}_{2}^{ + } = \left( {\begin{array}{*{20}{c}} { - i\mu \eta } \\ p \\ {{{\varphi }_{y}}{{\varphi }_{z}}} \\ 0 \end{array}} \right),$
а собственному значению –γ отвечают два других собственных вектора

(3.8)
$v_{1}^{ - } = \left( {\begin{array}{*{20}{c}} q \\ {i\varepsilon \eta } \\ 0 \\ {{{\varphi }_{y}}{{\varphi }_{z}}} \end{array}} \right),\quad v_{2}^{ - } = \left( {\begin{array}{*{20}{c}} {i\mu \eta } \\ p \\ {{{\varphi }_{y}}{{\varphi }_{z}}} \\ 0 \end{array}} \right).$

При ${{\varphi }_{y}}{{\varphi }_{z}} \ne 0$ они линейно независимы.

3.1.3. Редукция уравнений Максвелла – частный случай. Уравнения Максвелла в нулевом приближении (s = 0) асимптотического разложения при ${{\varphi }_{y}} \equiv 0$ редуцируются к системе обыкновенных дифференциальных уравнений первого порядка [2]:

(3.9)
$\frac{{\partial{ \vec {u}}}}{{\partial x}} + {{A}^{0}}(x,y,z)\vec {u} = \vec {0},$
и двум дополнительным соотношениям
(3.10)
$E_{x}^{0} = {{\varepsilon }^{{ - 1}}}{{\varphi }_{z}}H_{y}^{0},\quad H_{x}^{0} = - {{\mu }^{{ - 1}}}{{\varphi }_{z}}E_{y}^{0},$
где искомая вектор-функция $\vec {u}$ содержит величины $\vec {u}\left( {x;y,z} \right) = \mathop {(\begin{array}{*{20}{c}} {E_{y}^{0}}&{H_{z}^{0}}&{H_{y}^{0}}&{E_{z}^{0}} \end{array})}\nolimits^T $, которые описывают распределение по x соответствующей компоненты поля в каждой точке (y, z). Матрица A0 определена следующим образом

(3.11)
${{A}^{0}}\left( {x,y,z} \right) = i{{k}_{0}}\left( {\begin{array}{*{20}{c}} 0&\mu &0&0 \\ { - \frac{{\varphi _{z}^{2}}}{\mu } + \varepsilon }&0&0&0 \\ 0&0&0&{ - \varepsilon } \\ 0&0&{\frac{{\varphi _{z}^{2}}}{\varepsilon } - \mu }&0 \end{array}} \right).$

В слое с постоянными диэлектрической и магнитной проницаемостями $\varepsilon ,\mu $ решение системы дифференциальных уравнений (3.10) можно представить в следующем виде [2], который получен с помощью функции Eigenvectors пакета LinearAlgebra в Maple [14]:

(3.12)
$\begin{gathered} \vec {u}\left( {x;y,z} \right) = A\left( {\begin{array}{*{20}{c}} 0 \\ 0 \\ {i\varepsilon {{\eta }^{0}}} \\ {{{p}^{0}}} \end{array}} \right){{e}^{{{{\gamma }^{0}}x}}} + B\left( {\begin{array}{*{20}{c}} { - i\mu {{\eta }^{0}}} \\ {{{p}^{0}}} \\ 0 \\ 0 \end{array}} \right){{e}^{{{{\gamma }^{0}}x}}} + \\ \, + C\left( {\begin{array}{*{20}{c}} 0 \\ 0 \\ { - i\varepsilon {{\eta }^{0}}} \\ {{{p}^{0}}} \end{array}} \right){{e}^{{ - {{\gamma }^{0}}x}}} + D\left( {\begin{array}{*{20}{c}} {i\mu {{\eta }^{0}}} \\ {{{p}^{0}}} \\ 0 \\ 0 \end{array}} \right){{e}^{{ - {{\gamma }^{0}}x}}}, \\ \end{gathered} $
где ${{p}^{0}} = \varphi _{z}^{2} - \varepsilon \mu $, ${{\eta }^{0}} = \sqrt {\varphi _{z}^{2} - \varepsilon \mu } $, ${{\gamma }^{0}} = {{k}_{0}}{{\eta }^{0}}$, а A, B, C, D – неопределенные константы в каждой точке (y, z).

Замечание. Собственному значению γ0 = = ${{k}_{0}}\sqrt {\varphi _{z}^{2}\, - \,\varepsilon \mu } $ отвечают два собственных вектора

(3.13)
$v_{1}^{ + } = \left( {\begin{array}{*{20}{c}} 0 \\ 0 \\ {i\varepsilon {{\eta }^{0}}} \\ {{{p}^{0}}} \end{array}} \right),\quad v_{2}^{ + } = \left( {\begin{array}{*{20}{c}} { - i\mu {{\eta }^{0}}} \\ {{{p}^{0}}} \\ 0 \\ 0 \end{array}} \right),$
а собственному значению –γ0 отвечают два других собственных вектора

(3.14)
$v_{1}^{ - } = \left( {\begin{array}{*{20}{c}} 0 \\ 0 \\ { - i\varepsilon {{\eta }^{0}}} \\ {{{p}^{0}}} \end{array}} \right),\quad v_{2}^{ - } = \left( {\begin{array}{*{20}{c}} {i\mu {{\eta }^{0}}} \\ {{{p}^{0}}} \\ 0 \\ 0 \end{array}} \right).$

Причем, по сравнению с общим случаем

$\begin{gathered} {{p}^{0}} = \mathop {lim}\limits_{{{\varphi }_{y}} \to 0} p,\quad {{\eta }^{0}} = \mathop {lim}\limits_{{{\varphi }_{y}} \to 0} \eta , \\ {{\gamma }^{0}} = \mathop {lim}\limits_{{{\varphi }_{y}} \to 0} \gamma . \\ \end{gathered} $

Главная мотивация исследования. Неудобство использования приведенных выше символьных представлений собственных векторов состоит в том, что при ${{\varphi }_{y}} \to 0$ собственные векторы (3.7), (3.8) не переходят в (3.13), (3.14) и становятся линейно зависимыми. С физической точки зрения φy может равняться нулю, причем, так как φy искомая величина, заранее мы не знаем при каких y, z она обратится в ноль, и потому нам нужно универсальное представление собственных векторов. А именно такое, при котором собственные векторы в общем случае переходят в (3.13), (3.14) при ${{\varphi }_{y}} \to 0$ и, аналогично, при ${{\varphi }_{z}} \to 0$. С целью добиться такого соответствия проведем дополнительное исследование возможных символьных представлений собственных векторов, основанное на преобразованиях перестановок.

3.2. Основа алгоритма исследования – преобразование перестановок

Рассмотрим систему обыкновенных дифференциальных уравнений с постоянными коэффициентами

(3.15)
$\vec {y}{\kern 1pt} '\, - A\vec {y} = \vec {0},$
где $\vec {y} = \mathop {(\begin{array}{*{20}{c}} {{{y}_{1}}}&{{{y}_{2}}}& \ldots &{{{y}_{N}}} \end{array})}\nolimits^T $ – вектор искомых величин, A – квадратная невырожденная матрица коэффициентов. Рассмотрим эту же систему, но при другом порядке следования искомых величин $\vec {z} = P\vec {y}$, где матрица перестановок P определяет новый порядок следования искомых величин. Система (3.15) при новом порядке следования искомых величин $\vec {z}$ преобразуется следующим образом:

(3.16)
$\vec {z}{\kern 1pt} '\, - PA{{P}^{T}}\vec {z} = \vec {0}.$

При выводе (3.16) использовано свойство матрицы перестановок ${{P}^{{ - 1}}} = {{P}^{T}}$.

Будем рассматривать случай диагонализируемой матрицы A. В таком случае общее решение однородной системы (3.15) можно записать в следующем виде

(3.17)
$\vec {y}\left( x \right) = {{A}_{1}}{{\vec {v}}_{1}}{{e}^{{{{\lambda }_{1}}x}}} + {{A}_{2}}{{\vec {v}}_{2}}{{e}^{{{{\lambda }_{2}}x}}} + \ldots + {{A}_{N}}{{\vec {v}}_{N}}{{e}^{{{{\lambda }_{N}}x}}}$
где $\mathop {\left\{ {{{{\vec {v}}}_{j}}} \right\}}\nolimits_{j = 1}^N $ – собственные векторы матрицы A, отвечающие ее собственным значениям $\mathop {\left\{ {{{\lambda }_{j}}} \right\}}\nolimits_{j = 1}^N $. Нетрудно проверить, что собственные значения матрицы $PA{{P}^{T}}$ совпадают с собственными значениями матрицы A. Собственные векторы $\mathop {\left\{ {\mathop {\vec {w}}\nolimits_j } \right\}}\nolimits_{j = 1}^N $ матрицы $PA{{P}^{T}}$ связаны с собственными векторами матрицы A следующим образом
(3.18)
${{\vec {w}}_{j}} = P{{\vec {v}}_{j}},$
где нумерация векторов выбрана таким образом, чтобы при фиксированном  j векторы ${{\vec {w}}_{j}}$ и ${{\vec {v}}_{j}}$ отвечали одному и тому же собственному значению λj. Отсюда же следует связь между матрицами собственных векторов
(3.19)
$W = PV,$
где столбцами матриц W и V выступают векторы ${{\vec {w}}_{j}}$ и ${{\vec {v}}_{j}}$ соответственно.

Замечание. Существует множество различных невырожденных преобразований, фактически любая невырожденная матрица P определяет некоторое преобразование искомых величин исследуемой системы дифференциальных уравнений $\vec {z} = P\vec {y}$. В настоящей работе мы рассматриваем только преобразования перестановок.

Как показывает практика в случае, когда коэффициенты системы (3.15) заданы в символьном виде, в зависимости от преобразования порядка следования искомых величин $\vec {z} = P\vec {y}$ символьное представление собственных векторов различно. Другими словами, в зависимости от порядка следования искомых величин системы (3.15) могут получаться различные символьные представления собственных векторов.

Задача настоящего исследования состоит в том, чтобы исследовать все возможные различные символьные представления собственных векторов, получаемые при различных перестановках искомых величин, заданных различными матрицами перестановок P.

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

3.3. Алгоритм исследования

Ниже приведем формулировку метода, с помощью которого в рамках работы мы изучаем собственные векторы, определяющие решение системы (3.15).

1. Для изучения различных вариантов символьного представления собственных векторов сначала сгенерируем все возможные матрицы перестановок ${{P}_{j}} = P({{s}_{j}})$, где sj определяет текущую перестановку из n элементов, причем ${{s}_{1}} = (\begin{array}{*{20}{c}} 1&2& \ldots &n \end{array})$ – тождественная перестановка.

Генерацию всех возможных различных перестановок ${{s}_{j}},j = \overline {1,n!} $ в Maple осуществляем с помощью функций firstperm и nextperm пакета combinat [14].

Для каждой перестановки sj определяем матрицу перестановок Pj, элементы которой состоят из нулей и единиц: ${{P}_{j}}\left[ k \right][{{s}_{{jk}}}] = 1$, а остальные равны нулю.

2. Далее формируем матрицу ${{B}_{j}} = {{P}_{j}}AP_{j}^{T}$ и для нее символьно определяем матрицу собственных векторов Wj.

Символьное вычисление собственных векторов осуществляем в Maple с помощью функции Eigenvectors пакета LinearAlgebra [14].

Учитывая формулу (3.19) каждая матрица собственных векторов Wj определяет матрицу собственных векторов ${{V}_{j}} = {{P}_{j}}{{W}_{j}}$ исходной матрицы A. Собственные значения исходной матрицы A и любой матрицы ${{B}_{j}} = {{P}_{j}}AP_{j}^{T}$ совпадают, однако символьное представление собственных векторов могут отличаться.

3. После вычисления всех матриц собственных векторов ${{V}_{j}} = {{P}_{j}}{{W}_{j}}$, полученных всеми возможными различными перестановками sj, необходимо избавиться от всех линейно зависимых собственных векторов, отвечающих каждому собственному значению.

Если векторы $\vec {u}$ и $\vec {v}$ линейно зависимы, то $\left\| {\alpha \vec {u} - \vec {v}} \right\| = 0$ при $\alpha = {{\vec {u}}^{H}}\vec {v}{\text{/}}{{\vec {u}}^{H}}\vec {u}$ и из векторов $\vec {u}$ и $\vec {v}$ достаточно оставить только один любой вектор.

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

Применим описанный алгоритм для исследования собственных векторов, определяющих решения системы (3.9) в частном случае ${{\varphi }_{y}} = 0$ и системы (3.3) в общем случае.

3.4. Численный эксперимент

В качестве демонстрации применения полученных символьных результатов проведем численный эксперимент, в качестве структуры рассматриваем структуру из [2]. В рамках предположения линейной фазы из [2] ${{\varphi }_{z}} = \beta $, где β – коэффициент фазового замедления. Вычислим коэффициенты фазового замедления и соответствующие им направляемые волноводные моды для волноводной структуры из [2].

Структура представляет собой трехслойный волновод с двумя полубесконечными слоями – подложкой с ${{\varepsilon }_{s}} = n_{s}^{2}$, занимающей область x < 0 и покровным слоем с ${{\varepsilon }_{c}} = n_{c}^{2}$, занимающей область x > h. Между двумя полубесконечными слоями $0 < x < h$ расположен волноводный слой с ${{\varepsilon }_{f}} = n_{f}^{2}$. На границах между слоями x = 0 и x = h должны выполняться условия непрерывности компонент Ey, Hz, ${{H}_{y}},{{E}_{z}}$. Поля направляемых мод характеризуются также убыванием при бесконечном удалении от волноводного слоя, то есть для направляемых мод должны выполняться асимптотические условия вида $\left\| {\vec {E}} \right\|\;\xrightarrow[{x \to \pm \infty }]{}\;0$, $\left\| {\vec {H}} \right\|\;\xrightarrow[{x \to \pm \infty }]{}\;0$.

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

(3.20)
$\begin{gathered} \lambda = 0.55;\quad {{k}_{0}} = 2\pi {\text{/}}\lambda ;\quad h = 2\lambda \\ {{n}_{c}} = 1.0;\quad {{n}_{f}} = 1.565;\quad {{n}_{s}} = 1.47; \\ {{\mu }_{c}} = {{\mu }_{f}} = {{\mu }_{s}} = 1. \\ \end{gathered} $

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

1. первое представление соответствует случаю ${{\varphi }_{y}} = 0$;

2. второе представление соответствует случаю ${{\varphi }_{y}} = \delta \ne 0$, собственные векторы для которого выберем с помощью описанного в разделе 3 символьного алгоритма;

3. третье представление соответствует случаю ${{\varphi }_{y}} = \delta \ne 0$, собственные векторы для которого определены с помощью встроенной функции Eigenvectors (без использования описанного алгоритма).

Мы ожидаем, что при стремлении δ к нулю коэффициенты фазового замедления и поля соответствующих волноводных мод, полученные на основе собственных векторов с помощью предложенного алгоритма (второе представление), сходятся к аналогичным результатам для случая ${{\varphi }_{y}} = 0$ (первое представление). Исследуем также сходимость коэффициентов фазового замедления и полей соответствующих волноводных мод, построенных на основе собственных векторов, определенных с помощью встроенной функции Eigenvectors (третье представление) к аналогичным результатам для случая ${{\varphi }_{y}} = 0$ (первое представление).

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

(3.21)
$\vec {u}_{m}^{{}}\left( {x;y,z} \right) = \left\{ \begin{gathered} \vec {u}_{m}^{c}\left( {x;y,z} \right),\quad x > h \hfill \\ \vec {u}_{m}^{f}\left( {x;y,z} \right),\quad 0 < x < h \hfill \\ \vec {u}_{m}^{s}\left( {x;y,z} \right),\quad x < 0 \hfill \\ \end{gathered} \right.$
где m = 1, 2, 3. Для каждого из трех представлений система граничных уравнений формулируется из условий сопряжения при x = 0 и x = h:
(3.22)
$\vec {u}_{m}^{c}\left( {h;y,z} \right) - \vec {u}_{m}^{f}\left( {h;y,z} \right) = \vec {0},$
(3.23)
$\vec {u}_{m}^{f}\left( {0;y,z} \right) - \vec {u}_{m}^{s}\left( {0;y,z} \right) = \vec {0},$
где m = 1, 2, 3. В результате получаем для каждого m = 1, 2, 3 однородную систему линейных алгебраических уравнений относительно искомых величин $A_{j}^{\alpha }$ для m = 1, $B_{j}^{\alpha }$ для m = 2 и $D_{j}^{\alpha }$ для m = 3. Коэффициенты каждой системы зависят от $\beta = {{\varphi }_{z}}$ при фиксированном ${{\varphi }_{y}} = \delta $. Условие разрешимости однородной системы линейных алгебраических уравнений – это равенство нулю определителя матрицы коэффициентов системы ${\text{det}}{{M}_{m}}(\beta )$ = 0, m = 1, 2, 3.

Отыскиваем численно такие значения ${{\beta }_{{s,m}}}$ (индекс s определяет номер направляемой моды), при которых $det{{M}_{m}}\left( \beta \right) = 0$ при m = 1, 2, 3, а также для каждого из таких значений отыщем набор искомых величин, которые определяют поле соответствующей волноводной моды: $A_{j}^{\alpha }$ для m = 1, $B_{j}^{\alpha }$ для m = 2 и $D_{j}^{\alpha }$ для m = 3.

Для поиска нулей воспользуемся классическим методом деления отрезка пополам [15], для поиска решений однородных систем уравнений воспользуемся командой Eigenvectors и выберем тот собственный вектор, который соответствует минимальному по модулю собственному значению. Для почти вырожденной матрицы минимальное по модулю собственное значение близко к нулю, поэтому собственный вектор, отвечающий минимальному по модулю собственному значению, будет приближенным решением соответствующей однородной системы.

Определив искомые величины $A_{j}^{\alpha }$ для m = 1, $B_{j}^{\alpha }$ для m = 2 и $D_{j}^{\alpha }$ для m = 3, можно построить поля волноводных мод.

Замечание. Поля волноводных мод определяются с помощью посчитанных численно компонент собственного вектора, отвечающего почти нулевому собственному значению матрицы. Собственный вектор известен с точностью до мультипликативной константы, поэтому и соответствующие поля также известны с точностью до мультипликативной константы. Поэтому, для корректного сравнения полей необходимо их нормировать в каком-либо функциональном пространстве. Учитывая только непрерывность компонент поля на границах раздела и их убывание на бесконечности, каждую компоненту поля можно рассматривать как элемент пространства ${{L}_{2}}\left( \mathbb{R} \right)$ на всей вещественной оси. Поэтому рассматриваемые нами вектор-функции ${{\vec {u}}_{m}}\left( {x;y,z} \right)$, m = 1, 2, 3, состоящие из четырех компонент поля, могут рассматриваться как элементы пространства $H = \mathop {\left\{ {{{L}_{2}}\left( \mathbb{R} \right)} \right\}}\nolimits^4 $, состоящего из таких вектор-функций $\vec {u}$, для которых справедливо ${{u}_{j}} \in {{L}_{2}}\left( \mathbb{R} \right),$ $j = \overline {1..4} $. Скалярное произведение в пространстве H определено как

(3.24)
$\left( {\vec {u},\vec {v}} \right) = \int\limits_{ - \infty }^{ + \infty } \,\left\{ {\sum\limits_{j = 1}^4 \,{{u}_{j}}\left( x \right) \cdot \overline {{{{v}}_{j}}\left( x \right)} } \right\}dx,$
где черта сверху обозначает комплексное сопряжение. Норма элемента $\vec {u}$ определена следующим образом $\left\| {\vec {u}} \right\| = \sqrt {\left( {\vec {u},\vec {u}} \right)} $.

Вычисленные для каждого ${{\beta }_{{s,m}}}$ вектор-функции ${{\vec {u}}_{{s,m}}}\left( {x;y,z} \right)$

• во-первых нормируем, чтобы норма каждой функции была единичной, то есть $\vec {u}_{{s,m}}^{{norm}}$ = ${{\vec {u}}_{{s,m}}}{\text{/}}\left\| {{{{\vec {u}}}_{{s,m}}}} \right\|$ для m = 1, 2, 3;

• во-вторых подбираем такие комплексные числа ${{\tau }_{{s,2}}}$ и ${{\tau }_{{s,3}}}$, что $\left| {{{\tau }_{{s,2}}}} \right|$ = 1, $\left| {{{\tau }_{{s,3}}}} \right|$ = 1, чтобы $\left\| {\vec {u}_{{s,1}}^{{norm}} - {{\tau }_{{s,2}}} \cdot \vec {u}_{{s,2}}^{{norm}}} \right\|$ и $\left\| {\vec {u}_{{s,1}}^{{norm}} - {{\tau }_{{s,3}}} \cdot \vec {u}_{{s,3}}^{{norm}}} \right\|$ были минимальны (для этого воспользуемся алгоритмом, представленным ниже в разделе 3.5).

В рамках численных экспериментов оценим для различных значений δ и для различных направляемых мод (индекс s определяет номер направляемой моды) следующие величины

(3.25)
${{\Delta }_{2}}{{\beta }_{s}}\left( \delta \right) = \left| {\frac{{{{\beta }_{{s,1}}} - {{\beta }_{{s,2}}}}}{{{{\beta }_{{s,1}}}}}} \right|,$
(3.26)
${{\Delta }_{3}}{{\beta }_{s}}\left( \delta \right) = \left| {\frac{{{{\beta }_{{s,1}}} - {{\beta }_{{s,3}}}}}{{{{\beta }_{{s,1}}}}}} \right|,$
характеризующие близость вычисленных коэффициентов фазового замедления при ${{\varphi }_{y}} \ne 0$ с использованием приведенного алгоритма (m = 2) и без его использования (m = 3) к соответствующим коэффициентам фазового замедления при ${{\varphi }_{y}} = 0$ (m = 1).

Также оценим для различных значений $\delta $ и для различных направляемых мод (индекс s определяет номер направляемой моды) следующие величины

(3.27)
${{\Delta }_{2}}{{u}_{s}}\left( \delta \right) = \left\| {\vec {u}_{{s,1}}^{{norm}} - {{\tau }_{{s,2}}} \cdot \vec {u}_{{s,2}}^{{norm}}} \right\|,$
(3.28)
${{\Delta }_{3}}{{u}_{s}}\left( \delta \right) = \left\| {\vec {u}_{{s,1}}^{{norm}} - {{\tau }_{{s,3}}} \cdot \vec {u}_{{s,3}}^{{norm}}} \right\|,$
характеризующие близость вычисленных (нормированных) полей направляемых волноводных мод при ${{\varphi }_{y}} \ne 0$ с использованием приведенного алгоритма (m = 2) и без его использования (m = 3) к соответствующим (нормированным) полям направляемых мод при ${{\varphi }_{y}} = 0$ (m = 1).

3.5. Алгоритм подбора констант для сравнения полей

Рассмотрим абстрактное гильбертово пространство H [16, 17]. Пусть имеется два элемента $u,v \in H$ такие, что $\left\| u \right\| = \left\| v \right\| = 1$. Ищем такую комплексную константу $\alpha \in \mathbb{C}$, $\left| \alpha \right| = 1$, чтобы элементы ${{v}_{\alpha }} = \alpha v$ и u отличались минимально друг от друга в метрике рассматриваемого пространства H. Другими словами ищем такую комплексную константу $\alpha \in \mathbb{C}$, $\left| \alpha \right| = 1$ (при таком выборе норма элемента $\alpha v$ остается единичной, так как $\left\| {\alpha v} \right\| = \left| \alpha \right|\left\| v \right\| \equiv \left\| v \right\|$) чтобы $\mathop {\left\| {\alpha v - u} \right\|}\nolimits^2 \to min$. Для простоты рассмотрения введем параметр $t \in \left[ {0,2\pi } \right]$: α = ${{e}^{{it}}}$ и отыщем минимум неотрицательной функции

(3.29)
$f\left( t \right) = \mathop {\left\| {{{e}^{{it}}}v - u} \right\|}\nolimits^2 \geqslant 0$
на отрезке $t \in \left[ {0,2\pi } \right]$.

Правая часть f(t) в развернутом виде имеет следующий вид $f\left( t \right) = ({{e}^{{it}}}v - u,{{e}^{{it}}}v - u)$:

(3.30)
$f\left( t \right) = \mathop {\left\| u \right\|}\nolimits^2 + \mathop {\left\| v \right\|}\nolimits^2 - {{e}^{{it}}}\left( {v,u} \right) - {{e}^{{ - it}}}\left( {u,v} \right).$

Для поиска точек экстремума функции f(t) пользуемся необходимым условием экстремума [18]: находим ее производную

(3.31)
$f{\kern 1pt} '\left( t \right) = i({{e}^{{ - it}}}\left( {u,v} \right) - {{e}^{{it}}}\left( {v,u} \right)),$
и решаем уравнение $f{\kern 1pt} '\left( t \right) = 0$. Уравнение $f{\kern 1pt} '(t)$ = 0 тождественно выполнено, если ${{e}^{{2it}}} = \left( {u,v} \right){\text{/}}\left( {v,u} \right)$. Обозначим $c = \left( {u,v} \right)$, поэтому $\bar {c} = \left( {v,u} \right)$ и с учетом этих обозначений:

(3.32)
${{\alpha }_{ \pm }} = {{e}^{{i{{t}_{ \pm }}}}} \equiv \pm \frac{c}{{\left| c \right|}}.$

Замечание. В случае, когда элементы $u,v \in H$ ортогональны друг другу, c = 0. В таком случае задача не имеет смысла, так как  f(t) не имеет ни минимума, ни максимума, потому что $f\left( t \right) \equiv {\text{const}}$ при $t \in \left[ {0,2\pi } \right]$. Поэтому рассматриваем случай $c \ne 0$.

Одно из двух значений ${{\alpha }_{ \pm }} = {{e}^{{i{{t}_{ \pm }}}}}$ обеспечивает минимум функции f(t), а второе – максимум на отрезке $t \in \left[ {0,2\pi } \right]$. Для определения искомого значения константы α воспользуемся достаточным условием минимума [18]: находим вторую производную функции $f\left( t \right)$

(3.33)
$f{\kern 1pt} '{\kern 1pt} '\left( t \right) = {{e}^{{ - it}}}\left( {u,v} \right) + {{e}^{{it}}}\left( {v,u} \right)$
и проверяем в какой из точек ${{t}_{ \pm }}$ значение второй производной положительно.

Непосредственной подстановкой проверяется, что $f{\kern 1pt} '{\kern 1pt} '\left( {{{t}_{ + }}} \right) = 2\left| c \right| > 0$ при $c \ne 0$. Поэтому искомое значение константы α определено следующим образом

(3.34)
$\alpha = {{e}^{{i{{t}_{ + }}}}} \equiv \frac{c}{{\left| c \right|}}.$

Замечание. Для решения описанной задачи требуется вычислить только одно скалярное произведение для расчета константы c.

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

Применим алгоритм, описанный в предыдущем разделе, сначала для исследования собственных векторов системы (3.9), описывающей важный частный случай ${{\varphi }_{y}} = 0$.

4.1. Символьные результаты для частного случая

Исследуя этот частный случай при любых возможных преобразованиях перестановок, мы получаем тривиальный случай – каждому собственному значению $\lambda _{1}^{0} = \lambda _{2}^{0} = {{\gamma }^{0}}$, $\lambda _{3}^{0} = \lambda _{4}^{0} = - {{\gamma }^{0}}$ соответствует единственное символьное представление собственного вектора, а именно:

(4.1)
$\begin{gathered} \vec {v}_{1}^{0} = \left( {\begin{array}{*{20}{c}} { - i\mu } \\ {\sqrt {\varphi _{z}^{2} - \varepsilon \mu } } \\ 0 \\ 0 \end{array}} \right),\quad \vec {v}_{2}^{0} = \left( {\begin{array}{*{20}{c}} 0 \\ 0 \\ {i\varepsilon } \\ {\sqrt {\varphi _{z}^{2} - \varepsilon \mu } } \end{array}} \right), \\ \vec {v}_{3}^{0} = \left( {\begin{array}{*{20}{c}} {i\mu } \\ {\sqrt {\varphi _{z}^{2} - \varepsilon \mu } } \\ 0 \\ 0 \end{array}} \right),\quad \vec {v}_{4}^{0} = \left( {\begin{array}{*{20}{c}} 0 \\ 0 \\ { - i\varepsilon } \\ {\sqrt {\varphi _{z}^{2} - \varepsilon \mu } } \end{array}} \right), \\ \end{gathered} $
где ${{\gamma }^{0}} = {{k}_{0}}\sqrt {\varphi _{z}^{2} - \varepsilon \mu } $.

Замечание. В действительности мы получили по два ортогональных собственных вектора, отвечающих каждому из двух двукратных собственных значений $ \pm {{\gamma }^{0}}$. Однако, если за первое и второе собственные значения принять $\lambda _{1}^{0} = \lambda _{2}^{0} = {{\gamma }^{0}}$ и условиться за первый собственный вектор принять тот из собственных векторов, у которого первые две компоненты отличны от нуля, а за второй – тот, у которого две последних компоненты отличны от нуля, то в таком случае представление собственных векторов, отвечающих первым двум собственным значениям $\lambda _{1}^{0} = \lambda _{2}^{0} = {{\gamma }^{0}}$, единственно. Также единственны представления третьего и четвертого собственных векторов, принимая аналогичную условность для них, отвечающих третьему и четвертому собственным значениям –γ0.

4.2. Символьные результаты для общего случая

В общем случае (${{\varphi }_{y}} \ne 0,$ ${{\varphi }_{z}} \ne 0$) мы получаем нетривиальный случай – каждому собственному значению ${{\lambda }_{1}} = {{\lambda }_{2}} = \gamma $, ${{\lambda }_{3}} = {{\lambda }_{4}} = - \gamma $ соответствует по два различных представления каждого собственного вектора, а именно:

(4.2)
$\begin{gathered} {{{\vec {v}}}_{{1,1}}} = \left( {\begin{array}{*{20}{c}} { - i\mu \sqrt {\varphi _{y}^{2} + \varphi _{z}^{2} - \varepsilon \mu } } \\ {\varphi _{z}^{2} - \varepsilon \mu } \\ {{{\varphi }_{y}}{{\varphi }_{z}}} \\ 0 \end{array}} \right), \\ {{{\vec {v}}}_{{1,2}}} = \left( {\begin{array}{*{20}{c}} {\varphi _{y}^{2} - \varepsilon \mu } \\ { - i\varepsilon \sqrt {\varphi _{y}^{2} + \varphi _{z}^{2} - \varepsilon \mu } } \\ 0 \\ {{{\varphi }_{y}}{{\varphi }_{z}}} \end{array}} \right), \\ \end{gathered} $
(4.3)
$\begin{gathered} {{{\vec {v}}}_{{2,1}}} = \left( {\begin{array}{*{20}{c}} {{{\varphi }_{y}}{{\varphi }_{z}}} \\ 0 \\ {i\varepsilon \sqrt {\varphi _{y}^{2} + \varphi _{z}^{2} - \varepsilon \mu } } \\ {\varphi _{z}^{2} - \varepsilon \mu } \end{array}} \right), \\ {{{\vec {v}}}_{{2,2}}} = \left( {\begin{array}{*{20}{c}} 0 \\ {{{\varphi }_{y}}{{\varphi }_{z}}} \\ {\varphi _{y}^{2} - \varepsilon \mu } \\ {i\mu \sqrt {\varphi _{y}^{2} + \varphi _{z}^{2} - \varepsilon \mu } } \end{array}} \right), \\ \end{gathered} $
(4.4)
$\begin{gathered} {{{\vec {v}}}_{{3,1}}} = \left( {\begin{array}{*{20}{c}} {i\mu \sqrt {\varphi _{y}^{2} + \varphi _{z}^{2} - \varepsilon \mu } } \\ {\varphi _{z}^{2} - \varepsilon \mu } \\ {{{\varphi }_{y}}{{\varphi }_{z}}} \\ 0 \end{array}} \right), \\ {{{\vec {v}}}_{{3,2}}} = \left( {\begin{array}{*{20}{c}} {\varphi _{y}^{2} - \varepsilon \mu } \\ {i\varepsilon \sqrt {\varphi _{y}^{2} + \varphi _{z}^{2} - \varepsilon \mu } } \\ 0 \\ {{{\varphi }_{y}}{{\varphi }_{z}}} \end{array}} \right), \\ \end{gathered} $
(4.5)
$\begin{gathered} {{{\vec {v}}}_{{4,1}}} = \left( {\begin{array}{*{20}{c}} {{{\varphi }_{y}}{{\varphi }_{z}}} \\ 0 \\ { - i\varepsilon \sqrt {\varphi _{y}^{2} + \varphi _{z}^{2} - \varepsilon \mu } } \\ {\varphi _{z}^{2} - \varepsilon \mu } \end{array}} \right), \\ {{{\vec {v}}}_{{4,2}}} = \left( {\begin{array}{*{20}{c}} 0 \\ {{{\varphi }_{y}}{{\varphi }_{z}}} \\ {\varphi _{y}^{2} - \varepsilon \mu } \\ { - i\mu \sqrt {\varphi _{y}^{2} + \varphi _{z}^{2} - \varepsilon \mu } } \end{array}} \right), \\ \end{gathered} $
где $\gamma = {{k}_{0}}\sqrt {\varphi _{y}^{2} + \varphi _{z}^{2} - \varepsilon \mu } $.

Замечание 1. Нумерация собственных значений выбрана по аналогии с частным случаем, причем $\mathop {lim}\limits_{{{\varphi }_{y}} \to 0} {{\lambda }_{j}} = \lambda _{j}^{0},$ $j = \overline {1,4} $. Нумерация собственных векторов выбрана исходя из условия $\mathop {lim}\limits_{{{\varphi }_{y}} \to 0} {{\vec {v}}_{{j,1}}} = {{\alpha }_{{j,1}}}\vec {v}_{j}^{0},$ $\mathop {lim}\limits_{{{\varphi }_{y}} \to 0} {{\vec {v}}_{{j,2}}} = {{\alpha }_{{j,2}}}\vec {v}_{j}^{0},$ $j = \overline {1,4} $. Другими словами, нумерация выбрана таким образом, чтобы каждый собственный вектор в общем случае переходил при ${{\varphi }_{y}} \to 0$ в соответствующий собственный вектор при ${{\varphi }_{y}} = 0$.

Замечание 2. Фактически любая линейная комбинация ${{\vec {v}}_{j}} = A{{\vec {v}}_{{j,1}}} + B{{\vec {v}}_{{j,2}}}$ – собственный вектор, отвечающий собственному значению λj.

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

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

Первое представление ${{\vec {u}}_{1}}\left( {x;y,z} \right)$ характеризуется следующим видом решения для каждого из слоев:

(4.6)
$\vec {u}_{1}^{с}\left( {x;y,z} \right) = (A_{3}^{c} \cdot \vec {v}_{{3,c}}^{1} + A_{4}^{c} \cdot \vec {v}_{{4,c}}^{1}){{e}^{{ - \gamma _{c}^{1} \cdot \left( {x - h} \right)}}},$
(4.7)
$\begin{gathered} \vec {u}_{1}^{f}\left( {x;y,z} \right) = (A_{1}^{f} \cdot \vec {v}_{{1,f}}^{1} + A_{2}^{f} \cdot \vec {v}_{{2,f}}^{1}){{e}^{{\gamma _{f}^{1} \cdot x}}} + \\ \, + (A_{3}^{f} \cdot \vec {v}_{{3,f}}^{1} + A_{4}^{f} \cdot \vec {v}_{{4,f}}^{1}){{e}^{{ - \gamma _{f}^{1}tx}}}, \\ \end{gathered} $
(4.8)
$\vec {u}_{1}^{s}\left( {x;y,z} \right) = (A_{1}^{s} \cdot \vec {v}_{{1,s}}^{1} + A_{2}^{s} \cdot \vec {v}_{{2,s}}^{1}){{e}^{{\gamma _{s}^{1} \cdot x}}},$
где $\vec {v}_{{j,\alpha }}^{1} = {{\left. {\vec {v}_{j}^{1}} \right|}_{{\varepsilon = {{\varepsilon }_{\alpha }},\mu = {{\mu }_{\alpha }}}}}$, $j = \overline {1,4} $ и $\gamma _{\alpha }^{1} = \mathop {\left. {{{\gamma }^{0}}} \right|}\nolimits_{\varepsilon = {{\varepsilon }_{\alpha }},u = {{\mu }_{\alpha }}} $ при $\alpha = c,f,s$, причем $\vec {v}_{j}^{0}$ определены согласно формулам (4.1).

Второе представление ${{\vec {u}}_{2}}(x;y,z)$ характеризуется следующим видом решения для каждого из слоев:

(4.9)
$\vec {u}_{2}^{c}\left( {x;y,z} \right) = (B_{3}^{c} \cdot \vec {v}_{{3,c}}^{2} + B_{4}^{c} \cdot \vec {v}_{{4,c}}^{2}){{e}^{{ - \gamma _{c}^{2} \cdot \left( {x - h} \right)}}},$
(4.10)
$\begin{gathered} \vec {u}_{2}^{f}\left( {x;y,z} \right) = (B_{1}^{f} \cdot \vec {v}_{{1,f}}^{2} + B_{2}^{f} \cdot \vec {v}_{{2,f}}^{2}){{e}^{{\gamma _{f}^{2} \cdot x}}} + \\ \, + (B_{3}^{f} \cdot \vec {v}_{{3,f}}^{2} + B_{4}^{f} \cdot \vec {v}_{{4,f}}^{2}){{e}^{{ - \gamma _{f}^{2}tx}}}, \\ \end{gathered} $
(4.11)
$\vec {u}_{2}^{s}\left( {x;y,z} \right) = (B_{1}^{s} \cdot \vec {v}_{{1,s}}^{2} + B_{2}^{s} \cdot \vec {v}_{{2,s}}^{2}){{e}^{{\gamma _{s}^{2} \cdot x}}},$
где $\vec {v}_{{j,\alpha }}^{2} = {{\left. {\vec {v}_{{j,1}}^{{}}} \right|}_{{\varepsilon = {{\varepsilon }_{\alpha }},\mu = {{\mu }_{\alpha }}}}}$, $j = \overline {1,4} $ и $\gamma _{\alpha }^{2} = \mathop {\left. \gamma \right|}\nolimits_{\varepsilon = {{\varepsilon }_{\alpha }},\mu = {{\mu }_{\alpha }}} $ при $\alpha = c,f,s$, причем ${{\vec {v}}_{{j,1}}}$ определены согласно формулам (4.2)(4.5).

Третье представление ${{\vec {u}}_{3}}\left( {x;y,z} \right)$ характеризуется следующим видом решения для каждого из слоев:

(4.12)
$\vec {u}_{3}^{c}\left( {x;y,z} \right) = (D_{3}^{c} \cdot \vec {v}_{{3,c}}^{3} + D_{4}^{c} \cdot \vec {v}_{{4,c}}^{3}){{e}^{{ - \gamma _{c}^{3} \cdot \left( {x - h} \right)}}},$
(4.13)
$\begin{gathered} \mathop {\vec {u}}\nolimits_3^f \left( {x;y,z} \right) = (D_{1}^{f} \cdot \vec {v}_{{1,f}}^{3} + D_{2}^{f} \cdot \vec {v}_{{2,f}}^{3}){{e}^{{\gamma _{f}^{3} \cdot x}}} + \\ \, + (D_{3}^{f} \cdot \vec {v}_{{3,f}}^{3} + D_{4}^{f} \cdot \vec {v}_{{4,f}}^{3}){{e}^{{ - \gamma _{f}^{3}tx}}}, \\ \end{gathered} $
(4.14)
$\vec {u}_{3}^{s}\left( {x;y,z} \right) = (D_{1}^{s} \cdot \vec {v}_{{1,s}}^{3} + D_{2}^{s} \cdot \vec {v}_{{2,s}}^{3}){{e}^{{\gamma _{s}^{3} \cdot x}}},$
где $\vec {v}_{{1,\alpha }}^{3} = {{\left. {\vec {v}_{1}^{ + }} \right|}_{{\varepsilon = {{\varepsilon }_{\alpha }},\mu = {{\mu }_{\alpha }}}}}$, $\vec {v}_{{2,\alpha }}^{3} = {{\left. {\vec {v}_{2}^{ + }} \right|}_{{\varepsilon = {{\varepsilon }_{\alpha }},\mu = {{\mu }_{\alpha }}}}}$, $\vec {v}_{{3,\alpha }}^{3}$ = = ${{\left. {\vec {v}_{1}^{ - }} \right|}_{{\varepsilon = {{\varepsilon }_{\alpha }},\mu = {{\mu }_{\alpha }}}}}$, $\vec {v}_{{4,\alpha }}^{3} = {{\left. {\vec {v}_{2}^{ - }} \right|}_{{\varepsilon = {{\varepsilon }_{\alpha }},\mu = {{\mu }_{\alpha }}}}},$ и $\gamma _{\alpha }^{3} = \mathop {\left. \gamma \right|}\nolimits_{\varepsilon = {{\varepsilon }_{\alpha }},\mu = {{\mu }_{\alpha }}} $ при $\alpha = c,f,s$, причем $\vec {v}_{{1,2}}^{ \pm }$ определены согласно формулам (3.7), (3.8).

Для каждого из трех представлений полей согласно разделу 3.4 формируется система граничных уравнений (3.22), (3.23), после чего вычисляются коэффициенты фазового замедления и коэффициенты, определяющие поля волноводных мод. Полученные численно величины коэффициентов фазового замедления и полей соответствующих волноводных мод анализируются в разделе ниже.

4.3. Численные результаты

Ниже приведем вычисленные c точностью 10–14 коэффициенты фазового замедления для первого случая ${{\beta }_{{s,1}}}$, где $s = \overline {1,4} $:

(4.15)
$\begin{aligned} & {{\beta }_{{1,1}}} = 1.55149273806928903, \\ & {{\beta }_{{2,1}}} = 1.55018111589009942, \\ & {{\beta }_{{3,1}}} = 1.51175061453743748, \\ & {{\beta }_{{4,1}}} = 1.50727495127641732. \\ \end{aligned} $

Ниже приведены графики величин ${{\Delta }_{{2,3}}}{{\beta }_{s}}\left( \delta \right)$, определенных формулами (3.25), (3.26) для первой моды (s = 1).

Для старших мод $s = 2,3,4$ поведение величин ${{\Delta }_{{2,3}}}{{\beta }_{s}}\left( \delta \right)$ полностью идентично их поведению при s = 1, поэтому ограничимся только одним графиком.

Далее приведем графики величин ${{\Delta }_{{2,3}}}{{u}_{s}}\left( \delta \right)$, определенных формулами (3.27), (3.28) для четырех направляемых мод $s = \overline {1,4} $.

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

5.1. Символьные результаты

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

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

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

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

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

Настоящее исследование не претендует на исчерпывающее множество всех различных символьных представлений собственных векторов, а скорее дает минимальное множество, которое получается при преобразованиях перестановок. Это минимальное множество тем не менее позволило получить собственные векторы, которые обеспечивают корректный предельный переход при обнулении символьных констант модели. Однако существует значительно больше преобразований – например, преобразования поворота – которые могут дать большее количество различных символьных представлений собственных векторов и дать исследователю еще больше вариантов различных символьных представлений собственных векторов для дальнейшего численного и качественного анализа.

5.2. Численные результаты

Замечение. Все приведенные выше графики необходимо читать в направлении убывания δ – то есть справа налево.

В рамках численных расчетов демонстрируется практическая польза от разработанного символьного алгоритма (см. раздел 3). В рамках расчетов показано, что при любом выборе собственных векторов – как с ипользованием предложенного алгоритма, так и с использованием встроенной функции Eigenvectors – коэффициенты фазового замедления сходятся при уменьшении δ к значениям при δ = 0 (см. рис. 1).

Рис. 1.

Графики величин ${{\Delta }_{{2,3}}}{{\beta }_{s}}\left( \delta \right)$, определенных формулами (3.25), (3.26) для первой моды (s = 1).

Практическая польза от применения алгоритма видна при рассмотрении графиков сходимостей полей направляемых волноводных мод – см. рис. 2–5. Как видно из рисунка 2, погрешность вычисленных полей первой волноводной моды без использования предложенного алгоритма ${{\Delta }_{3}}{{u}_{1}}\left( \delta \right)$ на несколько порядков меньше, чем при использовании предложенного алгоритма ${{\Delta }_{2}}{{u}_{1}}\left( \delta \right)$, эти величины имеют порядок 10–5 и 10–12 соответственно.

Рис. 2.

Графики величин ${{\Delta }_{{2,3}}}{{u}_{s}}\left( \delta \right)$, определенных формулами (3.27), (3.28) для первой моды (s = 1).

Рис. 3.

Графики величин ${{\Delta }_{{2,3}}}{{u}_{s}}\left( \delta \right)$, определенных формулами (3.27), (3.28) для второй моды (s = 2).

Рис. 4.

Графики величин ${{\Delta }_{{2,3}}}{{u}_{s}}\left( \delta \right)$, определенных формулами (3.27), (3.28) для третьей моды (s = 3).

Рис. 5.

Графики величин ${{\Delta }_{{2,3}}}{{u}_{s}}\left( \delta \right)$, определенных формулами (3.27), (3.28) для четвертой моды (s = 4).

Рисунок 3 демонстрирует эффективность использования предложенного алгоритма: погрешность вычисленного поля второй волноводной моды без использования предложенного алгоритма ${{\Delta }_{3}}{{u}_{2}}\left( \delta \right)$ растет с уменьшением параметра δ и достигает значения порядка 1.1 при $\delta = {{10}^{{ - 17}}}$, причем сравниваемые поля имеют единичную норму (то есть относительная погрешность более 100%). При этом погрешность вычисленного поля второй волноводной моды с использованием предложенного алгоритма ${{\Delta }_{2}}{{u}_{2}}\left( \delta \right)$ имеет порядок 10–5 при $\delta = {{10}^{{ - 17}}}$.

Рисунок 4 показывает, что обе величины ${{\Delta }_{3}}{{u}_{3}}\left( \delta \right)$ и ${{\Delta }_{2}}{{u}_{3}}\left( \delta \right)$ достаточно малы для численных расчетов и имеют порядок 10–12 и 10–5 соответственно.

Рисунок 5 демонстрирует эффективность использования предложенного алгоритма: величина ${{\Delta }_{3}}{{u}_{4}}\left( \delta \right)$ растет с уменьшением параметра δ и достигает значения порядка 1.3, при том что сравниваемые поля имеют единичную норму (то есть относительная погрешность более 100%) при $\delta = {{10}^{{ - 17}}}$. При этом относительная погрешность вычисленного поля четвертой волноводной моды с использованием предложенного алгоритма ${{\Delta }_{2}}{{u}_{4}}\left( \delta \right)$ имеет порядок 10–5.

Другими словами, использование собственных векторов, полученных в символьном виде без использования предложенного в разделе 3 алгоритма, позволяет корректно вычислять коэффициенты фазового замедления направляемых волноводных мод при малых δ. Однако поля соответствующих волноводных мод могут вычисляться с недопустимо большими погрешностями, достигающими порядка 100% и более, и потому предложенный алгоритм необходим для проведения корректных практических расчетов.

Расчеты с использованием разработанного символьного алгоритма показывают невысокую точность: относительная погрешность вычисленных полей направляемых волноводных мод составляет 10–5. Однако, эта величина относительной погрешности сохраняется для всех вычисленных полей направляемых волноводных мод, в то время как вычисления без использования разработанного алгоритма могут характеризоваться относительными погрешностями более единицы.

Отдельного внимания также заслуживает вопрос неравномерной сходимости вычисленных полей. По мнению авторов численные артефакты, которые наблюдаются на рис. 2–5, связаны с тем, что при расчетах были использованы символьные представления собственных векторов без дополнительной численной нормировки. Однако вопрос разработки устойчивого численного алгоритма вычисления полей выходит за рамки настоящего исследования и будет отдельно рассмотрен в рамках будущих исследований.

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

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

Численные расчеты с применением предложенного алгоритма демонстрируют его эффективность на фоне аналогичных результатов, полученных без его использования. На примере задачи вычисления полей направляемых волноводных мод планарного трехслойного волновода показано, что поля волноводных мод, вычисленные без применения предложенного алгоритма, могут вычисляться с недопустимо большими погрешностями, достигающими порядка 100% и более.

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

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

  1. Sevastyanov L.A., Sevastyanov A.L., Tyutyunnik A.A. Analytical Calculations in Maple to Implement the Method of Adiabatic Modes for Modelling Smoothly Irregular Integrated Optical Waveguide Structures // Lecture Notes in Computer Science. 2014. V. 8660. P. 419–431. https://doi.org/10.1007/978-3-319-10515-4_30

  2. Divakov D.V., Sevastianov A.L. The Implementation of the Symbolic-Numerical Method for Finding the Adiabatic Waveguide Modes of Integrated Optical Waveguides in CAS Maple. // Lecture Notes in Computer Science. 2019. V. 11661. P. 107–121. https://doi.org/10.1007/978-3-030-26831-2_8

  3. Yee K. Numerical solution of initial boundary value problems involving maxwell’s equations in isotropic media // IEEE Transactions on Antennas and Propagation. 1966. V. 14 (3.1). P. 302–307. https://doi.org/10.1109/TAP.1966.1138693

  4. Taflove A. Application of the Finite-Difference Time-Domain Method to Sinusoidal Steady-State Electromagnetic-Penetration Problems // IEEE Transactions on Electromagnetic Compatibility. 1980. V. EMC-22 (3.1). P. 191–202. https://doi.org/10.1109/TEMC.1980.30387

  5. Joseph R., Goorjian P., Taflove A. Direct time integration of Maxwell’s equations in two-dimensional dielectric waveguides for propagation and scattering of femtosecond electromagnetic solitons // Optics Letters. 1993. V. 18 (3.5). P. 491–493. https://doi.org/10.1364/OL.18.000491

  6. Sveshnikov A.G. The incomplete Galerkin method. Dokl. Akad. Nauk SSSR. 1977. V. 236 (3.3). P. 1076–1079.

  7. Petukhov A.A. Joint application of the incomplete Galerkin method and scattering matrix method for modeling multilayer diffraction gratings. Mathematical Models and Computer Simulations. 2014. V. 6. P. 92–100. https://doi.org/10.1134/S2070048214010128

  8. 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 11661, 412–429 (2019). https://doi.org/10.1007/978-3-030-26831-2_27

  9. Kantorovich L.V., Krylov V.I. Approximate Methods of Higher Analysis. Wiley, New York, 1964.

  10. Gusev A.A., Chuluunbaatar O., Vinitsky S.I., Derbov V.L. Solution of the Boundary-Value Problem for a Systems of ODEs of Large Dimension: Benchmark Calculations in the Framework of Kantorovich Method // Discrete and Continuous Models and Applied Computational Science. 2016. № 3. P. 31–37.

  11. Bathe K.J. Finite Element Procedures in Engineering Analysis. Prentice Hall, Englewood Cliffs. 1982.

  12. Bogolyubov A.N., Mukhartova Yu.V., Gao J., Bogolyubov N.A. Mathematical modeling of plane chiral waveguide using mixed finite elements // Progress in Electromagnetics Research Symposium. 2012. 1216–1219.

  13. Babich V.M., Buldyrev V.S. Asymptotic Methods in Short-Wave Diffraction Problems. Nauka, Moscow, 1972. [English translation: Springer Series on Wave Phenomena 4. Springer, Berlin Heidelberg New York, 1991.]

  14. Mathematics-based software and services for education, engineering, and research https://www.maplesoft.com/

  15. Hamming R.W. Numerical Methods for Scientists and Engineers. Dover Publications; 2nd Revised ed. edition. 1987.

  16. Reed M., Simon B. Methods of Modern Mathematical Physics. I: Functional Analysis (v. 1). Academic Press; 1st edition. 1972.

  17. Lang S. Real and Functional Analysis. Springer-Verlag New York; 3rd edition, 1993.

  18. Hardy G.H. A course of pure mathematics. Cambridge, At the University Press; 3rd edition, 1921.

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