1. На главную
  2. Электронные версии
  3. Журнал вычислительной математики и математической физики
  4. T. 60, № 8, 2020

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

Вариационный метод расчета лучевых траекторий и фронтов волн цунами, порожденных локализованным источником

С. Ю. Доброхотов 1, 3, *, М. В. Клименко 2, И. А. Носиков 2, **, А. А. Толченников 1, 3

1 ИПМех
117526 Москва, пр-т Вернадского, 101, Россия

2 КФ ИЗМИРАН
236035 Калининград, ул. Пионерская, 61, Россия

3 МФТИ
141700 М.о., Долгопрудный, Институтский пер., 9, Россия

* E-mail: dobr@ipmnet.ru
** E-mail: igor.nosikov@gmail.com

Поступила в редакцию 11.11.2019
После доработки 15.02.2020
Принята к публикации 09.04.2020

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

Аннотация

Представлен вариационный подход к решению граничной задачи о расчете лучевых траекторий и фронтов океанических волн. В основе метода лежит принцип Ферма о стационарности времени распространения волны. Отличительной особенностью предложенного подхода является прямая оптимизация функционала Ферма без необходимости решения уравнения Эйлера–Лагранжа, причем положения источника и пункта регистрации волны зафиксированы. В работе предлагается решение проблемы многолучевости граничной задачи на основе поиска различных типов стационарных точек функционала Ферма. Верификация метода проведена на основе численных расчетов методом бихарактеристик с использованием аналитических моделей океанического дна. В работе представлены преимущества использования вариационного подхода и перспективы его дальнейшего развития для решения задач расчета океанических волн. Обсуждается вопрос о взаимосвязи различных типов стационарных точек функционала луча, каустик и фокусов. Библ. 26. Фиг. 6.

Ключевые слова: океанические волны, цунами, лучи, фронты, принцип Ферма, функционал, метод бихарактеристик.

1. ВВЕДЕНИЕ. ПОСТАНОВКА ЗАДАЧИ

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

В работах [6]–[8] был разработан алгоритм быстрого аналитико-численного моделирования распространения длинных волн (например, волн цунами) и вихрей в океане. В этих работах были получены асимптотические формулы для волн и вихрей, возбуждаемых локализованными источниками и описываемых либо линеаризованными уравнениями мелкой воды с переменной глубиной, либо потенциальной моделью (учитывающей эффекты дисперсии). Эти асимптотические формулы получены с помощью модифицированного канонического оператора Маслова, приспособленного для построения решений, локализованных в окрестности точек и кривых. Полученные формулы достаточно просты и используют минимум информации, необходимой для описания качественных и количественных характеристик волн. Первый шаг этого алгоритма заключается в нахождении в каждый момент времени $t$ волнового фронта (возможно, с точками самопересечения) – кривой, получаемой проекцией на координатную плоскость концов траекторий системы Гамильтона (системы характеристик) ${\mathbf{r}} = {\mathbf{R}}(\phi ,t),$ ${\mathbf{p}} = {\mathbf{P}}(\phi ,t)$:

(1)
${\mathbf{\dot {r}}} = {{H}_{{\mathbf{p}}}}({\mathbf{r}},{\mathbf{p}}) = c({\mathbf{r}})\frac{{\mathbf{p}}}{{\left| {\mathbf{p}} \right|}},\quad {\mathbf{\dot {p}}} = - {{H}_{{\mathbf{r}}}}({\mathbf{r}},{\mathbf{p}}) = - \nabla c({\mathbf{r}})\left| {\mathbf{p}} \right|$
(с гамильтонианом $H({\mathbf{r}},{\mathbf{p}}) = \left| {\mathbf{p}} \right|c({\mathbf{r}})$, где $c({\mathbf{r}}) = \sqrt {gD({\mathbf{r}})} $, $g$ – ускорение свободного падения, $D({\mathbf{r}})$ – глубина океана в точке с радиус-вектором ${\mathbf{r}}$), выпущенных из заданной точки ${\mathbf{r}}{\kern 1pt} {{{\text{|}}}_{{t = 0}}} = {{{\mathbf{r}}}_{0}} = ({{x}_{0}},{{y}_{0}})$ (точки, где локализовано начальное возмущение ${{\eta }^{0}}(({\mathbf{r}} - {{{\mathbf{r}}}_{0}}){\text{/}}l)$) с начальными импульсами, лежащими на единичной окружности ${\mathbf{p}}{\kern 1pt} {{{\text{|}}}_{{t = 0}}} = {\mathbf{n}}(\phi ): = (cos\phi ,sin\phi )$. В каждый момент времени фронт $\gamma (t)$ определяется концами траекторий $\gamma (t) = \{ {\mathbf{r}} = {\mathbf{R}}(\phi ,t),\phi \in [0,2\pi )\} $. Пример фронтов и проекций траекторий приведен на фиг. 1.

Фиг. 1.

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

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

$\eta ({\mathbf{r}},t) \approx \sqrt {\frac{l}{{\left| {{{{\mathbf{R}}}_{\phi }}} \right|}}} {\kern 1pt} \sqrt[4]{{\frac{{D({{{\mathbf{r}}}_{0}})}}{{D({\mathbf{R}}(\phi ,t))}}}}{{\left. {{\text{Re}}\left[ {{{e}^{{ - i\pi m(\phi ,t)/2}}}F\left( {\frac{{S({\mathbf{r}},t)}}{l},\phi } \right)} \right]} \right|}_{{\phi = \phi ({\mathbf{r}},t)}}},$
где $l$ – ширина начального возмущения ${{\eta }^{0}}(({\mathbf{r}} - {{{\mathbf{r}}}_{0}}){\text{/}}l)$, функция $\phi ({\mathbf{r}},t)$ определяется из уравнения $\left\langle {{\mathbf{r}} - {\mathbf{R}}(\phi ,t),{{{\mathbf{R}}}_{\phi }}(\phi ,t)} \right\rangle = 0$, фаза $S({\mathbf{r}},t) = \left\langle {{\mathbf{P}}(\phi ({\mathbf{r}},t),t),{\mathbf{r}} - {\mathbf{R}}(\phi ({\mathbf{r}},t),t)} \right\rangle $, целое число $m(\phi ,t)$ – индекс Маслова, совпадающий с индексом Морса траектории $({\mathbf{R}}(\phi ,t),{\mathbf{P}}(\phi ,t))$ и также определяющий число фокальных точек на траектории при $(0 < \tau < t)$. “Функция профиля” $F$ связана с начальным возмущением ${{\eta }^{0}}(({\mathbf{r}} - {{{\mathbf{r}}}_{0}}){\text{/}}l)$ следующим образом:
$F(z,\phi ) = \frac{{{{e}^{{ - i\pi /4}}}}}{{\sqrt {2\pi } }}\int\limits_0^\infty {{{{\tilde {\eta }}}^{0}}(\rho {\mathbf{n}}(\phi )){{e}^{{iz\rho }}}\sqrt \rho d\rho } ,$
и ${{\tilde {\eta }}^{0}}({\mathbf{p}})$ – преобразование Фурье ${{\eta }^{0}}({\mathbf{r}})$. Пример расчета профилей волны в регулярных точках фронта приведен на фиг. 2.

Фиг. 2.

Эскизы профилей волны в разных участках фронта, разделенных точками поворота фронта.

Заметим при этом, что формулы [6]–[8] носят “локальный характер” , распространение волн идет по траекториям (характеристикам) и их структура не связана с поведением волны в точках, далеких от выделенных траекторий. Поэтому если нас интересует поведение волнового поля и временная история в выделенной точке ${{{\mathbf{r}}}_{1}}$ (например, находящейся около берега), то нужно построить только траектории, которые приходят в ту выделенную точку, а остальные траектории, как и полный волновой фронт, строить не нужно. Время расчета отдельной траектории, разумеется, существенно меньше времени расчета всего волнового фронта.

Таким образом, возникает нелинейная краевая задача нахождения траекторий, выпущенных с неизвестным углом из заданной точки ${{{\mathbf{r}}}_{0}}$ и приходящих за некоторое (также неизвестное) время в заданную точку наблюдения ${{{\mathbf{r}}}_{1}}$:

${\mathbf{R}}(\phi ,t) = {{{\mathbf{r}}}_{1}}.$
Один из способов решения такой задачи – метод стрельбы (пристрелка), который в нелинейном случае зачастую неустойчив и требует значительных вычислительных затрат [9]. Альтернативным подходом к решению граничной задачи является использование прямого вариационного метода (основанного на принципе Ферма), который и представлен в данной работе.

Ранее были предложены различные формулировки вариационной задачи, основанные прежде всего на минимизации фазового пути [10]–[12]. Также известны работы, посвященные решению вариационного уравнения [12], [13]. Однако с практической точки зрения прямой метод оптимизации [10] является наиболее эффективным. Идея метода заключается в том, что некая первоначально заданная траектория последовательно трансформируется в оптимальную, причем ее концы на протяжении всего процесса оптимизации зафиксированы в соответствии с граничными условиями. Важным достоинством такого подхода является автоматическое выполнение граничных условий в виде закрепленных начальной и конечной точек. Варианты такого метода известны в различных областях науки: химии и физике твердого тела [14]–[16], физике магнитных систем [17], [18], сейсмологии [10], [11], радиофизике [12], [19], [20] и т.д. Идея об использовании вариационных методов в расчете траекторий в задаче о волнах цунами была отмечена в [1]. Она была реализована для модельной задачи распространения волнового фронта и его отражения от берега при условии, что дно линейно (без связи с формулами для волнового поля, определяемого этим фронтом). Также вариационный подход к решению граничной задачи применялся в сравнительно недавней работе [21] о распространении волн цунами.

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

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

2. ВАРИАЦИОННЫЙ МЕТОД

Законы геометрической оптики основаны на принципе Ферма, согласно которому траектория луча обеспечивает экстремум функционалу времени прихода волны:

(2)
$T[\gamma ] = \int\limits_A^B \,\frac{1}{{c(\vec {r})}}dl.$
Здесь интегрирование производится вдоль кривой $\gamma $, задающей траекторию луча, которая соединяет точки $A$ и $B$, $c\left( {\vec {r}} \right)$ – скорость волны в точке $\vec {r} = \left( {x,y} \right)$, лежащей на кривой $\gamma $, и $dl$ – элемент длины вдоль $\gamma $. Аппроксимируем интеграл (2), используя формулу трапеций, и обозначаем приближенное значение через $T({\mathbf{r}})$:
(3)
$T[\gamma ] \approx T({\mathbf{r}}) = \frac{1}{2}\sum\limits_{i = 0}^P \,\left( {\frac{1}{{{{c}_{{i + 1}}}}} + \frac{1}{{{{c}_{i}}}}} \right)\left| {{{{\vec {r}}}_{{i + 1}}} - {{{\vec {r}}}_{i}}} \right|,$
где ${\mathbf{r}} = \left( {{{{\vec {r}}}_{1}},{{{\vec {r}}}_{2}}, \ldots ,{{{\vec {r}}}_{P}}} \right)$ – многомерный вектор, компонентами которого являются векторы $\mathop {\vec {r}}\nolimits_i = ({{x}_{i}},{{y}_{i}})$ точек траектории; ${{\vec {r}}_{0}} = {{\vec {r}}_{A}}$; ${{\vec {r}}_{{P + 1}}} = {{\vec {r}}_{B}}$; ${{c}_{i}} = c({{\vec {r}}_{i}})$. Кривая $\gamma $ представляется в виде ломаной линии с $P + 2$ точками, где конечные точки зафиксированы, а положения $P$ промежуточных точек последовательно модифицируются в ходе итерационного процесса к оптимальной форме. Таким образом, граничная задача о расчете лучевых траекторий сводится к поиску экстремумов функционала (3). Ранее в работе [19] было показано, что лучи могут соответствовать двум типам стационарных решений: минимумам и седловым точкам первого порядка. Это позволило создать универсальный метод оптимизации для поиска различных типов экстремума, получивший название метод обобщенной силы [20]. Рассмотрим основные положения поиска минимумов и седловых точек методом обобщенной силы.

Поиск минимума основан на модификации градиента целевой функции:

(4)
${{{\mathbf{F}}}^{{{\text{min}}}}} = - \mathop {\left. {\nabla T({\mathbf{r}})} \right|}\nolimits_ \bot + {{{\mathbf{F}}}^{{\text{s}}}}.$
Здесь $\mathop {\left. {\nabla T({\mathbf{r}})} \right|}\nolimits_ \bot $ – поперечная компонента градиента многомерной функции, ${{{\mathbf{F}}}^{{\text{s}}}}$ – сила упругости, действующая между точками траектории и позволяющая контролировать плотность распределения точек вдоль траектории [22]. Обобщенная сила ${{{\mathbf{F}}}^{{{\text{min}}}}}$ формулируется в соответствии с методом упругой нити, применяемым для расчета оптимальных путей в различных областях науки [23], [24]. Поперечная компонента силы определяет смещение точек в направлении к искомой траектории, соответствующей решению краевой задачи, а замена продольной компоненты на упругие силы взаимодействия точек обеспечивает связаность цепочки точек и контролирует их распределение вдоль траектории [22].

Запишем выражение для компоненты силы упругости ${{{\mathbf{F}}}^{{\text{s}}}} = (\vec {F}_{1}^{{\text{s}}},\vec {F}_{2}^{{\text{s}}}, \ldots ,\vec {F}_{P}^{{\text{s}}})$, действующей на каждую $i$-ю точку траектории:

(5)
$\vec {F}_{i}^{{\text{s}}} = \kappa \left( {\left| {{{{\vec {r}}}_{{i + 1}}} - {{{\vec {r}}}_{i}}} \right| - \left| {{{{\vec {r}}}_{i}} - {{{\vec {r}}}_{{i - 1}}}} \right|} \right){{\vec {\tau }}_{i}},$
где $\kappa $ – константа силы упругости. Единичный вектор касательной ${{\vec {\tau }}_{i}}$ к траектории в каждой $i$-й точке может быть рассчитан, используя следующую формулу:

(6)
${{\vec {\tau }}_{i}} = \frac{{{{{\vec {r}}}_{{i + 1}}} - {{{\vec {r}}}_{{i - 1}}}}}{{\left| {{{{\vec {r}}}_{{i + 1}}} - {{{\vec {r}}}_{{i - 1}}}} \right|}}.$

Для поиска седловых точек обобщенная сила формулируется в соответствии с методом минимальной моды [16], [25]. Работа метода минимальной моды основана на модификации вектора силы таким образом, чтобы седловая точка преобразовывалась в локальный минимум. Для этого осуществляется расчет собственных значений матрицы Гессе (матрицы вторых производных) для функционала (3). Матрица Гессе определяется в соответствии с выражением:

(7)
${{H}_{{i,j}}} = \frac{{{{\partial }^{2}}T}}{{\partial {{\alpha }_{i}}\partial {{\alpha }_{j}}}},$
где ${{\alpha }_{i}}$ – компоненты многомерного вектора ${\mathbf{r}} = ({{\vec {r}}_{1}},{{\vec {r}}_{2}}, \ldots ,{{\vec {r}}_{P}}) = ({{\alpha }_{1}},{{\alpha }_{2}}, \ldots ,{{\alpha }_{M}})$, $M = P \times N$ – общее число компонент $P$ подвижных точек в $N$-мерном пространстве. Затем определяется минимальное собственное значение и соответствующий ему собственный вектор, представляющий собой минимальную моду. В соответствии с определением седловой точки инверсия минимальной моды приводит к выражению для модифицированной силы, соответствующей локальному минимуму:
(8)
${{{\mathbf{F}}}^{{{\text{saddle}}}}} = \left\{ \begin{gathered} (\mathop {\left. {\nabla T({\mathbf{r}})} \right|}\nolimits_ \bot \cdot {{{\mathbf{Q}}}_{\lambda }}){{{\mathbf{Q}}}_{\lambda }} + {{{\mathbf{F}}}^{{\text{s}}}},\quad {\text{если}}\quad {{\lambda }_{{{\text{min}}}}} \geqslant 0, \hfill \\ - \mathop {\left. {\nabla T({\mathbf{r}})} \right|}\nolimits_ \bot + 2(\mathop {\left. {\nabla T({\mathbf{r}})} \right|}\nolimits_ \bot \cdot {{{\mathbf{Q}}}_{\lambda }}){{{\mathbf{Q}}}_{\lambda }} + {{{\mathbf{F}}}^{{\text{s}}}},\quad {\text{если}}\quad {{\lambda }_{{{\text{min}}}}} < 0, \hfill \\ \end{gathered} \right.$
где ${{{\mathbf{Q}}}_{\lambda }}$ – нормированный собственный вектор, соответствующий минимальному собственному значению ${{\lambda }_{{\min }}}$, или минимальная мода. Выражения (4) и (8) составляют основу итерационной процедуры оптимизации, в ходе которой первоначальная конфигурация точек итерационно сходится к нулю обобщенной силы. Финальная конфигурация точек траектории есть дискретное представление лучевой траектории. Отметим, что процедура сходимости силы к нулю может быть реализована различными подходами, например, методом наискорейшего спуска, сопряженных градиентов и др. В нашей работе был выбран метод проецирования скорости, который является вариантом моделирования классической динамики с трением [26].

Представленный вариант вариационного подхода имеет ряд преимуществ перед существующими методами решения двухточечной задачи. Во-первых, фиксация граничных точек во время процесса оптимизации автоматически удовлетворяет краевым условиям, что делает его наиболее подходящим для решения двуточечной задачи, нежели применение численного решения уравнения эйконала и метода стрельбы. Во-вторых, впервые представлен вариант прямого вариационного метода, определяющий как минимумы, так и седловые точки, что предпочтительнее методов, основанных только на минимизации функционала [10], [12]. В-третьих, переходные свойства седловых точек позволили на основе метода обобщенной силы создать алгоритм систематического поиска лучей без необходимости подбора начальных приближений. Взаимосвязь стационарных точек функционала позволила решить проблему многолучевости последовательным поиском от предыдущего решения к последующему. Данный алгоритм получил название глобальная оптимизация и подробно описан в [20].

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

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

3.1. Модель подводной горы

Рассмотрим область океанического дна, задаваемая функцией глубины $D\left( {x,y} \right)$:

(9)
$D(x,y) = 1 - \frac{1}{2}{{(1 + {{(x - 2)}^{2}} + {{y}^{2}})}^{{ - 7/2}}}.$
Скорость распространения волны задается формулой Лагранжа $c = \sqrt {gD(x,y)} $, где в качестве упрощения $g = 1$.

Вначале воспользуемся вариационным подходом к расчету лучевых траекторий с заданными граничными условиями в данной аналитической модели среды. Зададим положения источника $A$ и приемника $B$ в точках с радиус-векторами ${{\vec {r}}_{A}} = (0,0)$ и ${{\vec {r}}_{B}} = (5,0.5)$. Инициализация метода осуществляется заданием начальных положений промежуточных точек траектории количеством $P = 13$ равномерно вдоль прямой линии, соединяющей точки $A$ и $B$. Как показано на фиг. 3, в результате итерационной процедуры минимизации начальное приближение (траектория 0) сходится к оптимальной конфигурации, соответствующей лучу 1 (локальный минимум). В данном случае принято, что процесс оптимизации считается завершенным, когда абсолютное значение обобщенной силы, действующей на каждую точку траектории, не превышает ${{10}^{{ - 9}}}$.

Фиг. 3.

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

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

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

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

Далее осуществим сравнение вариационного подхода и традиционного лучевого метода, основанного на решении уравнения эйконала. Для этого получим эталонное семейство лучей методом бихарактеристик (см. формулу (1)). На этапе инициализации зададим положение источника $A$ в точке ${{\vec {r}}_{A}} = (0,0)$ и последовательность начальных импульсов ${{\alpha }_{i}}$ в диапазоне –45–45° с постоянным шагом $\Delta \alpha = 1$. Решив гамильтонову систему уравнений методом Рунге–Кутты 4-го порядка с данными начальными условиями, получим семейство лучей, представленное на фиг. 4a. Численные результаты показывают, что лучевые траектории, преломляясь в области подводной горы, образуют фокус и простую каустику. Как следствие, рассчитанный по лучевому семейству фронт волны после прохождения подводной горы образует складку (см. фиг. 4a).

Фиг. 4.

Результаты расчетов методом бихарактеристик (а) и вариационным методом (б). Сплошные линии представляют лучевые траектории и фронты волны. Жирным выделены фронты в моменты времени ${{t}_{1}} = 1$, ${{t}_{2}} = 2.98$ и ${{t}_{3}} = 7$. В электронной версии статьи зеленые и фиолетовые линии представляют лучевые траектории, соответствующие минимумам и седловым точкам функционала (2), красные линии – фронты волны.

Воспроизведем полученные результаты с использованием вариационного метода. При этом важно учесть, что для вариационного подхода начальные углы излучения, в отличие от метода бихарактеристик, являются одними из искомых параметров. Поэтому для исследования выбранного сектора (см. фиг. 4a) последовательно закрепим точку $B$ на окружности радиуса $R = 8$ c заданным угловым шагом $\Delta \alpha = 1$. Применяя вариационный метод и алгоритм глобальной оптимизации для каждого положения точки $B$, получаем семейство лучевых траекторий и соответствующий фронт волны (см. фиг. 4б).

Как показано на фиг. 4, восстановленные фронты, полученные двумя различными методами, полностью согласуются. Выделим некоторые особенности применения рассмотренных методик. Семейство лучей, полученных методом бихарактеристик, равномерно покрывает исследуемое пространство в некоторой окрестности точки $A$. Однако из-за неоднородных свойств среды лучи, равномерно испускаемые в некотором угловом диапазоне, испытывают сгущение или разрежение вдали от источника. Как видно на фиг. 4a, плотность лучей вдоль волновых фронтов значительно меняется. В результате возможны случаи, когда количествa лучей, необходимых для детального исследования параметров волны в заданном пространстве, будет недостаточно. Для решения данной проблемы применяется метод стрельбы, который может приводить к значительным вычислительным затратам.

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

Отметим, что решение задачи методом бихарактеристик позволяет находить семейство лучей вне зависимости от их типа. Преимуществом вариационного метода является избирательность в определении типов лучей. Так, на фиг. 4б представлены два типа траекторий: минимумы и седловые точки первого порядка. Характер распространения семейства лучей, образованных минимумами и седловыми точками, различен. Так, только семейство лучей, соответствуюших седловым точкам первого порядка, формирует фокальную точку и каустику (см. фиг. 4).

Кроме того, применение вариационного подхода не требует больших временных затрат, что также является преимуществом. Среднее время расчета лучей для одной двуточечной задачи на персональном компьютере с процессором Intel Core i7-4770 составляет менее 1 сек. Для набора граничных условий (см. фиг. 4б) среднее время расчета не превышает 10 сек.

3.2. Модель подводного хребта

Рассмотрим другой пример морского дна в виде длинного подводного хребта, задаваемого функцией глубины:

(10)
$D(x,y) = 1 - \frac{1}{2}\mathop {\left( {1 + \mathop {\left( {\frac{{x - 15}}{{20}}} \right)}\nolimits^2 + {{y}^{2}}} \right)}\nolimits^{ - 2} .$
Аналогично предыдущему примеру, поместим источник в начало координат и последовательно воспользуемся двумя методиками. Лучевые траектории и фронт волны, представленные на фиг. 5а, получены методом бихарактеристик для начальных углов излучения ${{\alpha }_{i}}$ в диапазоне –45–45° с постоянным шагом $\Delta \alpha = 1$°. Особенностями распространения волны вдоль подводного хребта являются образование пространственно-временных каустик, сложной структуры фронта (см. фиг. 5a) и усложнение профиля волны.

Фиг. 5.

Результаты расчетов методом бихарактеристик (а) и вариационным методом (б). Обозначения те же, что и на фиг. 4. Жирным выделены фронты в моменты времени ${{t}_{1}} = 8$, ${{t}_{2}} = 12$ и ${{t}_{3}} = 15$.

Применение вариационного метода, как и в случае подводной горы, связано с последовательным закреплением конечной точки. В данном случае точка $B$ закреплялась вдоль прямой, проходящей по гребню подводной горы, с шагом $\Delta x = 0.5$ (см. фиг. 5б). Аналогично примеру, представленному в п. 3.1, для каждого положения точки $B$ начальное приближение задавалось вдоль прямой линии, после чего применялся систематический поиск всех решений. Лучевое семейство, полученное вариационным методом в ходе глобальной оптимизации (см. фиг. 5б), не полностью воспроизвело расчеты методом бихарактеристик. Были найдены только лучевые траектории, представленные на фиг. 5б, которые образовали волновые фронты, показанные сплошными линиями. На фиг. 5б пунктирными линиями отмечены участки волновых фронтов, которые не были восстановлены вследствие поиска минимумов и седловых точек первого порядка. При этом отметим, что ненайденные лучевые траектории, образующие каустические структуры, имеют характерную форму волноводного распространения вдоль максимума подводного хребта. Анализ собственных значений матрицы Гессе, выполненный для этих лучей, показал, что лучи, касающиеся каустик $n$ раз, соответствуют седловым точкам $n$ порядка. Следовательно, порядок седловой точки совпадает с индексом Маслова. Стоит отметить, что в данном примере луч, проходящий вдоль подводного хребта, является максимумом целевой функции, который был найден без труда прямой максимизацией функционала (2). Каждый тип стационарного решения образует соответствующий участок волнового фронта, который по мере распространения вдоль подводного хребта принимает сложную структуру.

В данной работе для поиска седловых точек высокого порядка впервые была апробирована модификация метода обобщенной силы, где преобразование седловых точек n-го порядка в локальный минимум (см. разд. 2) осуществлялось инверсией всех отрицательных мод гессиана целевой функции (2). При этом важной составляющей успешного поиска седел высшего порядка является правильное задание начального приближения. Для определения седловой точки $n$-го порядка начальное приближение необходимо задавать в окрестности найденного седла $(n - 1)$-го порядка. Это правило соответствует концепции глобальной оптимизации [20] и в данном случае имеет смысл определения полного числа решений, в результате последовательного поиска вариационным методом стационарных точек: минимум $ \to $ cедло I порядка $ \to $ cедло II порядка $ \to $ … $ \to $ максимум $ \to $ … $ \to $ cедло II порядка $ \to $ cедло I порядка $ \to $ минимум и т.д. Применение данного подхода позволило определить дополнительные решения граничной задачи. На фиг. 6 показаны все найденные лучи, среди которых впервые найдены волноводные траектории, соответствующие седловым точкам 2-, 3- и 4-го порядка. Успешный поиск новых решений граничной задачи позволил восстановить фронт волны и добиться полного согласия с методом бихарактеристик. Развитие глобальной оптимизации и метода обобщенной силы с учетом поиска всех типов решений граничной задачи является предметом дальнейших исследований.

Фиг. 6.

Результаты расчетов лучей вариационным методом с заданными граничными условиями. Градиентом выделена область подводного хребта. Цифрами 0–5 обозначены лучи, соответствующие минимумам, седловым точкам 1, 2, 3, 4-го порядка и максимуму функционала (3).

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

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

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

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

Авторы выражают благодарность П.Ф. Бессарабу за неоценимый вклад в развитие вариационного подхода и многолетнее плодотворное сотрудничество.

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

  1. Марчук А.Г., Чубаров Л.Б., Шокин Ю.И. Численное моделирование волн цунами. М.: Наука, 1983. 175 с.

  2. Пелиновский Е.Н. Гидродинамика волн цунами. Н. Новгород: ИПФ РАН, 1996. 276 с.

  3. Berry M.V. Focused tsunami waves // P. Roy. Soc. A-Math. Phy. 2007. T. 463. № 2087. P. 3055–3071.

  4. Satake K. Effects of bathymetry on tsunami propagation: Application of ray tracing to tsunamis // Pure Appl. Geophys. 1988. V. 126. № 1. P. 27–36.

  5. Ward S.N. Landslide tsunami // J. Geophys. Res. Solid Earth. 2001. V. 106. № B6. P. 11201–11215.

  6. Dobrokhotov S.Y., Sekerzh-Zenkovich S.Y., Tirozzi B., Tudorovskii T.Y. Description of tsunami propagation based on the Maslov canonical operator // Dokl. Math. 2006. V. 74. № 1. P. 592–596.

  7. Dobrokhotov S.Y., Shafarevich A.I., Tirozzi B. Localized wave and vortical solutions to linear hyperbolic systems and their application to linear shallow water equations // Russ. J. Math. Phys. 2008. V. 15. № 2. P. 192–221.

  8. Dobrokhotov S.Y., Nazaikinskii V.E. Punctured Lagrangian manifolds and asymptotic solutions of the linear water wave equations with localized initial conditions // Math. Notes. 2017. V. 101. № 5. C. 1053–1060.

  9. Калиткин Н.Н. Численные методы. 2 изд. БХВ-Петербург, 2011. 592 с.

  10. Um J., Thurber C. A fast algorithm for two-point seismic ray tracing // Bull. Seismol. Soc. Am. 1987. T. 77. № 3. P. 972–986.

  11. Moser T.J., Nolet G., Snieder R. Ray bending revisited // Bull. Seismol. Soc. Am. 1992. V. 82. № 1. P. 259–288.

  12. Coleman C.J. Point-to-point ionospheric ray tracing by a direct variational method // Radio Sci. 2011. V. 46. № 05. P. 1–7.

  13. Бензик А.В. Прогнозирование траекторных характеристик распространения ДКМ-радиоволн на основе принципа Ферма // Техн. радиосвязи. 2014. № 1. С. 32–41.

  14. Mills G., Jónsson H. Quantum and thermal effects in H 2 dissociative adsorption: Evaluation of free energy barriers in multidimensional quantum systems // Phys. Rev. Lett. 1994. T. 72. № 7. P. 1124.

  15. Jónsson H., Mills G., Jacobsen K.W. Nudged elastic band method for finding minimum energy paths of transitions / Classical and Quantum Dynamics in Condensed Phase Simulations. World Scientific, 1998. P. 385–404.

  16. Henkelman G., Jónsson H. A dimer method for finding saddle points on high dimensional potential surfaces using only first derivatives // J. Chem. Phys. 1999. V. 111. № 15. P. 7010–7022.

  17. Bessarab P.F., Uzdin V.M., Jónsson H. Size and shape dependence of thermal spin transitions in nanoislands // Phys. Rev. Lett. 2013. V. 110. № 2. P. 020604.

  18. Bessarab P.F., Uzdin V.M., Jónsson H. Method for finding mechanism and activation energy of magnetic transitions, applied to skyrmion and antivortex annihilation // Comp. Phys. Commun. 2015. V. 196. P. 335–347.

  19. Nosikov I.A., Klimenko M.V., Bessarab P.F., Zhbankov G.A. Application of the nudged elastic band method to the point-to-point radio wave ray tracing in IRI modeled ionosphere // Adv. Space Res. 2017. V. 60. № 2. P 491–497.

  20. Nosikov I.A., Klimenko M.V., Zhbankov G.A., Podlesnyi A.V., Ivanova V.A., Bessarab P.F. Generalized Force Approach to Point-to-Point Ionospheric Ray Tracing and Systematic Identification of High and Low Rays // IEEE Trans. Antennas Propag. 2020. V. 68. № 1. P. 455–467.

  21. Kabanikhin S., Hasanov A., Marinin I., Krivorotko O., Khidasheli D. A variational approach to reconstruction of an initial tsunami source perturbation // Appl. Numer. Math. 2014. V. 83. P. 22–37.

  22. Носиков И.А., Бессараб П.Ф., Клименко М.В. Применение метода поперечных смещений для расчета коротковолновых радиотрасс. Постановка задачи и предварительные результаты // Известия вузов. Радиофизика. 2016. Т. 59. № 1. С. 1–14.

  23. Einarsdóttir D.M., Arnaldsson A., Óskarsson F., Jónsson H. Path optimization with application to tunneling // International Workshop on Applied Parallel Computing. Berlin, Heidelberg, Springer: 2010. P. 45–55.

  24. Ásgeirsson V., Arnaldsson A., Jónsson H. Efficient evaluation of atom tunneling combined with electronic structure calculations // J. Chem. Phys. 2018. T. 148. № 10. P. 102334.

  25. Gutiérrez P.M., Argáez C., Jónsson H. Improved minimum mode following method for finding first order saddle points // J. Chem. Theory Comput. 2017. V. 13. № 1. P. 125–134.

  26. Andersen H.C. Molecular dynamics simulations at constant pressure and/or temperature // J. Chem. Phys. 1980. V. 72. № 4. P. 2384–2393.

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

Инструменты

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

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

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