CC BY
13
Вычислительные технологии, 2020, том 25, № 2, с. 63-79. © ИВТ СО РАН, 2020 Computational Technologies, 2020, vol. 25, no. 2, pp. 63-79. © ICT SB RAS, 2020
ISSN 1560-7534 eISSN 2313-691X
ВЫЧИСЛИТЕЛЬНЫЕ ТЕХНОЛОГИИ
D01:10.25743/ICT.2020.25.2.006
Дробно-рациональная аппроксимация в начально-краевых задачах с фронтами
С. В. Идимешев
Институт вычислительных технологий СО РАН, Новосибирск, Россия Контактный автор: Идимешев Семен В., e-mail: [email protected]
Поступила 18 .марта 2019 г., доработана 2 .марта 2020 г., принята в печать 16 марта 2020 г.
Предложен спектральный метод на основе дробно-рациональной аппроксимации. На примере решений уравнения Бюргерса с особенностями в виде фронтов показано, что дробно-рациональное приближение решений имеет существенные преимущества перед полиномиальным. Для эффективной реализации дробно-рациональной аппроксимации в работе использована барицентрическая интерполяционная формула Лагранжа, обеспечивающая быстроту вычислений и численную устойчивость. Для адаптации узлов интерполяции использован метод, основанный на аппроксимации положения особенности аналитического продолжения решения в комплексной плоскости. Предложено обобщение метода на случай нескольких особенностей. Описано построение спектрального метода и проведены расчеты на модельных задачах, в т. ч. с двумя фронтами.
Ключевые слова: дробно-рациональная интерполяция, полиномиальная интерполяция, спектральный метод, особенность в комплексной плоскости, барицентрическая интерполяционная формула Лагранжа, уравнение Бюргерса.
Цитирование: Идимешев С.В. Дробно-рациональная аппроксимация в начально-краевых задачах с фронтами. Вычислительные технологии. 2020. 25(2): 63-79.
Введение
Среди множества различных численных методов решения дифференциальных уравнений особое место занимают спектральные методы [1]. Спектральные методы — семейство численных методов, объединенных идеей представления решения одной глобальной функцией, определенной на всей расчетной области. На гладких решениях это позволяет получать экспоненциальный характер уменьшения погрешности О (e-W) при увеличении свободных параметров в представлении решения N. Традиционно в спектральных методах для аппроксимации решения используются алгебраические или тригонометрические полиномы. Но известно, что полиномиальная аппроксимация имеет существенные ограничения для быстроизменяющихся функций. По этой причине полиномиальные спектральные методы обычно не используются для решения уравнений с большими градиентми в виде фронтов, хотя идеологию спектрального подхода можно применить и для этого класса, если использовать другие базисные пространства функций.
Известно, что дробно-рациональная аппроксимация не уступает полиномиальной и принципиально превосходит ее для функций с особенностями. Сложность реализации
дробно-рациональной аппроксимации в спектральном методе заключается в нелинейности относительно неизвестных коэффициентов и неустойчивостей в виде "ложных" полюсов и дуплетов Фруассара [2].
В настоящей работе рассматривается специальный способ записи дробно-рациональной функции в виде барицентрической интерполяционной формулы Лагранжа. Такое представление позволяет построить эффективный метод интерполяции функций с особенностями в виде больших градиентов, например крутых фронтов. На основе этого представления можно построить псевдоспектральный метод — метод коллокаций, в котором для аппроксимации решения используется одна глобальная функция. Особенности в виде фронтов могут возникать при решении уравнения Бюргерса
иг = —иих + еихх, х Е [—1,1], Ь Е [0,Т], (1)
которое будет использовано для верификации метода. Преимущества и недостатки метода будут описаны в этой статье.
1. Барицентрическая интерполяционная формула Лагранжа
Традиционно для построения полиномиальной интерполяции используется интерполяционная формула Лагранжа (ФЛ). Несмотря на свою простоту и конструктивность, ФЛ обладает одним существенным недостатком. Как известно, для вычисления значения интерполянта в точке ФЛ требует 0(М2) операций (М — число узлов интерполяции). Существуют другие более экономичные формулы записи полинома. Например, во многих учебных пособиях альтернативной формулой записи полинома предлагается формула Ньютона. Она требует 0(М) операций, если предварительно рассчитана таблица разделенных разностей. Удобная при работе с небольшими значениями N эта формула проявляет численную неустойчивость с ростом числа узлов интерполяции [3]. Более удачной альтернативой является барицентрическая интерполяционная формула Лагранжа.
Барицентрическая интерполяционная формула Лагранжа (БФЛ) является модификацией классической ФЛ. Для ее построения выпишем интерполяционный полином в форме Лагранжа Р^ (ж)
N
N
Рм(х) = ^2 /(хп) Ьп(х), Ьп(х) = Д
п=0
т=0,т=п
X Хт
где Хг — попарно различные узлы интерполяции, £(х) — интерполируемая функция. Далее введем обозначения
N
^^По^ - хп), к
п=0
1
N
(
т=0,т=п
и перепишем (2) в виде
N
Рм (х) = ш(х)^2 / (хп)~
п=0
\г,
_
Будем называть Ап весами ФЛ. Очевидно, если /(х) дественно равен 1:
А„,
= 1, то сам полином также тож-
1 = ш(х)
N
I ^ гр _ гр
.ь .ь п
п=0
Если левую и правую части уравнения (4) разделить на левую и правую части (5) соответственно, то получим общий вид барицентрической интерполяционной формулы Лагранжа:
N Хп
Е _п 1(хп)
п=0х хп
Рм (х)
м А
Е -А^х-
п=0х хп
(6)
Без потери общности рассмотрим задачу интерполяции на отрезке [-1,1]. Для вычисления по формуле (6) требуется предварительно рассчитать параметры Ап, зависящие от узлов интерполяции (3). Если в качестве узлов интерполяции использовать корни или экстремумы многочлена Чебышёва, то их вычисление можно значительно упростить. Экстремумы многочлена Чебышёва удобны при решении краевых задач, так как среди них содержатся граничные значения расчетной области - 1 и 1 . Для экстремумов многочлена Чебышёва
хг.
/пп\
\ж)
сое — , п = 0,...,ЛТ
веса Ап имеют аналитический вид
1
Ап
N
(-1) $п, =
{Г
п = 0 или п = N, иначе.
Полученную в итоге форму записи интерполяционного многочлена Лагранжа назовем барицентрической интерполяционной формулой Лагранжа
Рм (х) =
N £'
п=0
(-1)г
х хГ,
N £'
п=0
(-1)?
,х 'х г.
где символ "указывает, что первый и последний члены суммы домножаются на 1/2. Барицентрическая интерполяционная формула Лагранжа для вычисления значения в точке требует О^) операций, и она обладает важным свойством прямой устойчивости [3]. Более подробную информацию о сравнении различных полиномиальных интерполяционных формул можно найти в [3].
Ограничения полиномиальной интерполяции. Как известно, интерполяционный процесс в корнях или экстремумах многочлена Чебышёва сходится для всех функций, удовлетворяющих условиям Дини — Липшица [2, 4]. Рассматриваемые далее решения уравнения Бюргерса (1) бесконечное число раз дифференцируемы, а следовательно, удовлетворяют условию Дини — Липшица. Но на практике сходимость может
10°
л
о 10-5
5
а
6
о С
П10
10
-10
e = 0.0001
V e = 0.001
\ £~ = 0.01 \
e = 0.1 \
10
10
n
10
10
Рис. 1. Моделирующая фронт функция F(х, е), где параметр е характеризует "крутизну" фронта (а); погрешность интерполяционного процесса (8) для различных е (б) Fig. 1. Function F(х,е) simulates front, with the parameter e characterizes the thickness of the front (a); the interpolation error (8) for various e (б)
a
быть катастрофически медленной. Покажем это на примере приближения функции, являющейся решением уравнения Бюргерса и моделирующей фронт (рис. 1):
' = ( -Г ) •
Параметр е характеризует "крутизну" фронта — чем меньше его значение, тем более узким является фронт. Также можно сказать, что функция Р(х,е) имеет внутренний погранслой.
Из рис. 1 видно, что точность интерполяции с использованием БФЛ высока лишь при относительно больших е. Как только фронт становится узким, т. е. при уменьшении £, интерполяционный процесс начинает медленно сходиться. Для случая е = 0.0001 даже при аппроксимации полиномом степени 1000 погрешность имеет порядок 0(1).
Рис. 2. Распределение узлов интерполяции для функции F(х,е) при е = 0.01. Экстремумы многочлена Чебышёва сгущаются вблизи концов отрезка и разрежены в центре Fig. 2. The distribution of the interpolation points for the function F(x,e) with e = 0.01. The Chebyshev points are concentrated near the ends of the segment and are sparse in the center
Чем уже фронт, тем меньше точек интерполяции на него попадают (рис. 2), что в данном примере усугубляется особенностью распределения узлов интерполяции. Экстремумы многочлена Чебышёва (как и его корни) сгущаются к краям отрезка и наиболее разрежены в его центре. Очевидно, что такое расположение узлов интерполяции для данных функций неудачное.
К сожалению, произвольное расположение узлов невозможно. Достаточно вспомнить пример Рунге [5], демонстрирующий недостатки равномерного распределения узлов. При полиномиальной интерполяции к расположению узлов предъявляются строгие требования, гарантирующие сходимость интерполяционного процесса. Необходимым условием сходимости является медленный рост константы Лебега [2], обеспечить его намного проще, если перейти от полиномиальной интерполяции к дробно-рациональной.
2. Дробно-рациональная аппроксимация и барицентрическая формула
В барицентрической формуле (6) трудно узнать полином. Функция (6) — это полином только потому, что веса Ап имеют явную зависимость (3) от узлов интерполяции. Но если их рассмотреть как независимые параметры, то получим общий вид дробно-рациональной интерполяции в барицентрической форме. Действительно, простой подстановкой можно показать, что при любых ненулевых значениях Ап формула (6) является дробно-рациональной функцией и удовлетворяет условиям интерполяции. Используя возникшие дополнительные степени свободы, можно построить устойчивый интерполяционный процесс на более широком классе сеток. Но каким образом выбрать веса Ап и узлы интерполяции?
Элегантный и практичный анализ погрешности интерполяции можно провести, если рассмотреть аналитическое продолжение приближаемой функции в комплексную плоскость. Для этого определим эллипс Бернштейна — эллипс в комплексной плоскости с фокусами на вещественной оси в точках ±1. Радиусом эллипса Бернштейна р назовем сумму его главных полуосей (р > 1). Тогда условия сходимости интерполяционного процесса определяет следующая теорема [6]:
Теорема. Пусть И1 и — области в С, включающие вещественные интервалы 3 = [— 1,1] и I соответственно. Пусть д : И1 ^ — конформное отображение, такое что д(3) = I. Если для функции f : ^ С отображение /од : И1 ^ С аналитично внутри эллипса Бернштейна радиуса р, то для дробно-рациональной интерполяции (ДР БФЛ)
Е т—г) /
(х — ¿п) '
Ям (х) = -, & = д(хп),
(-1)" (х - Сп)
п=0
где хп = cos(п^т/N) — экстремумы многочлена Чебышёва, верно
\Ян(х) - /(х)| = О(р).
Из этой оценки следует, что именно величина радиуса эллипса Бернштейна определяет точность интерполяционного процесса. Чем больше радиус, тем выше точность приближения. Из этой теоремы также следует, что при фиксированных весах Хп (7) в качестве узлов интерполяции можно использовать конформно отображенные экстремумы многочлена Чебышёва. При этом сходимость процесса будет определяться радиусом эллипса Бернштейна функции (/од), а не f, как в случае полиномиальной интерполяции (д = 1). Так, удачно выбранным будет отображение д(х), увеличивающее радиус эллипса Бернштейна по сравнению с исходной функцией.
Каким образом можно оценить радиус эллипса Бернштейна для конкретной функции f (х)? Это наибольший эллипс с фокусами в ±1, полностью вписанный в область аналитичности функции f (z), z Е C. Как известно, область аналитичности функции ограничивается положением особенностей в комплексной плоскости. Следовательно, эллипс Бернштейна будет ограничен первой особенностью в комплексной плоскости, которая встретится при его "раздувании", т. е. при увеличении радиуса р.
Способы адаптации сетки. Kosloff и Tal-Ezer [7] предложили д(х), отображающее экстремумы многочлена Чебышёва на почти равномерную сетку:
. . arcsin(ax)
д(х) =-tVV, « Е (0; 1).
arcsin(a)
При а = 0 экстремумы многочлена Чебышёва не перемещаются, а при а ^ 1 стремятся к равномерному распределению на отрезке [— 1,1]. Очевидно, что равномерная сетка более универсальна, чем любая другая. Но в этой работе мы рассматриваем функции с узколокализованными особенностями. Поэтому равномерная сетка, как и любая не адаптирующаяся под конкретную функцию, не позволит эффективно воспользоваться потенциалом ДР БФЛ. Обзор других вариантов отображения д(х) приведен в работе Jafari-Varzaneh и Hosseini [8].
В настоящей работе использовано отображение, предложенное Tee и Trefethen в [9]. Оно замечательно тем, что конструктивно направлено на увеличение радиуса эллипса Бернштейна. Увеличение радиуса происходит за счет переноса особенности, ограничивающей область аналитичности, на периферию отрезка [-1,1]:
д(х) = 8 + С sh
(a* (¥) - (^)) ^ + arsh (^)
:iq)
Отображение построено на предположении, что область аналитичности функции / ограничена в комплескной плоскости двумя лучами (8+8+тог) и (8 — £г; 8-тог). Если точки 8 ± лежат близко к отрезку [—1,1], то радиус эллипса Бернштейна будет близок к 1, т. е. к минимальному значению. Но отображение (10) устроено так, что радиус для функции / о д всегда будет ограничен снизу некоторой величиной, превышающей единицу, что обеспечит повышение точности приближения. Таким образом, насколько близко бы ни была особенность к отрезку [—1,1], отображение (10) отразит ее на достаточное расстояние от отрезка. Параметр 8 можно интерпретировать как положение точки на отрезке [—1,1], в окрестности которой требуется сгущение, а параметр £ — как степень сгущения (чем он меньше, тем сильнее сгущение).
Можно ли для конкретной функции / получить параметры 8 и £? Авторы отображения (10) предложили алгоритм определения этих параметров [9]. Он основан на гипотезе о том, что особенность можно аппроксимировать полюсом второго порядка.
Сравнение положения особенности, радиуса эллипса Бернштейна и погрешности интерполяции для функции F (х, е) при различных е Comparison of the singularity position, Bernstein ellipse radius and interpolation error for the
function F(x, e) for various e
F (x,e) S id P (F) P (F ° 9) Погрешность
полиномиальная дробно-рациональная
F(x, 0.1) 1.1e- 10 6.3e- 1 1.803 2.872 2.5e-15 2.5e-15
F(x, 0.01) 1.1e- 11 6.3e- 2 1.068 1.552 2.2e-3 2.52-14
F(x, 0.001) 1.4e- 13 6.3e- 3 1.0068 1.309 3.2e-1 8.12-11
F(x, 0.0001) 1.1e- 13 6.3e- 4 1.0006 1.213 4.7e-1 1.09-7
F(x, 0.00001) 2.0e- 14 6.3e- -5 1.0001 1.162 4.9e-1 1.75-6
Для этого функция / при помощи аппроксимации Паде локально приближается дробно-рациональной функцией с полиномом второго порядка в знаменателе. Его комплексно-сопряженные корни и определяют полюс. Вещественная часть полюса и абсолютное значение мнимой части есть и соответственно.
На примере функции Р(х, е) в таблице приведены рассчитанные параметры 8, £ и р для различных е. Погрешность интерполяции вычислялась при N = 100. Чем более крутой фронт, тем ближе особенность к вещественной оси (£ ^ 0) и меньше радиус эллипса Бернштейна (р(Р) ^ 1) для исходной функции Р. Но после преобразования (10), которое учитывает положение особенности, радиус р(Р о д) значительно увеличивается, что особенно важно при малых е. Все это отражается на точности интерполяции в соответствии с теоремой. Из таблицы видно, что точность дробно-рациональной интерполяции существенно выше, и даже для очень узких фронтов ( е = 0.00001) наблюдается высокая точность.
а б
100
л
5 10
о а
а
ш
6
С10"10
10
^£ = 0.01 ^-e = 0.001
— £ = 0.0001 "
■ e = 0.1 \ \
ЖШ1Ш1ШИШП08И1П
е = 0.01
100
101
102
103
n
Рис. 3. Погрешность дробно-рационального интерполяционного процесса (9) для функции F(х,е) — (а); распределение адаптированных узлов интерполяции для функции F(х,е) при е = 0.01 — (б)
Fig. 3. The rational interpolation error (9) for the function F(x,e) — (a); distribution of adapted interpolation points for e = 0.01 — (б)
Из рис. 3 видно, что интерполяция с адаптивной сеткой (10) в ДР БФЛ принципиально точнее исходной БФЛ (для сравнения см. рис. 1). Разница особенно чувствительна для узких фронтов. Если положение особенности в комплексной плоскости известно, отображение (10) эффективно сгущает узлы интерполяции в окрестности фронта (рис. 3, б).
Обоснование практичности вышеописанного подхода можно найти в том, что не требуется с высокой точностью определять параметры отображения д(х). Достаточно, чтобы отображение увеличивало радиус эллипса Бернштейна. Но если особенность очень сильная, т.е. ее положение стремится к отрезку [-1,1], то алгоритм становится чувствительным к параметрам 8 и £ .В этом случае особенность нужно приближать с высокой точностью или интерполяционный алгоритм может расходиться. Но практика показала, что описанный выше алгоритм адаптации сетки применим для решения уравнения Бюргерса с е > 0.0001.
Учет нескольких особенностей. На практике могут возникать решения с несколькими фронтами. Рассмотренный выше алгоритм учитывает лишь одну особенность. В этой работе предложен метод построения сетки для учета нескольких особенностей одновременно. Рассмотрим в качестве примера функцию
ч 1 п ( -х - 0.5\ 1 п ( -х + 0.3 \
с двумя фронтами различной "крутизны" £\ и е2 (рис. 4). Требуется разработать алгоритм адаптации сетки, учитывающий положение фронтов и их крутизну. Для определенности положим е± = 0.0001 и е2 = 0.01. Алгоритм построения сетки можно условно разделить на два этапа. На первом этапе локализуются особенности и определяются их параметры. На втором строится д<х), которое отображает экстремумы многочлена Чебышёва на новую сетку с учетом особенностей. В данной работе используется упрощенный подход к локализации особенностей. Отрезок [-1,1] разбивается на заданное число равных подотрезков, и на каждом из них используется алгоритм Tee и Trefethen для определения особенности.
Опишем алгоритм построения отображения д<х) на примере функции F2<x,£i,£2) с двумя фронтами. Если особенности определены, то известны отображения д\<х) и д2<х),
Рис. 4. Функция F2(x, £\,£2) с двумя фронтами / Fig. 4. The function F2(x, £\,£2) with two fronts
построенные на основе (10). Каждое из них учитывает экстремумы многочлена Чебы-шёва с учетом соответствующей особенности (см. верхние рис. 5). Нам необходимо построить отображение, которое переносит экстремумы многочлена Чебышёва на сетку с учетом двух особенностей. Идея построения функции д(х) заключается в суммировании обратных отображений д-1 и д-1 (см. нижний рис. 5). Тогда прямое отображение
есть
Обратные отображения д- (у) имеют явный вид
a+ — а +2 arsh
9i (У)
)
а+ + а
arsh
)
1, 2.
Следовательно, и среднее арифметическое обратных функций
9
- ^ 9\1 + 921 2
можно выписать в явном виде. К сожалению, обратную функцию к д-1 нельзя выписать в элементарных функциях. Но для практических целей не требуется явный вид функции д(х). Важно лишь, как отображение д(х) действует на экстремумы многочлена Чебышёва. Для этого нужно решить N +1 независимых уравнений вида
9~1(Уг)
Xi
0,...,N,
:i2)
где Хг — экстремумы многочлена Чебышёва, у^ — адаптированные узлы интерполяции. Так как функция д-1 строго монотонно возрастающая, численное решение уравнений (12) не представляет трудности. При реализации алгоритма автором использовался векторизованный код для метода бисекций, который позволяет решать все уравнения одновременно.
На рис. 6 представлена погрешность интерполяционного процесса на основе полиномиальной БФЛ и предложенного алгоритма сгущения сетки. Видно, что в первом
92
k^ + gi1)
Рис. 5. Действие отображений gi(x) на экстремумы многочлена Чебышёва Fig. 5. The mapping of gi(x) on the extremes of the Chebyshev polynomial
10L
л
н
g 10" а
а
10"
10
100
101
БФЛ
ДРБФЛ
102
103
N
Рис. 6. Погрешность полиномиального (8) и дробно-рационального (9) интерполяционных процессов для функции F2(x,e1,е2) — (а); распределение адаптированных узлов интерполяции — (б)
Fig. 6. The polynomial (8) and rational (9) interpolation error for the function F2(x,e1,e2) — (а); distribution of adapted interpolation points — (б)
а
случае сходимости нет, тогда как адаптация позволяет с хорошей точностью приблизить функцию, хотя из рисунка видно, что предложенный алгоритм численно менее устойчив, так как увеличением числа узлов нельзя добиться точности, близкой к 10-15.
Описанный выше алгоритм можно применить и для большего числа фронтов. В итоге отображение, учитывающее М особенностей, есть
9= ( + ^ + - + 9]-1 )-1, (13)
где g- 1 определены в (11)
3. Решение начально-краевых задач
Разработанный аппарат применим для решения начально-краевых задач. Опишем псевдоспектральный метод (метод коллокаций) с адаптивными узлами. Для реализации метода необходимо дифференцировать дробно-рациональную функцию (9). В полиномиальных спектральных методах с этой целью применяется эффективная техника, основанная на матрицах дифференцирования [11]. Она позволяет свести задачу дифференцирования к задаче линейной алгебры, а именно к умножению матрицы определенного вида на столбец значений функции в узлах коллокации. Подобная техника может быть использована и для дробно-рационального приближения на основе (9).
Пусть х — вектор-столбец длины N + 1, состоящий из узлов интерполяции, /(х) — вектор-столбец данных в узлах интерполяции, а (х) — вектор значений производных интерполирующей функции К^ (х) (9) в узлах х. В работе [10] для общего случая дробно-рациональной БФЛ показано, что
Л,
(х) = И /(х), И*
^ г =
Л1(Х% х^) N
- Е Щк, ^ ¿ = з, к=0, к=]
аналогично и для второй производном
(х) = И2 /(X), И?
2 п*(п* - ^ )
\ ^г ^J /
N
- Е Щ.
к
^ г = ^ ^ г = 3.
к=0, к=
Быстрые алгоритмы для построения матриц И1 и И2 для (9) можно найти в работе [9].
Продемонстрируем построение псевдоспектального метода на примере простой краевой задачи:
%х
= к(х), х е [-1,1],
и(-1) = Ь_, и(1)
Ь+.
'14)
:15)
Найдем приближенное решение задачи (14), (15) в узлах и = и(х^). Представим решение задачи в виде дробно-рациональной функции (9). Тогда систему линейных алгебраических уравнений (СЛАУ) для уравнения (14) можно записать в виде
И2 и = ВД.
:1б)
Поскольку полученная СЛАУ имеет неполный ранг, ее нужно дополнить краевыми условиями. Для этого удалим условия коллокаций для граничных точек и''(±1) = ^(±1) и добавим условия и(±1) = Ь±. Технически это значит удалить первое и последнее уравнение в (16) и заменить их на 1_ и и 1+ и. Здесь 1± — вектор-строка длины N +1 из нулей и с единицей в первом компоненте для 1_ и последнем для 1+. В итоге получим систему
I _\ / ъ_
И1и1 I и = I НХсгЛ)
1+/ V ь+
где — матрица И2 без первой и последней строки, хсу^ — вектор х без первого и последнего компонента. Для вычисления и необходимо решить полученную СЛАУ. Более подробно техника матриц дифференцирования описана в [11].
Реализация метода для уравнения Бюргерса. Применим этот подход для численного решения уравнения Бюргерса:
щ = -иих + еихх, х е [-1,1], ¿е [0,Т], и(х, 0) = и0(х), и(— 1, ¿) = и_(£), и(1, ¿) = и+ (¿).
:17)
Решению этого уравнения посвященно множество публикаций [16]. Линеаризуем уравнение (3) методом Ньютона:
щ = Ь(и), Ь(и) = -иих - ихи + £ихх + иих,
где и — значения на предыдущей итерации по нелинейности. С использованием техники матриц дифференцирования Ь(и) можно переписать в виде
£(и) = —diag (иж) и — diag (и) Б1 и + еВ2 и + и * . иж = = (—diag (иж) — diag (и) В1 + е В2) и + и * . иж,
где diag (и) — диагональная матрица с и на диагонали, *. — поэлементное умножение векторов.
Для дискретизации по времени использовались две схемы. Явный метод Рунге — Кутты 4-го порядка с переменным шагом по времени в реализации Дормана — Прин-са [13]. Подробное описание метода можно найти в [14]. Если для устойчивости требовался очень мелкий шаг по времени, применялась неявная четырехслойная формула дифференцирования назад 3-го порядка [12]
ип+3 + ——ип+2 + -ип+1 — 2ип = -тЬ(ип+3), 11 11 11 11 v
где ип = и(х, пт) — решение на п-м слое по времени, т — шаг сетки по времени. Во всех расчетах величина подбиралась таким образом, что ее уменьшение не влияло на погрешность численного решения.
Для адаптации сетки на каждом шаге по времени к численному решению применялся описанный выше алгоритм сгущения сетки.
Задача с известным решением. В качестве тестового примера рассмотрим задачу с известным решением
/ \ 1 1 , (—х +1/2 — 0.5\ г .
и(х, ¿) = - + - Ш --^ , ге [0,1], (18)
моделирующим фронт, движущийся вправо и сохраняющий постоянную форму (рис. 7). Здесь для начального условия и(х, 0) применялся алгоритм поиска особенности и адаптации сетки, описанный выше. На следующем шаге по времени строилась новая сетка.
а б
Рис. 7. Решение уравнения Бюргерса с начальными данными (18) — (а) (движущийся фронт сохраняет форму); погрешность численного решения для различных £ — (б) Fig. 7. Solution of the Burgers equation with initial data (18) — (a) (the moving front keeps its shape). The numerical solution error for various £ — (б)
При этом найденная ранее особенность использовалась в качестве начального приближения для локализации новой особенности. Таким образом, разработанный метод относится к численным методам с подвижной адаптивной сеткой.
Погрешность численного решения представлена на рис. 7. Для метода на основе полиномиальной БФЛ (8) сходимость решения возможна только при больших N, так как начальное условие требуется приближать с хорошей точностью. Напротив, адаптивная сетка ДР БФЛ позволяет уже при относительно небольшом числе узлов наблюдать сходимость. К сожалению, здесь не удается получать погрешности, близкие к машинной точности. Существует несколько источников накопления погрешности. Во-первых, это погрешность решения СЛАУ. Второй источник — погрешность метода дискретизации по времени. Она имеет место, несмотря на то, что мы использовали методы высокого порядка. Третья и самая важная причина — это погрешность при расчете особенности. Если при определении особенности начального приближения и(х, 0) алгоритм работает с точной функцией, то уже на следующем шаге по времени он имеет дело с приближением, т. е. дробно-рациональной функцией. Поэтому очень важно, чтобы решение на каждом шаге по времени достаточно хорошо приближало решение задачи.
Гауссова функция. Рассмотрим другой показательный пример решения уравнения Бюргерса (рис. 8):
и(х, 0) = е-8, ж(-М)= ж(М) = 0, £ е [0, 0.6]. (19)
В отличие от предыдущего примера, фронт при движении "набирает крутизну". Траектории узлов подвижной сетки представлены на рис. 8. В начальный момент времени сгущения не требуется, так как интерполяции в экстремумах многочлена Чебышёва достаточно для начального приближения (при £ = 0). Со временем сетка начинает сгущаться в окрестности образующегося фронта. На рис. 9, б представлены траектории движения особенностей в координатах (6,£). В начальный момент времени особенность несущественна, т. е. £ достаточно велика. Затем особенность начинает сдвигаться
а б
Адаптирующаяся сетка
Рис. 8. Решение уравнения Бюргерса с начальными данными (19) — а; траектории узлов подвижной сетки — б
Fig. 8. Solution of the Burgers equation with initial data (19) — a; trajectories of the moving mesh nodes — б
а б
Перемещение особенности во времени
Рис. 9. Картина перемещения особенности во времени при различных £ — а; сходимость численного решения при различных значениях £ (за эталонное принято решение с максимальным N) — б
Fig. 9. The trajectories of moving singularities in time for various £ — a; the convergence of the numerical solution for various e (the solution with the maximum N was taken as reference) — б
вправо (т. е. увеличивается 5) и усиливаться (т. е. £ уменьшается). Чем меньше значение параметра г, тем более сильна особенность. Важно отметить, что особенность в определенный момент выходит на асимптоту и больше не усиливается (т. е. £ ограничена снизу). Отсюда следует, что в этой задаче сетка имеет предельную степень сгущения.
Объединение фронтов. В качестве последнего примера рассмотрим задачу моделирования взаимодействия двух фронтов. В работе Huang и Russell [15] приведено решение уравнения Бюргерса, описывающее движение двух фронтов. Более быстрый фронт догоняет медленный, и образуется единый более крутой фронт (рис. 10). Явный вид решения:
а б
Адаптирующаяся сетка
Рис. 10. Решение (а) и траектории (б) узлов подвижной сетки для задачи (20) Fig. 10. Solution (а) and trajectories (б) of moving mesh nodes for the problem (20)
/-ж+0.Б-4.9Б^ {-х+0.5-0.75Ь\ /-Ж+0.37Б ^
( ) = 0.1е( 20, )+0.5е( 4, )+ е( 2, )
ЩХ,1) (-ж+0.Б-4.9Бг ) (-ж+0.Б-0.7Бг ) (-ж+0.375) . (20)
е( 20г ) + е( 4г ) + е( 2г )
Для решения задачи использовано предложенное в этой работе отображение для сгущения сетки с учетом нескольких особенностей (13). Из рис. 10 видно, что траектории узлов подвижной сетки адаптируются к положению и крутизне фронтов. Интересно, что предложенный алгоритм также позволил учесть объединение фронтов. При е = 0.0005 и N = 150 погрешность расчета имеет порядок 10-5. Более подробный анализ погрешности решения требует отдельной работы с функцией (20), так как она при малых е становится численно неустойчивой.
Заключение
Дробно-рациональная аппроксимация позволяет сделать существенный шаг в развитии спектральных методов. При этом важной задачей по-прежнему остается эффективная реализация дробно-рациональной аппроксимации. Из представленных расчетов видно, что на рассмотренном классе решений предложенный метод сохраняет отличительное свойство спектральных методов, а именно высокую точность. Для функций с особенностями сходимость может наблюдаться уже при относительно малом числе узлов сетки. Также возможность учесть несколько особенностей в решении расширяет область применимости метода.
Однако дробно-рациональная аппроксимация не избавляет нас от сложности самой дифференциальной задачи. В данной работе это перемещающаяся во времени особенность — крутой фронт. Для получения решений высокой точности, а иногда только для обеспечения сходимости требуется использовать малый шаг по времени даже для неявных схем. Такое ограничение возникает из-за необходимости адаптировать сетку на каждом шаге по времени. Чем меньше шаг, тем точнее численное решение передаст информацию с предыдущего шага. Очевидно, что существуют способы оптимизировать процедуру расчета особенности, например анализируя дифференциальное уравнение и его влияние на особенность. Эта работа будет продолжена в следующих исследованиях.
Благодарности. Работа выполнена при финансовой поддержке РФФИ (грант мол_а № 18-31-00202\18).
Список литературы
[1] Boyd J.P. Chebyshev and Fourier spectral methods: second revised edition. New York: Dover Publications; 2000: 611.
[2] Trefethen L.N. Approximation theory and approximation practice. Society for Industrial and Applied Mathematics; 2013: 295.
[3] Higham N.J. The numerical stability of barycentric Lagrange interpolation. IMA Journal of Numerical Analysis. 2004; 24(4):547-556.
[4] Мысовских И.П. Лекции по методам вычислений. СПб.: Изд-во Санкт-Петербургского ун-та; 1998: 784.
[5] Cheney W., Light W. Course in approximation theory. W. Light Brooks, Cole; 2000: 360.
[6] Baltensperger R., Berrut J.-P., Noel B. Exponential convergence of a linear rational interpolant between transformed Chebyshev points. Mathematics of Computation. 1999; 68(227):1109-1121.
[7] Kosloff D., Tal-Ezer H. A modified Chebyshev pseudospectral method with an O(N-1) time step restriction. Journal of Computational Physics. 1993; 104(2):457-469.
[8] Jafari-Varzaneh H.A., Hosseini S.M. A new map for the Chebyshev pseudospectral solution of differential equations with large gradients. Numerical Algorithms. 2015; 69(1):95-108.
[9] Tee T.W., Trefethen L.N. A rational spectral collocation method with adaptively transformed chebyshev grid point. SIAM Journal on Scientific Computing. 2006; 28(5):1798-1811.
[10] Schneider C., Werner W. Some new aspects of rational interpolation. Mathematics of Computation. 1986; 47(175):285-299.
[11] Trefethen L.N. Spectral methods in Matlab. Philadelphia: SIAM University City Center; 2000: 165.
[12] Холл Д., Уатт Д. Современные численные методы решения обыкновенных дифференциальных уравнений. М.: Мир; 1979: 312.
[13] Dormand J.R., Prince P.J. A family of embedded Runge —Kutta formulae. Journal of Computational and Applied Mathematics. 1980; 6(1):19-26.
[14] Shampine L.F., Reichelt M.W. The MATLAB ODE suite. SIAM Journal on Scientific Computing. 1997; 18(1):1-22.
[15] Huang W., Russell R.D. Adaptive moving mesh methods. AMS Book Sciences, vol. 174. New York: Springer; 2011: 432. DOI:10.1007/978-1-4419-7916-2
[16] Shapeev V.P., Vorozhtsov E.V. CAS application to the construction of the collocations and least residuals method for the solution of 3D Navier —Stokes equations. Lecture Notes in Computer Science. Computer Algebra in Scientific Computing. 2013; (8136):381-392.
Вычислительные технологии, 2020, том 25, № 2, с. 63-79. © ИВТ СО РАН, 2020 ISSN 1560-7534
Computational Technologies, 2020, vol. 25, no. 2, pp. 63-79. © ICT SB RAS, 2020 elSSN 2313-691X
COMPUTATIONAL TECHNOLOGIES
DOI:10.25743/ICT.2020.25.2.006 Rational approximation for initial boundary problems with the fronts
Idimeshev Semen V.
Institute of Computational Technologies SB RAS, 630090, Novosibirsk, Russia Corresponding author: Idimeshev Semen V., e-mail: [email protected]
Received March 18, 2019, revised March 2, 2020, accepted March 16, 2020
Abstract
A spectral method with adaptive rational approximation is proposed. In traditional spectral polynomial interpolation, the interpolation points are fixed, usually at the roots or extrema of orthogonal polynomials. Free selection of interpolation points is impossible due to the effect described in the Runge example. The key feature of rational interpolation is the free distribution of interpolation nodes without the occurrence of the Runge phenomenon. Nevertheless, in practice it is very important to implement rational approximation effectively. Here rational approximation is implemented using the barycentric Lagrange form. This leads to fast computations and numerical stability comparable with the polynomial interpolation. It is shown that rational interpolation has
significant advantages over polynomial on functions that have singularities in the form of fronts. The key idea is that rational interpolation allows adapting interpolation points according to function singularities. An effective method of grid adaptation that accounts for singularity location was used. Method was generalized to the case of several singularities, for example, for solutions with several fronts. For the solutions of the Burgers equation with singularities in the form of fronts, it is shown that rational interpolation has significant advantages over polynomial. The implementation of spectral method is described, and calculations results on model problems, including problems with two fronts, are presented.
Keywords: rational interpolation, polynomial interpolation, spectral method, singularity in the complex plane, barycentric Lagrange interpolation form, Burgers equation.
Citation: Idimeshev S.V. Rational approximation for initial boundary problems with the fronts. Computational Technologies. 2020; 25(2): 63-79. (In Russ.)
Acknowledgements. This research was supported by RFBR (grant mol_a No. 18-31-00202\18).
References
1. Boyd J.P. Chebyshev and Fourier spectral methods: second revised edition. New York: Dover Publications; 2000: 611.
2. Trefethen L.N. Approximation theory and approximation practice. Society for Industrial and Applied Mathematics; 2013: 295.
3. Higham N.J. The numerical stability of barycentric Lagrange interpolation. IMA Journal of Numerical Analysis. 2004; 24(4):547-556.
4. Mysovskih I.P. Lectures on numerical methods. Springer Netherlands; 1969: 344.
5. Cheney W., Light W. Course in approximation theory. W. Light Brooks, Cole; 2000: 360.
6. Baltensperger R., Berrut J.-P., Noel B. Exponential convergence of a linear rational interpolant between transformed Chebyshev points. Mathematics of Computation. 1999; 68(227):1109-1121.
7. Kosloff D., Tal-Ezer H. A modified Chebyshev pseudospectral method with an O(N-1) time step restriction. Journal of Computational Physics. 1993; 104(2):457-469.
8. Jafari-Varzaneh H.A., Hosseini S.M. A new map for the Chebyshev pseudospectral solution of differential equations with large gradients. Numerical Algorithms. 2015; 69(1):95-108.
9. Tee T.W., Trefethen L.N. A rational spectral collocation method with adaptively transformed chebyshev grid point. SIAM Journal on Scientific Computing. 2006; 28(5):1798-1811.
10. Schneider C., Werner W. Some new aspects of rational interpolation. Mathematics of Computation. 1986; 47(175):285-299.
11. Trefethen L.N. Spectral methods in Matlab. Philadelphia: SIAM University City Center; 2000: 165.
12. Hall G., Watt J. Modern numerical methods for ordinary differential equations. Oxford University Press; 1976: 350.
13. Dormand J.R., Prince P.J. A family of embedded Runge — Kutta formulae. Journal of Computational and Applied Mathematics. 1980; 6(1):19-26.
14. Shampine L.F., Reichelt M.W. The MATLAB ODE suite. SIAM Journal on Scientific Computing. 1997; 18(1):1-22.
15. Huang W., Russell R.D. Adaptive mesh movement in 1D. Adaptive Moving Mesh Methods. Part of the Applied Mathematical Sciences. 2010; (174):27-135.
DOI:https://doi.org/10.1007/978-1-4419-7916-2_2
16. Shapeev V.P., Vorozhtsov E.V. CAS application to the construction of the collocations and least residuals method for the solution of 3D Navier — Stokes equations. Lecture Notes in Computer Science. Computer Algebra in Scientific Computing. 2013; (8136):381-392.
