CC BY
11
Вычислительные технологии
Том 24, № 2, 2019
Численное моделирование газа солитонов в рамках уравнений типа Кортевега — де Вриза*
Е. Г. ДиДЕНКУЛОВА^ А. В. КокоринА, А. В. СлЮНЯЕВ Институт прикладной физики РАН, Нижний Новгород, Россия ^Контактный e-mail: eshurgalina@mail.ru
Приведены детали численной схемы и способа задания начальных условий для моделирования нерегулярной динамики ансамблей солитонов в рамках уравнений типа Кортевега — де Вриза на примере модифицированного уравнения Кортевега— де Вриза с фокусирующим типом нелинейности. Дано качественное описание эволюции статистических характеристик для ансамблей солитонов одной и разных полярностей. Обсуждаются результаты тестовых экспериментов по столкновению большого числа солитонов.
Ключевые слова: солитонный газ, численное моделирование, уравнение Кортевега—де Вриза, модифицированное уравнение Кортевега —де Вриза, уравнение Гарднера.
Библиографическая ссылка: Диденкулова Е.Г., Кокорина А.В., Слюняев А.В. Численное моделирование газа солитонов в рамках уравнений типа Кортевега — де Вриза // Вычислительные технологии. 2019. Т. 24, № 2. С. 52-66. DOI: 10.25743/ICT.2019.24.2.005.
Введение
Под газом солитонов подразумевают ансамбль солитонов со случайными параметрами. Свойство солитонов взаимодействовать упруго друг с другом порождает очевидную ассоциацию с газом упруго сталкивающихся частиц. Строго говоря, солитонный газ — это детерминированная динамическая система в силу интегрируемости уравнений, описывающих эволюцию волн (солитонов) [1]. Однако из-за сложности ее поведения (обусловливаемой большим числом участвующих солитонов и нелинейным характером их взаимодействий) динамика системы может рассматриваться как случайная и, соответственно, исследоваться методами, типичными для таких задач. По этой причине эволюцию солитонного газа рассматривают как объект солитонной или интегрируемой турбулентности. В частности, могут вычисляться функции распределения вероятностей параметров волн, статистические моменты и т. д.
В первую очередь, исследования солитонного газа представляют собой фундаментальную задачу. Поведение волн в интегрируемых системах демонстрирует неординарные эффекты, например классический парадокс Ферми — Пасты — Улама [2], а также недавно обнаруженное возвращение функций распределения к линейным законам
* Title translation and abstract in English can be found on page 66.
© ИВТ СО РАН, 2019.
на больших временах [3]. Кинетические уравнения для ансамблей солитонов были записаны практически одновременно с построением метода обратной задачи рассеяния [4, 5]. Они описывают транспорт спектральных данных ассоциированной задачи рассеяния, но не дают информации о фазах (полярности) солитонов и собственно о полях волн, потому непригодны для исследования их статистики. Прямое численное моделирование ансамблей волн здесь выступает в качестве удобной альтернативы.
Для исследования эволюции волн на очень больших временах, когда наблюдаются эффекты рекуррентности, точность алгоритмов расчета эволюции нелинейных волн имеет огромное значение [3, 6]. Однако для практических приложений большое значение имеет изучение связи между параметрами солитонов и соответствующими вероятностными свойствами солитонного газа [7]. Для такого класса задач (относительно небольшие времена расчета) требования к численному коду значительно мягче, большее значение имеет возможность накопления богатых статистических ансамблей.
Отметим, что описанные проблемы актуальны не только в рамках интегрируемых систем, но и для близких к ним, допускающих долгоживущие решения в виде "почти солитонов". Кроме того, они могут рассматриваться и в более широком смысле — применения подхода Монте-Карло для расчета стохастической волновой динамики. Так, численное моделирование ансамблей уединенных волн типа солитонов и солитонов огибающей на фоне других волн проводилось во многих работах (см. литературу в [8] и более свежие работы [3, 7, 9, 10]). Численное моделирование ансамблей локализованных групп типа одиночных и мультисолитонов в рамках интегрируемых и неинтегрируемых уравнений выполнялось для оптических волноводов в целях определения безопасных для линий связи кодирующих паттернов [11]. Интенсивное численное моделирование приближенных интегрируемых уравнений и исходных уравнений гидродинамики для потенциальных движений проведено в последние годы в рамках исследований так называемых волн-убийц в океане (см. обзор [12]). Для таких существенно нелинейных и негауссовых событий предположения кинетической теории нарушаются, и для расчета связи между частотно-угловыми спектрами волн и их вероятностными характеристиками производится прямое численное решение нестационарных уравнений гидродинамики для ансамблей нерегулярных волн [13-17].
В настоящей работе обсуждается специфика расчета ансамблей солитонов с целью определения вероятностных свойств в рамках модельного уравнения — модифицированного уравнения Кортевега — де Вриза (мКдВ) с фокусирующим типом нелинейности. Нужно отметить, что схожим образом авторами выполнено моделирование солитон-ного газа и в рамках других родственных интегрируемых уравнений — классического уравнения Кортевега — де Вриза (КдВ) [18, 19] и уравнения Гарднера, представляющего собой комбинацию этих двух уравнений [20]. Отметим, что из-за более высокой степени нелинейности численное моделирование уравнения мКдВ сложнее, чем расчет классического уравнения КдВ с квадратичной нелинейностью.
Следует добавить, что качественное представление о динамике волн в не слишком плотном солитонном газе можно получить уже на примере предельного случая столкновения двух солитонов. В силу наличия соответствующих аналитических решений и сохраняющихся интегралов для уравнений картина двухсолитонного взаимодействия может быть рассмотрена достаточно детально, включая получение некоторых аналитических соотношений. Столкновение солитонов как элементарный акт турбулентности было исследовано в рамках уравнения КдВ в работе [21], для его модифицированной версии — в публикациях [22, 23] и для уравнения Гарднера — в [24, 25]. Хотя взаи-
модействие множества солитонов также представляется формально известным точным решением, исследование которого из-за сложности с использованием аналитических подходов сильно ограничено. В плане исследования мультисолитонов наиболее ярким недавним результатом являются формулы для максимально достижимой амплитуды волн, представляющие собой линейную суперпозицию парциальных амплитуд солитонов и уровня фона для фокусирующего нелинейного уравнения Шрёдингера [26, 27], а также для уравнений мКдВ [28] и Гарднера [29] фокусирующих типов.
Настоящая работа имеет следующую структуру. В разд. 1 кратко описан численный алгоритм решения модифицированного уравнения Кортевега — де Вриза. Способ задания начальных условий в виде квазислучайного набора солитонов и основные результаты численного моделирования солитонного газа описаны в разд. 2. Возможности использованного кода описывать процессы столкновения большого числа солитонов, в том числе возникающие "волны-убийцы", исследованы в разд. 3. Далее приводятся заключительные замечания.
1. Алгоритм решения модифицированного уравнения Кортевега — де Вриза
Для численного интегрирования уравнений длинных нелинейных волн типа КдВ в литературе представлен широкий спектр численных подходов. Среди основных можно выделить использование явных и неявных конечно-разностных схем, псевдоспектральный подход, метод сплит-степ [30-33]. В настоящей работе моделируется начальная задача для модифицированного уравнения Кортевега — де Вриза, записанного в стандартной безразмерной форме
ди + би2 + ®3и _ о (1)
дЬ дх дх3
для действительной функции координаты и времени и(х,Ь), заданной в начальный момент времени и(х,Ь _ 0). Для численного интегрирования уравнения (1) используем псевдоспектральный метод [33], который позволяет вычислять с высокой точностью производные по пространственной переменной в пространстве Фурье, и конечно-разностный подход для интегрирования по времени.
Функция и(х,Ь) вычисляется на пространственном интервале длиной Ь и является периодической на нем. Она задается в N узлах координатной сетки, так что в реальном пространстве расстояние между точками одинаково и равно Ах _ Ь/Ы. Решение представимо в виде дискретного ряда Фурье по пространственной переменной х^, 3 ,
N-1
u(xj,t) = 2^ Ck(t) e k=0
L
где Ck (t) = F {и} — компоненты пространственного преобразования Фурье и к — целые числа, нумерующие массив волновых чисел. Использование представления Фурье позволяет аккуратно вычислять производные по координате; производная п-го порядка записывается следующим образом:
д пи
.............! и !
дхп
F-1{(ik)nF{u}} .
Нелинейное слагаемое в уравнении (1) удобно представить в виде ди3/дх, что позволяет сократить количество операций вычисления преобразования Фурье. С учетом (2) уравнение (1) может быть записано в пространстве Фурье в виде системы обыкновенных дифференциальных уравнений для функций С к:
ССк а( (-1\3
+ 6(г k)Qk + (г к)3Ск = 0, (3)
где Qk = F{u3/3}. Для интегрирования системы (3) по времени используется неявная трехслойная конечно-разностная схема (leapfrog) [33]
y(t+At) , n(t-At)
j rl(l+Al) ^ П(
U sy(t) _ ck + Ck
сИ к 2Аt
Для повышения устойчивости такой численной схемы в работе [34] предложено использовать неявное представление линейного дисперсионного слагаемого в виде полусуммы компонент Фурье на временных слоях (£ — А£) и (£ + А£) (схема Кранка — Николсон):
С(¿+А*) _ С(г-ы) (г к)3&? = (г к)3 С-—С-. (4)
Тогда с учетом (4) для интегрирования уравнения мКдВ (1) применяется разностная схема
= 1-А€) (1 + г к3 АЪ \ Ук Ск I 1 ¿и3**) I 1 _¡к3Аг
c(t+At) = c(t-At) (1 + jк3At \ (t) ( 2ikAt \ Ck =Ck \1 - ik3At) - 6Qk \1 -г k3At)
которая имеет второй порядок точности по времени и является безусловно устойчивой [33]. Для вычисления функции на среднем слое на первом шаге моделирования используется явная схема [32]
сЦ+А) = С® (1 + гк3А£) — %гкАtQк).
Для контроля аккуратности вычислений в экспериментах отслеживалось сохранение интегралов уравнения (1), условно соответствующих "массе" М и "энергии" Е:
М = udx, Е = u2dx.
Поскольку тестовые примеры не выявили необходимости применения деалиазинга (de-aliasing) [33], эта процедура не применялась.
2. Численное моделирование ансамбля солитонов
Начальное условие при t = 0 для компьютерного моделирования задавалось в виде набора солитонов, являющихся точными решениями уравнения мКдВ:
An , I 1
u(x, 0) / USol,n(x xn, °) , USol,n(x, t) 1/.1 / ЧО| \ \, Xn d [ П I .
cosh (An(x — Ant)) V 2/
n=l
В (6) функция и301,п(х, ¿) — точное односолитонное решение (1) с действительной амплитудой Ап (которая может быть как положительной, так и отрицательной) и скоростью АП > 0. Интегралы уравнения (5) для солитона (6) с амплитудой А определяются соотношениями
В каждой реализации задавалось N3 = 30 солитонов с одинаковым набором амплитуд, чьи абсолютные величины располагались эквидистантно в интервале [1, 3]. Таким образом, распределение абсолютных величин амплитуд в каждой реализации соответствовало равномерному в этом интервале. Полярность солитонов (знаки ^п) задавалась либо одинаковой (однополярный газ), либо число положительных и отрицательных со-литонов выбиралось равным, так что интеграл массы для всего поля был равен нулю, М = 0 (разнополярный газ). Координаты вершин соседних солитонов отличались на заданную величину й = 20, которая выбиралась достаточно большой, чтобы обеспечить хорошее спадание экспоненциальных "хвостов" солитонов (рис. 1). В то же время расстояние между солитонами не было настолько большим, чтобы газ не оказался слишком разреженным. Характерная плотность газа в расчетах составляла р = М3/Ь = 0.05 (Ь — размер вычислительной области) при критической плотности 2(Ап)/ж2 = 0.4 [35]. Порядок солитонов с заданными амплитудами и полярностями выбирался случайным образом.
Отметим, что заданное таким способом начальное условие уже содержит в себе погрешность из-за переналожения "хвостов" соседних солитонов; задаваемое поле не является точным ^-солитонным решением уравнения мКдВ. Погрешность в задании поля может быть оценена как равная А весЬ(А^/2), что для выбранного диапазона амплитуд солитонов дает величины от 10-4 до 10-13. В принципе, задавать начальное условие можно и с помощью точного мультисолитонного решения, находимого с помощью численных методов (тогда начальные координаты солитонов могут быть случайными). На сегодня предложены способы построения мультисолитонных решений, включающих порядка сотни солитонов [10, 11], но эти процедуры не слишком просты с вычислительной точки зрения, поскольку основаны на расчете вронскианов высоких порядков от функций с экспоненциальными зависимостями. Наконец, под солитонами обычно подразумевают локализованные решения на бесконечном интервале, в то время как нами моделируется динамика в ограниченной области с периодическими граничными условиями (в этом случае говорят о конечно-зонных решениях).
МаЫ = ж А, ЕаЫ = 2|А|.
3
2
-2
-3--1-1-1-1-1-
0 100 200 300 400 500 600
х
Рис. 1. Пример начального условия для солитонов разной полярности и(х,Ь = 0)
Невырожденное мультисолитонное решение предполагает, что у солитонов все амплитуды различны и область распространения солитонов бесконечна. В нашем случае способ задания амплитуд солитонов формально гарантирует невырожденность решения. Разумеется, возможно задание амплитуд солитонов, подчиняющихся другим функциям распределения. Не слишком малая разница между амплитудами (и, соответственно, скоростями) солитонов также важна для исключения ситуаций связывания солитонов разной полярности, как обсуждается далее в разд. 3. Задание амплитуд соли-тонов и их позиций фактически происходит детерминистским способом, но из-за разных скоростей солитонов и нелинейного характера их взаимодействия (в частности, нелинейных сдвигов координат, возникающих в процессах столкновений) вскоре картина движения солитонов приобретает хаотичный характер (рис. 2).
Для рассматриваемого класса задач перечисленные несоответствия между математически строгой постановкой задачи и фактически исследуемой не представляются значительными. В частности, в работе [36] выполнено численное моделирование ансамблей солитонов в рамках интегрируемого уравнения Кортевега — де Вриза и его неинтегри-руемой модификации, так называемой регуляризованной модели Бенджамина — Бона — Махони (Benjamin — Bona — Mahony). Здесь существует точное аналитическое решение в виде уединенной волны, близкое к солитонному решению уравнения КдВ, но такие солитоноподобные волны взаимодействуют не полностью упруго. Численное моделирование интегрируемого и неинтегрируемого уравнений не выявило заметных различий в поведении статистических характеристик волн.
В наших экспериментах расчет эволюции солитонов проводился до времени t = 100 в расчетной области размером L = 600 на N = 214 узлах координатной сетки с шагом по времени At = 2.5 • 10-5. Результат моделирования эволюции в рамках уравнения (1) сохранялся для последующей обработки каждые 0.05 единиц времени. При выбранных условиях типичная относительная ошибка сохранения интегралов движения М и Е есть 10-13 и 10-6. Траектории движения солитонов в одной реализации показаны на рис. 2. Видно, что в основном преобладают парные взаимодействия, хотя иногда в небольшой окрестности сходятся 3-4 солитона.
Рассматриваемые далее асимметрия волн А3 и эксцесс А4 задаются через статистические моменты тп следующим образом:
X
Рис. 2. Траектории движения солитонов в одной реализации
где
сю
тп^) = J уп4х, п = 1, 2,...
—ю
Напомним, что в силу законов сохранения уравнения мКдВ (5) первый момент (среднее значение и(х) при фиксированном времени, т\ = М) и второй момент (интеграл энергии, т2 = Е) с хорошей точностью остаются постоянными во времени. Статистические моменты более высоких порядков изменяются с течением времени.
Как можно было ожидать, параметр асимметрии после усреднения по большому числу реализаций примерно равен нулю в случае газа разнополярных солитонов и больше нуля в ситуации солитонов положительной полярности. Эксцесс положителен в обеих ситуациях — солитонов одного или разных знаков. Хотя начиная с какого-то момента времени значения Л3 и Л4 остаются примерно постоянными, по результатам численного моделирования обнаружена переходная область на стадии установления квазистационарных значений величин Л3 и Л4 (рис. 3 и 4). Начальное условие и(х,1 = 0) является вырожденным в том смысле, что в этот момент солитоны практически не взаимодействуют. По этой причине на рис. 3 и 4 хорошо видна переходная стадия в течение первых примерно 5 единиц времени от абсолютно изолированных к взаимодействующим солитонам. Наблюдаемые отличия квазистационарных значений асимметрии и эксцесса от начальных соотносятся с результатами анализа статистических моментов для примеров парных взаимодействий солитонов, который может быть легко выполнен
т
Рис. 3. Эволюция асимметрии для случая однополярных солитонов
15.5
12 -
11.5'-1-1-1-1-
0 20 40 60 80 100
?
Рис. 4. Эволюция эксцесса для случаев однополярных (сплошная кривая) и разнополярных (пунктирная кривая) солитонов
на основе аналитических двухсолитонных решений [21]. Асимметрия однополярных солитонов несколько меньше, чем в начальный момент времени (см. рис. 3), а эксцесс — меньше для однополярных солитонов, но больше для солитонов разной полярности (см. рис. 4).
Ситуация заметной динамики статистических моментов исключает возможность усреднения данных по времени и требует больших ансамблей солитонов. Поскольку необходимое на вычисление преобразования Фурье время растет с числом узлов быстрее, чем линейная функция, более выигрышным является наращивание объема статистических данных за счет большего числа реализаций, а не большей области расчета. Приведенные здесь результаты статистической обработки основываются на расчетах 50 реализаций солитонных полей (что включает в себя 1500 солитонов для каждых условий). После наступления стационарного состояния дополнительно может использоваться и усреднение по времени в предположении эргодичности процесса. Как упоминалось выше, быстрая стохастизация волн и достижение квазистационарного состояния, видимые на рис. 3 и 4, получаются в результате нелинейных сдвигов солитонов. Достигаемое ансамблем солитонов квазистационарное состояние может представляться "термодинамическим пределом" солитонной системы, к которому относятся получаемые статистические результаты [37, 38].
Эксцесс характеризует соотношение между вероятностями больших и малых значений поля и и для нормального случайного распределения равен А4 = 0. Таким образом, динамика солитонов разной полярности демонстрирует в статистическом смысле более экстремальное поведение, чем газ солитонов одной полярности с теми же характерными амплитудами по абсолютной величине (см. рис. 4). Эти выводы согласуются с распределениями вероятности амплитуд волн, построенными по данным численных экспериментов (рис. 5). Здесь под амплитудами понимались высоты локальных максимумов поля. Действительно, в ситуации солитонов одной полярности максимально наблюдавшееся значение и не превосходило 3 — максимальной амплитуды солитона в наборе. В поле солитонов разной полярности наблюдалось множество событий с большими амплитудами.
Обнаруженное качественное различие статистических характеристик газа солито-нов одного и разных знаков согласуется с прогнозом, сделанным на основе точного
10° 10"1
Ьн
ю-2 10~3
1 1.5 2 2.5 3 3.5 4 4.5 5 \А\
Рис. 5. Функции распределения вероятности амплитуд волн для случаев однополярных (мелкие точки) и разнополярных (крупные точки) солитонов. Штриховая линия — начальное равномерное распределение амплитуд солитонов
мультисолитонного решения уравнения мКдВ [28]. В точке столкновения нескольких солитонов максимально достижимая амплитуда равна сумме парциальных амплитуд солитонов (при правильно выбранных знаках солитонов). Учитывая диапазон амплитуд заданных солитонов (от 1 до 3), из рис. 5 можно сделать вывод об эффективных парных взаимодействиях солитонов в проведенных численных расчетах, несмотря на кажущиеся многосолитонные столкновения на рис. 2. Отсюда также следует заключение о редкости столкновений большого числа солитонов. Подобные события возможны при большем объеме статистических данных. Следы многосолитонных столкновений наблюдались в распределениях вероятности по результатам численного моделирования ансамблей солитонов огибающей нелинейного уравнения Шрёдингера в [10]. Возможность используемого нами алгоритма интегрирования уравнения (1) аккуратно описывать такие события исследуется в разд. 3.
3. Столкновение большого числа солитонов
Распространение одного солитона и столкновение пар солитонов являются стандартными тестами для проверки численных кодов интегрируемых уравнений. Как обсуждалось выше, в нашем случае важно аккуратное описание столкновения большого (больше чем два) числа солитонов. Наиболее резкая эволюция с формированием очень высоких волн происходит в случае так называемых синхронной или оптимальной фокусировки солитонов [28, 39]. При этом максимально высокая волна образуется в ситуации знакопеременного цуга фокусирующихся солитонов.
Для проведения теста по моделированию оптимальной фокусировки трех солитонов в качестве начального условия используем точное трехсолитонное решение уравнения мКдВ (1) в момент времени £ = -10 (рис. 6) [28]. Оно представляет собой последовательность солитонов чередующейся полярности, которые упорядочены по скорости: левый солитон на рис. 6 при £ = -10 — самый быстрый. С течением времени за счет разницы в скоростях солитоны фокусируются в момент £ = 0 в точке х = 0, а затем снова расходятся. На рис. 6 профили волн представлены в системе отсчета, сопровождающей самый быстрый солитон. За счет нелинейного сдвига фаз в процессе взаимодействий быстрые солитоны сдвигаются вправо, медленные — влево. На рис. 6 профили построены двумя линиями для двух разных шагов по времени. При выборе шага = 5 • 10—5 (широкая красная линия) результат численного расчета неотличим на глаз от соответствующего точного решения уравнения (1). Такой выбор дискретизации примерно соответствует условиям моделирования солитонных ансамблей, описанного в разд. 2. При некотором огрублении параметров эксперимента (увеличении шагов по времени или координате) процесс фокусировки солитонов качественно сохраняется: солитоны фокусируются с образованием очень высокой волны, но в процессе их расхождения численное решение отклоняется от аналитического. В первую очередь это проявляется для большего шага по времени ДЪ = 1 • 10—4 в изменении положения солитонов при сохранении ими амплитуды (рис. 6, тонкая черная линия).
При еще более грубой дискретизации может возникать другой эффект динамики разнополярных солитонов — их связывание в стабильную структуру переменной формы — бризер мКдВ. Видимо, впервые такой эффект наблюдался при численном моделировании уравнения мКдВ в [40], где было показано, что наличие малой численной вязкости способно приводить к образованию бризеров из пар солитонов разной полярности с близкими скоростями. На языке задачи рассеяния для мКдВ этот эффект соответ-
Рис. 6. Численное решение уравнения мКдВ для трех фокусирующихся солитонов с амплитудами 3, 2.85, 2.7 и шагами по времени ДЪ = 5 • 10-5 (широкая красная линия) и ДЪ = 1 • 10-4 (тонкая черная линия). Шаг по координате в обоих случаях Дх ~ 0.03. В первом случае результат неотличим от аналитического решения. Профили решения в различные моменты времени даны в системе отсчета, сопровождающей самый быстрый солитон со скоростью 9
ствует бифуркации двух близких дискретных собственных значений на оси комплексной спектральной плоскости в пару близких комплексно-сопряженных значений.
Связывание разнополярных солитонов в численном расчете может возникать, когда скорости солитонов становятся слишком близкими, так что ошибка накапливается за время очень долгого расхождения и приводит к формированию бризера. Такая ситуация показана на рис. 7 для тех же параметров численного эксперимента, что и для успешного расчета, приведенного на рис. 6 (Д£ = 5 • 10-5, Ах ~ 0.03), но при меньшей разнице между амплитудами солитонов. Два солитона с меньшими скоростями не расходятся после фокуса, £ > 0, а периодическим образом сталкиваются (см. момент £ ^ 9 на рис. 7), при этом полярности ведущей и следующей волн попеременно меняются. Таким образом, для выбранных шагов по координате и времени (при определенной амплитуде солитонов) процесс столкновения солитонов может описываться качественно неверно, если их амплитуды слишком близки.
Отметим, что изменение положения солитонов на рис. 6 для случая ДЪ =1 • 10-4 (тонкая черная линия) в сравнении с точным решением (которое совпадает с широкой красной линией на рис. 6) соответствует их тенденции к связыванию, но в этом случае солитоны успели разойтись и бризера не возникло. При расчете эволюции волн, представленных на рис. 6, с еще более крупным шагом (Д£ = 2 • 10-4) происходит связывание солитонов, подобно случаю на рис. 7. Способ задания амплитуд солитонов в начальных условиях для расчета ансамблей солитонов с определенной минимальной разницей, описанный в разд. 2, позволяет исключить такие паразитные эффекты.
С ростом числа сфокусированных солитонов амплитуда волны в фокусе увеличивается (равна сумме амплитуд всех солитонов, если они знакопеременны, и обладает полярностью, совпадающей с самым быстрым солитоном), а ее характерный размер по координате уменьшается. Соответственно, численное моделирование фокусировки большего числа солитонов усложняется. Для параметров успешного эксперимента на
Рис. 7. Численное моделирование распада точного трехсолитонного решения с амплитудами 3, 2.97, 2.94, приводящего к образованию связанного состояния солитонов — бризеру, шаг по времени ДЪ = 5 • 10-5, шаг по координате Дх ~ 0.03
рис. 6 заметное отклонение численного решения от аналитического проявляется уже при фокусировке четырех солитонов с амплитудами 3, 2.85, 2.7, 2.55 (пара солитонов разных знаков со средними скоростями связывается).
Заключение
Приведены описание кода и детали подхода для задания начальных условий и проведения численного моделирования ансамблей солитонов, которые были использованы для исследования статистических характеристик солитонного газа в интегрируемых уравнениях: Кортевега —де Вриза, модифицированного КдВ и уравнения Гарднера [18, 35, 41]. Модификации описанного в разд. 1 алгоритма расчета уравнения мКдВ на случаи квадратичной нелинейности или смешанной квадратичной и кубической нелинейности тривиальны. В работе приведены результаты моделирования наиболее интересных тестовых задач столкновения множества солитонов.
Подчеркнем исключительную пригодность использованного кода для решения рассмотренных задач. Несмотря на простоту алгоритма и относительно невысокую точность сохранения второго интеграла уравнения Е, примененная схема обеспечивает достаточную точность при описании даже сложных ситуаций фокусировки множества солитонов с формированием очень больших волн. Скорость вычисления, имеющая существенное значение при расчете больших ансамблей солитонов, оказывается также высокой. Пробные расчеты задачи столкновения трех солитонов с теми же амплитудами и той же степенью дискретизации, как в успешном эксперименте на рис. 6 (Д£ = 5 • 10-5, Дх ~ 0.03), с помощью кодов, использующих сплит-степ-Фурье [42] и метод Рунге — Кутты четвертого порядка, дали неудовлетворительные результаты (неверные позиции солитонов и их связывание при заметном увеличении компьютерного времени расчета).
Для расчета сложных ситуаций столкновения солитонов с близкими скоростями должны использоваться более сложные алгоритмы решения эволюционного уравнения, например использующие адаптивные дискретные сетки [38].
В заключение отметим, что схожие с обсуждаемыми в разд. 2 переходные эффекты в эволюции вероятностных характеристик из заданных изначально квазислучайных полей наблюдались и в численном моделировании других задач, инициированном работами по проблеме волн-убийц. Они присутствуют в ситуациях волновых полей общего вида (состоящих из солитонов и квазилинейных волн) в рамках интегрируемых уравнения КдВ и нелинейного уравнения Шрёдингера, а также при моделировании нерегулярных волн в рамках исходных уравнений гидродинамики и в лабораторных условиях, если начальные условия отвечают определенным критериям, благоприятным для формирования когерентных солитоноподобных волн или групп волн [8].
Благодарности. Результаты, представленные в разд. 2, получены при поддержке гранта РФФИ № 16-32-60012. Работы по разд. 1 и 3 выполнены в рамках проекта РФФИ № 18-02-00042 и при поддержке программы фундаментальных исследований президиума РАН "Нелинейная динамика: фундаментальные проблемы и приложения".
Список литературы / References
[1] Захаров В.Е., Манаков С.В., Новиков С.П., Питаевский Л.П. Теория солитонов. Метод обратной задачи. М.: Наука, 1980. 319 с.
Novikov, S.P., Manakov, S.V., Pitaevskiy, L.P., Zakharov, V.E. Theory of solitons: the inverse scattering method. N.Y.: Consultants Bureau, 1984. 319 p.
[2] Zabusky, N.J., Kruskal, M.D. Interaction of "solitons" in a collisionless plasma and the recurrence of initial states // Phys. Review Letters. 1965. Vol. 15, iss. 6. P. 240-243.
[3] Agafontsev, D.S., Zakharov, V.E. Intermittency in generalized NLS equation with focusing six-wave interactions // Phys. Letters A. 2015. Vol. 379(40-41). P. 2586-2590.
[4] Захаров В.Е. Кинетическое уравнение для солитонов // Журн. эксп. и теор. физики. 1971. Т. 60. С. 993-1000.
Zakharov, V.E. Kinetic equation for solitons //J. of Exp. and Theor. Physics. 1971. Vol. 33, No. 6. P. 538-541.
[5] El, G.A., Kamchatnov, A.M. Kinetic equation for a dense soliton gas // Phys. Review Letters. 2005. Vol. 95. P. 204101.
[6] Slunyaev, A., Dosaev, A. On the incomplete recurrence of modulationally unstable deep-water surface gravity waves // Communic. in Nonlinear Sci. and Numer. Simulation. 2019. Vol. 66. P. 167-182.
[7] Soto-Crespo, J.M., Devine, N., Akhmediev, N. Integrable turbulence and rogue waves: breathers or solitons? // Phys. Review Letters. 2016. Vol. 116. P. 103901.
[8] Слюняев А.В., Сергеева А.В. Стохастическое моделирование однонаправленных интенсивных волн на глубокой воде в приложении к аномальным морским волнам // Письма в Журн. эксп. и теор. физики. 2011. Т. 94, № 10. С. 850-858.
Slunyaev, A.V., Sergeeva, A.V. Stochastic simulation of unidirectional intense waves in deep water applied to rogue waves // JETP Letters. 2011. Vol. 94, No. 10. P. 779-786.
[9] Slunyaev, A., Sergeeva, A., Pelinovsky, E. Wave amplification in the framework of forced nonlinear Schrodinger equation: the rogue wave context // Physica D. 2015. Vol. 303. P. 18-27.
[10] Gelash, A., Agafontsev, D. Strongly interacting soliton gas and formation of rogue waves // Phys. Review E. 2018. Vol. 98, No. 4. P. 042210.
[11] Frumin, L.L., Gelash, A.A., Turitsyn, S.K. New approaches to coding information using inverse scattering transform // Phys. Review Letters. 2017. Vol. 118, No. 22. P. 223901.
[12] Слюняев А.В. Морские "волны-убийцы": прогноз возможен? // Вестн. МГУ. Сер. 3. Физика и астрономия. 2017. № 3. С. 33-47.
Slunyaev, A.V. Predicting rogue waves // Moscow Univ. Physics Bull. 2017. Vol. 72. P. 236-249.
[13] Xiao, W., Liu, Y., Wu, G., Yue, D.K.P. Rogue wave occurrence and dynamics by direct simulations of nonlinear wave-field evolution // J. of Fluid Mech. 2013. Vol. 720. P. 357-392.
[14] Sergeeva, A., Slunyaev, A. Rogue waves, rogue events and extreme wave kinematics in spatio-temporal fields of simulated sea states // Nat. Hazards and Earth Syst. Sci. 2013. Vol. 13. P. 1759-1771.
[15] Slunyaev, A., Sergeeva, A., Didenkulova, I. Rogue events in spatiotemporal numerical simulations of unidirectional waves in basins of different depth // Nat. Hazards. 2016. Vol. 84, No. 2. P. 549-565.
[16] Annenkov, S.Y., Shrira, V.I. Spectral evolution of weakly nonlinear random waves: kinetic description versus direct numerical simulations //J. of Fluid Mech. 2018. Vol. 844. P. 766-795.
[17] Chalikov, D. Numerical modeling of surface wave development under the action of wind // Ocean Sci. 2018. Vol. 14. P. 453-470.
[18] Pelinovsky, E., Shurgalina, E. KDV soliton gas: interactions and turbulence. Challenges in Complexity: Dynamics, Patterns, Cognition / I. Aronson, A. Pikovsky, N. Rulkov, L. Tsimring (Eds) // Nonlinear Systems and Complexit. 2017. Vol. 20. P. 295-306.
[19] Шургалина Е.Г., Пелиновский Е.Н., Горшков К.А. Эффект отрицательной скорости частиц в солитонном газе в рамках уравнений типа Кортевега — де Вриза // Вестн. МГУ. Сер. 3. Физика и астрономия. 2017. № 5. С. 10-16.
Shurgalina, E., Pelinovsky, E., Gorshkov, K. The effect of a negative particle velocity in a soliton gas within the Korteweg — de Vries-type equations // Moscow Univ. Physics Bull. 2017. Vol. 72(5). P. 441-448.
[20] ШШургалина Е.Г. Механизм образования волн-убийц в результате взаимодействия соли-тонов внутренних волн в стратифицированном водоеме // Изв. РАН. Механика жидкости и газа. 2018. Т. 1. С. 61-67.
Shurgalina, E.G. The mechanism of the formation of freak waves in the result of interaction of internal waves in stratified basin // Fluid Dynamics. 2018. Vol. 53, No. 1. P. 59-64.
[21] Pelinovsky, E.N., Shurgalina, E.G., Sergeeva, A.V. et al. Two-soliton interaction as an elementary act of soliton turbulence in integrable systems // Phys. Letters A. 2013. Vol. 377, No. 3-4. P. 272-275.
[22] Anco, S.C., Ngatat, N.T., Willoughby, M. Interaction properties of complex mKdV solitons // Physica D. 2011. Vol. 240. P. 1378-1394.
[23] Пелиновский Е.Н., ШШургалина Е.Г. Двухсолитонное взаимодействие в рамках модифицированного уравнения Кортевега —де Вриза // Изв. вузов. Радиофизика. 2014. Т. 57, № 10. С. 825-833.
Pelinovsky, E.N., Shurgalina, E.G. Two-soliton interaction in the frameworks of modified Korteweg —de Vries equation // Radiophysics and Quantum Electronics. 2015. Vol. 57, No. 1. P. 737-744.
[24] Шургалина Е.Г. Особенности двухсолитонного взаимодействия в рамках уравнения Гарднера // Изв. вузов. Радиофизика. 2017. Т. 60, № 9. С. 787-792.
Shurgalina, E.G. The features of the paired soliton interactions within the framework of the Gardner equation // Radiophysics and Quantum Electronics. 2018. Vol. 60, No. 9. P. 703-708.
[25] ШШургалина Е.Г. Механизм образования волн-убийц в результате взаимодействия солитонов внутренних волн в стратифицированном водоеме // Изв. РАН. Механика жидкости и газа. 2018. Т. 1. С. 61-67.
Shurgalina, E.G. The mechanism of the formation of freak waves in the result of interaction of internal waves in stratified basin // Fluid Dynamics. 2018. Vol. 53, No. 1. P. 59-64.
[26] Slunyaev, A. Nonlinear analysis and simulations of measured freak wave time series // Europ. J. of Mechanics B. Fluids. 2006. Vol. 25. P. 621-635.
[27] Wang, L., Yang, C., Wang, J. He, J. The height of an nth-order fundamental rogue wave for the nonlinear Schrodinger equation // Phys. Letters A. 2017. Vol. 381. P. 1714-1718.
[28] Slunyaev, A.V., Pelinovsky, E.N. The role of multiple soliton interactions in generation of rogue waves: the mKdV framework // Phys. Review Letters. 2016. Vol. 117. P. 214501(1-5).
[29] Slunyaev, A. On the optimal focusing of solitons and breathers in long wave models // Submitted to Studies in Appl. Mathematics. 2018. P. 1-29.
[30] Taha, R.T., Ablowittz, M.J. Analytical and numerical aspects of certain nonlinear evolution equations. I. Analytical // J. of Comput. Physics. 1984. Vol. 55, iss. 2. P. 192-202.
[31] Taha, R.T., Ablowittz, M.J. Analytical and numerical aspects of certain nonlinear evolution equations. III. Numerical, Korteweg —de Vries equation // J. of Comput. Physics. 1984. Vol. 55, iss. 2. P. 231-253.
[32] Taha, R.T., Ablowittz, M.J. Analytical and numerical aspects of certain nonlinear evolution equations. IV. Numerical, Modified Korteweg —de Vries equation // J. of Comput. Physics. 1988. Vol. 77, iss. 2. P. 540-548.
[33] Fronberg, B. A practical Guide to pseudospectral methods. Cambrige: Cambrige Univ. Press, 1998. 231 p.
[34] Chan, T.F., Kerkhoven, T. Fourier methods with extended stability intervals for the Korteweg —de Vries equation // SIAM J. of Numer. Analysis. 1985. Vol. 22, No. 3. P. 441-454.
[35] Shurgalina, E.G., Pelinovsky, E.N. Nonlinear dynamics of a soliton gas: modified Korteweg —de Vries equation framework // Phys. Letters A. 2016. Vol. 380, No. 24. P. 20492053.
[36] Dutykh, D., Pelinovsky, E. Numerical simulation of a solitonic gas in KdV and KdV-BBM equations // Phys. Letters A. 2014. Vol. 378, No. 42. P. 3102-3110.
[37] El, G.A., Krylov, A.L., Molchanov, S.A., Venakides, S. Soliton turbulence as a thermodynamic limit of stochastic soliton lattices // Physica D: Nonlinear Phenomena. 2001. Vol. 152. P. 653-664.
[38] Agafontsev, D.S., Zakharov, V.E. Integrable turbulence and formation of rogue waves // Nonlinearity. 2015. Vol. 28, No. 8. P. 2791.
[39] Sun, Y.-H. Soliton synchronization in the focusing nonlinear Schrodinger equation // Phys. Review E. 2016. Vol. 93. P. 052222.
[40] Иванычев Д.Н., Фрайман Г.М. Образование солитонных пар в нелинейных средах с малой диссипацией // Теор. и матем. физика. 1997. Т. 110. С. 254-271.
Ivanychev, D.N., Fraiman, G.M. Generation of soliton pairs in nonlinear media with weak dissipation // Theor. and Math. Physics. 1997. Vol. 110, No. 2. P. 199-213.
[41] Шургалина Е.Г., Пелиновский Е.Н. Динамика ансамблей нерегулярных волн в прибрежной зоне. Нижний Новгород: НГТУ, 2015. 179 с.
Shurgalina, E.G., Pelinovskiy, E.N. Dynamics of irregular wave ensembles in the coastal zone. Nizhny Novgorod: NGTU, 2015. 179 p.
[42] Lo, E., Mei, C.C. A numerical study of water-wave modulation based on a higher-order nonlinear Schrodinger equation // J. of Fluid Mech. 1985. Vol. 150. P. 395-416.
Поступила в 'редакцию 25 декабря 2018 г.
Numerical simulation of soliton gas within the Korteweg — de Vries type equations
Didenkulova, Ekaterina G.*, Kokorina, Anna V., Slyunyaev, Alexey V.
Institute of Applied Physics of the RAS, Nizhny Novgorod, 603950, Russia * Corresponding author: Didenkulova, Ekaterina G., e-mail: eshurgalina@mail.ru
The details of the numerical scheme and the method of specifying the initial conditions for the simulation of the irregular dynamics of soliton ensembles within the framework of equations of the Korteweg — de Vries type are given using the example of the modified Korteweg — de Vries equation with a focusing type of nonlinearity. The numerical algorithm is based on a pseudo-spectral method with implicit integration over time and uses the Crank-Nicholson scheme for improving the stability property. The aims of the research are to determine the relationship between the spectral composition of the waves (the Fourier spectrum or the spectrum of the associated scattering problem) and their probabilistic properties, to describe transient processes and the equilibrium states.
The paper gives a qualitative description of the evolution of statistical characteristics for ensembles of solitons of the same and different polarities, obtained as a result of numerical simulations; the probability distributions for wave amplitudes are also provided. The results of test experiments on the collision of a large number of solitons are discussed: the choice of optimal conditions and the manifestation of numerical artifacts caused by insufficient accuracy of the discretization. The numerical scheme used turned out to be extremely suitable for the class of the problems studied, since it ensures good accuracy in describing collisions of solitons with a short computation time.
Keywords: soliton gas, Korteweg — de Vries equation, modified Korteweg — de Vries equation, Gardner equation, numerical simulation.
Cite : Didenkulova, E.G., Kokorina, A.V., Slunyaev, A.V. Numerical simulation of soliton gas within the Korteweg — de Vries type equations // Computational Technologies. 2019. Vol. 24, No. 2. P. 52-66. (In Russ.) DOI: 10.25743/ICT.2019.24.2.005.
Acknowledgements. This research was partly supported by RFBR grant No. 16-3260012 (Section 2); and also by RFBR grant No. 18-02-00042, and by the Fundamental Research Programme of RAS "Nonlinear Dynamics: Fundamental problems and applications" (Sections 1 and 3).
Received December 25, 2018
© ICT SB RAS, 2019
