CC BY
68
Вычислительные технологии
Том 16, № 6, 2011
Компактная диссипативная схема для нелинейного уравнения Шредингера*
В. И. Паасонен, М. II. Федорук Институт вычислительных технологий СО РАН, Новосибирск, Россия Новосибирский государственный университет, Россия e-mail: paas@ict .nsc.ru, mif e@ict .nsc.ru
Для нелинейного уравнения Шредингера построена компактная разностная схема четвертого порядка точности по временной переменной. Кроме того, для компактной разностной схемы разработан механизм регулируемой диссипации. В серии численных экспериментов показано ее существенное преимущество в точности в сравнении с диссипативным аналогом схемы Кранка — Николсон.
Ключевые слова: уравнение Шредингера, компактная разностная схема, диссипативная схема, повышенный порядок точности, нелинейная волоконная оптика, волоконный лазер, волоконно-оптическая линия связи.
1. Постановка задачи
Рассматривается задача Дирихле или периодическая начально-краевая задача для одномерного нелинейного уравнения Шредингера [1, 2]
с комплексной искомой функцией А, зависящей от циклической переменной г и времени ¿. В левой части г — мнимая единица, коэффициенты Ь и в в правой части вещественны.
Данное уравнение является основным уравнением нелинейной волоконной оптики и активно используется, например, для математического моделирования волоконно-оптических линий связи, волоконных лазеров, различных оптических устройств [2, 3]. В настоящее время для численного решения этого уравнения, как правило, используется метод расщепления по физическим процессам (см., например, [2]). Этот метод позволяет получить приближенное решение в предположении, что при распространении оптического сигнала на малое расстояние д по эволюционной переменной г нелинейность и дисперсия действуют независимо. Относительно большая скорость реализации данного алгоритма N где N — количество точек по временной переменной)
по сравнению с большинством методов конечных разностей связана с использованием алгоритма быстрого преобразования Фурье (БПФ) для решения линейной части
*Работа выполнена при финансовой поддержке Министерства образования и науки РФ (Гос. контракт № 02.740.11.5129), РФФИ (грант № 11-01-00294-а) и Правительства РФ для государственной поддержки научных исследований, проводимых под руководством ведущих ученых (грант № 11.G34.31.0035).
исходного уравнения. Однако для адекватного численного решения задач в области волоконно-оптических линий связи или волоконных лазеров типичное число точек по временной переменной может составлять N = 106 — 107, что приводит к огромным затратам машинного времени и к невозможности решения подобных задач даже на самых мощных персональных компьютерах. Единственным реальным путем расчета этих задач является параллельная реализация численных алгоритмов и использование высокопроизводительных вычислительных систем. Вместе с тем хорошо известно, что алгоритмы БПФ обладают низкой эффективностью распараллеливания,
С другой стороны, конечно-разностные схемы легко допускают параллельную реализацию [4, 5], хотя значительно уступают методам, основанным на БПФ, в точности расчета на фиксированной сетке по временной переменной, В связи с этим весьма актуальной является задача построения компактных конечно-разностных схем повышенного порядка аппроксимации, допускающих эффективную параллельную реализацию. Цель данного исследования состоит в разработке более экономичной компактной схемы четвертого порядка точности по временной переменной и сравнении ее с модифицированной схемой Кранка — Николеон на типичных линейных и нелинейных тестах.
Хотя рассматриваемое уравнение не является параболическим, формально оно очень напоминает обычное уравнение теплопроводности, и при выборе схем для расчета воспользуемся классическим опытом построения экономичных схем [6] для параболического уравнения.
2. Разностные схемы
Пусть g и т суть шаги равномерной сетки по z и t соответственно, а Л — обычный разностный оператор, аппроксимирующий двойное дифференцирование по времени. Тогда для уравнения Шредингера разностная схема с весами имеет вид
An+1 _ An b
г-= -A(aAn+1 + (1 - a)An) - s(afn+1 + (1 - a)fn), f = |A|2A
q 2
Известно, что в параболическом случае (когда слева нет мнимой единицы) однородная схема абсолютно устойчива при a > 1/2, т.е. множители возрастания гармоник
1 + r(a — 1\фк 2bq 2 кт
Рк= 1 + гафк ' Г = —' ^ = smT
по модулю не превышают единицы при любом соотношении шагов, а при меньших весах схема, в том числе и явная, условно устойчива. При этом важно отметить, что устойчивость именно сильная, поскольку неравенства, ограничивающие рост гармоник, являются строгими,
В нашем случае из дисперсионного соотношения схемы легко заметить, что при весе a > 1/2 имеет место абсолютная устойчивость, так как множители
i + r(a — 1 )фк 2bq . . 2 кт
Рк =---1-, г = —, фк = sin —
i + гафк т2 2
по модулю не превышают единицы, но при меньших весах устойчивости нет вовсе, т, е, область условной устойчивости отсутствует. Заметим, что наибольший интерес, очевидно, представляют чисто неявная схема с единичным весом и схема Кранка —Николеон
с равными весами. Первая имеет наибольший запас устойчивости, но только первый порядок по циклической переменной, в то время как вторая имеет второй порядок по обеим переменным. Однако в случае а = 1/2 множитель возрастания гармоник по модулю равен точно единице, что означает отсутствие сильной устойчивости, В некоторых типичных случаях решений с большими градиентами при расчете на значительное число циклов слабой устойчивости может быть не достаточно для получения надежного результата, В целях устранения этого недостатка целесообразно модифицировать схему Кранка —Николсон, взяв вес чуть больше половины (а = 1/2 + сц, с > 0) с отличием па величину порядка шага по циклической переменной. Это обеспечивает сильную устойчивость при сохранении второго порядка точности. Если зафиксировать специальное значение веса
2 3т т2
и специальным же образом [6] сформировать правую часть разностной схемы, приходим к компактной схеме повышенной точности 0(ц2 + т4):
Ап+1 _ Ап Ь (¡п+1 + {
2
. 5 2Л+ (1 - а)Ап) - 8 ^ 2' ' + -Л/^ , / = \А\2А.
Для однородного уравнения теплопроводности (г = 1) гармонический анализ указывает на абсолютную и сильную устойчивость схемы повышенной точности, В этом случае множители возрастания гармоник
1 _ (1/3 + т/2)фк
Рк
1 _(1/3_ т/2)фк
строго меньше единицы.
Однако в нашем случае ситуация не такая благоприятная. Из дисперсионного соотношения для к-й гармоники получаем коэффициент возрастания в виде отношения комплексно-сопряженных величин
= (1 - -фк/3)г - г/2фк Рк (1 — фк/3) + т/2фк '
Следовательно, хоть устойчивость и абсолютна, но сильной устойчивости нет, как и в случае схемы Кранка —Николсон, Меняет ситуацию слабое увеличение веса схемы
1 г ц г = 26—, с>0.
2 3т т2
Это изменение эквивалентно добавлению в уравнение искусственного диссипативного слагаемого
2 дА ах
по порядку малости находящегося в пределах погрешности аппроксимации. Таким образом, получаем компактную схему точности 0(ц2 + т4) с сильной устойчивостью. Действительно, квадрат модуля для коэффициента возрастания гармоники
2 = 1 _ 2 сдгЩ к \г + гафк\2
оказывается строго меньшим единицы.
3. Некоторые результаты расчетов
В простейшем случае линейного уравнения Шредингера, когда 5 = 0, разностные схемы реализуются непосредственно трехточечной прогонкой на каждом шаге по г, В случае неоднородного уравнения (в = 0) приведенные схемы имеют кубичную нелинейность, присутствующую и на верхнем слое, что требует организации внутренних итераций по нелинейности. Представляется целесообразным применить простой способ линеаризации, подобный двухшаговой процедуре, применяемой в схемах предиктор-корректор. Для каждого фиксированного слоя п по эволюционной переменной определим последовательность у0,у1,...,ук,... приближений решения на следующем (п + 1)-м слое, В качестве начального приближения используем, например, решение на верхнем слое, полученное по явной схеме. От текущего приближения ук к следующему ^к+1 переходим с помощью линейной схемы, В случае схемы повышенной точности итерационный процесс имеет вид
гд х ^к+1 — Ап Ъд п \ьк¡V + \ЛП\2ЛП т2 .п.2 „Л А
?---Л---КАП + ^ ^—1-1—1--1--к\Ап\2Ап = 0.
4 ] д 2 V 2 12 1 1 у
Можно рассматривать также вариант представления функции f на верхнем слое в виде \ук\2ук+1, это несколько меняет коэффициенты трехдиагональной матрицы. Совершенно аналогично строится итерационный процесс для схемы с весами. На практике в отдельных вариантах снимались верхние ограничения по числу внутренних итераций, чтобы убедиться в их сходимости, а в основных расчетах достаточно ограничиться несколькими итерациями с выходом из итерационного цикла по сходимости. Сначала схемы тестировались на точном решении
= —.= ехр ( — ) , ги = Т02 — ъЬг
однородного уравнения (в = 0) при значениях пара метров Ъ = —1, с = 0.01 Т0 = 1.8, Задача решалась в области (0 < г < 10) х (—10 < £ < 10), Результаты расчетов приведены в табл. 1, Здесь в первых двух столбцах приводятся размеры сетки по временной и циклической переменной, третий и пятый столбцы содержат глобальную С-норму ошибки численного решения
5 = тах \Ап — А(гп,^)\
по схемам второго и четвертого порядка точности, а в четвертом и шестом столбцах приведены коэффициенты убывания ошибки при детализации сетки для рассматриваемых схем. Полученные значения коэффициентов свидетельствуют о совпадении теоретических и практически наблюдаемых порядков сходимости. Очевидно явное преимущество компактной схемы,
В табл. 1 приведены также результаты расчета нелинейного уравнения с параметром в = 1 в области (0 < г < 4п) х (—10 < £ < 10), Точное решение задачи
= ехр {Ц^ 8еск{Ь)
имеет вид солитона.
Таблица 1. Результаты расчета линейной и нелинейной задач на сгущающихся сетках по схемам второго и четвертого порядка точности
Nt Nz ¿2 к2 ¿4 кА
Линейная задача
20 20 8.7144-02 1.2963-02
40 80 2.2095-02 3.94 8.2559-04 15.69
80 320 5.4733-03 4.02 5.1446-05 16.05
160 1280 1.3645-03 4.02 3.2122-06 16.02
320 5120 3.4089-04 4.00 2.0072-07 16.05
640 20480 8.5210-05 4.00 1.2545-08 15.99
Нелинейная задача
20 20 1.4643-00 7.6688-01
40 80 7.1275-01 2.06 6.6585-02 11.52
80 320 2.2153-01 3.22 6.2449-03 10.66
160 1280 7.6367-02 2.90 5.5693-04 11.21
320 5120 2.6835-02 2.85 4.9901-05 11.16
640 20480 9.4732-03 2.83 4.3687-06 11.42
Таблица 2. Зависимость ошибки расчета от параметра диссипации компактной схемы
с 0.0 0.0001 0.001 0.01 0.1
Сетка 320 х 1280 Сетка 640 х 2360 0.0444 0.0105 0.0440 0.0103 0.0399 0.0087 0.0021 0.0081 0.4057 0.1763
В данном примере внутренние итерации проводились до различия между итерациями е = 10-8 с верхним ограничением до пяти итераций. Практически же для такого порога сходимости не требовалось более пяти итераций на шаге. Из результатов расчета следует, что коэффициенты убывания ошибки для обеих схем в нелинейном случае несколько ниже, чем в линейном, что объясняется влиянием ошибок округления. Анализ коэффициентов показывает, что практически наблюдаемые порядки сходимости для схем второго и четвертого порядка аппроксимации для нелинейной задачи составляют примерно величины 1.6 и 3.5 вместо 2 и 4,
С целью определения влияния на ошибку величины коэффициента диссипации с проведена серия расчетов этой же нелинейной задачи, но на большем интервале по циклической переменной (0 < г < 20п), Из результатов расчета (табл. 2) видно, что при возрастании параметра с ошибка 84 сначала падает, а затем возрастает, Следова-
с
умеренная искусственная диссипация в сравнении с полным ее отсутствием улучшает расчет. Экспериментально установлено, что во всех расчетах наилучшее значение коэффициента искусственной диссипации находится в пределах 0.001 ч—0.01,
Список литературы
[1] Кившарь Ю.С., Агравал Г.П. Оптические солитоны. От волоконных световодов к фотонным кристаллам. М.: Физматлит, 2005. 647 с.
[2] Agrawal G.P. Nonlinear Fiber Optics. N.Y.: Acad. Press, 2001. 446 c.
[31 Agrawal G.P. Aplications of Nonlinear Fiber Optics. N.Y.: Acad. Press, 2001. 458 p.
[4] Яненко H.H., Коновалов A.H., Бугров A.H., Шустов Г.В. Об организации параллельных вычислений и распараллеливании прогонки // Численные методы механики сплошной среды. 1978. Т. 9, № 7. С. 136-139.
[5] Паасонен В.И. Параллельный алгоритм для компактных схем в неоднородных областях // Вычисл. технологии. 2003. Т. 8, № 3. С. 98-106.
[6] Микеладзе Ш.Е. О численном интегрировании уравнений эллиптического и параболического типов // Изв. АН СССР. Математика. 1941. Т. 5, № 1. С. 57-74.
Поступила в редакцию 12 мая 2011 г.
