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

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

Том 17, № 3, 2012

Компактная безытерационная схема с искусственной диссипацией для нелинейного уравнения Шрёдингера*

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

Для уравнения Шрёдингера построена компактная разностная схема четвёртого порядка точности с искусственной диссипацией. Схема не требует итераций по нелинейности. Ее существенное преимущество в точности в сравнении с дисси-пативным аналогом схемы Кранка — Николсон доказано в серии численных экспериментов. Кроме того, значительно сокращено время счета по сравнению с итерационной компактной схемой. Отмечено также улучшение точности схемы при удачном выборе параметра диссипации в сравнении со схемой без диссипации.

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

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

Рассматривается задача Дирихле или периодическая начально-краевая задача для одномерного нелинейного уравнения Шрёдингера (см., например, [1, 2])

4А=§£ - * /=|АГА

с комплексной искомой функцией А, зависящей от эволюционной переменной г и времени £ (здесь г — мнимая единица, коэффициенты Ь и в — вещественны). Данное уравнение, являющееся основным уравнением нелинейной волоконной оптики, моделирует процессы, происходящие в волоконно-оптических линиях связи и в различных оптических устройствах (см., например, [2, 3]). Для численного решения этого уравнения и его многомерных версий используется метод расщепления по физическим процессам (см., например, [2]), опирающийся на предположение о независимом влиянии нелинейности и дисперсии при малом масштабе изменения эволюционной переменной г. Применение быстрого преобразования Фурье (БПФ) для решения линейной части исходного уравнения даёт хорошую скорость реализации алгоритма (~ N log2 где N — количество точек по временной переменной) в сравнении с большинством традиционных методов конечных разностей. Однако для решения с достаточной точностью типичных задач

* Работа выполнена при поддержке Министерства образования и науки РФ (гос. контракт № 11.519.11.4018), гранта Правительства РФ для государственной поддержки научных исследований, проводимых под руководством ведущих ученых, № 11.G34.31.0035 и РФФИ (грант № 11-01-00294^).

волоконно-оптических линий связи или волоконных лазеров число точек по временной переменной может составлять N = 106 —107, что делает невозможным их решение даже на самых мощных персональных компьютерах.

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

2. Разностные схемы второго и четвёртого порядка точности

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

АП+1 _ Ап -

г-= - Л(аАп+1 + (1 — а)Ап) — в(а/п+1 + (1 — а)/п), / = |А|РА.

д 2

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

Заметим, что из схем обычной точности наибольший интерес, очевидно, представляет схема Кранка — Николсон с равными весами, имеющая второй порядок по обеим переменным. Однако в случае а =1/2 множитель возрастания гармоник по модулю в точности равен единице, что означает отсутствие сильной устойчивости (естественной диссипации). В некоторых типичных случаях решений с большими градиентами при расчёте на значительное число циклов для получения надёжного результата слабой устойчивости может быть не достаточно. Для устранения этого недостатка целесообразно модифицировать схему Кранка — Николсон, выбрав вес несколько больше 1/2 (а = 1/2 + сд, с > 0) с отличием на величину порядка шага по циклической переменной. Это обеспечивает сильную устойчивость при сохранении второго порядка точности.

Рассмотрим компактную схему. Формально уравнение Шрёдингера подобно уравнению теплопроводности, поэтому все известные подходы построения компактных схем так или иначе восходят к схеме Микеладзе [8] для параболического уравнения, хотя многие авторы и не ссылаются на эту основополагающую работу. Если зафиксировать специальное значение веса

1 г оа (

а =о - V' т = 2Ь1.

2 3т т2

и специальным образом [8] сформировать правую часть разностной схемы, то получим компактную схему повышенной точности 0((2 + т4)

An+1 _ An b / fn+1 + fn T2 N

i-q-= 2 Л(аАп+1 + (1 - a)An) - в[* 2+ f + -Л/п^ , / = |A|pA,

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

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

г/2фк - (1 - А/3)г

Рк -

т/2фк + (1 - А/3)г'

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

1 i г оа (

а =о+ с( - > т = 2Ь^.

2 3т т2

Это эквивалентно добавлению в правую часть уравнения искусственного диссипатив-ного слагаемого

Ь 2 д 3А 2С( Ы2дг

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

. 2с(т2-фк < _

|рк 1 = ^ -|г + таАк12 < 1.

Эта идея построения компактной схемы с внедрением в схему слабой искусственной диссипации с регулируемым множителем предложена в работе [6].

3. Анализ неоднородной схемы в линейном приближении

Проведём гармонический анализ схемы для неоднородного уравнения, когда правая часть линейна относительно искомой функции, т. е. / = wA, где w = const представляет

собой "замороженный" коэффициент в|А|2. Тогда наша схема приобретает вид

Лп+1 _ А'п - / Ап+1 + А'п т2 \

гА д А = 2Л(аАп+1 + (1 — а)Ап) — ы 2+ А + р^ЛАп) ,

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

1 i г оа д

- + сд — p—, г = 2Ь~,

2 3т т2

Множитель возрастания к-й гармоники для рассматриваемой линейной схемы имеет вид

= —тфк/2 — wд/2 + стдфк + ры дфк/3 + (1 — рфк/3)г тфк/2 + ыд/2 + стдфк + (1 — рфк/3)г '

В параболическом случае (г = 1) характер оценки для спектра схемы должен зависеть от знака правой части, а именно, при положительном ы максимум модуля собственных чисел не превышает величины 1 — тд (т > 0), а при отрицательном — величины 1 + 0(д). Аналогичная ситуация имеет место и в случае уравнения Шрёдингера. Дисперсионный анализ позволяет сделать следующие выводы о необходимых условиях сильной устойчивости схем с ненулевой линейной правой частью:

1) для модифицированной схемы Кранка — Николсон (р = 0) сильная устойчивость имеет место при совпадении знаков т и ы и неотрицательности параметра с;

2) в случае компактной схемы (р =1) для сильной устойчивости при совпадении знаков коэффициентов т и ы на параметр с накладывается менее жёсткое ограничение, чем неотрицательность, а именно

ы

_ 6т

В этом легко убедиться, проанализировав условие |рк | < 1, которое эквивалентно неравенству

(2ст + ры/3) ((т/д — ры/3)фк + ы) > 0.

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

4. Построение безытерационной компактной схемы

В простейшем случае линейного уравнения Шрёдингера при в = 0 разностные схемы реализуются непосредственно трёхточечной прогонкой на каждом шаге по г. В случае неоднородного уравнения (в = 0) приведённые схемы имеют нелинейность, присутствующую и на верхнем слое, что требует организации внутренних итераций по нелинейности. В работе [6] итерационный процесс организован обычным способом. Для каждого фиксированного слоя п по эволюционной переменной определяется последовательность V0, V1,..., ьк,... приближений решения на следующем п + 1-м слое. В качестве начального

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

(Е _ ^Л) ^ - АП _ Ьда" + , ( 1 ^ 1 ^ + I А" 1 РАП + т!Л(|А"|рА")) = 0. 4 ( 2 V 2 12 I

Аналогично строится итерационный процесс для схемы с весами и, в частности, для диссипативной модификации схемы Кранка — Николсон.

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

Для построения безытерационной схемы воспользуемся подходом, аналогичным использованному в работе [7]. С этой целью перепишем исходную компактную схему симметрично по трем слоям с двойным шагом по циклической переменной

А"+1 — А"-1 Ь / /"+1 + /"-1 т2 \ г---= 2Л(аА"+1 + (1 - а)А"-1) - я ^-+ , / = |А|РА.

Следует учесть, что шаг здесь удвоенный, поэтому в схеме изменяется определение веса

1 г (

2 3т т2

Схема содержит нелинейность на верхнем шаге только в полусумме правой части. Несколько преобразуем ее без потери точности, аппроксимируя степень модуля |А|Р на среднем слое, а сомножитель А в виде полусуммы — на крайних слоях:

£П+1 I £П—1 Лп+1 \ ЛП—1

/ +/ = |Ап|р А +А + 0((2).

В результате получим трёхслойную безытерационную схему

А"+1 _ А"-1 Ь / А"+1 + А"-1 т2 г---= 2Л(аАп+1 + (1 - а)Ап-1) - Л |Ап|р-+-+ -Л(|Ап|рА'

Линейный дисперсионный анализ трёхслойной схемы даст тот же результат, что и для двухслойной схемы, с точностью до переобозначения ( на 2(, поэтому по устойчивости данная схема не хуже исследованной выше. Таким образом, построена трёхслойная компактная схема с регулируемой диссипацией, не требующая итераций по нелинейности. Аналогично можно поступить с модифицированной схемой Кранка — Николсон, построив трёхслойную безытерационную схему второго порядка с искусственной диссипацией. Так как схема является трёхслойной, то требуется задать решение не только на нулевом, но и на первом слое по эволюционной переменной, причём с адекватной точностью, что легко решается использованием на первом шаге схемы повышенной точности [6].

Заметим, что в отличие от предлагаемой в схеме [7] нет искусственной диссипации (с = 0) и полной симметрии, поскольку разностная производная по д в [7] аппроксимируется по верхнему и среднему слоям, отчего порядок погрешности аппроксимации относительно эволюционной переменной снижается до первого, тогда как в нашем случае он остается вторым.

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

Сначала схемы тестировались на точном решении

Т0 ( ¿2 \ _ _ гр2

А(г, £) = —= ехр--, ы = ТО — гЬг

уы \ 2ы у

однородного уравнения (в = 0) при значениях параметров Ь = —1, с = 0.01, Т0 = 1.8. Задача решалась в области (0 < г < 10) х (—10 < £ < 10) с помощью двухслойных схем. Результаты расчётов приведены в табл. 1. Здесь в первых двух столбцах приводятся размеры сетки по временной и циклической переменным, третий и пятый столбцы содержат глобальную С -норму ошибки численного решения

5 = тах |Ап — А(гп,^)|

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

В табл. 2 приведены результаты расчёта нелинейного уравнения с параметром в = 1 в области (0 < г < 4п) х (—10 < £ < 10). Точное решение задачи

А(г,£) = ехр ^ ——^ йвсЬ(£)

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

Для итерационных схем внутренние итерации в данном примере проводились до достижения различия между итерациями е = 10-8 с верхним ограничением до пяти

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

N N ¿2 К2 ¿4 К4

20 20 8.7144-02 1.2963-02

40 80 2.2095-02 3.944 8.2559-04 15.702

80 320 5.4733-03 4.037 5.1446-05 16.048

160 1280 1.3645-03 4.011 3.2122-06 16.016

320 5120 3.4089-04 4.003 2.0072-07 16.003

640 20480 8.5210-05 4.001 1.2545-08 16.000

итераций. Практически же для такого порога сходимости не требовалось более трёх итераций на шаге, как это видно из табл. 3, в которой приведено время центрального процессора на расчёт одного варианта в операционной среде Ма^аЬ для рассматриваемой последовательности из шести сеток, степень детальности которых возрастает по столбцам слева направо и соответствует данным табл. 2.

Коэффициенты убывания ошибки для обеих схем в нелинейном случае оказались несколько ниже, чем в линейном (см. табл. 2). Из анализа коэффициентов следует, что практически наблюдаемые порядки сходимости для схем второго и четвёртого порядка аппроксимации на мелких сетках составляют для нелинейной задачи примерно величины 1.6 и 3.5 вместо теоретических 2 и 4. По мнению авторов, это связано не с качеством схем, а с обычным влиянием ошибок округления, которые в условиях сильной нелинейности при достаточно малых шагах сетки становятся сравнимыми с аппроксимационной ошибкой. Для подтверждения этого были проведены расчёты той же задачи на более грубых сетках. Соответствующие результаты приведены в табл. 4, из которой видно, что при детализации сетки коэффициент убывания ошибки сначала приближается к теоретическим значениям 4 и 16, а при дальнейшем уменьшении шага снижается ввиду влияния ошибок округления. Это обычное явление для любых краевых задач, но в задачах с сильной нелинейностью, как видно из расчётов, оно наступает быстрее.

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

N N2 ¿2 К2 ¿4 К4 5 К

20 20 2.4197-00 — 8.7966-01 — 1.7320-00 —

40 80 7.1275-01 3.39 6.6587-02 13.21 2.3239-02 74.53

80 320 2.2153-01 3.22 6.2449-03 10.66 1.8797-02 12.36

160 1280 7.6367-02 2.90 5.5693-04 11.21 1.6632-03 11.30

320 5120 2.6835-02 2.85 4.9901-05 11.16 1.441-05 11.31

640 20480 9.4732-03 2.83 4.3687-06 11.42 1.3006-06 11.31

Таблица 3. Время счёта на сгущающихся сетках по схемам второго и четвёртого порядка точности

0(Н2), итерационная схема 0.14 0.19 0.92 4.51 24.10 170.52

0(Л,4), итерационная схема 0.03 0.19 0.94 4.54 24.57 170.76

0(Л,2), безытерационная схема 0.05 0.12 0.51 2.48 12.82 74.54

0(Л,4), безытерационная схема 0.05 0.14 0.53 2.44 12.78 74.62

Таблица 4. Результаты расчёта нелинейной задачи по безытерационным схемам второго и четвёртого порядка точности на грубых сетках

N N2 ¿2 К2 ¿4 К4

4 8 1.4541-00 — 1.1474-00 —

8 32 2.0410-00 0.712 1.5821-00 0.725

16 128 2.4684-00 0.827 1.9784-00 0.800

32 512 1.3302-00 1.856 1.1152-01 17.740

64 2048 3.1587-01 4.211 6.6918-03 16.665

128 8192 1.0722-01 2.946 5.7731-04 11.591

Таблица 5. Зависимость ошибки решения от параметра диссипации для безытерационных схем второго и четвёртого порядка точности

С2 0 -0.2 -0.4 -0.6 -0.8 -1.0

102 х ¿2 2.65 2.37 2.08 1.79 1.52 21.67

С4 0 0.008 0.016 0.024 0.032 0.04

105 X ¿4 26.29 14.88 7.62 9.32 20.54 32.01

С целью определения влияния коэффициента диссипации на точность расчёта проведена серия численных экспериментов при различных значениях c для рассматриваемой нелинейной задачи в области (0 < z < 40п) х (-10 < t < 10). В табл. 5 показана зависимость ошибки решения от параметра диссипации c для двух безытерационных схем на сетке Nt х Nz = 320 х 1280 — для модифицированной схемы Кранка — Николсон и компактной схемы. Здесь первые две строки относятся к трёхслойной схеме второго порядка точности, а последующие — к трёхслойной безытерационной компактной схеме. Как и следовало ожидать, масштабы ошибок существенно различаются: для первой схемы оптимальный параметр c лежит в отрицательной области, для второй — в положительной, причем масштабы наилучших значений c различны. Разница в знаках показывает, что нелинейная картина сложнее, чем результаты анализа в рамках линейных приближений. Аналогичные исследования были проведены для двух вариантов двухслойных итерационных схем. Для них тоже наблюдалось существенное улучшение результатов при удачном выборе параметра диссипации c в сравнении со схемами с нулевым параметром диссипации. Таким образом, приведённые результаты указывают на полезность введения в схемы искусственной диссипативной (или антидиссипативной) составляющей со свободно выбираемым параметром. Вместе с тем вопрос о теоретических рекомендациях к выбору оптимальных параметров остается пока открытым.

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

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

[2] Agrawal G.P. Nonlinear Fiber Optics. New York: Acad. Press, 2001. 446 p.

[3] Agrawal G.P. Aplications of Nonlinear Fiber Optics. New York: Acad. Press, 2001. 458 p.

[4] Яненко Н.Н., Коновалов А.Н., Бугров А.Н., Шустов Г.В. Об организации параллельных вычислений и распараллеливании прогонки // Численные методы механики сплошной среды. 1978. Т. 9, № 7. C. 136-139.

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

[6] Федорук М.П., Паасонен В.И. Компактная диссипативная схема для нелинейного уравнения Шрёдингера // Там же. 2011. Т. 16, № 6. C. 68-73.

[7] Shü-Sen Xie, Guang-Xing Li, Sucheol Yi Compact finite difference schemes with high accuracy for one-dimensional nonlinear Schrödinger equation b,2 // Comput. Methods Appl. Mech. Engrg. 2009. Vol. 198. P. 1052-1061.

[8] Микеладзе Ш.Е. О численном интегрировании уравнений эллиптического и параболического типов // Изв. АН СССР. Математика. 1941. T. 5, № 1. C. 57-74.

Поступила в 'редакцию 8 декабря 2011 г., с доработки —17 апреля 2012 г.