Программирование, 2022, № 2, стр. 23-32

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

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

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

В рамках чисто численных подходов после постановки задачи сразу же следует построение приближенного метода ее решения – либо разностного [3–5], либо методов, основанных на разложениях Галеркина [6–8] (реже – методом Ритца [9]), Канторовича [10, 11] или метода конечных элементов [12, 13]. Однако, даже для генерации разностных схем также применяются системы компьютерной алгебры [14, 15].

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

1.3. Задача текущего исследования

Одним из первых методов исследования спектральных характеристик плавно-нерегулярных волноводов является метод поперечных сечений, предложенный в работах Каценеленбаума Б.З. [16, 17] для закрытых волноводов и обобщенный Шевченко В.В. на открытые волноводы [18, 19].

Особенностью метода поперечных сечений является символьный вид решения, в котором в явном виде присутствует искомая функция $\gamma \left( z \right)$, определяющая коэффициент фазового замедления соответствующей волноводной моды в каждом поперечном сечении волновода $z = {\text{const}}$. Для определения электромагнитного поля данной волноводной моды необходимо решать спектральную задачу специального вида, которая в фиксированном поперечном сечении волновода представляет собой задачу Штурма–Лиувилля для дифференциального оператора второго порядка на оси.

Аналогичную спектральную задачу можно получить и в рамках других подходов, например, как это сделано в работе авторов [20]. В рамках настоящей работы рассматривается система функциональных уравнений, получаемых из спектральной задачи, возникающей из метода поперечных сечений. В работе эта система функциональных уравнений исследуется на символьном уровне, для исследования разработан символьный алгоритм, реализованный в системе компьютерной алгебры Maple [21]. Зависимость $\gamma \left( z \right)$ в работе отыскивается численно.

2.1. Модель распространения TE-моды в волноводе

Рассмотрим волноводное распространение TE-моды в однородном открытом диэлектрическом волноводном переходе типа “рупор”, изображенном на рисунке 1.

Рассматривать будем волноводный переход толщина которого линейно увеличивается от значения h₁ при $z = 0$ до значения ${{h}_{2}}$ при $z = d$: $h(z)$ = = h₁ + δz, где $\delta = \frac{{{{h}_{2}} - {{h}_{1}}}}{d} \ll 1$ – малый параметр. Ведущая компонента TE-моды E_y удовлетворяет уравнению Гельмгольца в области волноводного перехода, что соответствует $0 \leqslant z \leqslant d$ [22]:

где $u = u\left( {x,z} \right)$ – искомая величина, а $n\left( {x,z} \right)$ определена следующим образом:

Асимптотические условия, отвечающие направляемым модам определены следующим образом [22]: . На границах раздела сред для TE-моды должны выполняться условия сопряжения [22]:

где введено обозначение ${{\left. {\left[ f \right]} \right|}_{{x = c}}} = f(c - 0)$ – f(c + 0), а $\vec {n}$ – нормаль к поверхности, поэтому:

2.2. Разложение решения в асимптотический ряд

Так как в уравнении и граничных условиях присутствует малый параметр $\delta $ можно строить решение в виде ряда по степеням параметра $\delta $: $u = {{u}_{0}} + \delta {{u}_{1}} + {{\delta }^{2}}{{u}_{2}} + \; \ldots $. Для такого вида решения можно сформулировать уравнение:

В нулевом приближении уравнение, граничные и асимптотические условия принимают следующий вид:

2.3. Аналог разделения переменных в нулевом приближении

Классический подход к построению решения такой задачи состоит в попытке разыскать решение в виде произведения функций ${{u}_{0}}(x,z) = X(x)Z(z)$ и подстановке в уравнение из (1) с дальнейшим преобразованием уравнения в попытке разделить переменные. В данном случае аналогичная процедура приводит к следующему соотношению:

где левая часть зависит только от $z$, а правая зависит от $x$ и $z$, в результате обе части должны равняться некоторой, пока не определенной функции $ - k_{0}^{2}\gamma \left( z \right)$. В результате можно сформулировать уравнение

Замечание 1. Из уравнения (2) следует, что функция X должна зависеть от $x$ и $z$, что, на первый взгляд, противоречит предположению X = = X(x). Однако, если X зависит от $x$ и $z$, а величина $\left| {\frac{{\partial X}}{{\partial z}}} \right|$ имеет порядок ${{\delta }^{p}}$, где $p \geqslant 1$, то в нулевом приближении по $\delta $ выкладки вполне справедливы. Другими словами X зависит от $x$ и $z$, однако для справедливости уравнения (1.8) необходимо, чтобы зависимость X от $z$ была настолько “медленной”, чтобы выполнялось

Такую зависимость будем обозначать $X = X(x;z)$.

Замечание 2. Учитывая, что $\left| {\frac{{\partial X}}{{\partial z}}} \right|$ имеет порядок ${{\delta }^{p}}$, то в приближении порядка $p$ вклад производной $\frac{{\partial X}}{{\partial z}}$ необходимо будет учитывать.

2.4. Спектральная задача в нулевом приближении

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

Искомыми величинами в задаче (3) выступают функция двух переменных X(x, z) и функция одной переменной $\gamma \left( z \right)$. Общее решение уравнения из (3), удовлетворяющее асимптотическому условию имеет следующий вид:

где ${{p}_{c}}\left( z \right)$, ${{p}_{f}}\left( z \right)$ и ${{p}_{s}}\left( z \right)$ определены следующим образом:

Для выполнения условий сопряжения из (3) на границах $x = 0,\;h\left( z \right)$ необходимо чтобы выполнялись следующие функциональные уравнения

причем искомыми величинами выступают $\gamma \left( z \right)$ и ${{A}_{c}}\left( z \right)$, ${{A}_{f}}\left( z \right)$, ${{B}_{f}}\left( z \right)$, ${{A}_{s}}\left( z \right)$.

Замечание. Зависимость $\gamma \left( z \right)$ в рамках настоящей работы отыскивается численно методом деления отрезка пополам, поэтому отысканию этой величины не будет уделено особого внимания.

Систему (6) можно записать в следующем виде:

где коэффициенты матрицы $M\left( {\gamma \left( z \right)} \right)$ зависят от искомой величины $\gamma \left( z \right)$ нелинейно, но при этом система линейна относительно искомых величин $\vec {A}\left( z \right) = {{\left( {{{A}_{c}}\left( z \right),{{A}_{f}}\left( z \right),{{B}_{f}}\left( z \right),{{A}_{s}}\left( z \right)} \right)}^{{\text{т}}}}$.

При любом фиксированном $z = {{z}_{j}}$ система (7) принимает вид однородной системы уравнений вида

у которой всегда существует тривиальное решение ${{\vec {A}}_{j}} = \vec {0}$.

В случае, если $\det M\left( {{{\gamma }_{j}}} \right) \ne 0$, то тривиальное решение единственно, а если $\det M\left( {{{\gamma }_{j}}} \right) = 0$, то у системы (8) кроме тривиального существует хотя бы одно нетривиальное решение.

Численные алгоритмы решения системы (7) как правило основываются на введении сетки z_j, $1 \leqslant j \leqslant N$ и многократном численном решении однородной системы типа (8) при каждом фиксированном z_j. При этом численное решение системы (8) состоит из двух задач:

Задача 1. Приближенное отыскание всех значений $\gamma _{j}^{s}$, являющихся корнями нелинейного уравнения $\det M\left( {{{\gamma }_{j}}} \right) = 0$.

2. Приближенное отыскание нетривиального решения однородной системы линейных алгебраических уравнений $M(\gamma _{j}^{s})\vec {A}_{j}^{s} = \vec {0}$ для каждого s.

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

При анализе численных результатов авторами было замечено, что присутствует некая функциональная зависимость величин ${{A}_{c}}\left( z \right)$, ${{A}_{f}}\left( z \right)$, ${{B}_{f}}\left( z \right)$, ${{A}_{s}}\left( z \right)$ от ${{p}_{f}}\left( z \right)$, ${{p}_{s}}\left( z \right)$. Поэтому у авторов возникла гипотеза, что искомые величины системы (7) могут выражаться через ${{p}_{f}}\left( z \right)$ и ${{p}_{s}}\left( z \right)$ в виде простых символьных выражений. Другими словами покажем, что система (7) может иметь символьное решение в зависимости от $\gamma \left( z \right)$. Для проверки этой гипотезы был разработан символьный алгоритм исследования системы функциональных уравнений (7), который приведен в следующем разделе. Если бы у рассматриваемой системы существовало символьное решение, которое в явном виде зависит от $\gamma \left( z \right)$ и состоит из достаточно компактных символьных выражений, то оно бы избавило исследователей от отыскания решения некорректной задачи путем многократного решения системы (8) на сетке. Также наличие символьного решения упростило бы качественный анализ полученных результатов.

3.1. Численный метод отыскания $\gamma \left( z \right)$

В работе рассматривается система (7), матрица которой имеет следующий вид

где все ${{p}_{\alpha }} = {{p}_{\alpha }}\left( z \right)$, $\alpha = c,s,f$ зависят от $\gamma \left( z \right)$ и определены в (5). Отыскание $\gamma \left( z \right)$ осуществляется численно с помощью введения сетки z_j, $1 \leqslant j \leqslant N$ и многократного численного решения нелинейного уравнения $\det M\left( {{{\gamma }_{j}}} \right) = 0$ для каждого z_j методом деления отрезка пополам [24].

Замечание. Важно заметить, что фактически мы отыскиваем не $\gamma \left( z \right)$, а ${{\tilde {\gamma }}_{j}}$ такое, что $\left| {{{{\tilde {\gamma }}}_{j}}\, - \,\gamma ({{z}_{j}})} \right|\, < \,\varepsilon $, для $1 \leqslant j \leqslant N$, где $\varepsilon $ – погрешность приближенного решения. Однако, важно учесть, что

1. Используя числа повышенной разрядности можно добиться произвольно малой для практических расчетов погрешности $\varepsilon $;

2. Увеличивая количество узлов сетки $N$ можно добиться такой близости соседних узлов z_j и ${{z}_{{j + 1}}}$, что кусочно-линейная аппроксимация $\tilde {\gamma }\left( z \right)$ гладкой функции $\gamma \left( z \right)$ была сколь угодно близка к функции $\gamma \left( z \right)$ в любой метрике.

Учитывая п. 1–2 замечания можно сказать, что величину $\gamma \left( z \right)$ можно определить с любой наперед заданной точностью. Учитывая также, что при численных расчетах все величины округляются до чисел конечной разрядности – как правило до чисел двойной точности – можно добиться одинаковой точности и величины $\gamma \left( z \right)$ и всех остальных величин, принимающих участие в численных расчетах. Поэтому величину $\gamma \left( z \right)$ можно считать численно известной.

Замечание. Как правило уравнение ${\text{det}}M(\gamma (z))$ = = 0 имеет несколько решений ${{\gamma }_{1}}\left( z \right),{{\gamma }_{2}}\left( z \right),\; \ldots $ – тем больше, чем больше толщина волновода.

3.2. Символьный метод отыскания $\vec {A}\left( z \right)$

Рассмотрим символьный алгоритм исследования системы (6), предлагаемый авторами. Учитывая изложенное в предыдущем разделе, величину $\gamma \left( z \right)$ будем считать известной. Для простоты изложения в настоящем разделе опустим зависимость от z у величин, входящих в уравнения (6).

Итак, обратим внимание на первые два уравнения системы (6)

которые представляют собой систему двух однородных уравнений с тремя неизвестными. Умножим первое уравнение на $\left( { - {{p}_{s}}} \right)$ и сложим со вторым уравнением, в результате получим уравнение которое можно сразу решить где A – произвольная символьная константа. Учитывая (12) из (10), можно сразу определить и величину ${{A}_{s}}$ и далее, зная ${{A}_{f}}$, ${{B}_{f}}$ и ${{A}_{s}}$, из оставшихся уравнений системы (6) можно определить ${{A}_{c}}$.

Замечание 1. Если начать рассмотрение системы (6) с последних двух уравнений, то сначала найдем из этих уравнений ${{A}_{f}}$, ${{B}_{f}}$ и ${{A}_{c}}$, а затем из первых двух определим ${{A}_{s}}$. Другими словами, из подсистемы, состоящей из двух однородных уравнений с тремя неизвестными, можно определить три неизвестных величины – то есть решить эту подсистему. Учитывая это, при реализации символьного алгоритма реализован поиск подсистем из двух уравнений с тремя неизвестными среди рассматриваемых уравнений и поочередное решение этих подсистем.

Замечание 2. После отыскания компонент решения ${{A}_{f}}$, ${{B}_{f}}$, ${{A}_{s}}$ (или ${{A}_{f}}$, ${{B}_{f}}$, ${{A}_{c}}$) оставшаяся единственная неизвестная компонента ${{A}_{c}}$ (или ${{A}_{s}}$) определяется из двух оставшихся уравнений – каждое с одной неизвестной. Однако, для определения одной неизвестной достаточно одного уравнения. Поэтому из каждого из двух оставшихся уравнений единственная неопределенная величина может определяться по-разному. Эти возможности учтены при реализации символьного алгоритма.

Замечание 3. Описанный алгоритм обобщен на случай систем уравнений произвольной размерности.

Описанный выше алгоритм реализован в системе компьютерной алгебры Maple. Он позволяет отыскивать различные представления символьного решения однородной системы уравнений типа (6). Псевдокод алгоритма приведен ниже.

3.3. Численный метод отыскания $\vec {A}\left( z \right)$

Численное отыскание $\vec {A}\left( z \right)$ состоит в приближенном отыскании нетривиального решения однородной системы $M\left( {{{{\tilde {\gamma }}}_{j}}} \right){{\vec {A}}_{j}} = \vec {0}$ на сетке ${{z}_{j}}$, $1 \leqslant j \leqslant N$, на которой приближенно определена величина $\gamma \left( z \right)$. Нетривиальное решение системы $M\left( {{{{\tilde {\gamma }}}_{j}}} \right){{\vec {A}}_{j}} = \vec {0}$ есть собственный вектор, отвечающий почти нулевому собственному значению матрицы $M\left( {{{{\tilde {\gamma }}}_{j}}} \right)$. Поэтому численное отыскание нетривиального решения осуществляется с помощью встроенного в Maple метода Eigenvectors пакета LinearAlgebra, который используюет QR-алгоритм для решения задачи.

Замечание 1. Задача отыскания нетривиального решения однородной системы $M\vec {A} = \vec {0}$ некорректна. Действительно, если вектор $\vec {A}$ есть нетривиальное решение этой системы, то вектор $\vec {B} = c \cdot \vec {A}$ также будет нетривиальным ее решением для любой константы $c \ne 0$.

Замечание 2. Для решения задачи (7) необходимо многократно для каждого $j$ на сетке ${{z}_{j}}$, $1 \leqslant j \leqslant N$ решать систему $M\left( {{{{\tilde {\gamma }}}_{j}}} \right){{\vec {A}}_{j}} = \vec {0}$. Имея два решения в соседних узлах сетки ${{\vec {A}}_{j}}$ и ${{\vec {A}}_{{j + 1}}}$ ($\left\| {{{{\vec {A}}}_{j}}} \right\| = \left\| {{{{\vec {A}}}_{{j + 1}}}} \right\|$ = 1) для обеспечения непрерывности необходимо подобрать константу $\alpha \in \mathbb{C}$, $\left| \alpha \right| = 1$ такую, что $\left\| {\alpha {{{\vec {A}}}_{{j + 1}}} - {{{\vec {A}}}_{j}}} \right\| \to \min $.

Учитывая замечания 1, 2 кроме отыскания нетривиальных решений систем $M\left( {{{{\tilde {\gamma }}}_{j}}} \right){{\vec {A}}_{j}} = \vec {0}$ для каждого j необходимо также на каждом шаге подбирать соответствующую константу $\alpha $. Решение этой задачи подбора константы подробно описано в разделе 3.5 работы [25].

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

Применим разработанный символьный алгоритм к рассматриваемой однородной системе (6). В результате получим следующие различные символьные представления решения однородной системы (6):

Замечание. Решение однородной системы (6) можно отыскивать как собственный вектор матрицы коэффициентов системы, отвечающий нулевому ее собственному значению. Однако встроенный в Maple метод Eigenvectors пакета LinearAlgebra не позволяет отыскать символьные представления собственных векторов: после получаса работы метода результат вычислен не был. Тем не менее, собственные значения вычисляются методом Eigenvalues, однако символьные представления собственных значений превышают допустимый лимит длины отображения выражений в Maple. Другими словами, выражения для собственных значений слишком длинные для отображения, и, соответственно, для анализа.

Приведем также символьное выражение для определителя матрицы M, полученное с помощью встроенного в Maple метода Determinant пакета LinearAlgebra:

После подстановки каждого из полученных символьных решений ${{\vec {A}}_{j}}$ в систему уравнений получим следующие вектора невязок ${{\vec {\delta }}_{j}} = M{{\vec {A}}_{j}}$:

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

Вычислим сначала приближенно величину $\gamma (z)$, а затем вычислим $\vec {A}\left( z \right)$ на основе явных символьных выражений (13)–(16) и на основе численного алгоритма из раздела 3.3. Также численно сравним между собой величины, вычисленные по каждой из четырех формул (13)–(16).

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

Ниже приведены графики величин ${{\gamma }_{s}}\left( z \right)$, вычисленных на отрезке $z \in \left[ {0,d} \right]$, $d = 100\lambda $ на сетке z_j, $1 \leqslant j \leqslant N$ для $N = 64$ с точностью $\varepsilon {{ = 10}^{{ - 15}}}$.

Для каждого из вычисленных ${{\gamma }_{s}}\left( z \right)$ вычислим на той же сетке численно $\vec {A}\left( z \right)$ и сравним с ${{\vec {A}}_{m}}\left( z \right)$, рассчитанными по формулам (13)–(16). Для сравнения введем следующую метрику:

где ${{\left\| {\; \cdot \;.} \right\|}_{E}}$ – евклидова векторная норма. Для первой моды (s = 1) данные приведены ниже.

Для $s = 2$ (вторая мода) данные приведены ниже.

Кроме того, приведем результаты сравнения величин ${{\vec {A}}_{m}}\left( z \right)$ между собой. Введем матрицу ${{\Delta }^{1}}$ (для $s = 1$) с элементами $\Delta _{{m,l}}^{1} = \left\| {{{{\vec {A}}}_{m}}\left( z \right) - {{{\vec {A}}}_{l}}\left( z \right)} \right\|$:

и аналогичную матрицу Δ² для $s = 2$:

5.2. Обсуждение численных результатов

В рамках численных экспериментов мы дали ответ на вопрос о численной эквивалентности всех символьных представлений решения системы (6). Как видно из (23), (24) все четыре символьных представления решения (13)–(16) дают одинаковые с точностью до константы решения, отличия составляют не более $3.2 \times {{10}^{{ - 13}}}$. Поэтому, не имеет значение какое из полученных символьных решений (13)–(16) системы (6) выбирать для анализа и использования в численных расчетах – здесь исследователь может выбрать наиболее удобное с точки зрения дальнейших манипуляций: представляется логичным выбрать самое компактное представление как для символьных манипуляций, так и для численных расчетов.

Важно также заметить, что символьные представления решения (13)–(16) системы (6) дают примерно одинаковую с численным решением точность порядка $2.3 \times {{10}^{{ - 13}}}$ (см. (21), (22)). Однако, численный метод решения однородной системы (6) содержит специальные численные алгоритмы, обеспечивающие гладкость решения, так как задача отыскания решения системы (6) является некорректной.

Расчет решения по символьным формулам не требует использования специальных дополнительных шагов для обеспечения гладкости решения, ибо символьные выражения непрерывны и дифференцируемы столько раз, сколько дифференцируемы $\gamma (z)$ и h(z). Другими словами, полученные символьные решения (13)–(16), благодаря непрерывной зависимости от непрерывных величин $\gamma (z)$ и h(z), автоматически обладают свойством непрерывности, в отличие от решений, получаемых численно.

В работе рассмотрен простейший случай, соответствующий трехслойному волноводу с одной плоской границей x = 0 и одной криволинейной $x = h(z)$, где функция h(z) описывает переменную толщину волноводного слоя. Существуют также и более сложные конфигурации с большим количеством слоев, для которых получение символьного решения системы, аналогичной системе (6), имеет большую важность. Для таких конфигураций характерно наличие трех слоев для одного интервала значений $z$, а для другого интервала значений $z$ слоев будет уже четыре. Другими словами, для четырехслойных конфигураций характерно наличие одного слоя переменной толщины, причем для определенного интервала значений $z$ эта толщина равна нулю (то есть структура трехслойная), а для другого интервала значений $z$ толщина плавно изменяется от нуля до своего максимального значения. На границе между этими интервалами необходимо также обеспечивать непрерывность решения, для чего также необходимо использовать разработанный алгоритм. Здесь, в частности, выявится предпочтительность одних символьных представлений решения перед другими.