CC BY
4
Вычислительные технологии
Том 22, № 5, 2017
Полностью консервативная численная схема
для нелинейного уравнения Шрёдингера с высшими
нелинейностями
А. Д. Булыгин*, А. А. Землянов
Институт оптики атмосферы им. В.Е. Зуева СО РАН, Томск, Россия *Контактный e-mail: b.a.d@iao.ru
Анализируется выполнение условий полной консервативности различных численных схем, наиболее широко используемых при численном исследовании нелинейного уравнения Шрёдингера (НУШ) в задаче о филаментации. Рассматривается стационарная радиально-симметричная модель, которая имеет два интеграла движения — функцию Гамильтона и "число частиц". Определен явный вид дисба-лансных членов, возникающих в несимметричных схемах расщепления по физическим факторам. На основе численных расчетов установлено, что использование стандартных процедур численной реализации НУШ с высшими нелинейностями приводит к существенному нарушению закона сохранения функции Гамильтона в процессе филаментации. Найдены условия на численную сетку, при выполнении которых удается удовлетворить законам сохранения для полностью симметричной разностной схемы и схемы расщепления по физическим факторам в ее дискретно-разностной реализации.
Ключевые слова: филаментация, нелинейное уравнение Шрёдингера, законы сохранения, полностью консервативная схема.
Введение
Исследование самофокусировки и филаментации на основе нелинейного уравнения Шрёдингера (НУШ) [1] в различных нелинейных средах проводится преимущественно численными методами, и причина этого в том, что строгих аналитических результатов в этой области исследования крайне мало. Такое положение дел оставляет исследователя практически без какого-либо инструментария для проверки достоверности результатов численных экспериментов. Пожалуй, единственным строгим результатом в теории самофокусировки является набор законов сохранения, которые обсуждаются, как правило, только для керровской среды [1, 2]. Знание этих законов позволяет отслеживать корректность численных расчетов на этапе самофокусировки [3]. Однако существуют законы сохранения в среде более общего типа нелинейности [4]. Частным классом общего вида нелинейности является полиномиальный вид нелинейности, а уравнение с таким полиномиальным видом нелинейности мы будем называть уравнением Шрёдингера с высшими нелинейностями. Такой вид уравнений Шрёдингера используется при описании процесса филаментации [1, 5, 6]. Казалось бы естественным, чтобы численное моделирование процесса филаментации для таких сред с высшими
© ИВТ СО РАН, 2017
нелинейностями типа [4, 5] также проверялось на корректность с использованием знаний об этих интегралах. Однако в литературе мы не нашли информации о соответствии результатов численных расчетов по филаментации этим единственно известным и необходимым условиям достоверности, именно условиям выполнения законов сохранения в процессе филаментации.
В настоящей работе поставлен вопрос о выполнении необходимых условий достоверности в отношении наиболее употребляемых численных схем, используемых при численном решении НУШ. Прежде всего это метод расщепления по физическим факторам [7, 8], когда последовательно реализуются шаг дифракции и шаг нелинейности, причем шаг дифракции может быть реализован методами Кранка — Николсон [7, 8, 12, 13] и быстрого фурье-преобразования [3, 8] (в том числе в радиально симметричной схеме методом Ханкеля [9]). Кроме того, рассмотрена симметричная разностная схема с итерационным методом вычисления нелинейного вклада, один из вариантов которой рассмотрен в работе [6]. Для анализа выполнения условий полной консервативности различных разностных схем выписан явный вид дисбалансных членов в соответствии со стандартной методикой [10]. Проанализировано влияние этих дисбалансных членов на примере конкретных численных реализаций. Анализ этих численных экспериментов позволит выделить схемы и сеточные условия для них, при которых проявление дисбалансных членов пренебрежимо мало.
1. Уравнение динамики и законы сохранения
Рассмотрим НУШ в его наиболее простой стационарной форме [1, 7], которое в безразмерных координатах [1, 11] имеет вид1
4" - = 0. (1)
дх 8и* ' ;
Комплексная амплитуда светового поля и(х, х) [1] (нормированная на свое начальное максимальное значение) зависит от координат двумерного вектора х (нормированного на начальный размер пучка Д0), который принадлежит плоскости, трансверсо-нальной направлению распространения излучения (нормированного, в свою очередь, на длину пучка Рэлея Ьд = к0Щ,/2, где к0 = 2и/А0 — волновое число, А0 — длина волны излучения на основной частоте). Функция Гамильтона определяется соотношением
н = / ^ + миIV. (2)
1 Здесь и далее для индексов используем обозначения Эйнштейна, если нигде специальным образом не оговорено другое. Так, например, в рассматриваемой нами ситуации
ди ди* ди ди
д, ид1 и * = — —- +
3 дх\ дх\ ' дх2 дх2
Для интеграла по всему пространству для кратности используем обозначение
г г ^
/ /¿х = / fdxldx2.
'то ./то
Потенциал hp будем рассматривать в виде
, Lr 4 LR 2 m
h = -rT\U 1 + ™T \U1 •
lilLim
Здесь первое слагаемое соответствует кубичной (керровской) нелинейности, которая определяет эффект самофокусировки [1, 7]; Lk = (k0I0n2)-1 — длина самовоздействия, обусловленная керровской нелинейностью [1] и определяемая через коэффициент п2 при нелинейной добавке к показателю преломления среды п0 и начальную максимальную интенсивность лазерного излучения 10. Второе слагаемое учитывает эффект насыщения нелинейности [1] (Lm = (k0imnm)-1 — длина самовоздействия, обусловленная нелинейностью порядка m и определяемая, соответственно, коэффициентом nm) либо эффективно учитывает плазмообразование и дефокусировку в наведенной плазме [6]. В любом случае второе слагаемое обеспечивает остановку коллапса (достаточным условием его реализации является достижение критического значения мощности пучка, которое для гауссова пучка Рсг = У^/(2кп2щ)) и формирование высокоинтенсивного пространственно-локализованного состояния светового поля — филамента [1].
Обсудим два закона сохранения, которые соответствуют двум глобальным симмет-риям: внешней — сдвиг по эволюционной координате z ^ z + а и внутренней — сдвиг фазы комплексного поля U ^ Uem. Они имеют соответственно вид
H = J hdx = const, (3)
где h = hk + hp — плотность функции Гамильтона, имеющая кинетическую и потенциальную составляющие, через которые определяется безразмерная диэлектрическая проницаемость е = dhp/d \ U|2,
Е = J \U\2 dx = const. (4)
Эту величину будем отождествлять с мерой числа частиц (по аналогии с квантовой механикой).
Заметим, что лапласиан, который входит в (1) (и генерируется, соответственно, hk) часто называют оператором дифракции [3], однако с физической точки зрения это принципиально неправильно. Так, в эйкональном приближении действие лапласиана остается, а дифракции уже нет. Дифракционное слагаемое содержится в кинетической части
функции Гамильтона и легко выделяется при переходе к амплитудно-фазовому пред-
^j ^ j
ставлению комплексного поля hk = ^ 2 + hd [11]. Тем не менее мы будем время от
времени использовать эту терминологию, понимая, что это некоторый сленг.
Следствием интегральных соотношений (3) и (4) являются два уравнения дивер-
2
гентного вида
dh dQk d \U\2 dSk
dz dxk} dz dxk
2Есть еще интегралы движения, которые связаны с симметриями НУШ относительно сдвига пространственных координат и поворотами вокруг оси. Однако в настоящей работе мы будем рассматривать задачи в рамках аксиальной симметрии, а поэтому не будем интересоваться здесь этими интегралами.
Здесь введены следующие компоненты векторов:
„к идки* - и*дки , „к ^д1идки* - д'дМ*дки як =---, дк = еБк +—--——--.
Подчеркнем, что выполнение законов сохранения интегралов движения должно рассматриваться как необходимое условие достоверности численных расчетов. В связи этим проведем серию численных тестов для различных схем, используемых при решении НУШ.
1а. Метод расщепления по физическим факторам [7]. В нашем случае эта схема реализуется следующим образом: сначала осуществляется шаг линейного распространения посредством прямого и обратного преобразования Фурье [8]. Затем выполняется шаг прохождения фазового экрана, который наведен нелинейностью. В случае радиальной симметрии вместо преобразования Фурье перейдем к преобразованию Ханкеля [9].
1Ь. Метод расщепления по физическим факторам, когда реализация шага линейного распространения осуществляется посредством разностной схемы, являющейся консервативной только для оператора дифракции [7].
2. Полностью консервативная разностная схема, реализуемая посредством ее симметричного варианта, которая в декартовых координатах использована в работе [6]. В радиальной же постановке задачи нам незнакома работа, в которой она была бы описана и реализована в приложении к задачам о самофокусировке и филаментации. Поэтому рассмотрим ее более подробно.
Запишем уравнения (1) в радиальных координатах
иг = 4 (игг + 1и,) + /(Г, г), / = е(|и|2)и. (5)
Симметричная разностная реализация определяет следующую аппроксимацию конечными разностями :
и = иГ -и и = (и]+1-иП) + (и+ -иЬ)
и,
4 ¿г
( и]+1 + и-1 - 2и]+1) + (и]+1 + и]-1 - 2и])
2 2
(е(|и]+1|2)и]+1) + (е(|иЛ2)и]) ; 2 . Приведенная система является нелинейной. Решать ее будем итерационным методом Ньютона [10]. В качестве нулевой итерации выбираем значение функций на предыдущем шаге. На каждом итерационном шаге I решается матричное уравнение
А,и1+1 -си^1 + ви+1 = -М, (6)
3Заметим, что выбранная нами схема конечно-разностной аппроксимации оператора Лапласа, как можно показать путем несложных арифметических процедур, полностью эквивалентна дивергентной численной реализации:
(и„ + ) = И (__Ч) = 4. и1/2М! - ЧА-^\ + 0«1_).
\ г ) _ _ V _ / г^а_\ аг аг )
Рис. 1. Распределение величины ЬИ = sgn(H)log(1 + |И|) вдоль трассы распространения для различных численных схем: расщепление по физическим факторам с прохождением шага дифракции методом быстрого преобразования Фурье на декартовой сетке (кривая I), методом преобразования Ханкеля для радиальной реализации (2), с помощью симметричной разностной схемы (3) и схемы Кранка — Николсон (4)
где матрицы А, В, С, Г легко находятся из уравнения (5) по стандартной процедуре, которая описана, например, в [10]. Уравнение (6) решалось методом прогонки. Приведенная здесь схема является в некотором смысле аксиальным аналогом симметричной схемы, рассмотренной в работе [6].
Надо отметить, что в работе [6] совершенно правомерно замечено, что в случае использования метода расщепления по физическим факторам вопрос о полной консервативности требует отдельного исследования. Этому отчасти и посвящена данная работа. Однако и в случае реализации полностью консервативной схемы необходимо следить за параметрами сетки, при которых с достаточной степенью точности выполняются законы сохранения. Авторы работы [6] приводят неправильный вид функции Гамильтона для уравнения динамики, что могло бы рассматриваться как опечатка, однако дальнейший анализ, проводимый в статье, говорит об обратном. В действительности автор выписывает функцию Гамильтона только для среды керровского типа.
Для тестового примера рассматривался гауссов пучок миллиметрового радиуса с длиной волны 800 нм и мощностью Р0 = 7РСГ (для воздуха Рсг = 3.18). Значение параметра т = 8, величина Ьт подбиралось так, чтобы максимальное значение интенсивности достигало типичного уровня при филаментации в воздухе 1т & 1014 Вт/см2 [1]. Это вполне соответствует работам [5, 6]. Размер поперечного шага сетки принят равным 10, а продольного шага ¿г выбран по адаптивной схеме ¿г = 0.05/1тах, где 1тах — максимальное значение величины \и|2 [7, 8].
Как видно из графиков (рис. 1), закон сохранения энергии хорошо выполняется для всех схем, но лишь в среде керровского типа, т. е. до начала филаментации (понимаемой в нашем контексте как включение высших нелинейностей). Это абсолютно согласуется с результатами работ [3, 7, 8, 11]. При этом закон сохранения первого интеграла движения — "числа частиц" выполняется практические идеально для всех схем. Важно обратить внимание лишь на то, что ни одна из схем не обеспечивает закона сохранения второго интеграла с момента начала филаментации.
2. Анализ дисбалансных членов в различных схемах и условия их подавления
Как известно, причина нарушения законов сохранения на конечно-разностных схемах обусловлена появлением фиктивных источников энергии, вызванных формированием численных дисбалансов [10]. Для того чтобы проанализировать, при каких параметрах сетки роль этих фиктивных источников будет пренебрежимо мала, нам представляется полезным выписать их в явной форме. Выделение дисбаланса членов основано на следующей форме записи для приращения плотности функции Гамильтона:
hl+1 - hl = dzd,Qj + öh(1)dz2 + o(dz2).
Здесь öh(1 — коэффициент при дисбалансном вкладе первого порядка. В схеме расщепления по физическим факторам несложно вычислить коэффициент при дисбалансном вкладе первого порядка:
(1) = |U |2 djt dj б + д,е dzSj + д36вЦсУ 5eU = 4 + 2 + 4
Здесь 5ßU — дисбалансное слагаемое, порождаемое шагом дифракции. Перейдем к рассмотрению конкретных реализаций.
1a. Случай фурье-преобразования характеризуется тем, что даже шаг линейного распространения, реализованный по явной схеме, формально порождает дисбаланс первого порядка. Однако в линейном случае, поскольку ряд экспоненты сходится к решению задачи, фактическое суммарное влияние дисбалансов всех порядков обнуляется. Это замечание поучительно в том смысле, что схема с дисбалансами первого порядка может быть лучше, чем так называемые полностью консервативные схемы, к которым относится и рассматриваемая нами симметричная схема. Однако для более сложной, рассматриваемой нами ситуации экспоненциальный ряд не сходится к решению. Поэтому, вообще говоря, априорно мы отдаем предпочтение полностью консервативным схемам. Итак, при реализации шага дифракции методом фурье-преобразования последнее слагаемое, входящее в формулу (7), примет вид
SeU = F-1( kß)F-1(k3Ü),
где F — оператор преобразования Фурье; kj — координата в импульсном пространстве; U = FU — фурье-образ поля U.
1b. В случае аксиальной симметрии, когда используется преобразование Ханкеля, все построения такие же, как для преобразования Фурье на декартовой сетке, но с соответствующей модификацией кинетического слагаемого:
11 5eU = Fh (J kjUh)Fh (JkjUh).
Здесь введен оператор преобразования Ханкеля Fh [9], действующий аналогично фурье-преобразованию Uh = FhU; Jkj — координата в импульсном пространстве, связанном с преобразованием Ханкеля.
1c. В случае прохождения шага линейного распространения методом Кранка — Ни-колсон дисбаланс, порожденный кинетическим вкладом, обнуляется в силу симметричности метода.
(7)
той ее вариант, поддающийся оценкам, — это схема расщепления с локальными выражениями для дифракционного слагаемого, т.е. схема Кранка — Николсон. Здесь понадобится, в свою очередь, оценка для "силы дифракции [11]" д^На. Значение силы дифракции должно быть того же порядка, что и максимальное значение силы рефракции, которую часто оценивают из балансных соображений [1]. С учетом сказанного имеем
Отсюда следует оценка сверху для шага, который может обеспечить хорошую численную сходимость к постоянному значению функции Гамильтона:
Здесь Г/ — радиус филамента, имеющий для воздуха значение порядка 50-100 мкм.
Таким образом, полученные оценки дают условия на построение адаптивного шага ¿х, существенно отличающиеся от стандартно использующихся, отраженных в литературе по численному моделированию филаментации [6-8, 11]. Итак, только в случае полностью симметричной разностной схемы исчезают линейные по ¿х слагаемые в выражениях для дисбалансов. Однако ошибка, порожденная линейным дисбалансом, не всегда будет меньше ошибки, порожденной квадратичным по ¿х дисбалансом. Для того чтобы посмотреть, как ведут себя интегралы движения, надо уменьшать шаг по ¿х с учетом вышевыписанной оценки до тех пор, пока не будут получены приемлемые результаты. Как показывают численные эксперименты, можно использовать более эффективный способ выбора адаптивного шага. На рис. 2 приведены примеры расчета для схем при ¿г = ^/('Ц \и\4) с различными значениями параметра
Удивительно, но численные расчеты, выполненные с использованием схемы расщепления Кранка — Николсон, показывают сопоставимые и даже лучшие результаты, чем с применением симметричной схемы, при соответствующем выборе продольного и поперечного ¿г шагов, причем из-за отсутствия итерационного цикла схема реализуется быстрее. Этот факт, по-видимому, объясняется тем, что шаг нелинейности реализуется в схеме расщепления через действие экспоненты, которая удачно экстраполирует истинный результат.
В отношении схем с расщеплением, включающих преобразование Фурье, можно сказать, что нелокальный вид дисбалансов проявился в отсутствие сколько-нибудь существенного улучшения результатов для малого шага по йг и ¿г по сравнению с грубой сеткой (рис. 3), кроме того, его реализация требует времени на два и более порядков больше. В связи с этим для ускорения счета в настоящей работе использован распараллеленный вариант быстрого преобразования Фурье [8]. Для его реализации мы воспользовались вычислительным ресурсом кластера ССКЦ СО РАН. Тем не менее неочевидно, что при более сложной конфигурации поля, реализующейся в произвольном случае, мы всегда достигнем столь замечательного результата, как при расчетах с использованием схемы Кранка — Николсон, т. е. при конкретных расчетах всегда нужно отслеживать выполнение полной консервативности и в случае необходимости пользоваться полностью симметричной схемой, когда заведомо можно подобрать параметр
дj(дгд1 \и\\и\) &\и\2д-е.
■1-0-. н
г/ьЕ
-2.0-1—|-1-1-'-1-'-1-'-1—
0.00 0.25 0.50 0.75 1.00
Рис. 2. Пример распределения величины И = sgn(H) log(1 + |И|), полученной с использованием симметричной разностной схемы: кривые 1а (^ = 0.2), 1Ъ (^ = 0.1), 1с (^ = 0.05) при ¿г = 10; схема расщепления Кранка — Николсон (^ = 0.05) при ¿г = 10 — кривая 2 и ¿г = 1 — кривая 3
о.о-, т
0.0 0.2 0.4 0.6 0.8 г/Ьп
Рис. 3. Распределение плотности функции Гамильтона, полученной при использовании схемы расщепления с параметрами сетки ^ = 0.05 и ¿г = 1: 1 — фурье-преобразование на декартовой сетке; 2 — преобразование Ханкеля
сетки нужным образом. Чтобы более подробно проанализировать, как проявляются дисбалансы, выведем пространственные распределения таких физических величин, как к, кр, I = |2 на дистанции ^ = Ь^.
Как видно из графиков (рис. 4 и 5), несмотря на существенное отклонение в распределении для к, которое достигает двух порядков, неконсервативность схемы не приводит в данном случае к существенному искажению профиля I, за исключением малых низкоинтенсивных областей. Соответственно, область, где реализуются отклонения для величины потенциальной энергии кр, несущественна по сравнению с отклонениями для величины к. Это говорит о том, что в неконсервативных схемах наибольшей ошибке подвержено кинетическое слагаемое.
г_п-1-1—......i-1-1—......i-1-1—1111
0.01 0.1 1 г, мм
Рис. 4. Распределение величин Ьк = sgn(к) ^(1 + \к\) и Ькр = sgn(кp) log(1 + \кр|) от расстояния до центра соответственно для схемы Кранка — Николсон (кривые 1 и 2 — консервативный и неконсервативный варианты) и схемы с преобразованием Ханкеля (кривые 3 и 4 — консервативный и неконсервативный варианты)
Рис. 5. Распределение величины I для консервативного варианта схемы Кранка — Николсон (1) и неконсервативного варианта схемы с преобразованием Ханкеля (2)
Заключение
Рассмотрено выполнение условия полной консервативности при реализации различных численных схем, использующихся для решения НУШ, моделирующего филаментацию. На основе этого анализа установлены следующие факты.
1. Стандартные схемы и способы разбиения сеткой пространственных координат, используемые при численном решении НУШ с высшими нелинейностями, не удовлетворяют необходимому условию достоверности сохранения интегралов движения.
2. Однако из этих численных схем можно выбрать те, которые при существенном ужесточении требований к параметрам сетки удовлетворяют законам сохранения ин-
тегралов движения НУШ с высшими нелинейностями и обеспечивают достаточно эффективную численную реализацию. Это схемы расщепления по физическим факторам, в том числе схема Кранка — Николсон при прохождении дифракционного слоя, и абсолютно симметричная разностная схема.
Следует отметить, что дискретно-разностные схемы более предпочтительны с точки зрения не только возможностей распараллеливания и выбора адаптивных сеток по пространственным координатам [6], но, самое главное, достоверности численных расчетов, условия для выполнения которых в приложении к задаче о филаментации фемто-секундного лазерного импульса впервые сформулированы в данной работе.
Благодарности. Работа выполнена при финансовой поддержке РНФ (соглашение № 16-17-10128).
Список литературы / References
[1] Boyd, R.W., Lukishova, S.G., Shen, Y.R. Self-focusing: Past and present: Fundamentals and prospects. Topics in Applied Physics. Vol. 114. New York: Springer Verlag, 2009. 605 p.
[2] Власов С.Н., Таланов В.И. Самофокусировка волн. Н. Новгород: ИПФ РАН, 1997. 220 с.
Vlasov, S.N., Talanov, V.I. Self-focusing of waves. Nizhniy Novgorod: IPF RAN, 1997. 220 p. (In Russ.)
[3] Кандидов В.П., Косарева О.Г., Колтун А.А. Нелинейно-оптическая трансформация мощного фемтосекундного лазерного импульса в воздухе // Квантовая электроника. 2003. Т. 33, № 1. С. 69-75.
Kandidov, V.P., Kosareva, O.G., Koltun, A.A. Nonlinear optical transformation of a high-power femtosecond laser pulse in air // Quantum Electronics. 2003. Vol. 33, No. 1. P. 69-75. (In Russ.)
[4] Алексеенко В.Н. Об интегралах движения нелинейных уравнений типа Шрёдингера // Диф. уравнения. 1976. Т. 12, № 12. С. 1121-1122.
Alekseenko, V.N. The integrals of the motion of nonlinear equations of Schrodinger type // Dif. Uravneniya. 1976. Vol. 33, No. 1. P. 1121—1122. (In Russ.)
[5] Hui Gao, Weiwei Liu, See Leang Chin. Post-filamentation multiple light channel formation in air // Laser Physics. 2014. Vol. 24, No. 5. P. 1581-1596.
[6] Балашов А.Д., Пергамент А.Х. Математическое моделирование распространения фемтосекундного импульса // Матем. моделирование. 2006. Т. 18, № 4. С. 3-18. Balashov, A.D., Pergament, A.Kh. Mathematical modeling of femtosecond pulse propagation // Math. Models and Comput. Simulations. 2006. Vol. 18, No. 4. P. 3-18. (In Russ.)
[7] Chesnokov, S.S., Egorov, K.D., Kandidov, V.P., Vysloukh, V.A. The finite element method in problems of nonlinear opticsа //J. Numer. Meth. Eng. 1979. Vol. 14, No. 1. P. 1581-1596.
[8] Дергачев А.А. Формирование и характеристики плазменных каналов при филамента-ции фемтосекундного лазерного излучения в воздухе: Дис. ... канд. физ.-мат. наук. М.: МГУ, 2014. 146 с.
Dergachev, A.A. Formation and characteristics of plasma channels under filamentation of femtosecond laser radiation in air: Dis. ... kand. fiz.-mat. nauk. Moscow: MGU, 2014. 146 p. (In Russ.)
[9] Pritchett, T.M. Fast Hankel transform algorithms for optical beam propagation. Army Research Laboratory (ARL-TR-2492), 2001. 29 p. Available at: http://printfu.org/hankel
[10] Самарский А.А., Попов Ю.П. Разностные методы решения задач газовой динамики М.: Наука, 1980. 351 с.
Samarskiy, A.A., Popov, Yu.P. Finite difference methods for solving gas dynamics problems. Moscow: Nauka, 1980. 351 p. (In Russ.)
[11] Землянов А.А., Булыгин А.Д., Гейнц Ю.Э. Дифракционная оптика светового фила-мента, образованного при самофокусировке фемтосекундного лазерного импульса в воздухе // Оптика атмосферы и океана. 2011. Т. 24, № 10. С. 839-847.
Zemlyanov, A.A., Bulygin, A.D., Geints, Yu.E. Diffraction optics of a light filament generated during self-focusing of femtosecond laser pulse in air // Optika Atmosfery i Okeana. 2011. Vol. 24, No. 10. P. 839-847. (In Russ.)
[12] Cunliang Ma, Wenbin Lin. Parallel simulation for the ultra-short laser pulses' propagation in air. Available at: http://arxiv.org/pdf/1507.05988.pdf[physics.optics] (accessed 23.06.2017).
[13] Couairon, A., Brambilla, E., Corti, T., Majus, D., Ramrez-Gongora, O. de J., Kolesik, M. Practitioner's guide to laser pulse propagation models and simulation. Numerical implementation and practical usage of modern pulse propagation models // Europ. Phys. J. 2011. Vol. 199, No. 3. P. 5-76.
Поступила в 'редакцию 15 мая 2017 г., с доработки — 2 июля 2017 г.
Fully conservative numerical scheme for nonlinear Schrodinger equation with higher nonlinearities
Bulygin, Andrey D.*, Zemlyanov, Aleksandr A.
V.E. Zuev Institute of Atmospheric Optics SB RAS, Tomsk, 634055, Russia Corresponding author: Bulygin, Andrey D., e-mail: b.a.d@iao.ru
The paper examines the performance conditions to satisfy the complete conservatism property for most widely used various numerical schemes in the context of filamentation problem of the nonlinear Schrodinger equation (NLSE). Stationary radial model has two integrals of motion which are the Hamilton function and the "number of particles".
We set the explicit form of the unbalanced terms arising in asymmetrical schemes of splitting into physical factors. On the basis of numerical calculations it is understood that the standard scheme constructing the adaptive meshes with filamentation leads to a catastrophic violation of the conservation laws. We found the conditions for the numerical mesh for which it is possible to satisfy the conservation law with any acceptable accuracy for both fully symmetric scheme and the scheme with splitting on physical factors that includes the Crank —Nicolson scheme.
Keywords: filamentation, nonlinear Schrodinger equation, conservation laws, fully conservative scheme.
Acknowledgements. The work was supported by the Russian Science Foundation (grant № 16-17-10128).
Received 15 May 2017 Received in revised form 2 July 2017
© ICT SB RAS, 2017
