Теоретическое исследование трёхмерных диссипативных топологических солитонов в лазерах. тема диссертации и автореферата по ВАК РФ 01.04.05, кандидат наук Веретенов Николай Александрович

• Веретенов Николай Александрович
• кандидат наук
• 2019

Оглавление диссертации кандидат наук Веретенов Николай Александрович

Введение

Глава 1. Комплексы лазерных пуль

1.1 Физическая модель

1.2 Стационарная сферически симметричная задача

1.3 Одиночные лазерные пули

1.4 Слабая связь лазерных пуль

Глава 2. Солитоны в интерферометре

2.1 Тёмные одномерные солитоны в интерферометре с дефокусирующей нелинейностью

2.2 Двумерные светлые солитоны

2.2.1 Физическая модель

2.2.2 Одиночные и связанные солитоны

2.2.3 Движение солитонных комплексов

Глава 3. Топологические солитоны в лазерах

3.1 Физическая модель

3.2 Измеряемые величины

3.3 Вихревой тор, "прецессон"

3.4 Солитоны с несколькими вихревыми линиями

3.4.1 "Яблоко"

3.4.2 Анализ траекторий вектора Пойнтинга для симметричного "яблока"

3.4.3 Простой двойной тор

3.4.4 Тривиальный узел (unknot)

3.4.5 Трилистник (trefoil)

3.4.6 Зацепление Соломона

3.4.7 Распад зацепления Соломона

3.4.8 Зацепление Хопфа

Стр.

Заключение

Словарь терминов

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

Список рисунков

Приложение А. Используемые численные методы и программы

А.1 От физической модели к компьютерной программе

А.2 Численное интегрирование многомерной задачи, программа Lbecon 113 А.3 Решение краевой стационарной задачи методом Ньютона

Введение диссертации (часть автореферата) на тему «Теоретическое исследование трёхмерных диссипативных топологических солитонов в лазерах.»

Введение

Локализованные состояния в различных нелинейных системах привлекают внимание как теоретиков, так и экспериментаторов. Во-первых, своей универсальностью, благодаря которой структуры совершенно различной природы могут быть очень похожими. Во-вторых, своей устойчивостью к шумам и флуктуаци-ям параметров схемы. Любая устойчивая локализованная структура — это запись информации на некотором однородном фоне. Её можно считать, можно стереть, а при определённых условиях она может исчезнуть или существенно измениться после взаимодействия с другой подобной структурой. Благодаря этим свойствам можно рассчитывать на их применение в обработке, передаче и хранении информации. К таким структурам относятся солитоны, у которых можно выделить некоторый центр, вдали от которого состояние системы не будет отличаться от однородного [1]. Солитоны также могут быть различной природы и размерности, но их можно разделить на два класса: консервативные и диссипативные. Первые существуют в системах без притока энергии, вторые — в системах с потерями и усилением. Классический консервативный солитон описывается нелинейным уравнением Шрёдингера (НУШ). Для существования солитонов НУШ необходим баланс между дисперсией среды, уширяющей импульс, и компенсирующей её нелинейностью (для временных солитонов) или аналогично между дифракцией и нелинейностью (для пространственных).

В настоящей работе исследуются диссипативные солитоны в широко-апертурных лазерах класса с безинерционно насыщающимися поглощением и усилением и бистабильных интерферометрах с когерентной накачкой, заполненных средой с керровской нелинейностью. Диссипативные структуры являются сильными аттракторами и отличаются большой устойчивостью к возмущениям. Лазерная схема описывается обобщённым комплексным уравнением Гинзбурга-Ландау. Для описания интерферометра подходит уравнение типа НУШ с двумя пространственными измерениями. Оба уравнения описывают медленную динамику огибающей электрического поля. Существование трёхмерных консервативных стационарных солитонов в однородной среде с насыщающейся нелинейностью было предсказано в работе [2]. В оптике трёхмерные локализованные в пространстве и времени импульсы известны как "лазерные пули" [3].

Работа является продолжением исследований диссипативных структур, которые были найдены в модели с двумя пространственными измерениями. Знания свойств двумерных солитонов использовались для построения трёхмерных фундаментальных (без фазовых дислокаций) и топологических солитонов. Многие из предсказанных ранее двумерных решений обнаружены экспериментально, например, солитоны с дислокациями фазы в поверхностно-излучающих лазерах с вертикальным резонатором (VCSEL). Трёхмерные топологические солитоны, в том числе, достаточно сложные, содержащие несколько вихревых линий, исследуются численно в различных моделях достаточно давно [4; 5], причём во многих случаях их устойчивость обеспечивается удерживающим потенциалом или нелокальностью отклика среды (см., например, обзор [6]), неоднородной дефокусирующей нелинейностью, стабилизирующей трёхмерные солитоны с топологическим зарядом [7-9] либо периодической модуляцией коэффициента преломления и коэффициентов усиления/потерь [10]. В данной работе рассматривается случай однородного и изотропного пространства (либо со слабой анизотропией) и устойчивость солитонов не связана с дополнительными стабилизирующими факторами.

Один из главных результатов работы — обнаружение устойчивых заузлен-ных лазерных солитонов со сложной внутренней структурой. Такие солитоны обладают набором вихревых линий, на которых интенсивность излучения обращается в нуль, а фаза при обходе линии по замкнутому контуру изменяется на величину 2пт, где целое число т называют топологическим зарядом. Благодаря перекрывающимся областям устойчивости у разных типов заузленных солито-нов, их можно использовать в качестве элементов оптической памяти с большим числом состояний (битностью). Различные типы солитонов — движущиеся, взаимодействующие и превращающиеся друг в друга, — образуют внутри оптической системы своего рода "химию", с помощью которой можно записывать и обрабатывать информацию. Непрекращающиеся исследования узловых конфигураций в физике [11-13], химии и биологии [14-16] заставляют задуматься о возможных приложениях найденных солитонов в других областях, не только в оптике.

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

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

1. Создать комплекс программ для эффективного решения дифференциальных уравнений и визуализации данных.

2. Найти параметры оптической схемы, при которых возможно существование новых типов солитонов.

3. Построить начальное распределение, способное релаксировать к устойчивому солитону при этих параметрах.

4. Всесторонне исследовать найденные решения, используя численные и полуаналитические методы.

Научная новизна:

1. Впервые в численном эксперименте были получены устойчивые комплексы лазерных пуль и смоделированы их столкновения.

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

3. Впервые в численном эксперименте был получен и исследован новый класс лазерных заузленных солитонов.

Практическая значимость: обработка и хранение информации оптическими методами

Mетодология и методы исследования. Большинство представленных результатов получены в численных экспериментах. В случае лазерной схемы моделировался лазер с безинерционно насыщающимися поглощением и усилением и нулевыми частотными расстройками. Использовалось приближение медленно меняющейся огибающей электрического поля, для описания динамики которой подходит обобщённое комплексное уравнение Гинзбурга-Ландау. Для прямого численного интегрирования нелинейных дифференциальных уравнений в частных производных использовались метод конечных разностей, метод расщепления, а для нахождения стационарных решений — метод релаксации и метод стрельбы. Для получения асимптотик использовались аналитические методы решения линеаризованных дифференциальных уравнений. Для проверки устойчивости стационарных симметричных решений использовался полуаналитический метод: численное решение линеаризованной возмущённой задачи.

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

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

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

Достоверность результатов, полученных численными методами, предварительно подтверждается проверками, теоретическими оценками и расчётом с использованием альтернативной численной схемы. Двумерные диссипативные пространственные солитоны, в том числе с дислокациями фазы, предсказанные ранее (см., например [17-19]), уже достаточно давно наблюдаются в экспериментах [20-22]. Трёхмерные полностью локализованные решения такого типа на сегодняшний день исследуются численно, хотя известны успешные эксперименты по локализации света в трёх измерениях, например, в работе [23] наблюдались довольно необычные структуры, в виде цепочек из пространственно-светлых временно-тёмных солитонов.

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

Как уже отмечалось выше, значительная часть работы является дальнейшим развитием исследований двумерных диссипативных структур, предсказанных в ГОИ ранее. Численные методы, схемы и программы, используемые в данной работе, достаточно точно воспроизводят эти результаты в случае одного и двух пространственных измерений. При исследовании солитонов в ин-

терферометре с когерентной накачкой подтвердились результаты, полученные ранее другими авторами, а именно, область устойчивости одиночного светлого солитона [24]. В случае трёх пространственных измерений применяемый автором метод использовался для моделирования процессов разделения фаз и образования локализованных структур в вырожденном оптическом параметрическом генераторе, описываемом системой уравнений Свифта-Хоенберга [25; 26]. Этому исследованию предшествовало изучение двумерного случая, поэтому вначале были подтверждены имеющиеся в литературе результаты. Для случая трёх измерений многие устойчивые решения обладали симметрией, поэтому имелась возможность сравнения результатов моделирования динамической задачи с решениями стационарной краевой задачи, полученными методом Ньютона. Хорошее совпадение результатов, полученных принципиально разными методами, подтверждает точность интегрирования управляющего уравнения. Метод успешно применялся для исследования периодических структур с различной симметрией, возникающих в метаматериалах с положительной и отрицательной дисперсией [27; 28]. Кроме оптических схем моделировался также атомный конденсат Бозе-Эйнштейна [29].

Для некоторых стационарных структур, обладающих симметрией, можно с высокой точностью решать краевую задачу методом релаксации. Относительную ошибку (невязку) можно значительно уменьшить, если использовать в программе квадратичную точность либо прибегнуть к помощи специализированных пакетов математических программ, способных считать с произвольной точностью. Поэтому можно оценить точность решения, полученного прямым интегрированием, сравнив его с результатом работы метода Ньютона. Также можно сравнивать области устойчивости, полученные двумя способами: длительными расчётами, добавлением шума в динамической задаче и анализом спектра линеаризованной задачи для возмущённого решения стационарной задачи. В некоторых случаях мы можем сравнивать не только границы устойчивости, но и пространственную моду возмущения, разрушающего симметричное решение [30].

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

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

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

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

Личный вклад. Автор разрабатывал программы и проводил численные эксперименты.

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

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

В численных экспериментах обнаружен и исследован новый класс трёхмерных солитонов в виде вихревых заузленных солитонов с различной топологией.

Самым простым из них является солитон с одной вихревой линией, а самым сложным — солитон с двумя заузленными замкнутыми вихревыми линиями, внутри которых расположены две незамкнутые вихревые линии, уходящие на периферию.

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

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

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

Для интерферометра с накачкой: построены комплексы светлых солитонов с различной симметрией и продемонстрированы различные типы движения — прямолинейное, вращение и "лунное" движение.

Публикации. Основные результаты по теме диссертации изложены в 13 печатных изданиях, 10 из которых изданы в журналах, в рецензируемых журналах (рекомендованных ВАК или входящих в базу Scopus и Web of science), 3 — в тезисах докладов.

Объем и структура работы. Диссертация состоит из введения, трёх глав, заключения и приложения. Полный объём диссертации составляет 123 страницы, включая 75 рисунков. Список литературы содержит 72 наименования.

Глава 1. Комплексы лазерных пуль

1.1 Физическая модель

т-» • • О ЧУ

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

дЕ ]Лд2Е . ]Лд2Е ... _|2._ „ _

^ = (г + йх) ^ + + йу ) ^ + + ^ ^ + 1(\Е \)Е- С1-1-1)

Нелинейная функция /(\Е\2) записывается следующим образом:

1 (Е 1 ) = -1 + (ТТЖЖ+Щу - (1 + гд„)(1 +\ Е\2)' (1-1-2)

где ^ — продольная координата, х,у — поперечные координаты, т = £ — ^/УдГ — время в системе отчета, движущейся с групповой скоростью УдГ, йх,йу,йт — коэффициенты диффузии, а0 и д0 — коэффициенты линейного

т 1 дат ^

поглощения и усиления, и 1д = --отношение интенсивностей насыщения

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

й\\, = й\\. Будем считать фиксированными коэффициенты: 1д = 10 и ао = 2. Параметр д0 будет нашим бифуркационным параметром. Рассмотрим стационарный однородный режим, для которого из уравнения 1.1.1 следует

— 1 +-^---ОЦ. = 0 (1.1.3)

1 + \ Е \ 2/1д 1 + \ Е \2 К)

или

I? + (1 + ао + 1д(1 - 9о))1а - 1д/0 = 0,/0 = -1 - ао + до. (1.1.4)

Зависимость 1е(д0) при 1д = 10, а0 = 2 показана на рисунке 1.1. В интервале от дтп & 1.9486 до дтах = 3. одному значению д0 соответствует два возможных устойчивых состояния: безгенерационный режим I = 0 и режим с ненулевой интенсивностью поля, а также одно неустойчивое состояние (пунктирная линия). Левее диапазона бистабильности при любой начальной интенсивности поля со временем установится безгенерационный режим, а при д0 > дтах безгенерационный режим неустойчив, как и любое локализованное решение с нулевой асимптотикой. Солитоны могут существовать только в диапазоне бистабильно-сти.

20 15

4 10

5 0

Рисунок 1.1 — Зависимость интенсивности однородного стационарного режима от д0 при 1д = 10, а0 = 2. Сплошные линии — устойчивые ветви, пунктирная —

неустойчивая ветвь.

Трёхмерная задача является обобщением двумерной задачи (подробно исследованной, например, в [30] и [19]) на случай трёх измерений. Дополнительное измерение даёт нам значительно больше возможностей для строительства

§0

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

Рассматриваемая модель подходит для описания лазера с кольцевым резонатором, заполненным нелинейной средой с указанными свойствами.

1.2 Стационарная сферически симметричная задача

Независящее от эволюционной переменной решение можно получить, решив краевую задачу с правильными граничными условиями. Запишем уравнение 1.1.1 в сферической системе координат. В комплексной амплитуде выделим функцию, зависящую только от радиальной координаты т, и периодическую зависимость от г: Е(т,г) = А(т)е-гуг, где V — нелинейный сдвиг волнового числа. Симметричное решение не зависит от ф и 6. Для решения стационарной задачи, приравняем производную по г к нулю. Сдвиг волнового числа V заранее неизвестен, он будет одним из неизвестных параметров при решении методом стрельбы.

В итоге мы получаем следующее одномерное уравнение для функции А(г):

^ + + гЪ+/ (А(г)12)]А(г) = 0. а*«

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

убывает А(г) ^ 0,т ^ о. Введём комплексную переменную р = т* выбрав ветвь Яв(р) > 0. Тогда уравнение 1.2.1 запишется в виде:

^ + /0

г + 1

д2А(р) 2дА(р) л, , ,

-¡—т+р-^г+А(р)=0' (1'2'2)

Требуемая асимптотика определяется следующим выражением:

С

А = — ехр(-р). (1.2.3)

Поскольку уравнение 1.2.1 инвариантно относительно фазы, можно считать константу С вещественной.

Из-за того, что в операторе Лапласа присутствует особенность 1/г, мы будем искать решение на интервале (х0,хтах) с достаточно малым х0 & 0.001. После того, как решение будет найдено, это значение можно уменьшить на несколько порядков, до 10—8 и менее. Комплексное дифференциальное уравнение второго

порядка 1.2.1 преобразуется к системе трех уравнений первого порядка для трех

1 да (г) д ф(г)

переменных : а(г) = \ А(г) \, р(г) = —— _ , д(г) = —. Благодаря тому,

а(г) дг дг

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

В итоге получаем три вещественных уравнения на функции а(а), р(г) и д(г) и одно дополнительное для фазы:

— = — р + д

йр йг йд

йг йа

йг йф

йг

2_ 2р_ 1 (а2)й + V

—2рд — 2- — г

г й2 + 1 д уй — 1 (а2) ' й2 + 1 '

(1.2.4)

= ар,

= д.

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

XI

х2

Х3 х4

А \2, дА 2

дг — ( \ А \ 2 дг

(1.2.5)

= -г А

( дг^)

А

дА г

Переменные связаны соотношением: 4x1 х2 = + х2. В этом случае получаем следующую систему уравнений:

dx

i

j = xs, dr

dr2 =--x2 + ^ixs - Vx4, (1.2.6)

dr r

dxs _ 2

—— = 2Vx + 2x2 - xs-, dr r

dx4 2

—— = -2^2xi - x±~. dr r

/(|A(r)|2)d - v /(|A(r)|2) + vd Здесь введены обозначения Vi = —— 70-, V2 = —— 70-.

1 + d2 1 + d2 Для симметричного солитона должно выполняться xi(rmin) = x2(rmin) =

xs(rmin) = x4(rmin) = 0. Эти задачи решились методом стрельбы с учётом асимптотик на бесконечности и условий в нуле. Суть метода стрельбы — серия выстрелов с разными значениями параметров от одной границы до другой, либо стрельба с крайних точек и попытки сшивания решения в некой промежуточной точке. Практически любой итерационный метод на больших расстояниях приводит к накоплению ошибок, поэтому второй способ предпочтительнее. Будет ли стрельба успешной зависит от того, какие именно параметры менялись при стрельбе и на какие переменные были наложены условия.

В своё время был использован хитроумный метод поиска решений способом стрельбы из бесконечности в ноль, реализованный на Фортран 77 (32 бит под DOS extender). Стартуя из достаточного большого r, там где интенсивность мала, мы движемся к нулю. Дойдя до точки r0, получаем значения a(r0), p(r0) и q(r0). Из соображений симметрии ясно, что для солитона должны выполняться следующие условия: p(r0) = q(r0) = 0, a(r0) > 0. В асимптотике присутствует свободный параметр C, кроме того неизвестен сдвиг волнового числа v. Алгоритм поиска нужного решения следующий. Берём достаточно большую область этих параметрам, сканируем её с большим шагом, отмечаем на ней границы областей, на которых p(r0) и q(r0) меняют знак. Точки их пересечения могут соответствовать искомым параметрам. Когда-то эти построения производились вручную на листе миллиметровки. Разумеется, для быстрого поиска решений такой метод не годится.

Рассмотрим програмные реализации метода стрельбы. На языке Wolfram Mathematica достаточно воспользоваться несколькими функциями для решения однородных дифференциальных уравнений (ОДУ) и поиска корней системы

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

Т Т VJ VJ W

Но основной проблемой, стоящей на пути к точному решению граничной задачи, является нахождение начального приближения для нескольких параметров, заранее неизвестных. Аналитические приближения тут практически бессильны, так как из асимптотик мы не можем найти значения свободных параметров, а именно — сдвига волнового числа и амплитуды поля на бесконечности. Если сдвиг предполагается малым, то амплитуда зависит от выбора точки rmax, откуда производится выстрел. Выбор этой точки зависит от точности метода интегрирования, чем она выше, тем дальше можно отойти от нуля, туда, где асимптотическое приближение 1.2.3 будет точнее. Для метода Рунге-Кутта 4-го порядка с контролем шага (5-й порядок точности) эта точка выбиралась в диапазоне 20..50. Уравнение, записанное в цилиндрической системе координат, имеет особенность в нуле, поэтому интегрировать его до нуля нельзя. Необходимо остановиться в некоторой точке rmin ^ 1. Желательно подобрать минимальное значения rmin, на практике оно должно быть много меньше величины шага Ar.

Интерактивная программа, написанная специально для решения задач такого типа, позволяла искать решения методом стрельбы, уточняя параметры с помощью простейшего графического интерфейса, затем найденные вручную параметры уточнялись методом Ньютона. Плоскость, на которой строились кривые p(rmin;pari,par2) = 0 и q(rmin;pari,par2) = 0, содержит особые точки, там где функция q терпит разрыв 2-го рода. Параметры pari и par2 это нелинейный сдвиг волнового числа v и коэффициент C в асимптотике 1.2.3. Теоретически, проблему можно было решить, устранив особенности в уравнении либо усложнить стрельбу: например, интегрировать из нуля и бесконечности навстречу друг другу и сшивать два решения в промежуточной точке. Но поскольку в программе присутствовала графика, можно было увеличивать область на плоскости параметров, постепенно выявляя область параметров, в которой обе функции p и q меняли знак. Уточнение корней после этого производилось методом Ньютона. Минус такого подхода — трудоёмкость поиска, плюсы в том, что на плоскости были видны все возможные решения, отвечающие необходимой симметрии, асимптотикам на бесконечности и условиям p(0) = 0, q(0) = 0 в центре структуры.

Список литературы диссертационного исследования кандидат наук Веретенов Николай Александрович, 2019 год

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

1. Dissipative Solitons / Ed. by Nail Akhmediev, Adrian Ankiewicz. — Springer Berlin Heidelberg, 2005. — URL: https://doi.org/10.1007%2Fb11728.

2. Vakhitov N.G., Kolokolov Aleksandr A. Stationary solutions of the wave equation in a medium with nonlinearity saturation // Radiophysics and Quantum Electronics. — 1973. — Vol. 16, no. 7. — Pp. 783-789.

3. Silberberg Yaron. Collapse of optical pulses // Optics letters. — 1990. — Vol. 15, no. 22. — Pp. 1282-1284.

4. Three-dimensional spatiotemporal optical solitons in nonlocal nonlinear media / Dumitru Mihalache, D Mazilu, F Lederer et al. // Physical Review E. — 2006. — Vol. 73, no. 2. — P. 025601.

5. Burlak Gennadiy, Malomed Boris A. Interactions of three-dimensional solitons in the cubic-quintic model // Chaos: An Interdisciplinary Journal of Nonlinear Science. — 2018. — Vol. 28, no. 6. — P. 063121.

6. Malomed Boris A. Multidimensional solitons: Well-established results and novel findings // The European Physical Journal Special Topics. — 2016. — Vol. 225, no. 13-14. — Pp. 2507-2532.

7. Three-dimensional hybrid vortex solitons / Rodislav Driben, Yaroslav V Kartashov, Boris A Malomed et al. // New Journal of Physics. — 2014. — Vol. 16, no. 6. — P. 063035.

8. Soliton gyroscopes in media with spatially growing repulsive nonlinearity / Rodislav Driben, Yaroslav V Kartashov, Boris A Malomed et al. // Physical review letters. — 2014. — Vol. 112, no. 2. — P. 020404.

9. Twisted toroidal vortex solitons in inhomogeneous media with repulsive nonlinear-ity / Yaroslav V Kartashov, Boris A Malomed, Yasha Shnir, Lluis Torner // Physical review letters. — 2014. — Vol. 113, no. 26. — P. 264101.

10. Three-dimensional topological solitons in PT-symmetric optical lattices / Yaroslav V Kartashov, Chao Hang, Guoxiang Huang, Lluis Torner // Optica. — 2016. — Vol. 3, no. 10. — Pp. 1048-1055.

11. Raymer D. M., Smith D. E. Spontaneous knotting of an agitated string // Proceedings of the National Academy of Sciences. — 2007. — oct. — Vol. 104, no. 42. — Pp. 16432-16437. — URL: https://doi.org/10.1073%2Fpnas.0611320104.

12. Spontaneous knotting of self-trapped waves / Anton S. Desyatnikov, Daniel Buc-coliero, Mark R. Dennis, Yuri S. Kivshar // Scientific Reports. — 2012. — oct. — Vol. 2, no. 1. — URL: https://doi.org/10.1038%2Fsrep00771.

13. Tying quantum knots / D. S. Hall, M. W. Ray, K. Tiurev et al. // Nature Physics.

— 2016. — jan. — Vol. 12, no. 5. — Pp. 478-483. — URL: https://doi.org/10. 1038%2Fnphys3624.

14. Molecular Trefoil Knot from a Trimeric Circular Helicate / Liang Zhang, David P August, Jiankang Zhong et al. // Journal of the American Chemical Society. — 2018. — Vol. 140, no. 15. — Pp. 4982-4985.

15. Fielden Stephen D. P., Leigh David A., Woltering Steffen L. Molecular Knots // Angewandte Chemie International Edition. — 2017. — aug. — Vol. 56, no. 37. — Pp. 11166-11194. — URL: https://doi.org/10.1002%2Fanie.201702531.

16. Translation and untying of DNA knots in extensional fields / Alexander Klotz, Beatrice Soh, Vivek Narsimhan, Patrick Doyle // Bulletin of the American Physical Society. — 2018.

17. Rozanov N.N., Khodova G.V. Autosolitons in bistable interferometers // Optics and Spectroscopy. — 1988. — Vol. 65. — Pp. 449-450.

18. Розанов Н.Н., Федоров С.В. Дифракционные волны переключения и автосо-литоны в лазере с насыщающимся поглощением // Оптика и спектроскопия.

— 1992. — Т. 72, № 6. — С. 1394-1399.

19. Розанов КН., Федоров С.В., Шацев А.Н. Структура энергетических потоков и её бифуркации для двумерных лазерных солитонов // Журнал экспериментальной и теоретической физики. — 2004. — Т. 125, № 3. — С. 486-498.

20. Observation of laser vortex solitons in a self-focusing semiconductor laser / J Jimenez, Y Noblet, P V Paulau et al. // Journal of Optics. — 2013. — Vol. 15, no. 4. — P. 044011. — URL: http://stacks.iop.org/2040-8986/15/i=4/a=044011.

21. Mihalache Dumitru. Localized optical structures: An overview of recent theoretical and experimental developments // Proceedings of the Romanian Academy Series A

— Mathematics Physics Technical Sciences Information Science. — 2015. — 01. — Vol. 16. — Pp. 62-69.

22. Barbay S., Kuszelewicz R., Tredicce J. R. Cavity Solitons in VCSEL Devices // Advances in Optical Technologies. — 2011. — Vol. 2011. — Pp. 1-23. — URL: https://doi.org/10.1155%2F2011%2F628761.

23. Three-Dimensional Spatiotemporal Pulse-Train Solitons / Oren Lahav, Ofer Kfir, Pavel Sidorenko et al. // Phys. Rev. X. — 2017. — Nov. — Vol. 7. — P. 041051.

— URL: https://link.aps.org/doi/10.1103/PhysRevX.7.041051.

24. Nesterov L.A, Veretenov N.A., Rosanov N.N. Quantum fluctuations of one-and two-dimensional spatial dissipative solitons in a nonlinear interferometer: II. Two-dimensional light solitons // Optics and Spectroscopy. — 2015. — Vol. 118, no. 5.

— Pp. 794-802.

25. Veretenov Nikolay, Tlidi Mustapha. Dissipative light bullets in an optical parametric oscillator // Physical Review A. — 2009. — Vol. 80, no. 2. — P. 023822.

26. Veretenov Nikolay. Degenarated Optical Parametric Oscillator in 3D, Youtube. — URL: https://www.youtube.com/watch?v=2fW1YE0p_aE.

27. Three-dimensional structures in nonlinear cavities containing left-handed materials / Philippe Tassin, Guy Van der Sande, Nikolay Veretenov et al. // Optics express.

— 2006. — Vol. 14, no. 20. — Pp. 9338-9343.

28. Veretenov Nikolay. Modulational instability in 3D, Youtube. — URL: youtube.com/watch?v=-s2uRdRtlhY.

29. Interferometric precision measurements with Bose-Einstein condensate solitons formed by an optical lattice / N. Veretenov, Yu. Rozhdestvensky, N. Rosanov et al. // The European Physical Journal D. — 2007. — mar. — Vol. 42, no. 3. — Pp. 455460. — URL: https://doi.org/10.1140%2Fepjd%2Fe2007-00129-2.

30. Topologically multicharged and multihumped rotating solitons in wide-aperture lasers with a saturable absorber / Sergey V Fedorov, Nikolay N Rosanov, Ana-toly N Shatsev et al. // IEEE Journal of Quantum Electronics. — 2003. — Vol. 39, no. 2. — Pp. 197-205.

31. Розанов Н.Н. Диссипативные оптические солитоны от микро- к нано- и атто.

— Москва: Физматлит, 2011. — 526 с.

32. Conditions for the existence of laser bullets / NA Veretenov, AG Vladimirov, NA Kaliteevskii et al. // Optics and spectroscopy. — 2000. — Vol. 89, no. 3.

— Pp. 380-383.

33. Розанов Н.Н., Федоров С.В., Шацев А.Н. Движение комплексов слабо связанных двумерных лазерных солитонов // Журнал экспериментальной и теоретической физики. — 2006. — Т. 129, № 4. — С. 625-635.

34. Veretenov NA, Rosanov NN, Fedorov SV. Motion of complexes of 3D-laser solitons // Optical and Quantum Electronics. — 2008. — Vol. 40, no. 2-4. — Pp. 253-262.

35. Veretenov NA, Rozanov NN, Fedorov SV. Clusters of three-dimensional laser solitons and their collisions // Optics and Spectroscopy. — 2008. — Vol. 104, no. 4.

— Pp. 563-566.

36. Veretenov Nikolay. 7-soliton complex spiral moving, Youtube. — URL: https:// www.youtube.com/watch?v=t4io3ulnDPI.

37. Lugiato L. A, Lefever R. Spatial Dissipative Structures in Passive Optical Systems // Physical Review Letters. — 1987. — may. — Vol. 58, no. 21. — Pp. 2209-2211. — URL: https://doi.org/10.1103%2Fphysrevlett.58.2209.

38. Розанов Н.Н. Оптическая бистабильность и гистерезис в распределенных нелинейных системах. — Наука М., 1997. — 336 с.

39. Veretenov Nikolay. Dark Soliton Birth, Youtube. — URL: https://www.youtube. com/watch?v=hAhJUAk3yQo.

40. Nesterov L.A, Veretenov N.A., Rosanov N.N. Quantum fluctuations of one-and two-dimensional spatial dissipative solitons in a nonlinear interferometer: I. One-dimensional dark solitons // Optics and Spectroscopy. — 2015. — Vol. 118, no. 5.

— Pp. 781-793.

41. Ackemann Thorsten, Firth William J, Oppo Gian-Luca. Fundamentals and applications of spatial dissipative solitons in photonic devices // Advances in atomic, molecular, and optical physics. — 2009. — Vol. 57. — Pp. 323-421.

42. Розанов Н.Н. О поступательном и вращательном движении нелинейно-оптических структур в целом // Оптика и спектроскопия. — 2004. — Т. 96, № 3. — С. 451-454.

43. Colet P, Gomila D, Jacobo A, Matias MA. Dissipative Solitons: From Optics to Biology and Medicine (N. Akhmediev and A. Ankiewicz, eds.). — 2008.

44. Motion of dissipative soliton complexes in a nonlinear interferometer with driving radiation / NA Veretenov, NN Rosanov, SV Fedorov, AN Shatsev // Optics and Spectroscopy. — 2010. — Vol. 109, no. 5. — Pp. 753-759.

45. Розанов Н.Н., Семенов В.Е., Ходова Г.В. Поперечная структура поля в нелинейных бистабильных интерферометрах. I. Волны переключения и стационарные профили // Квантовая электроника. — 1982. — Т. 9, № 2. — С. 354-360.

46. Veretenov Nikolay. Soliton's party, Youtube. — URL: https://www.youtube.com/ watch?v=I0PBTPuVkec.

47. Розанов Н.Н. Феноменологические уравнения движения диссипативных оптических солитонов // Оптика и спектроскопия. — 2007. — Т. 102, № 5. — С. 800-804.

48. Квантовые флуктуации пространственных диссипативных солитонов в нелинейном интерферометре / Л.А. Нестеров, Ал. С. Киселев, Ан. С. Киселев, Н.Н. Розанов // Оптика и спектроскопия. — 2009. — Т. 106, № 4. — С. 639657.

49. Rosanov N.N., Fedorov S.V., Shatsev A.N. Curvilinear motion of multivortex lasersoliton complexes with strong and weak coupling // Physical review letters. — 2005. — Vol. 95, no. 5. — P. 053903.

50. Rosanov N.N., Fedorov S.V., Shatsev A.N. Complexes of in-phase two-dimensional laser solitons // Quantum Electronics. — 2008. —jan. — Vol. 38, no. 1. — Pp. 4145. — URL: https://doi.org/10.1070%2Fqe2008v038n01abeh013562.

51. Adams Colin Conrad. The knot book: an elementary introduction to the mathematical theory of knots. — American Mathematical Soc., 2004. — 323 pp.

52. Veretenov Nikolay. 3D soliton precession, Youtube. — URL: https://www.youtube. com/watch?v=qFfZNeUxAdc.

53. Veretenov NA, Rosanov NN, Fedorov SV. Rotating and Precessing Dissipative-Optical-Topological-3D Solitons // Physical review letters. — 2016. — Vol. 117, no. 18. — P. 183901.

54. Veretenov Nikolay. Chaotic behaviour of 3D vortex dissipative Soliton, Youtube.

— URL: https://www.youtube.com/watch?v=YuLonUFzb5c.

55. Rosanov N.N., Fedorov S.V., Shatsev A.N. Two-dimensional laser soliton complexes with weak, strong, and mixed coupling // Applied Physics B. — 2005.

— nov. — Vol. 81, no. 7. — Pp. 937-943. — URL: https://doi.org/10.1007% 2Fs00340-005-1981-4.

56. Veretenov NA, Fedorov SV, Rosanov NN. Topological vortex and knotted dissipative optical 3D solitons generated by 2D vortex solitons // Physical review letters. — 2017. — Vol. 119, no. 26. — P. 263901.

57. Veretenov NA, Fedorov SV, Rosanov NN. Topological three-dimensional dissipative optical solitons // Phil. Trans. R. Soc. A. — 2018. — Vol. 376, no. 2124. — P. 20170367.

58. Fedorov SV, Rosanov NN, Veretenov NA. Structure of Energy Fluxes in Topological Three-Dimensional Dissipative Solitons // JETP Letters. — 2018. — Vol. 107, no. 5. — Pp. 327-331.

59. Veretenov Nikolay. Strange twisted torus m=1, Youtube. — URL: https://www. youtube.com/watch?v=FNY2LAe7QY4.

60. Veretenov Nikolay. Long Apple vortex soliton, Youtube. — URL: https://www. youtube.com/watch?v=zOTpholulk8.

61. Качественная теория динамических систем второго порядка / А.А. Андронов, Е.А. Леонтович, И.И. Гордон, А.Г. Майер. — М.: Наука, 1966. — 568 с.

62. Veretenov Nikolay. Double vortex ring dissipative soliton with three unclosed vortex lines, Youtube. — URL: https://www.youtube.com/watch?v=vjG5aJEF8sc.

63. Veretenov Nikolay. Trivial knot — vortex dissipative soliton, Youtube. — URL: https://www.youtube.com/watch?v=84ofpQlKC0c.

64. Veretenov Nikolay. Trefoil knotted soliton, Youtube. — URL: https://www. youtube.com/watch?v=2QlVDnf87OE.

65. Veretenov Nikolay. Stable asymmetrical solomon knot, Youtube. — URL: https: //www.youtube.com/watch?v=KrUVdW9nCNk.

66. Murasugi Kunio. Knot theory and its applications. — Springer Science & Business Media, 2007. — 348 pp.

67. Veretenov Nikolay. Wild decay of an unstable Hopf link, Youtube. — URL: https: //www.youtube.com/watch?v=ZPwCK0EUb0A.

68. Veretenov Nikolay. Stable Hopf Soliton do exist!, Youtube. — URL: https://youtu. be/LhlBGFrWFl0.

69. Теория солитонов: метод обратной задачи / В.Е. Захаров, С.П. Новиков, Л.П. Питаевский, С.В. Манаков. — Наука, 1980. — 320 с.

70. Lenhard Johannes. Computer simulation: The cooperation between experimenting and modeling // Philosophy of Science. — 2007. — Vol. 74, no. 2. — Pp. 176-194.

71. Ортега Джеймс М. Введение в параллельные и векторные методы решения линейных систем: Пер. с англ. — Мир, 1991.

72. Veretenov Nikolay. Asymmetrical Solomon's Knot in OpenDX, Youtube. — URL: https://www.youtube.com/watch?v=tQJKICGgd_o.

Список рисунков

1.1 Зависимость интенсивности однородного стационарного режима от д0 при 1д = 10, а0 = 2. Сплошные линии — устойчивые ветви, пунктирная — неустойчивая ветвь.....................12

1.2 Бифуркационные кривые для одиночных сферически симметричных солитонов: зависимости нелинейного сдвига волнового числа от коэффициента усиления слабого сигнала для трёх типов лазерных пуль, профили амплитуды которых изображены на врезках. Жирная линия — устойчивые решения.......................18

1.3 Область устойчивости сферических лазерных пуль на плоскости частотных расстроек Дд,Да для д0 = 2.156, а0 = 2,1д = 10, d = 0. ... 19

1.4 Образование движущегося устойчивого треугольника из двух синфазных и одного противофазного солитонов. (а) — начальное состояние, образование слабой связей, ^ = 10; (Ь) - неустойчивая цепочка складывается в треугольник, ^ = 400, (с) — образование слабой связи между солитонами в основании треугольника, начало движения, ^ = 580, — положение структуры при ^ = 840...... 21

1.5 а) — зависимость угла в между осью пары и осью ^ для

d± = 0.08, d\\ = 0.04 (сплошная линия) и d± = 0.04, d\\ = 0.08

(пунктирная); б) — траектория и ориентация пирамиды;

до = 2.169, dL = 0.08, dy = 0.04.......................22

1.6 Пары синфазных (а) и противофазных (Ь) солитонов, тетраэдр (с) из четырёх синфазных солитонов , пирамида с противофазным солитоном в вершине и тремя синфазными в основании; направление движения показано стрелкой. Поверхности солитонов показаны для уровня интенсивности 0.001 от максимума, шкала фазы приведена

слева; до = 2.165, d = 0.06..........................22

1.7 Изоповерхности интенсивности и вихревые линии у пирамиды из трёх синфазных и одного противофазного солитонов. На правом рисунке вид сверху..............................23

1.8 Нецентральное столкновение двух пирамидальных комплексов с обменом боковыми солитонами в основании пирамид, д0 = 2.169,

б = 0.07....................................25

1.9 Столкновение с образованием двух движущихся треугольников (из двух синфазных и одного противофазного солитона) и двух одиночных солитонов.............................26

1.10 Столкновение с образованием двух противофазных и двух

синфазных пар................................27

1.11 Столкновение двух пирамид и рождение "мельницы"...........28

1.12 Икосаэдры, состоящие из 12-ти (слева) и 13-ти солитонов........29

1.13 Установление межсолитонного расстояния I в икосаэдре после неточной установки солитонов на свои места...............29

1.14 7-солитонный асимметричный комплекс и траектория движения его центра.....................................29

2.1 Схемы интерферометра. М — зеркала, Ж — нелинейная среда,

Еп — поле когерентной накачки......................31

2.2 Профили интенсивности (сплошная линия) и фазы фундаментального (узкого) солитона при = 2.796, 6 = —2.65, соответствующее значение амплитуды накачки Еп = 1.69........32

2.3 Профили интенсивности (сплошная линия) и фазы модифицированного солитона 15=3.044, 6 = —2.95, соответствующее значение амплитуды накачки Еп = 1.777.................. 33

2.4 Профили интенсивности солитонов: модифицированного (1), с двумя минимумами (2), с тремя (3) и с шестью (4), 1а = 3.044, 6 = —2.95. . . 33

2.5 Область устойчивости солитонов в зависимости от 13 и 6: 1 — внешние границы этой области и 2 — граница области существования модифицированного солитона. На вставке представлены области существования солитонов с одной (3), двумя (4) осцилляциями и линия Максвелла (5). Кривая 6 соответствует

верхней границе области бистабильности.................34

2.6 Профили интенсивности (слева) и фазы одиночного солитона......39

2.7 Профили интенсивности I и фазы Ф для двух пар связанных солитонов с различным расстоянием между центрами солитонов. ... 39

2.8 Зависимость логарифма скорости движения центра тормозящегося солитона от времени.............................41

2.9 Поперечное распределение интенсивности (шкала справа) для комплекса пяти солитонов с одной осью симметрии в два момента времени....................................42

2.10 Поперечное распределение интенсивности (шкала справа) для треугольного комплекса с одной осью симметрии и различающимися межсолитонными расстояниями ( и (2...................43

2.11 График координаты центра комплекса, изображённого на рисунке 2.10. 43

2.12 Распределения интенсивности для вращающегося комплекса из 8 солитонов с межсолитонным расстоянием ( с центральной симметрией в моменты времени 2000 (слева), 50000 (посередине) и 80000 (справа). ............................... 43

2.13 Распределения интенсивности для вращающегося комплекса из 6 солитонов с межсолитонным расстоянием ( и (2 с симметрией к повороту на угол 2п/3............................44

2.14 Распределения интенсивности комплекса из 7 солитонов

с межсолитонным расстроянием совершающим "лунное"

движение, отмечена траектория центра комплекса.............44

3.1 (а) Изоповерхность по уровню 0.27 от Imax твердотельного солитона в виде деформированного тора. Распределение интенсивности вращается с постоянной угловой Q = \Q\ = 0.092 (период Trot = 68). Угловая скорость, определяемая через момент импульса, Qs = \ = 0.0081 (период Ts = 772). Главные оси J 1,2,3 и момент M вращаются вокруг Q. Абсолютные значения моментов инерции J1 = 33, J2 = 22.1, J3 = 21.7, а углового момента M = \M\ = 0.244. Фаза обозначена цветом, шкала для фазы справа. (b) Две вложенные друг в друга тороидальные поверхности, на которых дивергенция вектора Пойнтинга обращается в ноль. Они отделяют области с источниками от областей, где энергия излучения уходит в среду. (с) Солитон в разрезе и (d) — решение, полученное зеркальным отображением x ^ —ж, но приводимое к исходному поворотом (они идентичны). Параметры: g0 = 2.135, d = 0.06...............49

3.2 Область устойчивости твердотельных солитонов (а), график энергии W(г) при б = 0.06 (Ь), в заштрихованной области устойчивых решений не существует, установление интегральных характеристик солитона — W и максимального момента инерции .1\ (с^) при

до = 2.135, б = 0.06..............................50

3.3 Зависимость угла в от коэффициента диффузии бу при фиксированном значении = 0.07. При бу < 0.004 — устойчивых решений нет; При 0.004 < бу < 0.01 — переходные режимы, возможно хаотичное движение; при 0.01 < бу < 0.038 — режим прецессии; 0.038 < бу < 0.07 — колеблющийся прецессон, устойчивое вращение вокруг любой оси, лежащей в плоскости ху; бу = 0.07 — твердотельный прецессон, устойчивое вращение под любым углом к оси т; бу > 0.07 — твердотельный прецессон, устойчивое вращение только вокруг оси т (два возможных положения). 51

3.4 Зависимость главных моментов инерции от г при д0 = 2.137,

= 0.07 и (а) бу = 0.01 (хаос и перемежаемость), (Ь) бу = 0.03 (прецессия), (с) бу = 0.04 (периодические колебания), бу = 0.075

("твердотельный" режим, вращение)....................53

3.5 Схема построения трёхмерного начального распределения из двумерного решения: (а) образующая структура в виде двух сильносвязанных вихревых солитонов с единичным топологическим зарядом, (Ь) построение трёхмерного распределения поворотом

вокруг оси т с одновременным вращением вокруг своей оси на 3п. . . 54

3.6 Вихревые линии различных типов солитонов, стрелками показана ориентация линии...............................56

3.7 Расщепление солитона с т = 3 на элементарные вихри..........57

3.8 Типичная зависимость разности двух главных моментов инерции от

г для метастабильных вихрей........................57

3.9 Распределение интенсивности в сечении перпендикулярном оси т и проходящем через замкнутую вихревую линию (а) и в сечении, проходящем параллельно этой оси и проходящем через незамкнутую вихревую линию (Ь). Распределение фазы в этих сечениях, соответственно (с^).............................59

3.10 Изоповерхности интенсивности по уровню 0.5 асимметричного "яблока"; до = 2.128, б = 0.037....................... 60

3.11 Области устойчивости разных типов узловых солитонов и различные решения типа "яблока" (внизу). Зелёные кружки соответствуют параметрам, при которых получен устойчивый двойной тор (г), синий

кружок — устойчивое зацепление Хопфа (^................62

3.12 Графики зависимости энергии W(д0) (вверху) и максимального момента инерции (до) (внизу) при медленном изменении д0......63

3.13 Гистерезисная зависимость средних значений максимального момента инерции (наверху) и энергии (внизу) от коэффициента усиления д0 при ( = 0.037. 1 — асимметричное решение, 2 —

асимметричное решение второго типа с ещё более закрученными вихревыми линиями. На врезках показаны изоповерхности обоих

решений в области бистабильности..................... 65

3.14 Зависимость энергии и момента инерции ,1тах для двух устойчивых асимметричных решений от х при одинаковых параметрах. 1 — асимметричное решение, 2 — асимметричное решение второго типа.

При х « 60000 коэффициент д0 был изменён с 2.1288 до 2.129......66

3.15 Кажущийся гистерезис, наблюдающийся при медленном изменении д0 для возбуждённых решений между 6-угольным и асимметричным состоянием, ( = 0.037............................ 67

3.16 Зависимость модуля главного момента инерции (вверху) и энергии W от эволюционной переменной х для метастабильного 6-угольного состояния при д0 = 2.1271, ( = 0.037. После скачка при х « 60000 солитон становится асимметричным (появляется модуляция

вихревой линии с 5 периодами).......................67

3.17 Области источников (У£ > 0, красный цвет) энергии для осесимметричного (а) и асимметричного (Ь) "яблока", и "прецессона" (с), Замкнутые траектории на (а) и/или (с): вихревые линии 1, седловая 2, отталкивающие, неустойчивые 3, 4, и притягивающая, устойчивая 5 предельные линии. На — вихревая линия 1, предельные линии 2-4 и спираль на сепаратрисной поверхности, разделяющей области притяжения вихревой линии и удаляющихся от нее линий потока энергии. Стрелки показывают направление тангенциальной составляющей потока энергии вблизи вихревой линии 1...............................69

3.18 Особые элементы фазовой плоскости (а) и "рядовые" линии потоков энергии (Ь), демонстрирующие характер этих элементов для симметричного "яблока". Вертикальная пунктирная прямая р = 0 совпадает с незамкнутой вихревой линией. На ней расположены три вырожденные особые точки: устойчивый узел 6 и два седла 7 и 8. Близко расположенные сепаратрисы 9, выходящие из неустойчивого предельного цикла 10 и идущие в седла 2, 4 и 7, отделяют ограниченную область с траекториями, локализованными внутри солитона, от области с траекториями, уходящими на периферию. 1 — устойчивый фокус, 3 —седло, 5 — устойчивый узел. Неограниченная сепаратриса 11 отделяет траектории, уходящие на периферию вблизи незамкнутой вихревой линии, от приближающихся к "притягивающей" линии 12.........................70

3.19 Вращая полуплоскость (а) (смотри также рисунок 3.18(а)) вокруг оси т, получаем трёхмерное фазовое пространство (Ь). В разрезе (Ь) показаны особые поверхности, разделяющие линии потока энергии с качественно различающимся поведением. Они получены вращением вокруг оси т особых линий, изображённых на рисунке 3.18: предельного цикла 10 и сепаратрис, как приходящих в седла 2, 4, 7 из предельного цикла, так и соединяющих седла 3, 8, седла и узлы 2, 6, 2, 5, 3, 4, 4, 5, и уходящих на периферию сепаратрисы 11 и линии 12.

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

3.20 Двойной тор при небольших г и он же, потерявший симметрию

(справа, цвет показывает фазу). ...................... 73

3.21 Фазовый портрет симметричного (а) и асимметричного (Ь) двойного тора; до = 2.12, б = 0.06...........................73

3.22 Динамика устойчивого тривиального узла.................75

3.23 Изоповерхность тривиального узла, фаза изображена цветом......75

3.24 Трилистник с т = 3, фаза изображена цветом, изоповерхность интенсивности по уровню 0.5........................77

3.25 Главный момент инерции поворачивается на угол п/2 вдоль оси с наибольшим значением б..........................78

3.26 Поворот главного момента инерции в плоскость с большими значениями й.................................78

3.27 Траектория центра метастабильного зацепления Соломона с т = 3

(три незамкнутые линии)...........................80

3.28 Установление метастабильного узла Соломона с т = 3, слева график

3\(г), справа — W(г)............................80

3.29 Устойчивое асимметричное осциллирующее решение, т = 2. Параметры: д0 = 2.116, й = 0.06. Изоповерхность по уровню I = 0.7, цвет соответствует фазе...........................81

3.30 График зависимости д0(г) (а), энергия поля Wf (Ь) и отклонение энергии среды от фонового значения ЬWm (с) при увеличении и уменьшении д0................................84

3.31 Диаграммы начального (а, узел Соломона), конечного (^ два яблока) и основных промежуточных состояний. Пересечения пронумерованы, направление стрелок у вихревых линий определено по правилу буравчика. Внизу представлены зависимости топологических основных характеристик от г: (]*) линия 1 — минимальное число пересечений Сг, линия 2 — количество супервитков и, (к) линия 3 — коэффициент зацепления; линия 4 на обоих графиках соответствует обратному ходу, т.е. уменьшению д0

после формирования двух устойчивых яблок................86

3.32 Траектории вектора Пойнтинга для изолированных колец 2 и 3, отделившихся от незамкнутой линии 1. На схеме (а) цифрами 4-9 обозначены предельные циклы: неустойчивые циклы: 6 для линии 3, 4 и 5 для линии 2; сёдла: 7-9. Одинарные стрелки на линиях 1-3 показывают ориентацию линий (топологический заряд). На (Ь) показаны дополнительные спиральные траектории на поверхности, разделяющей области притяжения потоков колец 2 и 3. Двойные

стрелки показывают направление потоков энергии............87

3.33 Изоповерхность интенсивности для неустойчивого зацепления

Хопфа с т = 2................................87

3.34 Взаимодействие между замкнутой и незамкнутой вихревыми линиями, после реакции остаётся одна замкнутая и две незамкнутые линии......................................88

3.35 Две незамкнутые линии образуют узел "стремя"..............88

3.36 Рождение двух изолированных вихревых колец при расхождении образовавшихся "яблок"...........................89

3.37 Разлетающиеся "яблоки"...........................89

3.38 Окончательное разделение двух устойчивых "яблок", коллапс перемычки................................... 89

3.39 Зависимость энергии поля W от д0 при переходе от потерявшего устойчивость трилистника к зацеплению Хопфа; б = 0.05, д0 = 2.1223. 91

3.40 Изоповерхность интенсивности зацепления Хопфа с т = 3;

б = 0.05, до = 2.12...............................91

3.41 Графики установления постоянных значений энергии (слева) и моментов инерции (справа) для твердотельного зацепления Хопфа;

до = 2.122, б = 0.065............................. 92

3.42 Изоповерхность интенсивности твердотельного симметричного зацепления Хопфа с т = 3; б = 0.063, до = 2.12..............92

3.43 Фазовый портрет симметричного твердотельного зацепления Хопфа;

б = 0.07, до = 2.124..............................93

А.1 Один шаг интегрирования по г.......................116

А.2 Траектория движения 7-солитонного комплекса (слева), участок

траектории в месте её замыкания в кольцо (справа). Небольшой зазор возникает из-за численных ошибок. Синим цветом выделена область счёта......................................118

А.3 Предельные циклы твердотельного симметричного зацепления Хопфа, полученные при расчётах с различной точностью. Параметры: б = 0.07, до = 2.124; 1: бх = бу = бт = 0.658, бг = 0.1;

2: dx = dy = dT = 0.584, dz = 0.04; 3: dx = dy = dT = 0.584, dz = 0.01; 4: dx = dy = dT = 0.519, dz = 0.01;

5: dx = dy = dT = 0.459, dz = 0.01.....................119

А.4 Работа с программой lbecon в Gentoo Linux................121

Приложение А Используемые численные методы и программы

А.1 От физической модели к компьютерной программе

Решение дифференциальных уравнений в частных производных — одна из самых распространённых задач, стоящих перед теоретиками. Для многих нелинейных уравнений на сегодняшний день не существует аналитического решения, для некоторых существуют методы получения точного решения (например, метод обратной задачи рассеяния [69]). Если аналитического решения не существует, то у исследователя есть два пути. Можно искать аналитическое решение приближённого уравнения, например, с помощью теории возмущений, вводя дополнительные ограничения на параметры или функции. А можно искать приближённое решение исходного уравнения, заменяя производные конечными разностями и превращая его в дискретное, а затем интегрируя его численно. Возможны и комбинации обоих методов, когда, например, граничные условия для уравнения вначале вычисляются аналитически с помощью разложений, а само уравнение решается численно. В обоих случаях мы имеем приближённое решение, которое может довольно точно описывать реальность, а может лишь отдалённо напоминать её, причём далеко не всегда можно сказать, какой способ лучше. В первом случае мы ещё больше сужаем границы применимости физической модели, поэтому математически точное решение далеко не всегда отражает реальные физические процессы. В то же время, мы можем найти неустойчивые решения, которые не так-то просто получить вторым способом. Во втором случае мы не наносим вреда физической модели, но фактически работаем с конечно-разностными уравнениями и из-за этого получаем массу проблем из-за влияния сетки на устойчивость и динамику решений. То есть всё-таки портим — получаем дискретную модель. В каком случае мы оказываемся ближе к реальности? Некоторые примеры из истории компьютерных симуляций, например, первые попытки моделирования земной атмосферы, показывают, что порой с помощью неточных численных методов, прибегая к различным трюкам и даже нарушая законы физики в борьбе с численными ошибками, можно с успехом воссоздавать природные

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

Ниже будут рассмотрены эффективные и в то же время легко реализуемые способы численного интегрирования дифференциальных уравнений в частных производных. Обсуждаются некоторые проблемы, возникающие при поиске локализованных решений со сложной динамикой. Речь пойдёт об уравнениях параболического типа с двумя и тремя пространственными измерениями. Например, обобщённое комплексное уравнение Гинзбурга-Ландау, нелинейное уравнение Шрёдингера и Свифта-Хоенберга. Будут даны практические рекомендации по написанию собственных программ, поиску решений, проверке устойчивости и оценке ошибок. Это будет сделано без строгих математических доказательств, основываясь лишь на успешном опыте применения этих программ.

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

А.2 Численное интегрирование многомерной задачи, программа Lbecon

Эта программа на Fortran 90 подходит для решения параболических уравнений с различным числом пространственных измерений и использует OpenMP, то есть эффективно работает на многопроцессорных системах. Первый более или менее работающий код был написан в 2003 году. В основе лежит конечно-разностный метод Кранка-Николсона и split-step метод. Рассмотрим задачу 3+1 измерений: три пространственных измерения и одно эволюционное. Определившись вначале с количеством точек, величиной шага и фактическим размером области, мы должны подготовить трёхмерное начальное распределение и поместить его в счётную область в виде куба, содержащего N3 точек. Когда-то программа могла работать с разным количеством точек по всем измерениям, но с кубической областью работать оказалось намного проще. После этого надо делать шаги по эволюционной переменной и наблюдать как меняется решение. Как же сделать этот шаг в трёхмерной задаче с нелинейностью? Шаг разбивается на два — линейный и нелинейный (split-step). Для линейного шага годится метод Кранка-Николсона, нелинейный можно считать методом Эйлера или Рунге-Кутта. Метод Кранка-Николсона в одном измерении использует 6 точек, три на предыдущем слое и три на следующем. Его главное преимущество состоит в том, что он сохраняет энергию при линейном шаге, не давая затухания, но при этом очень устойчив. Расписывая вторые производные в операторе Лапласа по этой схеме, мы получаем две трёхдиагональные матрицы, одну из которых нам необходимо обратить. Можно расписать "честно" производные в двух и более измерениях, но сложность задачи увеличивается. А можно схитрить. Основная идея состоит в последовательном интегрировании одномерных задач, а нелинейный шаг, в свою очередь, разбивается на три, при этом каждая треть совмещается с каждым из линейных шагов. Таким образом, мы последовательно интегрируем наш куб вдоль x, вдоль y и затем вдоль z. Каждый раз куб разбивается на N2 одномерных задач, которые можно интегрировать независимо друг от друга — это легко распараллеливается. Желательно, чтобы число точек N нацело делилось на число потоков.

Рассмотрим процедуру перехода к конечно-разностному уравнению подробнее. Пусть E(x,y,z; t) наша неизвестная функция и линейная часть уравнения содержит такой оператор Лапласа: dE + dE + дЕ • Также в уравнении имеется нелинейный член, зависящий от интенсивности \E|2. Пусть Em — значение функ-

ции в 1-й точке по оси х в некий момент времени при фиксированных координатах ] по оси у и к по г. При этом вектор Ет содержит значения функции Е вдоль оси х. Тогда вторая производная по х запишется в виде

1( Ет++ - 2Ет+1 + Етх1 Ет 1 - 2Ет + Ет-г 2( дХ2 + дХ2 ) (А.2.1)

Аналогичным образом расписываются производные по оставшимся координатам. Метод можно назвать полунеявным — половина производной расписана на предыдущем т слое, половина на последующем т + 1 слое. В общем случае производные можно брать с весами, отличными от 1/2, например, можно

увеличить "неявность", а это может помочь при стабилизации начального расЕ т+1_е т

пределения. Производная по времени по двум точкам записывается как 1 м 1 . Линейные члены зависящие от Ет, если таковые имеются в уравнении, можно вписать в матрицы.

Для того, чтобы сделать линейный шаг по любому из измерений надо решить такую задачу:

Ет+1 = д-1рЕт. (А.2.2)

Здесь Ет — решение на предыдущей итерации, а Ет+1 — на следующей. Напомню, что речь идёт о векторе из значений функции вдоль одной из координатных осей при фиксированных двух других координатах. Матрицы Р и Ц размерности N х N могут содержать в себе зависимости от координат, если в уравнении есть линейные члены, например, параболический потенциал. Матрицы трёхдиа-гональные, если на границе функция (или её производная) обращается в ноль.

г + БсЯ -С 0 .. .0 0 0

-С г + 2Я -С . . .0 0 0

0 -С г + 2Я .. .0 0 0

0 0 0 ... -С г + 2Я -С

0 0 0 ... 0 С г + БСЯ

Q =

i-BM

C

0

C i - 2R C 0 C i- 2R

00 00 00

0 0

0 0

0 0

• •• C i- 2R

0 0 0

C

0

C i-BrR

(А.2.4)

Здесь R = и коэффициент Bc принимает два значения: 1 — если на границах обращается в ноль производная функции E и 2 — если сама функция. На каждом шаге приходится обращать матрицу Q. Делается это методом LU-разложения [71]. В детали этого метода лучше вникнуть самому, но практические реализации (желательно в виде исходных текстов, а не библиотек) можно взять готовые, например Linpack с сайта netlib.org. Исходные тексты нужных нам программ содержат несколько функций (zgbfa(), zgbsl() и т. д.), которым необходимо передать нужные комплексные массивы. Описание всех переменных есть в самих исходных текстах (в процессе работы я немного изменил их под свои нужды, уменьшив количество вызовов функций и увеличив быстродействие). Для хранения трёхдиагональной матрицы при этом требуется массив N х 4. Матрица Q, которую обращать не нужно, хранится в массиве N х 3.

На рисунке А.1 схематически изображён шаг по z. Его необходимо сделать для всех точек в плоскости xy, то есть для всех индексов i и j. Очевидно, интегрирование при этом происходит независимо, поэтому каждый шаг распараллеливается очевидным образом. В программе это реализовано в виде двух вложенных циклов по i и j , внешний из которых распараллеливается с помощью директивы $OMP PARALLEL DO. Для трёхмерного массива, в котором хранится функция E, надо не забыть указать атрибут SHARED, потому что куб у нас общий для всех потоков, а для временных переменных (уникальных в каждом из потоков) — атрибут PRIVATE. Для одного шага при этом потребуется порядка 3N3 действий.

Аналогичным образом делаются шаги по x и y. Нелинейный шаг делится между тремя линейными шагами. В принципе, его можно делать один раз в конце, но, распределяя его между тремя линейными шагами, мы в меньшей степени нарушаем изотропию пространства и заодно улучшаем производительность, потому что нелинейный шаг Эйлера или Рунге-Кутта, выполненный в виде отдельного цикла по всем элементам массива, потребует ещё N3 действий.

_ т+1 Ек-1

Е

т+1

Е

сЬ

т+1 к+1

С1

т

Ет-1/

1 1+С

Рисунок А.1 — Один шаг интегрирования по г.

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

к

При работе с программой довольно часто возникает ситуация, когда одиночное локализованное решение подходит к краю счётной области. Если граничные условия периодические, то решение просто переходит на другую сторону области. Но реализация периодических граничных условий требует дополнительных вычислительных затрат при обращении почти диагональных матриц с "уголками". Проще воспользоваться методом скользящего окна, в нашем случае — куба. При превышении на какой-либо границе порогового значения интенсивности решение сдвигается в центр, при это часть счётной области заполняется нулями. Пороговое значение определяется экспериментально, обычно оно на порядки больше величины фона (шумов с амплитудой 10"16). Если взять его слишком большим, то решение может взаимодействовать со стенками, менять траекторию движения и деформироваться. Если слишком мало — сдвиги будут слишком частые. В любом случае решение при сдвиге испытывает небольшое возмущение, иногда заметное на графиках, но в целом это работает почти в два раза быстрее периодических границ.

Как и любой split-step метод, этот работает хорошо и достаточно точно при условии, что шаги по пространству и времени выбраны не слишком большие, причём эти значения могут зависеть от вида решения и его динамики. Для быстро движущихся, осциллирующих решений очень важно отличать численные неустойчивости, а также факторы, влияющие на движение структуры при долгих расчётах. При движении неизбежны колебания всех параметров решения из-за его движения в сетке. В случае быстрого движения такие колебания могут исказить профиль солитона или вызвать его распад. Другим неприятным эффектом является искажение траектории движения солитона и дрожание вектора угловой скорости. В таком случае надо уменьшить шаг по пространству. Как правило, в расчётах ошибки (амплитуда колебаний) составляли доли процента. А вот искажения траектории движения центра в продолжительных расчётах были заметны даже на максимальной точности. Поэтому при изучении движущихся солитонов необходимо выбирать параметры так, чтобы физические эффекты были сильнее (заметнее) численных. Например, для анизотропного пространства разница между измерениями должна быть не очень маленькой. Если речь идёт о коэффициентах диффузии, то они должны отличаться хотя бы на несколько процентов.

Точность метода легко проверить длительным расчётом, результат которого можно предсказать заранее. Во второй главе на рисунке 2.14 показана траектория движения асимметричного солитонного комплекса в интерферометре. Из теоре-

тического анализа, связывающего тип движения комплекса с его симметрией, следует, что траектория его центра представляет собой окружность. Из-за медленного движения комплекса время разворота измеряется сотнями тысяч единиц. В итоге после миллионов итераций неизбежно накапливаются ошибки, которые становятся заметными при увеличении на рисунке А.2. С уменьшением величины шагов точность увеличивается, однако растёт и количество итераций, по-

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

10 0 10 20 30 40 50 60

X X

Рисунок А.2 — Траектория движения 7-солитонного комплекса (слева), участок траектории в месте её замыкания в кольцо (справа). Небольшой зазор возникает из-за численных ошибок. Синим цветом выделена область счёта.

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

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

563.6

563.4

12

563.2

563

Рисунок А.3 — Предельные циклы твердотельного симметричного зацепления Хопфа, полученные при расчётах с различной точностью. Параметры:

( = 0.07, до = 2.124; 1: (х = (у = (т = 0.658, (г = 0.1; 2: (х = (у = (т = 0.584, (г = 0.04; 3: (х = (у = (т = 0.584, (г = 0.01; 4: (х = (у = (т = 0.519, (г = 0.01; 5: (х = (у = (т = 0.459, (г = 0.01.

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

Как уже было отмечено ранее, в численном эксперименте очень важно графическое представление решения. Без него легко упустить искомое, проскочить интересную динамику, проглядеть развитие численной неустойчивости. При желании, "автоматизация" реализуется в виде скриптов или простой логики, но человеческий контроль, на мой взгляд, необходим. В конце концов, на необычные, красивые решения хочется посмотреть со всех сторон. Для визуализации данных использовалась программа IBM openDX, в которой есть удобная возможность управлять всем с помощью скриптов, из командной строки [72]. С помощью скриптов можно создавать различные графические файлы, например графики в формате PNG или трёхмерные STL. Для двумерных и простых трёхмерных графиков использовался Gnuplot, на вид простая, но очень мощная система, которая тоже может работать в неинтерактивном режиме. Используя эти программы, можно сделать простой web-интерфейс к своей виртуальной лаборатории и следить за результатами, собранными с удалённых машин и кластеров. Уравлять заданиями лучше всё-таки по ssh, причём с графикой неплохо справляется xpra.

Для сборки программы в настоящий момент требуется Intel Fortran Compiler, компилятор Фортрана под ОС Линукс. До 2015 года его можно было без проблем загрузить с официального сайта компании вместе с лицензией для некоммерческого использования. Эти лицензии (permanent uncounted) работают и сейчас, однако Интел поменял правила и теперь для получения лицензии необходимо каким-то образом подтвердить то, что вы являетесь студентом. В принципе, после некоторой доработки, программа может быть собрана в любом Linux-

ü verna@hitman:~/work/2018/S

verna@hitman ~/work/20l8/SAT_3 $

verna@hitman ~/work/20l8/SAT_3 $ ./lbecon refresh=l©. write_sequence=no

parsing command line arguments... ok! PROGRAM STARTED 18/10/2018 at 19:39:25. loading config file [lbecon.conf]... ok ?>time to calculate?

verna@hitman ~/work/20l8/SAT_3 $ ./rundx.sh 0 && pqiv dx.png

Starting render from 0 to 0 ...

Starting DX executive

Open Visualization Data Explorer

More Info at www.research.ibm.com/dx

and www.opendx.org

Version - 4.4.4

Memory cache will use 718 MB (44 for small items, 674 for large)

Starting 2 threads... allocating 89.057 Mb... ok! call matrix... ok! Calculate initilal distribution, brick: 40 x 4© x 40 92 : 60 : 60 # 0

5© : 100 : 100 # 2

Adding 7 bricks...

X/Y/Z step 0.559 0.559 0.

Output every 10.000

Hit "e1 to enter console

| No | t | time |

0: worker here [8684] 0: opening file '/home/veri 0: opening file '/home/veri 0: Colormap_l: begin 0: Colormap_l: huemap={[0 0.59253335 ] [0.16666667 0.5! 43333 ] [0.33333334 0.444399! ] [0.5 0.33329999 ] [0.55555!

0.22219998 ] [0.72222221 0.: 110002 ] [0.8888889 0.074066i 0: Colormap_l: satmap={[0 ] [0.22222222 1 ] [0.2777777!

[0.5 1 ] [0.55555558 1 ] [0 7779 1 ] [0.83333331 1 ] [0.1 0: Colormap_l: valmap={[0

na/work/20l8/SAT_3/data/0000O.im.dat' na/work/20l8/SAT_3/data/0000O.re.dat'

irface at I = 0.01

0 0.00 19:39:52 168812.999910179846

_>enter command, "help" or "q" to quit pump=? _> 2.17100000000000 flag О

time=?

_>current time = 9.40000000000000 flag 8

q

>Current time = >lnterrupt?[y/n]

$ import -window root lbecon.png

— verna@loginl:~/work/2018/sat... — verna@hitman:~

j verna@hitman:--/work/2018/S... verna@hitman:--/work/2018/S... xpqiv: dx.png (640x505) 100.00...

Рисунок А.4 — Работа с программой lbecon в Gentoo Linux.

дистрибутиве с помощью открытого компилятора gfortran, входящего в GCC. Одну из первоначальных версий программы я выложил под лицензией GPL, она доступна по адресу https://sourceforge.net/projects/lbecon. К сожалению, подробной документации для неё до сих пор не существует, а без неё работа с программой затруднительна. Эта версия для комплексного уравнения Гинзбурга-Ландау размерности 3+1. На сайте есть минимальное описание программы и примеры конфигурационных файлов для создания комплексов лазерных пуль.

А.3 Решение краевой стационарной задачи методом Ньютона

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

Приравнивая к нулю производную по времени в 2.1.4, получаем уравнение

{г + дХ - \<х)\2 - 0}еХ = Ein- (А.3.1)

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

= ^ (А.3.2)

дх

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

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

Выбрав шаг сетки, ширину области и число точек N, расписываем вторые производные в уравнении А.3.1 по трём точкам аналогично А.2.1. Дифференциальное уравнение сводится к системе нелинейных алгебраических уравнений,

которые уже можно решать методом Ньютона. После необходимых преобразований должно получиться следующее уравнение на вектор e, состоящий из значений неизвестной функции в узлах сетки e(xi) = 6i :

F (e) = 0 (А.3.3)

Метод Ньютона состоит из последовательных итераций, которые надо совершать до тех пор, пока норма вектора на n-й итерации F(en) не уменьшится до желаемой величины:

en+1 = en - J(en)-1F(en), (А.3.4)

где J — якобиан системы. Поскольку уравнение А.3.1 — комплексное, удобно превратить его в систему двух вещественных уравнений для Re(ei) и Im(ei). Тогда вектор e имеет размерность 2N, а якобиан представляет собой ленточную матрицу 2N х 2N. Уравнение А.3.1 достаточно простое, поэтому все производные вычисляются аналитически. Основная задача программиста сводится к безошибочному написанию элементов матриц. Для обращения якобиана используются указанные в предыдущей главе подпрограммы Linpack, в данном случае, для вещественных матриц.

Данный метод был реализован на Fortran 90 и использовался, в частности, для вычислений с квадратичной точностью (с невязкой порядка 10-60) при уточнении профиля тёмного солитона с последующим анализом квантовых флуктуаций. Также метод применялся для нахождения стационарных локализованных структур в вырожденном оптическом параметрическом осцилляторе.

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