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

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

Том 20, № 2, 2015

Компактные схемы третьего порядка точности на неравномерных адаптивных сетках

В.И. ПААСОНЕН

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

Для одномерного нелинейного уравнения Шредингера и уравнения теплопроводности рассматриваются компактные разностные схемы третьего порядка точности на неравномерных сетках с аппроксимацией второго порядка по эволюционной переменной. При численной реализации схем применялась динамическая адаптация сетки к решению с интерполяцией решения на новую сетку на каждом шаге. Проведено численное исследование созданных алгоритмов в сравнении с известными методами.

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

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

Рассматривается задача Дирихле для нелинейного уравнения Шредингера (см., напри-

где £ — линейная (эволюционная) переменная, а х — временная переменная, с = г — мнимая единица, а f = т\и\ри — правая часть специального вида, и — комплексная искомая функция.

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

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

2])

c dt а dx2 f,

dU д 2U

(1)

© ИВТ СО РАН, 2015

Другим инструментом улучшения эффективности численного метода может быть повышение порядка точности расчетов путем использования компактных схем. В предположении равномерности сетки компактные схемы четвертого порядка точности известны как для уравнения теплопроводности [3, 4], так и для уравнения Шредингера [5]. В работах [6, 7] компактные схемы для уравнения Шредингера усовершенствованы путем введения в них регулируемой диссипации. В сравнении с традиционными схемами компактные схемы даже на равномерных сетках обнаружили существенное преимущество в точности. С другой стороны, в работе [8] на примере уравнения Пуассона и уравнений Навье — Стокса продемонстрирована возможность построения компактных схем также и на неравномерных сетках. А в связи с вышесказанными соображениями возможность обобщения компактных схем для уравнения Шредингера на случай неравномерных сеток представляется дополнительным резервом для улучшения экономичности алгоритмов.

Рассматриваемые процессы являются существенно нестационарными, в связи с чем высокоточную технологию расчета на неравномерных сетках необходимо снабдить механизмом динамической адаптации сетки к решению, позволяющим обеспечить приблизительную равномерность локальной погрешности схемы на всем промежутке изменения переменной. Известны два основных подхода к решению проблемы адаптации сетки. Первый заключается в подборе подходящей, часто нелинейной, замены пространственной переменной, позволяющей благоприятным образом изменить линейные масштабы в зоне быстрого изменения функции, используя при этом в новых переменных равномерную сетку. Этот способ используется в различных вариациях уже около полувека, в основном в задачах численного моделирования пограничного слоя (см., например, [9]). Второй способ предполагает динамическую перестройку сетки по профилю решения на каждом временном слое с последующей интерполяцией решения на новую сетку. Алгоритмы построения новой сетки могут быть самыми разнообразными, например, по длине дуги кривой [10], по норме погрешности [11], по скорости изменения решения [12]. В случае компактных схем адаптация сетки, опирающаяся на замену переменной, вызывает затруднения, так как нелинейность в дифференциальном операторе является существенным препятствием для построения высокоточной схемы на компактном шаблоне. Поэтому предпочтение отдано построению на каждом эволюционном шаге новой сетки, квазиравномерной по приращению решения, по методу, близкому к [12], в сочетании с локальной интерполяцией сеточной функции с достаточным порядком точности.

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

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

2. Компактная разностная схема третьего порядка точности на неравномерной сетке

Зафиксируем произвольно узел неравномерной сетки по переменной х и обозначим через к+ и к- локальные значения шагов сетки справа и слева от этого узла, а через Д+ и Д- — соответствующие простейшие разделенные разности "вперед" и "назад". Введем также обозначения для разности, суммы и произведения соседних шагов к+ и к-:

( = к+ — к-, в = к+ + к-, р = к+к-.

Тогда для разностных аппроксимаций

Д= к-Л+ + ■■ -, Л

к-Д+ + к+Д- Д+ — Д_

в/2

операторов одно- и двукратного дифференцирования по х справедливы разложения

Дш = ш' + рш''' + О (к33), 6

Ль = ш'' + (ш''' + Ьш'''' + О(к3), Ь = (-+р,

3 12

где к — максимальный шаг.

Представим формально наше уравнение (1) в виде ОДУ

п92и = Р Р = 7 + г9и (2)

П~дХ? = Р Р = 7 + (2)

не обращая пока внимания на зависимость функций от £ и на конкретный вид правой части Р, зависящей от производной. Тогда по аналогии с рассуждениями работы [8] и в силу выписанных выше разложений легко построить трехточечную компактную разностную схему третьего порядка аппроксимации по переменной х:

(

пЛи = ЕР, Е = Е + -Д + ЬЛ.

' 3

Подставляя сюда выражение для Р из (2), а затем проводя дискретизацию эволюционной переменной £ и осредняя функции симметрично по двум соседним слоям, получим компактную разностную схему

сЕ и"+' — и" = пЛ и"+' + и" — 7", 7= = Е 7"+' + 7", (3)

т 2 2

где т — шаг сетки по который может быть как постоянным, так и меняющимся. По построению схема (3) аппроксимирует уравнение (1) с погрешностью О(т2 + к3). Из схемы (3) огрублением оператора Е = Е + (/3Д + ЬЛ получаются схемы первого (Е = Е) и второго (Е = Е + (/3Д) порядков точности относительно к. Именно эти три схемы сравнивались численно между собой на тестах с известными точными решениями.

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

7" = 7 2+7 +( 3Д + ьл] 7".

При нелинейности в правой части выражение !п+1 = f (Un+1) зависит от искомой функции на верхнем слое, поэтому схема (3) требует итераций по нелинейности на каждом шаге. Однако есть возможность избежать этого. С этой целью, следуя [5, 7], запишем схему (3) с удвоенным шагом по причем нелинейную правую часть возьмем со среднего слоя:

и^1 - иn-1 хип+1 + иn-1

сЕ---= аЛ-^--(4)

В схеме (4) нет нелинейности на верхнем слое, поэтому она реализуется на каждом шаге обычной трехточечной прогонкой без итераций. Однако ввиду трехслойности схемы для старта вычислений требуется знание и начального условия, и решения на первом слое. Поэтому для итерационного вычисления решения на первом слое необходимо предварительно использовать какую-нибудь линеаризацию схемы (3), например самую простую:

сЕУ- и« = аЛ- Е!!П±Ш, у« = и», к = 1,2,... (5)

т 2 2

Формально схемы (3)-(5) аппроксимируют и уравнение теплопроводности, и уравнение Шредингера. Отличие, как и в исходном уравнении (1), определяется параметром с, равным вещественной или мнимой единице. Следует заметить также, что в частном случае равномерной сетки

d = к+ — к_ = 0, 5 = к+ ± к_ = 2к, р = к+к- = к2,

и мы приходим в случае уравнения теплопроводности к классической схеме Микеладзе [3] четвертого порядка точности, а в случае уравнения Шредингера с точностью до слагаемых, ответственных за искусственную диссипацию, — к ранее исследованным устойчивым компактным схемам [6, 7] той же точности.

3. Построение адаптивной сетки и организация счета

По функции и на текущем слое эволюционной переменной обычным способом строится возрастающая сеточная функция переменной хг(г = 0, • • • , N)

тс = £ \и — и-_1|~ [

3=1 3

хг

ди

д£

После этого добавлением линейной функции строится обезразмеренная строго возрастающая функция

7 (Л ^ — х° ^ П

Zj = (1 — а)——± а-, а > 0,

хм — х0

где N — число шагов сетки.

Затем на интервале (0,1) изменения функции Z строится равномерная сетка

¿г = (г^ ^ (г = 0, ••• ^)

с тем же числом шагов, и тогда обратная функция х^) определяет однозначно новую сетку хг = х(^7г) с равномерным изменением функции Z и с почти равномерным изменением решения и на шаге. Линейная функция в Z добавляется для того, чтобы в случае

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

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

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

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

Для обоих дифференциальных уравнений сравнительные расчеты проводились на неравномерной адаптивной сетке по трем трехслойным схемам вида (4) точности О(т2 + к9), где д = 1, 2, 3.

Схемы тестировались на точном решении уравнения теплопроводности (1). Решение было выбрано в виде нормированной финитной тепловой волны, сосредоточенной на отрезке |ш| < 1, где

х — ьЬ

ш =-,

ц

V — скорость, а ц — полуширина волны. Форма волны имеет вид

и = (ш2 — 1)4 , |ш| < 1,

а правая часть

ц2 ц

7 = 8 (^ + ь-ш) у/и р = ш2 — 1, м< 1.

ц2 ц

Вне отрезка |ш| < 1 как решение и, так и правая часть / уравнения (1) равны нулю.

Задача решалась в области (0 < х < 1) х (0.4 < Ь < 1.6). Скорость волны V = 0.5, ее полуширина ц = 0.1. За указанный промежуток времени центр волны пробегает от х = 0.2 до х = 0.8. Вес а линейной части управляющей функции Z равен 0.2.

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

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

Для схем точности 0(т2 ± к4) интерполяция решения на новую сетку проводилась по д ±1 узлам. Увеличение порядка интерполяционной формулы на единицу практически не меняло картины. Существенное преимущество компактной схемы третьего порядка во всех расчетах было очевидно. На одинаково детальной сетке (последняя строка табл. 2) результат расчета по компактной схеме приблизительно в 15 раз точнее, чем по обычным схемам. Пятипроцентный порог точности достигается для компактной схемы на сетке 80 х 160, в 8 раз грубее по каждой переменной, чем для других схем.

Схемы тестировались также на точном решении уравнения Шредингера

и = ехр ^ — —^ зееЬ(х)

при нелинейной правой части ! = \и\2 и. В табл. 2 приведены результаты расчета на последовательности сеток в области (0 < £ < 10) х (—20 < х < 20) при а = 0.2. Здесь также левая часть таблицы соответствует одновременной детализации обеих сеток, а результаты, помещенные в правой части, получены при фиксированном мелком шаге по эволюционной переменной.

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

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

Пх П д = 1 д = 2 д = 3 Пх П д = 1 д = 2 д = 3

20 40 0.450 0.740 0.526 20 1280 0.316 0.383 0.277

40 80 0.344 0.321 0.116 40 1280 0.279 0.318 0.123

80 160 0.192 0.196 0.050 80 1280 0.185 0.195 0.051

160 320 0.100 0.101 0.013 160 1280 0.098 0.102 0.013

320 640 0.046 0.045 0.003 320 1280 0.045 0.045 0.003

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

Пх П д = 1 д = 2 д = 3 Пх П д = 1 д = 2 д = 3

20 40 0.193 0.277 0.1701 20 1280 0.199 0.246 0.1807

40 80 0.038 0.028 0.0144 40 1280 0.026 0.016 0.0016

80 160 0.141 0.011 0.0028 80 1280 0.371 0.114 0.0126

160 320 0.022 0.004 0.0005 160 1280 0.090 0.015 0.0031

320 640 0.051 0.002 0.0001 320 1280 0.109 0.003 0.0004

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

nx nt q = 2 q = 4 nx nt q = 2 q = 4

80 160 0.311 0.0195 80 1280 0.304 0.0159

160 320 0.066 0.0017 160 1280 0.064 0.0009

320 640 0.016 0.0003 320 1280 0.016 0.0001

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

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

Благодарности. Работа выполнена при финансовой поддержке РФФИ (грант № 1401-00191).

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

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

Kivshar, Yu.S., Agraval, G.P. Optical Solitons. From Fiber Light Guides to Photon Crystals. М.: Fizmatlit, 2005. 647 p. (in Russ.)

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

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

[4] Самарский А.А. Схемы повышенного порядка точности для многомерного уравнения теплопроводности // Журн. вычисл. математики и мат. физики. 1963. Т. 3, № 5. C. 812-840.

Samarskiy, А.А. Schemes of high-order accuracy for the multi-dimensional heat conduction equation // Zh. Vychisl. Matematiki i Mat. Fiziki. 1963. Vol. 3, No. 5. P. 812-840. (in Russ.)

[5] Shu-Sen Xie, Guang-Xing Li, Sucheol Yi. Compact finite difference schemes with high accuracy for one-dimensional nonlinear Schrodinger equation b,2 // Comput. Methods Appl. Mech. Engrg. 2009. No. 198. P. 1052-1061.

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

[7] Паасонен В.И., Федорук М.П. Компактная безытерационная схема с искусственной диссипацией для нелинейного уравнения Шредингера // Вычисл. технологии. 2012. Т. 17, № 3. С. 83-90.

Paasonen, V.I., Fedoruk, M.P. A compact noniterative scheme with artificial dissipation for nonlinear Schrodinger equation // Comput. Technologies. 2000. Vol. 17, No. 3. P. 83-90. (in Russ.)

[8] Паасонен В.И. Схема третьего порядка аппроксимации на неравномерной сетке для уравнений Навье —Стокса // Вычисл. технологии. 2000. Т. 5, № 5. С. 78-85. Paasonen, V.I. A third-order approximation scheme on a non-unoform grid for Navier — Stokes equations // Comput. Technologies. 2000. Vol. 5, No. 5. P. 78-85. (in Russ.)

[9] Петухов И.В. Преобразование уравнений пространственного пограничного слоя для численного расчета // Уч. зап. ЦАГИ. 1982. Т. 8, № 5. C. 69-78.

Petukhov, I.V. Transformation of the spatial boundary layer for the numerical calculation // Uch. Zap. TsAGI. 1982. Vol. 8, No. 5. P. 69-78. (in Russ.)

[10] Thompson, J.F. Grid generation techniques in computational fluid dynamics // AIAA J. 1984. Vol. 22, No. 11. P. 1505-1523.

[11] Rai, M.M., Anderson, D.A. Application of adaptive grids to fluid-flow problems with asymptotic solutions // AIAA J. 1982. Vol. 20, No. 4. P. 496-502.

[12] Dwyer, H.A. Grid adaptive for problem in fluid dynamics // AIAA J. 1984. Vol. 22, No. 12. P. 1705-1712.

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

Compact third-order accuracy schemes on nonuniform adaptive grids

Paasonen, Viktor I.

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

Compact difference schemes of the third order of accuracy (with the second order of approximation on an evolutionary variable) on non-uniform grids are created for both one-dimensional nonlinear Schrodinger equation and the heat conductivity equation.

© ICT SB RAS, 2015

64

EM. naacoHeH

The purpose of this study aims at the development of such difference method for the considered equations in which some approaches for the increase of the efficiency of calculations are combined. These approaches include a non-uniform grid, a high order of approximation, lack of iterations on nonlinearity, dynamic adaptation of the grid to the solution and rather exact interpolation of the solution on a new grid.

Accurate calculation of solutions with the concentrated waves in case of uniform grids requires quite small step. However, the small step is not required in zones with moderate values of gradient. Therefore, in this research we prefer to use the non-uniform grids. Another useful tool for improvement of the quality of calculations that we apply is the higher order of accuracy. For this purpose we implement a generalization of known compact schemes on a case of non-uniform grids in this research.

The high-accurate technology of calculation is supplemented with the mechanism of adaptation of a grid to the solution. For this purpose, the evolutionary variable is reconstructed on each step, and the solution is interpolated on the new grid with sufficient accuracy. Formally, here we apply the two-layer schemes with double step on three layers. The center layer is used only for the approximation of non-linear part in the right hand side of equations. It allows us to exclude the need for iterations.

Results of comparative calculations for compact schemes and known usual schemes are given in the final part of research. All carried-out calculations have shown the essential advantage of compact schemes in the accuracy and economy of resources.

Keywords: Shrodinger equation, compact difference scheme, scheme on nonuniform grids, high-order accuracy scheme, nonlinear fiber optics.

Acknowledgements. This research was supported by Russian Foundation for Basic Research (grant No. 14-01-00191).

Received 28 November 2014 Received in revised form 13 March 2015