Журнал вычислительной математики и математической физики, 2020, T. 60, № 10, стр. 1697-1710

Расчет частот собственных колебаний акустической среды внутри вытянутого сфероида модифицированным методом Абрамова

Т. В. Левитина^1, *

¹ Макс-Планк институт по изучению Солнечной системы
37077 Геттинген, Юстус-фон-Либиг-Вег 3, Германия

Поступила в редакцию 21.01.2020
После доработки 15.04.2020
Принята к публикации 03.06.2020

DOI: 10.31857/S0044466920100105

Аннотация

Представленный и исследованный в [1], [2] метод решения самосопряженных многопараметрических спектральных задач для слабо связанных систем обыкновенных дифференциальных уравнений основан на методе продолжения по параметру, вводимому в задачу. Несмотря на то что метод формально применим к системам обыкновенных дифференциальных уравнений с особенностями, непосредственное его использование для численного решения задачи, указанной в названии, ограничено. Предлагаемая ниже модификация метода позволяет проводить вычисления различных, в том числе высокочастотных, акустических колебаний внутри сфероидов как близких к сфере, так и сильно вытянутых. Библ. 15. Фиг. 3.

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

ВВЕДЕНИЕ

При моделировании акустических колебаний внутри вытянутого сфероида разделение переменных в трехмерном уравнении Гельмгольца в сфероидальных координатах приводит к двухпараметрической задаче Штурма–Лиувилля: разделенные угловое и радиальное волновые сфероидальные уравнения оказываются связаны парой спектральных параметров – константой разделения $\lambda $ и безразмерным волновым числом $c$, определяющим частоту колебательной моды. Каждая мода характеризуется мультииндексом $(l,n,m)$, $l,n,m = 0,1,2,...$, т. е. азимутальным индексом $m$ и числом внутренних нулей решений углового ($l$) и радиального ($n$) уравнений. Разрешимость задачи, несмотря на ее сингулярность, доказана для любого мультииндекса, но численное решение сопряжено с определенными трудностями.

Представленный в [1], [2] метод решения многопараметрических спектральных задач соединял ньютоновские итерации и движение по параметру, связующему более простую задачу с “замороженными” коэффициентами при спектральных параметрах с исходной задачей. В [1], [2] были доказаны разрешимость всех промежуточных задач и возможность вычислений ньютоновских итераций на каждом шаге на пути от “замороженной” задачи к исходной, среди прочего обратимость соответствующих матриц Якоби. Была получена оценка невязки расчетов в зависимости от максимальной длины промежуточных шагов. Тем не менее оказалось, что метод [1], [2] в его исходной форме неприменим к расчету частот мод типа “шепчущих галерей” (МШГ), важных для многочисленных практических приложений [3]–[5]. Прежде всего это связано с сингулярностью уравнений, возникающих в результате разделения переменных. Перенос условий ограниченности решений этих уравнений из особой точки в регулярную многократно обсуждался ранее [6]–[9]. Однако в отношении метода [1], [2] эта процедура требует дополнительных усилий. Также при численном моделировании высокочастотных мод, как например МШГ, приходится учитывать, что правые части вспомогательных уравнений могут принимать чрезвычайно большие значения, что приводит к резкому ступенчатому изменению решений этих уравнений. Введение так называемых сглаживающих функций помогает преодолеть связанные с этим осложнения.

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

1. РАЗДЕЛЕНИЕ ПЕРЕМЕННЫХ В УРАВНЕНИИ ГЕЛЬМГОЛЬЦА В ВЫТЯНУТЫХ СФЕРОИДАЛЬНЫХ КООРДИНАТАХ

Разделение переменных в уравнении Гельмгольца в вытянутых сфероидальных координатах многократно обсуждалось ранее [10], [11]. Ниже мы приводим лишь сведения, необходимые для дальнейшего изложения, используя обозначения [6]–[8]. Сами вытянутые сфероидальные координаты определим через их связь с декартовыми:

$x = \frac{d}{2}\sqrt {({{\xi }^{2}} - 1)(1 - {{\eta }^{2}})} cos\varphi ,\quad y = \frac{d}{2}\sqrt {({{\xi }^{2}} - 1)(1 - {{\eta }^{2}})} sin\varphi ,\quad z = \frac{d}{2}\xi \eta .$

Через $\varphi \in [0,2\pi )$ обозначен азимутальный угол, $\eta \in ( - 1,1)$ и $\xi \in (1,\infty )$ обозначают соответственно угловую и радиальную сфероидальные координаты. Координатные поверхности $\xi = {{\xi }_{0}}$ образуют софокусное семейство вытянутых сфероидов, а параметр $d$ определяет половину межфокусного расстояния.

В том случае, если на поверхности сфероидального резонатора $\mathcal{S}$: $\xi = {{\xi }_{\mathcal{S}}}$ заданы условия Дирихле или Неймана, разделение переменных в уравнении Гельмгольца

(1.1)

$ - \Delta W({\mathbf{r}}) = {{k}^{2}}W({\mathbf{r}}),\quad {\mathbf{r}} = (\varphi ,\eta ,\xi ),\quad \xi < {{\xi }_{\mathcal{S}}},$

в сфероидальных координатах представляется целесообразным. Частные решения уравнения (1.1) записываются в виде произведения

$W({\mathbf{r}}) = S(\eta )R(\xi )exp( \pm im\varphi ),\quad m = 0,1, \ldots \;.$

Функции $S(\eta )$ и $R(\xi )$ – это угловая и радиальная волновые вытянутые cфероидальные функции (ВВСФ). Они являются компонентами собственной функции двухпараметрической самосопряженной сингулярной задачи Штурма–Лиувилля (ЗШЛ), заданной уравнениями

(1.2)

$\frac{d}{{d\eta }}(1 - {{\eta }^{2}})\frac{d}{{d\eta }}S + {{Q}_{a}}(\eta ,\lambda ,{{c}^{2}})S = 0,\quad - {\kern 1pt} 1 < \eta < 1,$

(1.3)

$\frac{d}{{d\xi }}({{\xi }^{2}} - 1)\frac{d}{{d\xi }}R + {{Q}_{r}}(\eta ,\lambda ,{{c}^{2}})R = 0,\quad 1 < \xi < {{\xi }_{\mathcal{S}}}.$

Здесь $m$ – азимутальный индекс, $m = 0,1, \ldots $, константа разделения $\lambda $ и квадрат безразмерного волнового числa $c = kd{\text{/}}2$ суть компоненты двухпараметрического собственного значения (СЗ) $(\lambda ,{{c}^{2}})$,

(1.4)

${{Q}_{a}}(\eta ,\lambda ,{{c}^{2}}) = \lambda + {{c}^{2}}(1 - {{\eta }^{2}}) - \frac{{{{m}^{2}}}}{{1 - {{\eta }^{2}}}},$

(1.5)

${{Q}_{r}}(\xi ,\lambda ,{{c}^{2}}) = {{c}^{2}}({{\xi }^{2}} - 1) - \lambda - \frac{{{{m}^{2}}}}{{{{\xi }^{2}} - 1}}.$

Помимо уравнений, ВВСФ удовлетворяют следующим граничным условиям:

(1.6)

${\text{или}}\;\;S{\kern 1pt} '(0) = 0,\quad {\text{или}}\;\;S(0) = 0,$

(1.7)

${{\beta }_{a}}(\eta ) = (1 - {{\eta }^{2}})S{\kern 1pt} '(\eta ){\text{/}}S(\eta ) \to - m,\quad \eta \to {{1}_{{ - 0}}},$

(1.8)

${{\beta }_{r}}(\xi ) = ({{\xi }^{2}} - 1)R{\kern 1pt} '(\xi ){\text{/}}R(\xi ) \to m,\quad \xi \to {{1}_{{ + 0}}}.$

(1.9)

${\text{или}}\;\;R({{\xi }_{\mathcal{S}}}) = 0,\quad {\text{или}}\;\;R{\kern 1pt} '({{\xi }_{\mathcal{S}}}) = 0.$

Последнее условие непосредственно следует из граничного условия на поверхности сфероида $\mathcal{S}$. Условие (1.6) является следствием симметрии задачи: угловые ВВСФ должны быть либо четными, либо нечетными. Это же обстоятельство позволяет в дальнейшем ограничиться рассмотрением уравнения (1.2) на полуинтервале $0 \leqslant \eta < 1$.

Условия (1.7) и (1.8) выражают тот факт, что физический смысл имеют лишь решения уравнений (1.2), (1.3), ограниченные в особых точках $\eta = 1$ и $\xi = 1$, а именно те, чье поведение в окрестности особой точки определяется асимптотикой [10]

$S(\eta ) \sim {{(1 - {{\eta }^{2}})}^{{m/2}}},\quad \eta \to 1,\quad R(\xi ) \sim {{({{\xi }^{2}} - 1)}^{{m/2}}},\quad \xi \to 1.$

Заметим, что случай $m = 0$ следует рассматривать отдельно. Однако все последующие формулы справедливы для произвольного значения азимутального индекса.

Непосредственной проверкой можно подтвердить, что введенные условиями (1.7), (1.8) функции ${{\beta }_{a}}(\eta )$ и ${{\beta }_{r}}(\xi )$, т.е. логарифмические производные ограниченных в особой точке решений уравнений (1.2) и (1.3), удовлетворяют уравнениям Риккати

(1.10)

(1.11)

В свою очередь, подстановкой решений задач

(1.12)

$\begin{gathered} \theta _{a}^{'} = \frac{{si{{n}^{2}}({{\theta }_{a}})}}{{(1 - {{\eta }^{2}})}} + {{Q}_{a}}(\eta ,\lambda ,{{c}^{2}})co{{s}^{2}}({{\theta }_{a}}),\quad \eta \in (0,1), \\ \mathop {lim}\limits_{\eta \to {{1}_{{ - 0}}}} {{\theta }_{a}}(\eta ) = arctan(m), \\ \end{gathered} $

(1.13)

$\begin{gathered} \theta _{r}^{'} = \frac{{si{{n}^{2}}({{\theta }_{r}})}}{{({{\xi }^{2}} - 1)}} + {{Q}_{r}}(\xi ,\lambda ,{{c}^{2}})co{{s}^{2}}({{\theta }_{r}}),\quad \xi \in (1,{{\xi }_{\mathcal{S}}}), \\ \mathop {lim}\limits_{\xi \to {{1}_{{ + 0}}}} {{\theta }_{r}}(\xi ) = - \arctan (m), \\ \end{gathered} $

в уравнения (1.10), (1.11) легко показать, что

(1.14)

$tan{{\theta }_{a}}(\eta ) = - {{\beta }_{a}}(\eta ),\quad {\text{и}}\quad tan{{\theta }_{r}}(\xi ) = - {{\beta }_{r}}(\xi ).$

Логарифмические производные ${{\beta }_{a}}(\eta )$, ${{\beta }_{r}}(\xi )$ и связанные с ними фазовые функции, иначе углы Прюфера, ${{\theta }_{a}}(\eta )$ и ${{\theta }_{r}}(\xi )$, играют ключевую роль в дальнейшем изложении.

В отличие от логарифмических производных, неограниченно растущих в окрестности точек, где решения уравнений (1.2), (1.3) обращаются в ноль, углы Прюфера ${{\theta }_{a}}(\eta )$, ${{\theta }_{r}}(\xi )$ остаются непрерывными на всем интервале определения. Таким образом, (1.12), (1.13) переносят условие ограниченности из особых точек $\eta = 1$, $\xi = 1$ в произвольную регулярную точку интервалов $[0,1)$ и $(1,{{\xi }_{\mathcal{S}}}]$ соответственно.

Очевидно, пара $(\lambda ,{{c}^{2}})$ образует двухпараметрическое СЗ задачи (1.2)–(1.9) в том и только том случае, когда функции ${{\theta }_{a}}(\eta )$, ${{\theta }_{r}}(\xi )$ отвечают граничным условиям (1.6) и (1.9), т.е., с точностью до слагаемого, кратного $\pi $,

$\begin{gathered} {\text{или}}\;\;{{\theta }_{a}}(0) = 0,\quad {\text{или}}\;\;{{\theta }_{a}}(0) = \pi {\text{/}}2, \\ {\text{или}}\;\;{{\theta }_{r}}(0) = \pi {\text{/}}2,\quad {\text{или}}\;\;{{\theta }_{r}}(0) = 0. \\ \end{gathered} $

Подобно классическому однопараметрическому случаю, каждому двухпараметрическому СЗ можно приписать мультииндекс $(l,n)$, где $l$ и $n$ показывают, сколько раз соответствующие угловая и радиальная ВВСФ обращаются в ноль на интервалах $( - 1,1)$ и $(1,{{\xi }_{\mathcal{S}}})$, причем четность индекса $l$ совпадает с четностью угловой ВВСФ.

Отметим, что всякий раз, когда решения уравнений (1.2), (1.3) обращаются в ноль, отвечающий им угол Прюфера пересекает горизонталь $\theta = \pi k + \pi {\text{/}}2$, $k \in \mathbb{Z}$, возрастая при этом. Следовательно, число горизонталей, которое пересекает угол Прюфера ${{\theta }_{a}}(\eta )$ на полуинтервале $(0,1)$, совпадает с числом нулей, которое имеется там у ограниченного решения уравнения (1.2). Также число горизонталей вида $\theta = \pi k + \pi {\text{/}}2$, которое пересекает ${{\theta }_{r}}(\xi )$ на интервале $(1,\xi )$, совпадает с числом нулей ограниченного решения уравнения (1.3) на этом интервале.

Таким образом, поиск СЗ $({{\lambda }_{{l,n}}},c_{{l,n}}^{2})$ задачи (1.2)–(1.9) сводится к системе алгебраических уравнений

(1.15)

${{\varphi }_{a}}(\lambda ,{{c}^{2}}) = 0,\quad {{\varphi }_{r}}(\lambda ,{{c}^{2}}) = 0,$

где

(1.16)

${{\varphi }_{a}}(\lambda ,{{c}^{2}}) = - {{\theta }_{a}}(0) - \frac{{l\pi }}{2},$

и, в зависимости от заданного на поверхности сфероида $\mathcal{S}$ граничного условия,

(1.17)

${{\varphi }_{r}}(\lambda ,{{c}^{2}}) = {{\theta }_{r}}({{\xi }_{\mathcal{S}}}) - \left\{ \begin{gathered} n\pi \quad ({\text{условие}}\;{\text{Неймана}}), \hfill \\ \left( {n + \frac{1}{2}} \right)\pi \quad ({\text{условие}}\;{\text{Дирихле}}). \hfill \\ \end{gathered} \right.$

Алгебраическая формулировка (1.15) позволяет применить метод Ньютона для поиска СЗ $({{\lambda }_{{l,n}}},c_{{l,n}}^{2})$.

2. ПОИСК СОБСТВЕННОГО ЗНАЧЕНИЯ: ДВИЖЕНИЕ ПО ПАРАМЕТРУ

Существование СЗ $({{\lambda }_{{l,n}}},c_{{l,n}}^{2})$ для любой пары индексов $(l,n)$ следует из общей многопараметрической теории [12], поскольку

$det\left( {\begin{array}{*{20}{c}} 1&{(1 - {{\eta }^{2}})} \\ { - 1}&{({{\xi }^{2}} - 1)} \end{array}} \right) = {{\xi }^{2}} - {{\eta }^{2}} > 0\quad \forall (\eta ,\xi ) \in [0,1) \times (1,{{\xi }_{\mathcal{S}}}],$

см. также [2].

Тем самым задача (1.15) разрешима. Применительно к ней метод [1] состоит из следующих этапов.

Прежде всего вычисляется СЗ упрощенной задачи, в которой множители при спектральном параметре ${{c}^{2}}$ постоянны и равны $(1 - \eta _{f}^{2})$, $(\xi _{f}^{2} - 1)$, где ${{\eta }_{f}}$, ${{\xi }_{f}}$ – какие-либо точки интервалов $(0,1)$ и $(1,{{\xi }_{\mathcal{S}}})$ соответственно. Граничные условия (1.6)–(1.9) при этом остаются неизменными.

Очевидно, двухпараметрическая задача с замороженными коэффициентами при ${{c}^{2}}$ распадается на две стандартные (сингулярные) однопараметрические ЗШЛ относительно спектральных параметров

$\mu = \lambda + {{c}^{2}}(1 - \eta _{f}^{2}),\quad \chi = {{c}^{2}}(\xi _{f}^{2} - 1) - \lambda .$

Решение этих задач, пара $\mu = {{\mu }_{l}}$ и $\chi = {{\chi }_{n}}$, дает нам отправную точку $(\lambda _{{l,n}}^{{(0)}},с_{{l,n}}^{{2(0)}})$ для движения в сторону искомого СЗ $({{\lambda }_{{l,n}}},c_{{l,n}}^{2})$,

(2.1)

$\left( {\begin{array}{*{20}{c}} {\lambda _{{l,n}}^{{(0)}}} \\ {\mathop {{{c}^{2}}}\nolimits_{l,n}^{(0)} } \end{array}} \right) = \frac{1}{{\xi _{f}^{2} - \eta _{f}^{2}}}\left( {\begin{array}{*{20}{c}} {\xi _{f}^{2} - 1}&{\eta _{f}^{2} - 1} \\ 1&1 \end{array}} \right)\left( {\begin{array}{*{20}{c}} {{{\mu }_{l}}} \\ {{{\chi }_{n}}} \end{array}} \right).$

Введем параметр $\tau \in [0,1]$ и заменим функции ${{Q}_{a}}(\eta ,\lambda ,{{c}^{2}})$ и ${{Q}_{r}}(\xi ,\lambda ,{{c}^{2}})$ в уравнениях (1.2), (1.3) на

$\begin{gathered} {{Q}_{{a\tau }}}(\eta ,\lambda ,{{c}^{2}}) = \lambda + {{c}^{2}}[(1 - \eta _{f}^{2})(1 - \tau ) + (1 - {{\eta }^{2}})\tau ] - \frac{{{{m}^{2}}}}{{1 - {{\eta }^{2}}}}, \\ {{Q}_{{r\tau }}}(\xi ,\lambda ,{{c}^{2}}) = {{c}^{2}}[(\xi _{f}^{2} - 1)(1 - \tau ) + ({{\xi }^{2}} - 1)\tau ] - \lambda - \frac{{{{m}^{2}}}}{{{{\xi }^{2}} - 1}}. \\ \end{gathered} $

Нетрудно показать, что в окрестности особых точек $\eta = 1$ и $\xi = 1$ асимптотика решений полученных уравнений

(2.2)

$\frac{d}{{d\eta }}(1 - {{\eta }^{2}})\frac{d}{{d\eta }}S + {{Q}_{{a\tau }}}(\eta ,\lambda ,{{c}^{2}})S = 0,\quad - 1 < \eta < 1,$

(2.3)

$\frac{d}{{d\xi }}({{\xi }^{2}} - 1)\frac{d}{{d\xi }}R + {{Q}_{{r\tau }}}(\xi ,\lambda ,{{c}^{2}})R = 0,\quad 1 < \xi < {{\xi }_{\mathcal{S}}},$

остается прежней, а следовательно, сохраняются и граничные условия (1.7), (1.8), выражающие ограниченность решений уравнений (2.2), (2.3). Дополним постановку задачи граничными условиями (1.6), (1.9).

По мере движения $\tau $ от $\tau = 0$ к $\tau = 1$, задача с “замороженными” коэффициентами трансформируется в исходную задачу (1.2)–(1.9). Как и исходная, каждая промежуточная задача (2.2), (2.3), (1.6)–(1.9) разрешима [1], [2]. Kаковы бы ни были значения $l,n = 0,1,2, \ldots $ и $\tau \in [0,1]$, можно найти такие ${{\lambda }_{{l,n}}}(\tau )$, $с_{{l,n}}^{2}(\tau )$, что ограниченные в особых точках $\eta = 1$, $\xi = 1$ решения уравнений (2.2), (2.3) удовлетворяют граничным условиям (1.6), (1.9) и обращаются в ноль $l$ раз на интервале $( - 1,1)$ и $n$ раз на интервале $(1,{{\xi }_{\mathcal{S}}})$. При этом каждая из промежуточных задач имеет эквивалентную алгебраическую формулировку вида (1.15), где значения углов Прюфера ${{\theta }_{{a\tau }}}(\eta )$ и ${{\theta }_{{r\tau }}}(\xi )$ вычисляются по уравнениям (1.12), (1.13) с заменой функций ${{Q}_{a}}$ и ${{Q}_{r}}$ на ${{Q}_{{a\tau }}}$ и ${{Q}_{{r\tau }}}$:

(2.4)

$\begin{gathered} \theta _{{a\tau }}^{'} = \frac{{si{{n}^{2}}({{\theta }_{{a\tau }}})}}{{(1 - {{\eta }^{2}})}} + {{Q}_{{a\tau }}}(\eta ,\lambda ,{{c}^{2}})co{{s}^{2}}({{\theta }_{{a\tau }}}),\quad \eta \in (0,1), \\ \mathop {lim}\limits_{\eta \to {{1}_{{ - 0}}}} {{\theta }_{{a\tau }}}(\eta ) = arctan(m), \\ \end{gathered} $

(2.5)

$\begin{gathered} \theta _{{r\tau }}^{'} = \frac{{si{{n}^{2}}({{\theta }_{{r\tau }}})}}{{({{\xi }^{2}} - 1)}} + {{Q}_{{r\tau }}}(\xi ,\lambda ,{{c}^{2}})co{{s}^{2}}({{\theta }_{{r\tau }}}),\quad \xi \in (1,{{\xi }_{\mathcal{S}}}), \\ \mathop {lim}\limits_{\xi \to {{1}_{{ + 0}}}} {{\theta }_{{r\tau }}}(\xi ) = - arctan(m). \\ \end{gathered} $

СЗ промежуточных задач (2.2), (2.3), (1.6)–(1.9) образуют кривую $(\lambda ,{{c}^{2}}) = ({{\lambda }_{{l,n}}}(\tau ),c_{{l,n}}^{2}(\tau ))$, соединяющую точку $(\lambda _{{l,n}}^{{(0)}},c_{{l,n}}^{{2(0)}}) = ({{\lambda }_{{l,n}}}(0),c_{{l,n}}^{2}(0))$ и искомое СЗ исходной задачи $({{\lambda }_{{l,n}}},c_{{l,n}}^{2}) = ({{\lambda }_{{l,n}}}(1),c_{{l,n}}^{2}(1))$.

Положим $0 = {{\tau }_{0}} \leqslant {{\tau }_{1}} \leqslant \ldots \leqslant {{\tau }_{P}} = 1$. При движении по параметру $\tau $ от ${{\tau }_{0}}$ к ${{\tau }_{P}}$ шаг $p = 1,2, \ldots P$ состоит в том, что выполняется итерация метода Ньютона для расчета $({{\lambda }_{{l,n}}}(p),c_{{l,n}}^{2}(p))$ – приближения к СЗ $({{\lambda }_{{l,n}}}({{\tau }_{p}}),c_{{l,n}}^{2}({{\tau }_{p}}))$ – с использованием такового, полученного на шаге $p - 1$, в качестве предыдущей итерации:

(2.6)

$\left( {\begin{array}{*{20}{c}} {\lambda _{{l,n}}^{{(p)}}} \\ {\mathop {{{c}^{2}}}\nolimits_{l,n}^{(p)} } \end{array}} \right) = \left( {\begin{array}{*{20}{c}} {\lambda _{{l,n}}^{{(p - 1)}}} \\ {\mathop {{{c}^{2}}}\nolimits_{l,n}^{(p - 1)} } \end{array}} \right) - \Psi _{p}^{{ - 1}}\left( {\begin{array}{*{20}{c}} {{{\varphi }_{{ap}}}(\lambda _{{l,n}}^{{(p - 1)}},\mathop {c_{{l,n}}^{{2(p - 1)}}}\nolimits_{}^{} )} \\ {{{\varphi }_{{rp}}}(\lambda _{{l,n}}^{{(p - 1)}},\mathop {c_{{l,n}}^{{2(p - 1)}}}\nolimits_{}^{} )} \end{array}} \right).$

Здесь

$\begin{gathered} {{\varphi }_{{ap}}}(\lambda ,{{c}^{2}}) = - {{\theta }_{{a{{\tau }_{p}}}}}(0) - \frac{{l\pi }}{2}, \\ {{\varphi }_{{rp}}}(\lambda ,{{c}^{2}}) = {{\theta }_{{r{{\tau }_{p}}}}}({{\xi }_{\mathcal{S}}}) - \left\{ \begin{gathered} n\pi \quad ({\text{условие}}\;{\text{Неймана}}), \hfill \\ \left( {n + \frac{1}{2}} \right)\pi \quad ({\text{условие}}\;{\text{Дирихле}}), \hfill \\ \end{gathered} \right. \\ \end{gathered} $

а элементы матрицы Якоби ${{\Psi }_{p}}$,

${{\Psi }_{p}} = \left( {\begin{array}{*{20}{c}} {{{\Psi }_{{p11}}}}&{{{\Psi }_{{p12}}}} \\ {{{\Psi }_{{p21}}}}&{{{\Psi }_{{p22}}}} \end{array}} \right) = \mathop {\left( {\begin{array}{*{20}{c}} {\frac{{\partial {{\varphi }_{{ap}}}}}{{\partial \lambda }}}&{\frac{{\partial {{\varphi }_{{ap}}}}}{{\partial {{c}^{2}}}}} \\ {\frac{{\partial {{\varphi }_{{rp}}}}}{{\partial \lambda }}}&{\frac{{\partial {{\varphi }_{{rp}}}}}{{\partial {{c}^{2}}}}} \end{array}} \right)}\nolimits_{\begin{array}{*{20}{c}} {\lambda = \lambda _{{l,n}}^{{(p - 1)}},} \\ {{{c}^{2}} = с_{{l,n}}^{{2(p - 1)}}} \end{array}} ,$

вычисляются интегрированием вспомогательных уравнений относительно функций ${{\kappa }_{{i,j}}}$, $i,j = 1,2$:

$\begin{gathered} \frac{{\partial {{\kappa }_{{11}}}}}{{\partial \eta }} = 2sin{{\theta }_{{a{{\tau }_{p}}}}}cos{{\theta }_{{a{{\tau }_{p}}}}}\left[ {\frac{1}{{1 - {{\eta }^{2}}}} - {{Q}_{{a{{\tau }_{p}}}}}(\eta ,\lambda ,{{c}^{2}})} \right]{{\kappa }_{{11}}} + co{{s}^{2}}{{\theta }_{{a{{\tau }_{p}}}}}, \\ \frac{{\partial {{\kappa }_{{12}}}}}{{\partial \eta }} = 2sin{{\theta }_{{a{{\tau }_{p}}}}}cos{{\theta }_{{a{{\tau }_{p}}}}}\left[ {\frac{1}{{1 - {{\eta }^{2}}}} - {{Q}_{{a{{\tau }_{p}}}}}(\eta ,\lambda ,{{c}^{2}})} \right]{{\kappa }_{{12}}} + \\ \end{gathered} $

(2.7)

$\begin{gathered} \, + [(1 - \eta _{f}^{2})(1 - {{\tau }_{p}}) + (1 - {{\eta }^{2}}){{\tau }_{p}}]co{{s}^{2}}{{\theta }_{{a{{\tau }_{p}}}}}, \\ \frac{{\partial {{\kappa }_{{21}}}}}{{\partial \xi }} = 2sin{{\theta }_{{r{{\tau }_{p}}}}}cos{{\theta }_{{r{{\tau }_{p}}}}}\left[ {\frac{1}{{{{\xi }^{2}} - 1}} - {{Q}_{{r{{\tau }_{p}}}}}(\xi ,\lambda ,{{c}^{2}})} \right]{{\kappa }_{{21}}} - co{{s}^{2}}{{\theta }_{{r{{\tau }_{p}}}}}, \\ \end{gathered} $

$\begin{gathered} \frac{{\partial {{\kappa }_{{22}}}}}{{\partial \xi }} = 2sin{{\theta }_{{r{{\tau }_{p}}}}}cos{{\theta }_{{r{{\tau }_{p}}}}}\left[ {\frac{1}{{{{\xi }^{2}} - 1}} - {{Q}_{{r{{\tau }_{p}}}}}(\xi ,\lambda ,{{c}^{2}})} \right]{{\kappa }_{{22}}} + \\ \, + [(\xi _{f}^{2} - 1)(1 - \tau ) + ({{\xi }^{2}} - 1){{\tau }_{p}}]co{{s}^{2}}{{\theta }_{{r{{\tau }_{p}}}}} \\ \end{gathered} $

при условии, что

(2.8)

$\mathop {lim}\limits_{\eta \to {{1}_{{ - 0}}}} {{\kappa }_{{11}}}(\eta ) = \mathop {lim}\limits_{\eta \to {{1}_{{ - 0}}}} {{\kappa }_{{12}}}(\eta ) = \mathop {lim}\limits_{\xi \to {{1}_{{ + 0}}}} {{\kappa }_{{21}}}(\xi ) = \mathop {lim}\limits_{\xi \to {{1}_{{ + 0}}}} {{\kappa }_{{22}}}(\xi ) = 0.$

Окончательно,

${{\Psi }_{{11}}} = {{\kappa }_{{11}}}(0),\quad {{\Psi }_{{12}}} = {{\kappa }_{{12}}}(0),\quad {{\Psi }_{{21}}} = {{\kappa }_{{21}}}({{\xi }_{\mathcal{S}}}),\quad {{\Psi }_{{22}}} = {{\kappa }_{{22}}}({{\xi }_{\mathcal{S}}}).$

Обратимость матрицы Якоби на каждом шаге $p = 1,2, \ldots ,P$ доказана в [1]. Полученные по выполнении $P$ шагов значения $\lambda _{{l,n}}^{{(P)}}$, $c_{{l,n}}^{{2(P)}}$ дают искомое приближение к СЗ $({{\lambda }_{{l,n}}},с_{{l,n}}^{2})$.

Для регулярных многопараметрических ЗШЛ в [1] получена оценка погрешности вычисления СЗ. Обозначим $h = max({{\tau }_{p}} - {{\tau }_{{p - 1}}})$, $p = 1,2, \ldots ,P$, тогда

$\left| {{{\lambda }_{{l,n}}} - \lambda _{{l,n}}^{{(P)}}} \right| + \left| {с_{{l,n}}^{2} - с_{{l,n}}^{{2(P)}}} \right| < C{{h}^{2}},$

где константа $C$ не зависит от разбиения $\mathop {\{ {{\tau }_{p}}\} }\nolimits_{p = 0,1, \ldots ,P} $.

3. ПЕРЕНОС ГРАНИЧНЫХ УСЛОВИЙ ИЗ ОСОБОЙ ТОЧКИ

Численная реализация алгоритма, описанного в предыдущем разделе, затрудняется тем, что в особых точках $\eta = 1$, $\xi = 1$ правые части уравнений (2.4), (2.5) и (2.7) обращаются в бесконечность. Следуя [6]–[8], воспользуемся для отхода от особых точек тем обстоятельством, что обе логарифмические производные

${{\beta }_{{a\tau }}}(\eta ) = (1 - {{\eta }^{2}})S{\kern 1pt} '(\eta ){\text{/}}S(\eta ),\quad {{\beta }_{{r\tau }}}(\xi ) = ({{\xi }^{2}} - 1)R{\kern 1pt} '(\xi ){\text{/}}R(\xi ),$

отвечающие ограниченным решениям уравнений (2.2), (2.3), разлагаются вблизи этих точек в ряды

(3.1)

${{\beta }_{{a\tau }}}(\eta ) = \sum\limits_{k = 0}^\infty \,{{\beta }_{{a\tau k}}}{{(1 - \eta )}^{k}},$

(3.2)

${{\beta }_{{r\tau }}}(\xi ) = \sum\limits_{k = 0}^\infty \,{{\beta }_{{r\tau k}}}{{(\xi - 1)}^{k}}.$

Коэффициенты ${{\beta }_{{a\tau k}}}$, ${{\beta }_{{r\tau k}}}$ получаются подстановкой (3.1), (3.2) в уравнения, полученные из (1.10), (1.11) заменой функций ${{Q}_{a}}$ и ${{Q}_{r}}$ на ${{Q}_{{a\tau }}}$ и ${{Q}_{{r\tau }}}$ соответственно.

$\begin{gathered} {{\beta }_{{a\tau 0}}} = - m, \\ {{\beta }_{{a\tau 1}}} = \frac{{{{\rho }_{{a\tau }}}}}{{1 + m}}, \\ {{\beta }_{{a\tau 2}}} = \frac{{{{\beta }_{{a\tau 1}}} - {{\rho }_{{a\tau }}} + \beta _{{a\tau 1}}^{2} + 4{{c}^{2}}{{\tau }_{p}}}}{{2(2 + m)}}, \\ \end{gathered} $

$\begin{gathered} {{\beta }_{{a\tau 3}}} = \frac{{{{\beta }_{{a\tau 2}}} - 2{{c}^{2}}{{\tau }_{p}} + {{\beta }_{{a\tau 1}}}{{\beta }_{{a\tau 2}}}}}{{(3 + m)}}, \\ {{\beta }_{{a\tau 4}}} = \frac{{3{{\beta }_{{a\tau 3}}} + {{c}^{2}}{{\tau }_{p}} + 2{{\beta }_{{a\tau 1}}}{{\beta }_{{a\tau 3}}} + \beta _{{a\tau 2}}^{2}}}{{2(4 + m)}}, \\ {{\beta }_{{a\tau k}}} = \frac{{(k - 1){{\beta }_{{a\tau k - 1}}} + \sum\limits_{s = 1}^{k - 1} {{{\beta }_{{a\tau s}}}{{\beta }_{{a\tau k - s}}}} }}{{2(k + m)}},\quad k = 5,6, \ldots \,{\kern 1pt} , \\ \end{gathered} $

$\begin{gathered} {{\beta }_{{r\tau 0}}} = m, \\ {{\beta }_{{r\tau 1}}} = - \frac{{{{\rho }_{{r\tau }}}}}{{1 + m}}, \\ {{\beta }_{{r\tau 2}}} = - \frac{{{{\beta }_{{r\tau 1}}} + {{\rho }_{{r\tau }}} + \beta _{{r\tau 1}}^{2} + 4{{c}^{2}}{{\tau }_{p}}}}{{2(2 + m)}}, \\ \end{gathered} $

$\begin{gathered} {{\beta }_{{r\tau 3}}} = - \frac{{{{\beta }_{{r\tau 2}}} + 2{{c}^{2}}{{\tau }_{p}} + {{\beta }_{{r\tau 1}}}{{\beta }_{{r\tau 2}}}}}{{(3 + m)}}, \\ {{\beta }_{{r\tau 4}}} = - \frac{{3{{\beta }_{{r\tau 3}}} + {{c}^{2}}{{\tau }_{p}} + 2{{\beta }_{{r\tau 1}}}{{\beta }_{{r\tau 3}}} + \beta _{{r\tau 2}}^{2}}}{{2(4 + m)}}, \\ {{\beta }_{{r\tau k}}} = - \frac{{(k - 1){{\beta }_{{r\tau k - 1}}} + \sum\limits_{s = 1}^{k - 1} {{{\beta }_{{r\tau s}}}{{\beta }_{{r\tau k - s}}}} }}{{2(k + m)}},\quad k = 5,6, \ldots {\kern 1pt} \,, \\ \end{gathered} $

где ${{\rho }_{{a\tau }}} = \lambda + {{c}^{2}}(1 - \tau )(1 - \eta _{f}^{2})$ и ${{\rho }_{{r\tau }}} = {{c}^{2}}(1 - \tau )(\xi _{f}^{2} - 1) - \lambda $.

Как показано в [6], [8], радиус сходимости рядов (3.1), (3.2) определяется коэффициентами ${{\beta }_{{a\tau k}}}$, ${{\beta }_{{r\tau k}}}$, $k = 0,1, \ldots ,5$. Положим

${{\mathcal{B}}_{{a\tau }}} = \mathop {max}\limits_{k = 1,2,3,4,5} \left\{ {\sqrt[k]{{\left| {{{\beta }_{{a\tau k}}}} \right|}}} \right\},\quad {{\mathcal{B}}_{{r\tau }}} = \mathop {max}\limits_{k = 1,2,3,4,5} \left\{ {\sqrt[k]{{\left| {{{\beta }_{{r\tau k}}}} \right|}}} \right\}.$

При условии, что

(3.3)

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

Почленное дифференцирование рядов (3.1), (3.2) по спектральным параметрам $\lambda $, ${{c}^{2}}$ позволяет представить функции ${{\kappa }_{{ij}}}$, $i,j = 1,2$, в виде

(3.4)

$\begin{gathered} {{\kappa }_{{1j}}}(\eta ) = - \sum\limits_{k = 1}^\infty \,\gamma _{{1j}}^{{(k)}}{{(1 - \eta )}^{k}}{\text{/}}(1 + \beta _{{a\tau }}^{2}(\eta )), \\ {{\kappa }_{{2j}}}(\xi ) = - \sum\limits_{k = 1}^\infty \,\gamma _{{2j}}^{{(k)}}{{(\xi - 1)}^{k}}{\text{/}}(1 + \beta _{{r\tau }}^{2}(\xi )). \\ \end{gathered} $

Здесь первые четыре коэффициента вычисляются по формулам

$\begin{gathered} \gamma _{{11}}^{{(1)}} = \frac{1}{{1 + m}}, \\ \gamma _{{11}}^{{(2)}} = \frac{{\gamma _{{11}}^{{(1)}} - 1 + 2{{\beta }_{{a\tau 1}}}\gamma _{{11}}^{{(1)}}}}{{2(2 + m)}}, \\ \end{gathered} $

$\begin{gathered} \gamma _{{11}}^{{(3)}} = \frac{{\gamma _{{11}}^{{(2)}} + {{\beta }_{{a\tau 1}}}\gamma _{{11}}^{{(2)}} + {{\beta }_{{a\tau 2}}}\gamma _{{11}}^{{(1)}}}}{{(3 + m)}}, \\ \gamma _{{11}}^{{(4)}} = \frac{{3\gamma _{{11}}^{{(3)}} + 2{{\beta }_{{a\tau 1}}}\gamma _{{11}}^{{(3)}} + 2{{\beta }_{{a\tau 2}}}\gamma _{{11}}^{{(2)}} + 2{{\beta }_{{a\tau 3}}}\gamma _{{11}}^{{(1)}}}}{{2(4 + m)}}, \\ \end{gathered} $

$\begin{gathered} \gamma _{{12}}^{{(1)}} = \frac{{(1 - \tau )(1 - \eta _{f}^{2})}}{{1 + m}}, \\ \gamma _{{12}}^{{(2)}} = \frac{{\gamma _{{12}}^{{(1)}} - (1 - \tau )(1 - \eta _{f}^{2}) + 2{{\beta }_{{a\tau 1}}}\gamma _{{12}}^{{(1)}} + 4\tau }}{{2(2 + m)}}, \\ \end{gathered} $

$\begin{gathered} \gamma _{{12}}^{{(3)}} = \frac{{\gamma _{{12}}^{{(2)}} - 2\tau + {{\beta }_{{a\tau 1}}}\gamma _{{12}}^{{(2)}} + {{\beta }_{{a\tau 2}}}\gamma _{{12}}^{{(1)}}}}{{(3 + m)}}, \\ \gamma _{{12}}^{{(4)}} = \frac{{3\gamma _{{12}}^{{(3)}} + \tau + 2{{\beta }_{{a\tau 1}}}\gamma _{{12}}^{{(3)}} + 2{{\beta }_{{a\tau 2}}}\gamma _{{12}}^{{(2)}} + 2{{\beta }_{{a\tau 3}}}\gamma _{{12}}^{{(1)}}}}{{2(4 + m)}}, \\ \end{gathered} $

$\begin{gathered} \gamma _{{21}}^{{(1)}} = \frac{1}{{1 + m}}, \\ \gamma _{{21}}^{{(2)}} = - \frac{{\gamma _{{21}}^{{(1)}} - 1 + 2{{\beta }_{{r\tau 1}}}\gamma _{{21}}^{{(1)}}}}{{2(2 + m)}}, \\ \end{gathered} $

$\begin{gathered} \gamma _{{21}}^{{(3)}} = - \frac{{\gamma _{{21}}^{{(2)}} + {{\beta }_{{r\tau 1}}}\gamma _{{21}}^{{(2)}} + {{\beta }_{{r\tau 2}}}\gamma _{{21}}^{{(1)}}}}{{(3 + m)}}, \\ \gamma _{{21}}^{{(4)}} = - \frac{{3\gamma _{{21}}^{{(3)}} + 2{{\beta }_{{r\tau 1}}}\gamma _{{21}}^{{(3)}} + 2{{\beta }_{{r\tau 2}}}\gamma _{{21}}^{{(2)}} + 2{{\beta }_{{r\tau 3}}}\gamma _{{21}}^{{(1)}}}}{{2(4 + m)}}, \\ \end{gathered} $

$\begin{gathered} \gamma _{{22}}^{{(1)}} = - \frac{{(1 - \tau )(\xi _{f}^{2} - 1)}}{{1 + m}}, \\ \gamma _{{22}}^{{(2)}} = - \frac{{\gamma _{{21}}^{{(1)}} + (1 - \tau )(\xi _{f}^{2} - 1) + 2{{\beta }_{{r\tau 1}}}\gamma _{{21}}^{{(1)}} + 4\tau }}{{2(2 + m)}}, \\ \end{gathered} $

$\begin{gathered} \gamma _{{22}}^{{(3)}} = - \frac{{\gamma _{{22}}^{{(2)}} + 2\tau + {{\beta }_{{r\tau 1}}}\gamma _{{22}}^{{(2)}} + {{\beta }_{{r\tau 2}}}\gamma _{{22}}^{{(1)}}}}{{(3 + m)}}, \\ \gamma _{{22}}^{{(4)}} = - \frac{{3\gamma _{{22}}^{{(3)}} + \tau + 2{{\beta }_{{r\tau 3}}}\gamma _{{22}}^{{(1)}} + 2{{\beta }_{{r\tau 2}}}\gamma _{{22}}^{{(2)}} + 2{{\beta }_{{r\tau 1}}}\gamma _{{22}}^{{(3)}}}}{{2(4 + m)}}. \\ \end{gathered} $

Начиная с $k = 5$,

$\begin{gathered} \gamma _{{1j}}^{{(k)}} = \frac{{(k - 1)\gamma _{{1j}}^{{(k - 1)}} + 2\sum\limits_{s = 1}^{k - 1} {{{\beta }_{{a\tau s}}}\gamma _{{1j}}^{{(k - s)}}} }}{{2(k + m)}}, \\ \gamma _{{2j}}^{{(k)}} = - \frac{{(k - 1)\gamma _{{1j}}^{{(k - 1)}} + 2\sum\limits_{s = 1}^{k - 1} {{{\beta }_{{a\tau s}}}\gamma _{{2j}}^{{(k - s)}}} }}{{2(k + m)}},\quad k = 5,6, \ldots ,j = 1,2. \\ \end{gathered} $

Радиус сходимости рядов (3.4) определяется сходным образом. Обозначим

${{\Gamma }_{{ij}}} = \mathop {max}\limits_{k = 1,2,3,4,5} \left\{ {\sqrt[k]{{\left| {{{\gamma }_{{ijk}}}} \right|}}} \right\},\quad i,j = 1,2.$

Ряды (3.4) сходятся абсолютно и равномерно для всех $\eta $, $\xi $ таких, что в дополнение к условиям (3.3) выполняются неравенства

(3.5)

${{\Gamma }_{{1j}}}(1 - \eta ) < 1\quad {\text{и}}\quad {{\Gamma }_{{2j}}}(\xi - 1) < 1,\quad j = 1,2.$

Внутри интервалов сходимости рядов (3.1), (3.2) и (3.4) выберем точки ${{\eta }_{b}}$, ${{\xi }_{b}}$, в которые ряды (3.1), (3.2) и (3.4) переносят начальные условия (1.7), (1.8), (2.8) из особых точек $\eta = 1$, $\xi = 1$.

Как нетрудно показать [8], погрешность усечения любого из вышеперечисленных рядов удовлетворяет неравенству

(3.6)

где $N$ – число удержанных членов ряда, а $\mathcal{C}$ принимает значения

, или ${{\Gamma }_{{1j}}}(1 - {{\eta }_{b}})$, ${{\Gamma }_{{2j}}}({{\xi }_{b}} - 1)$, $j = 1,2$.

Таким образом, выбором точек отхода ${{\eta }_{b}}$, ${{\xi }_{b}}$ и числа $N$ определяется точность начальных условий

$\begin{gathered} {{\theta }_{{a\tau }}}({{\eta }_{b}}) = arctan\left[ {\sum\limits_{k = 0}^N \,{{\beta }_{{a\tau k}}}{{{(1 - {{\eta }_{b}})}}^{k}}} \right], \\ {{\theta }_{{r\tau }}}({{\xi }_{b}}) = arctan\left[ {\sum\limits_{k = 0}^N \,{{\beta }_{{r\tau k}}}{{{({{\xi }_{b}} - 1)}}^{k}}} \right], \\ \end{gathered} $

$\begin{gathered} {{\kappa }_{{1j}}}({{\eta }_{b}}) = - \sum\limits_{k = 1}^N \,\gamma _{{1j}}^{{(k)}}{{(1 - {{\eta }_{b}})}^{k}}\mathop {\left[ {1 + \mathop {\left( {\sum\limits_{k = 0}^N \,{{\beta }_{{a\tau k}}}{{{(1 - {{\eta }_{b}})}}^{k}}} \right)}\nolimits^2 } \right]}\nolimits^{ - 1} , \\ {{\kappa }_{{2j}}}({{\xi }_{b}}) = - \sum\limits_{k = 1}^N \,\gamma _{{2j}}^{{(k)}}{{({{\xi }_{b}} - 1)}^{k}}\mathop {\left[ {1 + \mathop {\left( {\sum\limits_{k = 0}^N \,{{\beta }_{{r\tau k}}}{{{({{\xi }_{b}} - 1)}}^{k}}} \right)}\nolimits^2 } \right]}\nolimits^{ - 1} ,\quad j = 1,2, \\ \end{gathered} $

для последующего интегрирования уравнений (2.4), (2.5) и (2.7).

4. СГЛАЖИВАЮЩИЕ ФУНКЦИИ

Перенос граничных условий из особой точки в близкую регулярную не полностью устраняет трудности, связанные с численным решением исходной задачи. Радиус сходимости рядов (3.1), (3.2), (3.4) уменьшается с ростом параметров $m$, $\lambda $ и c², в то время как правые части уравнений (2.4), (2.5), (2.7) растут как вблизи особых точек $\eta = 1$, $\xi = 1$, так и внутри интервалов $[0,1)$, $(1,{{\xi }_{\mathcal{S}}}]$. Все это существенно осложняет численное интегрирование уравнений (2.4), (2.5), (2.7) и приводит к резкому ступенчатому изменению их решений. Следуя рекомендациям, предложенным в [6], [7], введением специальной “сглаживающей” (масштабирующей) функции удается добиться менее резкого поведения функций ${{\theta }_{{a\tau }}}$, ${{\theta }_{{r\tau }}}$ и существенно расширить область применения описанных там численных методов (см. также [8], [13]).

Введем сглаживающие функции ${{\nu }_{{a\tau }}}(\eta ) = {{\nu }_{{a\tau }}}(\eta ,\lambda ,{{c}^{2}})$ и ${{\nu }_{{r\tau }}}(\xi ) = {{\nu }_{{r\tau }}}(\xi ,\lambda ,{{c}^{2}})$. Заменой

(4.1)

${{\beta }_{{a\tau }}}(\eta ) = - {{\nu }_{{a\tau }}}(\eta )tan{{\zeta }_{{a\tau }}}(\eta ),\quad {\text{и}}\quad {{\beta }_{{r\tau }}}(\xi ) = - {{\nu }_{{r\tau }}}(\xi )tan{{\zeta }_{{r\tau }}}(\xi )$

определим видоизмененные углы Прюфера.

Как нетрудно проверить, функции ${{\zeta }_{{a\tau }}}$, ${{\zeta }_{{r\tau }}}$ являются решениями уравнений

(4.2)

$\zeta _{{a\tau }}^{'} = - \frac{{\nu _{{a\tau }}^{'}(\eta )}}{{{{\nu }_{{a\tau }}}(\eta )}}sin({{\zeta }_{{a\tau }}})cos({{\zeta }_{{a\tau }}}) + \frac{{{{\nu }_{{a\tau }}}(\eta )si{{n}^{2}}({{\zeta }_{{a\tau }}})}}{{(1 - {{\eta }^{2}})}} + \frac{{{{\mathcal{Q}}_{{{\text{a}}\tau }}}(\eta )co{{s}^{2}}({{\zeta }_{{{\text{a}}\tau }}})}}{{{{\nu }_{{a\tau }}}(\eta )}},\quad \eta \in [0,1),$

(4.3)

$\zeta _{{r\tau }}^{'} = - \frac{{\nu _{{r\tau }}^{'}(\xi )}}{{{{\nu }_{{r\tau }}}(\xi )}}sin({{\zeta }_{{r\tau }}})cos({{\zeta }_{{r\tau }}}) + \frac{{{{\nu }_{{r\tau }}}(\xi )si{{n}^{2}}({{\zeta }_{{r\tau }}})}}{{({{\xi }^{2}} - 1)}} + \frac{{{{\mathcal{Q}}_{{{\text{r}}\tau }}}(\xi )co{{s}^{2}}({{\zeta }_{{{\text{r}}\tau }}})}}{{{{\nu }_{{r\tau }}}(\xi )}},\quad \xi \in (1,{{\xi }_{\mathcal{S}}}],$

причем с точностью, гарантированной неравенством (3.6), в точках ${{\eta }_{b}}$, ${{\xi }_{b}}$ выполняются условия

(4.4)

${{\zeta }_{{a\tau }}}({{\eta }_{b}}) = - arctan\left[ {\frac{{\sum\limits_{k = 0}^N {{{\beta }_{{a\tau k}}}{{{(1 - {{\eta }_{b}})}}^{k}}} }}{{{{\nu }_{{a\tau }}}({{\eta }_{b}})}}} \right],\quad {{\zeta }_{{r\tau }}}({{\xi }_{b}}) = - arctan\left[ {\frac{{\sum\limits_{k = 0}^N {{{\beta }_{{r\tau k}}}{{{({{\xi }_{b}} - 1)}}^{k}}} }}{{{{\nu }_{{r\tau }}}({{\xi }_{b}})}}} \right].$

Очевидно, формулировка условий на противоположных концах интервалов $[0,1)$, $(1,{{\xi }_{\mathcal{S}}}]$, в точках $\eta = 0$ и $\xi = {{\xi }_{\mathcal{S}}}$, остается прежней:

$\begin{gathered} 0 = {{\varphi }_{{a\tau }}}(\lambda ,{{c}^{2}}) = - {{\zeta }_{{a\tau }}}(0) - \frac{{l\pi }}{2}, \\ 0 = {{\varphi }_{{r\tau }}}(\lambda ,{{c}^{2}}) = {{\zeta }_{{r\tau }}}({{\xi }_{\mathcal{S}}}) - \left\{ \begin{gathered} n\pi \quad ({\text{условие}}\;{\text{Неймана}}), \hfill \\ \left( {n + \frac{1}{2}} \right)\pi \quad ({\text{условие}}\;{\text{Дирихле}}). \hfill \\ \end{gathered} \right. \\ \end{gathered} $

При условии, что сглаживающие функции положительны на всем своем интервале определения, сохраняется связь между числом внутренних нулей решений уравнений (2.2), (2.3) и числом горизонталей $\zeta = \pi k + \pi {\text{/}}2$, $k \in \mathbb{Z}$, которые пересекают модифицированные углы Прюфера на этих интервалах.

Как видно из уравнений (4.2), (4.3), при $m \ne 0$ для больших значений $\left| {{{Q}_{{a\tau }}}(\eta )} \right|$ оптимальный выбор сглаживающей функции – это ${{\nu }_{{a\tau }}}(\eta ) \sim \sqrt {\left| {{{Q}_{{a\tau }}}(\eta )} \right|(1 - {{\eta }^{2}})} $. Также для $\left| {{{Q}_{{r\tau }}}(\xi )} \right| \gg 1$ наилучшим образом подходит сглаживающая функция ${{\nu }_{{r\tau }}}(\xi ) \sim \sqrt {\left| {{{Q}_{{r\tau }}}(\xi )} \right|({{\xi }^{2}} - 1)} $.

Как и ранее, уравнения и начальные условия для вычисления элементов матрицы Якоби получаются дифференцированием соответственно уравнений (4.2), (4.3) и начальных условий (4.4) по спектральным параметрам $\lambda $, c², но ввиду громоздкости правых частей здесь не приводятся в общем виде. Вместо этого рассмотрим отдельные важные для приложений частные случаи.

4.1. Моделирование колебаний типа “шепчущей галереи”

Разделение переменных облегчает формулировку условий, выделяющих МШГ из общего числа собственных колебаний внутри сфероидального резонатора $\mathcal{S}$. А именно, МШГ характеризуются чрезвычайно большим азимутальным индексом $m$, $m \sim 100{\kern 1pt} - {\kern 1pt} 10\,000$, по сравнению со значениями $l$, $n$, $l,n \sim 0{\kern 1pt} - {\kern 1pt} 50$.

Казалось бы, гигантские значения индекса $m$ требуют введения “оптимальной” сглаживающей функции. Однако для расчета МШГ сфероидов, как сильно вытянутых, так и близких к шару, оказывается достаточно постоянного множителя ${{\nu }_{{a(r)}}} = m$.

Такая “сглаживающая” функция требует минимальных очевидных изменений в уравнениях и начальных условиях для вычисления модифицированных углов Прюфера и элементов матрицы Якоби.

Примеры расчетов частот МГШ для $m \sim 1000{\kern 1pt} - {\kern 1pt} 10\,000$, выполненные при помощи такой сглаживающей функции, представлены таблицей в [14].

На фиг. 1 изображены собственная функция сфероида ${{\xi }_{\mathcal{S}}} = 1.01$, отвечающая мультииндексу $(l,n) = (0,0)$ и $m = 100$, и ее компоненты, угловая и радиальная ВВСФ. На врезке (слева) возмущенная область показана в увеличенном масштабе. В используемой на трехмерных графиках цветовой схеме серый цвет соответствует отсутствию возбуждения в среде.

Фиг. 1.

Собственная функция сфероида ${{\xi }_{\mathcal{S}}} = 1.01$, отвечающая мультииндексу $(l,n) = (0,0)$ и $m = 100$.

Для вычисления соответствующего СЗ $({{\lambda }_{{0,0}}},c_{{0,0}}^{2}) = ( - 525688.569,\;536427.039)$ при точности интегрирования вспомогательных задач $\varepsilon = {{10}^{{ - 5}}}$ оказалось достаточно Р = 500 равноотстоящих шагов по параметру $\tau $. Здесь и далее для вычисления компонент собственной функции были использованы методы, представленные в [6]–[8].

Тот факт, что возмущение акустической среды на фиг. 1 сосредоточено на экваторе резонатора, объясняется поведением угловой и радиальной функций. Последнее, как известно (см., например, [15]), определяется положением точек поворота уравнений (1.2), (1.3), т.е. таких значений ${{\eta }_{{tp}}} \in [0,1)$, ${{\xi }_{{tp}}} \in (1,{{\xi }_{\mathcal{S}}}]$, при которых функции ${{Q}_{a}}(\eta ,{{\lambda }_{{0,0}}},c_{{0,0}}^{2})$ и ${{Q}_{r}}(\xi ,{{\lambda }_{{0,0}}},c_{{0,0}}^{2})$ меняют знак. Так, справа от ${{\eta }_{{tp}}} \approx 0.0368$, на интервале $({{\eta }_{{tp}}},1)$, ограниченное решение углового уравнения экспоненциально убывает к нулю, поскольку ${{Q}_{a}}(\eta ,{{\lambda }_{{0,0}}},c_{{0,0}}^{2}) < 0$. Сходным образом ведет себя радиальная функция слева от ${{\xi }_{{tp}}} \approx 1.0092$, где ${{Q}_{r}}(\xi ,{{\lambda }_{{0,0}}},c_{{0,0}}^{2}) < 0$. Таким образом, область, в которой возмущение существенно отлично от нуля, заключена между сфероидами $\xi = {{\xi }_{\mathcal{S}}}$, $\xi = {{\xi }_{{tp}}}$ и ограничена полостями гиперболоида $\eta = {{\eta }_{{tp}}}$.

Добавим, что при увеличении азимутального индекса $m$ точки поворота ${{\eta }_{{tp}}}$, ${{\xi }_{{tp}}}$ смещаются влево и вправо соответственно, что усиливает эффект “шепчущей галереи”, сжимая область, в которой сосредоточено возмущение, в кольцо, охватывающее экватор.

4.2. Сглаживающая функция $\nu = \sqrt {\left| \lambda \right|} $

Альтернативой такому выбору ${{\nu }_{{a(r)}}}$ служит другая постоянная функция, пригодная также и для вычисления колебательных мод, отвечающих небольшим по сравнению с компонентами мультииндекса $(l,n)$ величинам $m$, а именно ${{\nu }_{{a(r)}}}({\kern 1pt} \cdot {\kern 1pt} ) = \sqrt {\left| \lambda \right|} $, где $\lambda $ – текущее значение этого параметра. Заметим, что такой выбор ${{\nu }_{{a(r)}}}$ требует некоторых несложных изменений. Так, функции ${{\zeta }_{{a\tau }}}$, ${{\kappa }_{{1j}}}(\eta )$, $j = 1,2$, получаются интегрированием системы ОДУ

$\frac{{\partial {{\zeta }_{{a\tau }}}}}{{\partial \eta }} = \frac{{\sqrt {\left| \lambda \right|} si{{n}^{2}}({{\zeta }_{{a\tau }}})}}{{(1 - {{\eta }^{2}})}} + \frac{{{{{\text{Q}}}_{{{\text{a}}\tau }}}(\eta )co{{s}^{2}}({{\zeta }_{{{\text{a}}\tau }}})}}{{\sqrt {\left| \lambda \right|} }},$

$\begin{gathered} \frac{{\partial {{\kappa }_{{11}}}}}{{\partial \eta }} = 2sin({{\zeta }_{{a\tau }}})cos({{\zeta }_{{a\tau }}})\left[ {\frac{{\sqrt {\left| \lambda \right|} }}{{1 - {{\eta }^{2}}}} - \frac{{{{Q}_{{a\tau }}}(\eta ,\lambda ,{{c}^{2}})}}{{\sqrt {\left| \lambda \right|} }}} \right]{{\kappa }_{{11}}} + \frac{{co{{s}^{2}}({{\zeta }_{{a\tau }}})}}{{\sqrt {\left| \lambda \right|} }} + \\ \, + \left[ {\frac{{si{{n}^{2}}({{\zeta }_{{a\tau }}})}}{{(1 - {{\eta }^{2}})}} - \frac{{{{{\text{Q}}}_{{{\text{a}}\tau }}}(\eta )co{{s}^{2}}({{\zeta }_{{{\text{a}}\tau }}})}}{{\left| \lambda \right|}}} \right]\frac{{\operatorname{sign} (\lambda )}}{{2\sqrt {\left| \lambda \right|} }}, \\ \end{gathered} $

$\frac{{\partial {{\kappa }_{{12}}}}}{{\partial \eta }} = 2sin({{\zeta }_{{a\tau }}})cos({{\zeta }_{{a\tau }}})\left[ {\frac{{\sqrt {\left| \lambda \right|} }}{{1 - {{\eta }^{2}}}} - \frac{{{{Q}_{{a\tau }}}(\eta ,\lambda ,{{c}^{2}})}}{{\sqrt {\left| \lambda \right|} }}} \right]{{\kappa }_{{12}}} + [(1 - \eta _{f}^{2})(1 - \tau ) + (1 - {{\eta }^{2}})\tau ]\frac{{co{{s}^{2}}({{\zeta }_{{a\tau }}})}}{{\sqrt {\left| \lambda \right|} }},$

причем

$\begin{gathered} {{\zeta }_{{a\tau }}}({{\eta }_{b}}) = - arctan\left[ {\frac{1}{{\sqrt {\left| \lambda \right|} }}\sum\limits_{k = 0}^N \,{{\beta }_{{a\tau k}}}{{{(1 - {{\eta }_{b}})}}^{k}}} \right], \\ {{\kappa }_{{11}}}({{\eta }_{b}}) = \frac{1}{{\sqrt {\left| \lambda \right|} }}\left\{ {\frac{m}{{2\lambda }} + \sum\limits_{k = 1}^N \,\left[ {\frac{{{{\beta }_{{a\tau k}}}}}{{2\lambda }} - \gamma _{{11}}^{{(k)}}} \right]{{{(1 - {{\eta }_{b}})}}^{k}}} \right\}\mathop {\left[ {1 + \frac{1}{{\left| \lambda \right|}}\mathop {\left( {\sum\limits_{k = 0}^N \,{{\beta }_{{a\tau k}}}{{{(1 - {{\eta }_{b}})}}^{k}}} \right)}\nolimits^2 } \right]}\nolimits^{ - 1} , \\ {{\kappa }_{{12}}}({{\eta }_{b}}) = - \frac{1}{{\sqrt {\left| \lambda \right|} }}\sum\limits_{k = 1}^N \,\gamma _{{12}}^{{(k)}}{{(1 - {{\eta }_{b}})}^{k}}\mathop {\left[ {1 + \frac{1}{{\left| \lambda \right|}}\mathop {\left( {\sum\limits_{k = 0}^N \,{{\beta }_{{a\tau k}}}{{{(1 - {{\eta }_{b}})}}^{k}}} \right)}\nolimits^2 } \right]}\nolimits^{ - 1} . \\ \end{gathered} $

Введение такой сглаживающей функции осмысленно, очевидно, при условии, что $\left| \lambda \right| \gg 1$.

Частоты колебательных мод, изображенных на фиг. 2, 3 , были вычислены с использованием сглаживающей функции ${{\nu }_{{a(r)}}}({\kern 1pt} \cdot {\kern 1pt} ) = \sqrt {\left| \lambda \right|} $. В обоих случаях точность расчета вспомогательных функций с учетом погрешности усечения рядов (3.1), (3.2) и (3.4) $\varepsilon = {{10}^{{ - 5}}}$. Число шагов по параметру $\tau $, использованное для получения обоих СЗ, $P = 2000$. Структура фиг. 2, 3 такая же, как в случае фиг. 1.

Фиг. 2.

Собственная функция сфероида ${{\xi }_{\mathcal{S}}} = 2.0$, отвечающая мультииндексу $(l,n) = (10,100)$ и $m = 2$. Искомое СЗ $({{\lambda }_{{10,100}}},c_{{10,100}}^{2}) = ( - 32172.271,36110.535)$.

Фиг. 3.

Собственная функция сфероида ${{\xi }_{\mathcal{S}}} = 2.0$, отвечающая мультииндексу $(l,n) = (100,10)$ и $m = 10$. Здесь $({{\lambda }_{{100,10}}},c_{{100,10}}^{2}) = (8304.560,8058.216)$.

Точки поворота углового и радиального уравнений, отвечающих фиг. 2, ${{\eta }_{{tp}}} \approx 0.3301$, ${{\xi }_{{tp}}} \approx 1.00006$, и возмущение внутри резонатора распространяется в области внешней по отношению к сфероиду $\xi = {{\xi }_{{tp}}}$ и ограничено полостями гиперболоида $\eta = {{\eta }_{{tp}}}$.

Отметим, что с ростом индекса $n$ точка поворота ${{\xi }_{{tp}}}$ сдвигается в сторону особенности, и возбужденной оказывается вся область между полостями гиперболоида $\eta = {{\eta }_{{tp}}}$.

Увеличение индекса $l$ при фиксированных индексах $m$ и $n$ смещает точку поворота в угловом уравнении к особенности. Точка поворота радиального уравнения также смещается вправо, так что возмущенным оказывается тонкий приповерхностный слой резонатора. При этом амплитуда колебаний вблизи $\eta = 0$, что на трехмерном графике соответствует области экватора, оказывается существеннно меньше, чем вблизи особой точки, т.е. у полюсов.

Приведенная на фиг. 3 собственная функция сфероида ${{\xi }_{\mathcal{S}}} = 2.0$, а также ее угловая и радиальная компоненты отвечают мультииндексу $(l,n) = (100,10)$ и $m = 10$. Точки поворота уравнений (1.2), (1.3) в этом случае ${{\eta }_{{tp}}} \approx 0.994$, ${{\xi }_{{tp}}} = 1.4292$.

На всех представленных фигурах на поверхности резонатора выполняется условие Неймана.

ЗАКЛЮЧЕНИЕ

Модификация метода [1], [2], использующая разработанную ранее [6]–[8] процедуру переноса граничных условий в окрестности регулярной особой точки, позволила существенно расширить область изменения входных параметров, к примеру осуществить вычисление МШГ, отвечающих гигантским значениям азимутального индекса $m \sim 1000{\kern 1pt} - {\kern 1pt} 10\,000$, в том числе и для сильно вытянутых сфероидов [14].

Следует отметить, что мелкость разбиения интервала изменения параметра $\tau $, h = max(τ_p – ‒ ${{\tau }_{{p - 1}}})$, $p = 1,2, \ldots ,P$, требуемая для достижения заданной точности, существенным образом зависит от выбора точек ${{\eta }_{f}}$, ${{\xi }_{f}}$, где на первоначальном этапе “замораживаются” множители $(1 - {{\eta }^{2}})$, $({{\xi }^{2}} - 1)$. Объяснением этой зависимости служит соотношение (2.1). Как видно, сдвигая ${{\eta }_{f}}$, ${{\xi }_{f}}$ в сторону особенностей, можно добиться сколь угодно большого отклонения первоначальных значений $\lambda _{{l,n}}^{{(0)}}$, $\mathop {c_{{l,n}}^{{2(0)}}}\nolimits_{}^{} $ от искомых СЗ. Вычисления к фиг. 1–3 были выполнены для ${{\eta }_{f}} = 0.01$, ${{\xi }_{f}} = 1.009$.

Доказанная в [1], [2] разрешимость промежуточных задач при движении от “замороженной” к исходной задаче открыла более широкие возможности для применения метода Ньютона, к примеру, на базе высокоточных разностных схем. Такая работа была проведена автором совместно с П. Амодио и Дж. Сеттанни из университета г. Бари (Италия), а также А. Арнольдом и Е. Вайнмюллер из технического университета Вены (Австрия).

В [14] представлено к публикации подробное описание возможных модификаций метода [1], [2], проиллюстрированное результатами вычислений МШГ. Сходным образом идеи, заложенные в работах [1], [2], можно применять для вычисления МШГ и других собственных колебаний внутри сплюснутых сфероидов и прочих резонаторов, геометрия которых допускает разделение переменных в уравнении Гельмгольца.

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

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

Абрамов А.А., Ульянова В.И. Один метод решения самосопряженных многопараметрических спектральных задач для слабо связанных систем обыкновенных дифференциальных уравнений // Ж. вычисл. матем. и матем. физ. 1997. Т. 37. № 5. С. 566–571.
Абрамов А.А., Ульянова В.И. Один метод решения самосопряженных многопараметрических спектральных задач для систем уравнений с особенностями // Ж. вычисл. матем. и матем. физ. 1998. Т. 38. № 10. С. 1636–1640.
Matsko A.B., Savchenkov A.A., Strekalov D., Ilchenko V.S., Maleki L. Review of Applications of Whispering-Gallery Mode Resonators in Photonics and Nonlinear Optics // IPN Progress Report. 2005. V. 42. № 162. P. 1–51.
Matsko A.B. and Ilchenko V.S. Optical resonators with whispering gallery nodes – Part I: basics // IEEE J. Sel. Topics Quantum Electron. 2006. V. 12. P. 3–14.
Городецкий М.Л. Основы теории оптических микрорезонаторов. М.: Физический факультет МГУ им. М.В. Ломоносова, 2010.
Абрамов А.А., Дышко А.Л., Конюхова Н.Б., Пак Т.В., Парийский Б.С. Вычисление вытянутых сфероидальных функций решением соответствующих дифференциальных уравнений // Ж. вычисл. матем. и матем. физ. 1984. Т. 24. № 1. С. 3–18.
Абрамов А.А., Дышко А.Л., Конюхова Н.Б., Левитина Т.В. Вычисление радиальных волновых функций для сфероидов и трехосных эллипсоидов модифицированным методом фазовых функций // Ж. вычисл. матем. и матем. физ. 1991. Т. 31. № 6. С. 212–234.
Levitina T.V. and Brändas E.J. Computational techniques for Prolate Spheroidal Wave Functions in Signal Processing // J. Comp. Meth. Sci. & Engng. 2001. V. 1. P. 287–313.
Amodio P., Levitina T., Settanni G., Weinmüller E.B. On the calculation of the finite Hankel transform eigenfunctions // J. of Applied Mathematics and Computing. 2013. V. 43. № 1–2. P. 151–173.
Комаров И.В., Пономарев Л.И., Славянов С.Ю. Сфероидальные и кулоновские сфероидальные функции. М.: Наука, 1976.
Morse P.M. and Feshbach H. Methods of Theoretical Physics. McGrawHill, New York, 1953.
Volkmer H. Multiparameter Eigenvalue Problems and Expansion Theorems. Lecture Notes in Mathematics. V. 1356. Berlin: Springer-Verlag, 1988.
Конюхова Н.Б., Масалович С.Е., Староверова И.Б. О вычислении быстроосциллирующих собственных функций непрерывного спектра и несобственных интегралов от них // Ж. вычисл. матем. и матем. физ. 1995. Т. 35. № 3. С. 360–379.
Amodio P., Levitina T., Settanni G., Weinmüller E. On the Abramov approach for the numerical simulation of the whispering gallery modes in prolate spheroids // представлено к публикации в Applied Mathematics and Computations.
Федорюк М.В. Асимптотические методы для линейных обыкновенных дифференциальных уравнений. М.: Наука, 1983.

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

Инструменты

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

Журнал вычислительной математики и математической физики

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