Трехслойная безытерационная схема повышенного порядка точности для уравнения Гинзбурга-Ландау Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 20, № 3, 2015

Трехслойная безытерационная схема повышенного порядка точности для уравнения Гинзбурга — Ландау

В. И. ПААСОНЕН1'2'*, М.П. ФЕДОРУК1'2

1Институт вычислительных технологий СО РАН, Новосибирск, Россия 2 Новосибирский государственный университет, Россия *Контактный e-mail: paas@ict.nsc.ru

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

Ключевые слова: уравнение Шрёдингера, уравнение Гинзбурга —Ландау, компактная разностная схема, повышенный порядок точности, нелинейная волоконная оптика, волоконный лазер.

1. Постановка задачи

Рассматривается задача Дирихле или периодическая начально-краевая задача для одномерного уравнения Гинзбурга — Ландау [1, 2]. В нелинейной оптике оно обычно используется в форме

ЯТТ ПЯ2 тт Я2тт

г— + - — + |и12и = Ш + ге1и12и + + (^ - *^ (1)

с комплексной искомой функцией и, зависящей от эволюционной переменной ¿, которая имеет смысл нормированной длины распространения, и "медленного" времени х. Здесь г — мнимая единица, коэффициент И = 1 или И = —1 в зависимости от того, нормальна или аномальна дисперсия, а коэффициенты 8, е, в правой части веще-

ственны, причем для корректности задачи требуется неотрицательность коэффициента [3. В частном случае нулевых коэффициентов в правой части уравнение (1) переходит в нелинейное уравнение Шрёдингера.

Уравнения Шрёдингера и уравнения Гинзбурга — Ландау моделируют процессы, происходящие в волоконно-оптических линях связи, волоконных лазерах и других оптических устройствах (см., например, [2-4]). Типичные решения таких уравнений имеют вид функции с обширным фоном плавного изменения и отдельных сигналов в виде со-литонов, сосредоточенных в узких зонах. Ввиду наличия зон больших градиентов для

ИВТ СО РАН, 2015

успешного применения традиционных разностных схем необходима настолько детальная сетка, что процесс вычислений решения даже на современных компьютерах требует больших вычислительных затрат. Для решения таких задач чаще, чем разностные схемы, используются спектральные методы, основанные на применении быстрого преобразования Фурье [5, 6]. Они позволяют получить более высокую точность, однако в случае больших задач возникают серьезные затруднения при попытках параллельной реализации. Разностные схемы, напротив, очень перспективны именно в отношении распараллеливания [7], но имеют значительно более низкую точность по сравнению с алгоритмами, основанными на спектральных методах. Поэтому одним из актуальных путей для преодоления затруднений при использовании разностных схем и для повышения эффективности алгоритмов расчета представляется повышение порядка точности разностных схем. Исследования компактных разностных схем [8-10] применительно к уравнению Шрёдингера убеждают в перспективности этого направления.

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

Настоящая работа посвящена обобщению безытерационной компактной разностной схемы [10] на случай уравнения Гинзбурга — Ландау, оценке устойчивости и численному исследованию метода в сравнении со схемой Кранка — Николсон второго порядка аппроксимации.

2. Разностные схемы

Для удобства уравнение (1) представим в канонической форме, умножив его на —г и объединив слагаемые правой части. В результате получим уравнение

-тт -2тт

-^ = + №), Ци) = д(и )и, д(и ) = 6 + г\и |2 + ш\и |4, (2)

где

о

а = р + — г, г = £ + г, ш = ^ + иг.

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

Пусть к иг — шаги равномерной сетки по ж и Ь соответственно, а Л — обычный разностный оператор, аппроксимирующий двойное дифференцирование по х. Тогда для уравнения (2) разностная схема с весами имеет вид

Общая формулировка (3) включает в себя различные схемы. Например, при равных весах (а = 0.5) имеем схему Кранка — Николсон, а при специальном значении веса схемы

11 аТ /А\

" = 2- ш- к = F■ (4)

и специальной правой части

fn+1 i fn ь/2

Fn = í_—1__l _ Л f n

2 12 J

получаем классическую схему [12], имеющую погрешность аппроксимации 0(т2 + h4). Но правая часть зависит от решения, поэтому схема для своей реализации требует итераций по нелинейности на каждом шаге. Свободна от этого недостатка трехслойная схема

uп+1 _ ттп—1

--—U-= аЛ(аип+1 + (1 -a)Un-1) + Fn, (5)

2

которая получается из (3) удвоением шага т и изменением способа аппроксимации Fn: правая часть или полностью берется со среднего слоя

Fn = г + 12 ЛГ, f(U) = (S + z\U |2 + w\U |4 )U,

или из нее выделяется линейный сомножитель U, который записывается как полусумма на верхнем и нижнем слоях, чем также достигается и линейность схемы относительно решения на верхнем слое, и второй порядок аппроксимации по т. При этом в выражении (4) для веса схемы т следует также заменить на 2т, так как в схеме (5) шаг по эволюционной переменной удвоенный. Аналогично модифицируем и схему Кранка — Николсон, превращая ее в трехслойную.

В работе [10] в трехслойной безытерационной схеме для усиления устойчивости к весу а добавлялось слагаемое ст:

1 1 т^ ат

а =2+ - Ш- К = h2-

со свободным параметром , ответственным за искусственную диссипацию и позволяющим дополнительно улучшить результаты расчетов без понижения порядка точности. Это делалось по той причине, что коэффициент а при второй производной в случае уравнения Шрёдингера был чисто мнимый, и это порождало недиссипативную схему. В случае уравнения Гинзбурга — Ландау в этом мало смысла, поскольку Re а = [3 > 0 обеспечивает диссипативность. Но при [3 = 0 введение искусственной диссипации остается полезным.

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

v"+1-U0 = + (1 - ,)U») + ЛМГ!) + h2A/(U

т 2 12

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

3. Замечания об устойчивости схем

Исследуем устойчивость двухслойной схемы в линейном приближении. Результат в равной степени относится также и к трехслойной схеме, так как она в отсутствие нелинейной правой части представляет собой запись той же двухслойной схемы, но с двойным шагом. Дисперсионный анализ устойчивости приводит к выражению для множителя возрастания гармоники ехр(гкх):

1 — 4К(1 — а)Л 2 кк ат

р = 1 + 4к.Л' • Л = 8т2"- К = к- (6)

р

В случае уравнения Гинзбурга — Ландау а = Р + г—, и это отличает его от уравнения теплопроводности, в котором коэффициент теплопроводности вещественный, и от уравнения Шрёдингера, в котором коэффициент при второй производной чисто мнимый. Кроме того, в компактной схеме вес а также является комплексным, так как он связан с безразмерным параметром К алгебраическим соотношением (4).

Поясним схематично способ исследования неравенства \ р\ < 1 и приведем результат, опуская громоздкие выкладки. Записав выражения для а и К в виде комбинации вещественной и мнимой частей и выполнив операции умножения в числителе и знаменателе дроби (6), представим его в виде отношения комплексных чисел

г + г д

Р

R + iQ

Тогда исследуемое неравенство эквивалентно тому, что квадрат модуля знаменателя не меньше квадрата модуля числителя, т. е.

( R - r)(R + г) + (Q - q)(Q + q) > 0.

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

2(2 Rea - 1)\К|2A + ReK > 0. (7)

Когда к пробегает все номера гармоник, A пробегает, согласно его определению в (6), дискретное множество точек отрезка (0,1). Но неравенство (7) линейно по A, поэтому условие невозрастания всех гармоник требует выполнения неравенства (7) лишь на концах отрезка, т. е. в нуле и единице. Отсюда ввиду положительности [3 следует условие устойчивости

2(2 Re a - 1)\К\2 + Re К > 0. (8)

В случае компактной схемы в силу равенства (4) в (8) необходимо подставить выражение

1 1 КеК Пеа =2 — 12 .

Легко проверить, что в результате подстановки неравенство оказывается истинным, поэтому компактная схема для уравнения Гинзбурга — Ландау также абсолютно устойчива в линейном приближении.

4. Результаты численных экспериментов

Все точные решения уравнения Гинзбурга—Ландау, использованные в данном параграфе в качестве материала для тестирования разностных методов, взяты из обширной коллекции аналитических решений, построенных в [1].

Тестовая задача 1. Сначала схемы тестировались на точном решении уравнения (1), не содержащего слагаемых пятой степени, т.е. при р = V = 0. Согласно [1] точное решение имеет вид

и(х, Ь) = а(х) ехр \г(ф(х) — шЬ)], (9)

где ш — вещественная постоянная, а и ф — вещественные функции, связанные равенством ф(х) = (11п (а(х)), а d — так называемый параметр чирпа. Условия совместности двух обыкновенных дифференциальных уравнений для вещественной и мнимой компонент решения удовлетворяются [1] при следующих соотношениях между коэффициентами уравнения (1):

уДЕ^+М — 3 5

2 Н

5(1 — А2 + 4Р ¿)

ш = —- ,

2 К '

где

5 = 1 + 2е Р, Н = 2р — е, К = А — р + р А2. При этом функция а( х) однозначно определяется в явном виде

а(х) = ВС 8веЬ(Бх), (10)

где

3

С.,;м' + |р■), В

2 Н '

Задача решалась в области (—2 < х < 2) х (0 < £ < 1) при значениях параметров

Б = 1, 5=1, р = 1, е = 0.25

на сгущающихся вложенных сетках, начиная с сетки 100 х 400, с уменьшением шага по к вдвое и по г в четыре раза, с тем чтобы обеспечить постоянство отношения т/к2. Здесь и ниже на рисунках всюду N означает число шагов по х, а кривая N = 0 представляет точное решение.

На рис. 1 приведены результаты расчетов тестовой задачи 1 по трехслойной модификации схемы Кранка — Николсон и по компактной схеме (5) соответственно на конечный момент эволюционной переменной Ь = 1. Из рис. 1, а видно, что при детализации сетки

на интервале (—2 < х < 2) с числом шагов N = 200, N = 400 и т. д. численное решение по схеме Кранка — Николсон сходится к точному довольно медленно. На рис. 1, б приведены расчеты по компактной схеме только на сетках с N = 100 и N = 200, так как результаты расчета на более детальных сетках уже практически сливаются с точным решением.

В табл. 1 для обеих схем приведены нормы относительных ошибок в процентах. Видно, что для достижения точности около 3 % достаточно детальной для схемы второго порядка оказалась сетка с N = 3200 шагами. В противоположность ей компактная схе-

а

б

Таблица 1. Относительная ошибка в серии расчетов тестовой задачи 1, %

Схема Число шагов N

100 200 400 800 1600 3200

0{т2 + к2) 162.8 86.26 133.43 37.41 9.47 2.72

0(т2 + к4) 85.25 5.52 0.35 2.16Е-2 1.35Е-3 8.41Е-5

ма дает превосходные результаты с погрешностью менее 0.4 % уже на существенно более грубой сетке с N = 400 шагами. Заметим также, что соотношение погрешностей при каждом дроблении сетки (приблизительно 2Ч для схем точности 0(т2 + Ьч)) подтверждает соответствие результатов теоретически ожидаемым порядкам аппроксимации. Это в равной степени относится также и к последующим тестовым задачам.

Тестовая задача 2. Точное решение в предположении 8 = 0 при отсутствии слагаемых с нелинейностью пятой степени имеет тот же явный вид (9), (10), но входные параметры связаны иными соотношениями [1]:

„=— * 1±4?!Я2, с- .!*^^

/

2(3 V 2е

* = — 1 е = в' ' ' 32 — 1 (11) * 2(3 ' £ 3 4+18(32 ' (11)

где В — свободный параметр.

Задача решалась в области (—5 < х < 5) х (0 < £ < 20) при значениях параметров

в = 1, в = 1, ( = 1.

На рис. 2 приведены результаты расчетов по схемам 0(т2 + Ь2) и 0(т2 + Ь4), полученные на различных сетках, начиная с сетки 40 х 200. Из рис. 2, а видно, что при детализации сетки численное решение, полученное по модифицированной схеме Кран-ка — Николсон, сходится к точному так же медленно, как в первом тесте. Компактная же схема (рис. 2, б) благодаря значительному преимуществу в точности показывает на сетке 80 х 800 результаты, которые по схеме Кранка — Николсон не достигаются еще на сетке 320 х 12 800. В табл. 2 по результатам расчетов данного теста приведены нормы относительных ошибок в процентах.

Тестовая задача 3. В этой задаче с нелинейностью пятой степени точное решение [1], как и во второй задаче, определяется формулой (9) при значениях параметров

1 ± 4 32

6 = 0, ш = +4( В, 2 3

Таблица 2. Относительная ошибка в серии расчетов тестовой задачи 2, %

Схема Число шагов N

40 80 160 320 640

0{т2 + к2) 181.95 109.47 32.33 8.37 2.11

0(т 2к4) 94.22 4.57 0.28 1.76Е-2 1.10Е-3

а с1, Р и £ удовлетворяют соотношениям (11), при этом V и р не являются свободными, они связаны равенством

2 V р

8р<! — в? + 3 3Р — 2(1 — рй2 При сформулированных условиях функция а имеет вид

а(х) =

/

3^(1 + 4р 2)В

2р — £ + С сс8Ь(2л/Вх)'

........ N = 320

------N = 160

+ N = 80

* N = 40

- N = 0

********* * *

-5

а

2

0

5

где

' 18 д? и(1 + 4р 2)В

/

С =х(2р — е)2 +

8р(1 — в? + 3

а В — свободный параметр.

Задача решалась в области (—10 < х < 10) х (0 < Ь < 40) при значениях параметров:

Б = 1, р = 0.1, и = 0.5, В = 1.

а

б

-1 0 0 10 -10 0 10

Таблица 3. Относительная ошибка в серии расчетов тестовой задачи 3, %

Схема Число шагов N

40 80 160 320 640

0(т2 + h2) 101.5 17.48 4.72 1.21 0.305

0(т 2 + h4) 72.68 3.58 2.21E-1 1.38E-2 8.65E-4

Сетка, начиная с размера 40 х 100, многократно дробилась вдвое по ж и вчетверо по t. На рис. 3 приведены результаты расчетов по схемам 0(т2 + h2) и 0(т2 + h4), полученные на различных сетках. Соответствующие относительные ошибки численных решений представлены в табл. 3.

Из табл. 3, в частности, следует, что по схеме Кранка — Николсон на сетке 640 х 20 480 и по компактной схеме на сетке 160 х 1280 получены приблизительно одинаковые результаты, что убедительно свидетельствует о предпочтительности компактной схемы.

Таким образом, аналитически с помощью линейного дисперсионного анализа и численно в сериях численных экспериментов исследована компактная разностная схема точности 0(т2 + h4) для уравнения Гинзбурга — Ландау, аналогичная схеме, ранее предложенной авторами для нелинейного уравнения Шрёдингера. На ряде типичных существенно нелинейных задач с известными точными решениями в виде солитонов компактная схема обнаружила высокую эффективность, значительное преимущество по точности в сравнении со схемой второго порядка аппроксимации и соответствие практически наблюдаемого и теоретического порядка точности.

Благодарности. Работа выполнена при поддержке Российского научного фонда (Соглашение № 14-21-00110).

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

[1] Akhmediev, N.N., Afanasiev, V.V. Singularities and special soliton solutions of the cubic-quintic complex Ginsburg —Landau equation // Physical Review E. 1996. Vol. 53, No. 1. P. 1190-1201.

[2] Кившарь Ю.С., Агравал Г.П. Оптические солитоны. От волоконных световодов к фотонным кристаллам. М.: Физматлит, 2005. 647 c.

Kiwshar, Yu.S., Agrawal, G.P. Optical solitons. From fiber light guides to photon crystals. Moscow: Fizmatlit, 2005. 647 p. (in Russ.)

[3] Agrawal, G.P. Nonlinear fiber optics. N.Y.: Academic Press, 2001. 446 p.

[4] Agrawal, G.P. Aplications of nonlinear fiber optics. N.Y.: Academic Press, 2001. 458 p.

[5] Degond, P., Jin, S., Tang, M. On the time splitting spectral method for the complex ginzburg-landau equation in the large time and space scale limit // Siam Journal on Scientific Computing. 2008. Vol. 30. P. 2466-2487.

[7] Паасонен В.И. Параллельный алгоритм для компактных схем в неоднородных областях // Вычисл. технологии. 2003. T. 8, № 3. C. 98-106.

Paasonen, V.I. The parallel algorithm for the compact schemes in inhomogeneous domains // Computational Technologies. 2003. Vol. 8, No 3. P. 98—106. (in Russ.)

[8] Shu Sen Xie, Guang Xing Li, Sucheol Yi Compact finite difference schemes with high accuracy for one-dimensional nonlinear Schrodinger equation b,2 // Computer Methods in Applied Mechanics and Engineering. 2009. Vol. 198. P. 1052-1061.

[9] Паасонен В.И., Федорук М.П. Компактная диссипативная схема для нелинейного уравнения Шрёдингера // Вычисл. технологии. 2011. T 16, № 6. C. 68-73. Paasonen, V.I., Fedoruk, M.P. The compact dissipative scheme for the nonlinear Schrodinger equation // Computational Technologies. 2011. Vol. 16, No. 6. P. 68-73. (in Russ.)

[10] Паасонен В.И., Федорук М.П. Компактная безытерационная схема для нелинейного уравнения Шрёдингера // Вычисл. технологии. 2012. T. 17, № 3. C. 83-90. Paasonen, V.I., Fedoruk, M.P. The compact noniterative scheme with artificial dissipation for nonlinear Schrodinger equation // Computational Technologies. 2012. Vol. 17, No. 3. P. 8390. (in Russ.)

[11] Остапенко В.В. О сходимости разностных схем за фронтом нестационарной ударной волны // Журн. вычисл. математики и мат. физики. 1997. T. 37, № 10. C. 1201-1212. Ostapenko, V.V. On the convergence of difference schemes behind non-stationary shock wave front // Computational Mathematics and Mathematical Physics. 1997. Vol. 37, No 10. P. 1201-1212. (in Russ.)

[12] Микеладзе Ш!.Е. О численном интегрировании уравнений эллиптического и параболического типов // Известия АН СССР. Серия матем. 1941. Т. 5, № 1. С. 57-74. Mikeladze, Sh.E. On the numerical integration of the equations of elliptic and parabolic types // Izvestiya AN SSSR, Seriya mathematika. 1941. Vol. 5, No. 1. P. 57-74. (in Russ.)

Поступила в 'редакцию 20 апреля 2015 г., с доработки — 22 мая 2015 г.

Three-level non-iterative high accuracy scheme for Ginzburg — Landau equation

Paasonen, Viktor I.1,2'*, Fedoruk, Mikhail P.1,2

institute of Computational Technologies SB RAS, Novosibirsk, 630090, Russia 2Novosibirsk State University, Novosibirsk, 630090, Russia * Corresponding author: Paasonen, Viktor I., e-mail: paas@ict.nsc.ru

This paper presents a generalization of compact difference scheme previously developed and investigated by the authors for one-dimensional nonlinear Schrodinger equation to the case of the Ginzburg —Landau equation. We apply the increased order of accuracy as a very useful tool for improvement of the quality of calculations. The scheme approximates the Ginzburg —Landau equation with the second-order for the evolutionary variable and with the fourth order for the "slow time". The scheme is essentially of a two-level, but the scheme uses a double step at three levels, with the non-linear approximation on the right side at the middle level. This approach do not

require the need to iterate on the nonlinearity at all steps, except the first, and thus saves computing resources. To compute the solution in the first step we proposed to use a two-level iterative scheme providing the same order of approximation as for the main scheme. Stability of difference schemes was examined in the linear approximation using the analysis of variance of the behavior harmonics. The paper presents the results of calculations on the sequence of the test problems which use the condensing grids with known exact solutions obtained earlier by Ahmediev and Afanasiev in their famous work. A comparison with the results obtained by the three-layered modified Crank — Nicolson scheme with approximation of the second order. The above graphic and tabular material evidence of significant advantages of a compact difference scheme.

Keywords: Shrodinger equation, Ginzburg —Landau equation, compact difference scheme, high order accuracy, nonlinear fiber optics, fiber laser.

Acknowledgements. The research is supported by the Russian Science Foundation. Agreement No. 14-21-00110.

Received 20 April 2015 Received in revised form 22 May 2015