CC BY
127
Вычислительные технологии
Том 14, № 4, 2009
ENO- и WENO-алгоритмы сплайновой интерполяции
В. И. Пинчуков Учреждение Российской академии наук Институт вычислительных технологий СО РАН, Новосибирск, Россия
e-mail: pinchvi@ict.nsc.ru
Предложены ENO- и WENO-модификации нелокального сплайна третьей степени, которые позволяют уменьшить осцилляции сплайна возле разрывов интерполируемой функции или ее производной. Приводятся численные примеры решения тестовых задач.
Ключевые слова: интерполяция, кубический сплайн, ENO- и WENO-алгоритмы, монотонность.
Введение
Недостаток классического нелокального, т. е. определяемого посредством решения сеточной системы уравнений, сплайна третьей степени [1] — в появлении осцилляций в случае разрыва интерполируемой функции или ее производной. Поэтому проводятся поиски сплайнов других типов [2-7]. Так, проблема монотонности решена в локальных сплайнах [2, 3], обладающих свойствами изогеометрии. Эффектом подавления осцилляций за счет механизма натяжения обладают гиперболические сплайны [4, 5]. Применение алгоритмов монотонизации, ранее употреблявшихся при построении численных методов решения гиперболических уравнений первого порядка, в частности, уравнений Эйлера, позволило построить неосциллирующие сплайны (см. публикации [6, 7]). Алгоритмы этого типа развиваются и в настоящей работе.
Методы построения неосциллирующих схем высоких порядков для численного решения гиперболических уравнений (в частности, уравнений газовой динамики) можно разделить на две группы. Методы первой группы (см. обзор [8]) основаны на лимитировании слагаемых, отвечающих за высокий порядок этих методов, которое срабатывает, как правило, возле разрывов. Некоторые из алгоритмов этой группы (например, TVD-схема Хартена, см. [8]) гарантируют получение монотонных численных решений. Проблемой такого подхода является тенденция к снижению порядка аппроксимации и, соответственно, к ухудшению точности вблизи экстремумов решения. Аналогичный недостаток свойствен нелокальным монотонным сплайнам [7]. Методы второй группы (см. обзор [8]) основаны на алгоритмах ENO (Essentially NonOscillating — существенно неосциллирующие), заключающихся в локальном выборе аппроксимационной формулы с наименьшими осцилляционными свойствами из набора формул одного и того же порядка. Теоретически монотонность методов из этой обширной группы в большинстве случаев не доказана и на практике осцилляции возможны, однако, в соответствии с названием, лишь относительно небольшие. Предметом данной работы является развитие
© ИВТ СО РАН, 2009.
подобного подхода в конструировании кубических сплайнов, ранее предложенного в [6] для случая равномерной сетки.
Применение ENO- и WENO- (т. е. Weighted ENO) алгоритмов удобно, если использовать формулировку кубического сплайна s (x), основанную на значениях его первой производной s'(xj) = v, 0 < i < I, в узлах интерполяции 0 < xj < a (напомним, в этих узлах выполнены условия s(xj) = / и, при 0 < i < I, s'(xj+0) = s'(xj—0), s"(xj+0) = s"(xj—0)). Поскольку первая производная сплайна является кусочно-квадратичной непрерывной функцией, она может быть представлена на [xi-i, xj] формулой
s'(x) = Vj—i (1 — £) + Vj£ + c(1 — £)£, £ = (x — xj-i)/hj-i/2, hj-1/2 = xj — xj-i,
где c — произвольная константа, которую определим после интегрирования этой формулы, с учетом условий интерполяции s(xj—i) = /j—i, s(xj) = /j. В результате получаем
s(x) = /j—1(1 — 3£2 + 2£3) + /j(3C2 — 2£3) + [Vj—i (£ — 2£2 + С3) + 3 — £2)]hj—1/2- (1)
Это выражение гарантирует непрерывность первой производной сплайна. Продифференцировав его дважды, подставим полученную формулу в условие непрерывности второй производной s''(xj + 0) = s''(xj — 0). В результате можно получить уравнение
(Vj—i + 2Vj)/hj—1/2 + (2Vj + Vj+i)/hj+i/2 = 3£j+i/2/hj+i/2 + 35j—i/2/hj—1/2, (2)
¿j+1/2 = (/j+1 — /"j) /hj +1/2, ^j —1/2 = (/j — Л— 1)/hj—1/2-
В качестве граничных условий используем равенство нулю второй производной, эквивалентное соотношениям
2vo + Vi = 3ii/2, 2v/ + V/—1 = 1/2. (3)
Итак, получена замкнутая линейно независимая система уравнений с диагональным преобладанием для определения величин Vj, 0 < i < I на неравномерной сетке, реализующая классический сплайн третьего порядка.
Перенос алгоритмов решения гиперболических уравнений, основанных на автоматической замене возле разрывов аппроксимаций высоких порядков монотонными односторонними аппроксимациями первого порядка, позволяет получить уравнение для производных сплайна:
V — ip/hj —1/2 + (3 — pj )zjVj + V+ipj/hj+1/2 =
= 3MM(7^j+i/2Zj, ¿j—1/2/hj—1/2 + ¿j+1 /2/hj+1 /2, Y^j—i/2Zj), Y = л/2, (4)
MM (a, b, c) = max[—d, min(b, d)], d = min(| a |, | c |), zj = 1/h —1/2 + 1/hj+1/2,
pj = min[1, ymin(| ¿+1/2 |, | ¿j—1/2 |)z/(| ¿j+1/2 | /hj+i/2+ | ¿j —1/2 | /hj—1/2)]- (5)
В качестве граничных условий используются соотношения (3). Отметим, что функции min и max осуществляют лимитирование членов, отвечающих за высокий порядок аппроксимации кубического сплайна. В TVD-схемах решения уравнений Эйлера [8] аналогичное лимитрование осуществляется с помощью известной функции minmod.
Можно показать [7], что построенный сплайн является невозрастающим (неубывающим), если этим свойством обладают входные данные. Именно, верна
Теорема. Если значения интерполируемой функции в узлах интерполяции х, 0 < г < I, образуют невозрастающую (неубывающую) последовательность, то кубический сплайн, построенный по формулам (3)-(5), является непрерывной невозрастающей (неубывающей) функцией с непрерывной первой производной.
Для гладких монотонных интерполируемых функций этот сплайн тождествен исходному (2). Однако вблизи экстремума даже для гладкой интерполируемой функции отличие может иметь место, что повлечет снижение точности. Это очевидный недостаток подхода [7].
1. Е^-сплайн
Поэтому попытаемся построить обобщения ЕМО- и ШЕМО-алгоритмов, применяемых для решения уравнений Эйлера, и используем их для интерполяции функций. Ранее кубический ЕМО-сплайн на равномерной сетке предложен в [6]. Здесь этот подход применяется в другой версии (идеи ЕМО допускают простор реализаций, см. обзор [8]), а также строятся ШЕМО-алгоритмы. Рассмотрим уравнение, в котором снова используются коэффициенты р в левой части (см. уравнение (4)), однако правую часть, которую обозначим как Л, будем вычислять на основе идей ЕМО:
+ (3 - + р/Л-г+1/2 = Л. (6)
Определив посредине интервалов в узлах х+1/2 = (х + х+1 )/2 дискретную функцию .¿+1/2 = 6(/+1 — /¿)/(хг+1 — х)2, правую часть уравнения (2) можно записать в виде Л = (.£¿+1/2 + .¿-:1/2)/2. Очевидно, что это выражение может рассматриваться как результат интерполяции функции . с помощью полинома Ньютона первой степени: Л = .¿-1/2 + ^¿(ж — Х-1/2), где
N = (.¿+1/2 — .¿-1/2)/(Ж+1/2 — Ж-1/2) (7)
— разделенная разность первого порядка, посредством подстановки в него аргумента ж = ж* = (ж+1/2 + Ж-1/2)/2 = (ж-1 + 2ж + ж+1)/4.
Как первый шаг необходимо определить выражения, приближающие правую часть уравнения (2) на смещенных шаблонах, т. е. на шаблонах, отличных от шаблона формулы (7): (г — 1/2, г + 1/2), с тем же порядком аппроксимации, что и выражение Л = (.¿+1/2 + .¿-1/2)/2. Таким образом, следует построить линейные интерполяции для шаблонов (г — 3/2, г — 1/2) и (г + 1/2, г + 3/2). Применяя для вычисления функции . в точке ж = ж* полиномы Ньютона, построенные на смещенных шаблонах, можно получить выражения
Л- = .¿-1/2 + (ж* — ж¿—1/2 Л+ = .¿+1/2 + (ж* — ж¿+l/2)N¿+l. (8)
В соответствии с ЕМО-подходом, для аппроксимации правой части уравнения (6) следует выбирать то выражение, которое соответствует минимальной разделенной разности в интерполирующих полиномах Ньютона. Учитывая, что для "хорошей" дискретной функции / (без разрывов и изломов) должен реализоваться стандартный кубический сплайн, т. е. правая часть уравнения (6) должна принимать исходный вид, следует
ввести коэффициент сжатия Ь > 1 (в ЕМО-схемах решения уравнений Эйлера обычно используется Ь =2). Обозначив как Л± то из выражений Л—, Л+, которое соответствует минимальной из разделенных разностей N¿—1, N¿+1,
Л± = Л—, если | N¿—1 |<| N¿+1 Л± = Л+, если | N¿—1 |>| N¿+1 (9)
где Л—, Л+ определены формулами (8), получаем в итоге
Л = Л±, если Ь | N |> ш1п(| N¿—1 | N¿+1 |), (10)
Л = 1/2 + ^+1/0/2, если Ь | N |< ш1п(| N¿+1 | N¿—1 |), (11)
что замыкает процедуру определения ЕМО-сплайна.
2. WENO-сплайн
При решении гиперболических уравнений более популярным, чем ЕМО-подход, со временем стал ШЕМО-подход. Это связано с тем, что результат вычисления по ЕМО-алгоритмам не обладает непрерывностью от фигурирующих в них величин. Этим недостатком страдают также и формулы (9)—(11). В частности, при незначительном изменении значения функции в каком-либо узле возможно переключение на другую формулу и, соответственно, заметное изменение формы сплайна. Поэтому вместо жесткого переключения используем двухэтапный весовой подход. Для вычисления весов введем величины п = | N |, е > 0, — малая константа, предназначенная для предотвращения ситуаций с делением на 0. Как и в ЕМО-подходе, соотношение величин 1, п^, будет определять выбор формулы для правой части уравнения (6). Однако переход от одного варианта к другому должен осуществляться непрерывным образом. Поэтому предварительно введем весовой аналог формул (9):
Л± = (Л—Пг+1 + Л+ пг—1)/(пг+1 + пг—1). (12)
Снова вводим коэффициент сжатия Ь > 1, а также переходную зону с границами (Ь + 1)/2, (3Ь — 1)/2, шириной Ь — 1, с центром в точке Ь. Потребуем, чтобы при ш1п(п^+1, > пД3Ь — 1)/2 правая часть уравнения (6) вычислялась по формуле
Лг = (^—1/2 + ^¿+1/2)/2, при ш1п(пг—1, ) < пДЬ + 1)/2 следует подставить в правую часть уравнения (6) величину Л±, вычисленную по формуле (12), в противном случае необходимо применять весовую формулу, плавно сопрягающую эти два случая. Для численной реализации описываемого подхода определим весовой параметр д:
д = шах{0, ш1п[1, (ш1п(пг_1, п^+^/п — Ь/2 — 1/2)/(Ь — 1)]}. (13)
Итоговая формула ШЕМО-подхода, объединяющая все варианты, выглядит следующим образом:
Л = (£¿—1/2 + ^+1/2)д/2 + Л± (1 — д). (14)
Формулы (12)-(14) осуществляют мягкое переключение вариантов вычисления правой части уравнения (6).
Численный алгоритм построения сплайна замыкается соотношениями на границах. Наряду с условиями (3) рассмотрим граничные процедуры более высоких порядков аппроксимации. Продифференцировав трижды выражение (1), имеем
в''' (ж) = (V + ^¿+1)6/^2+1/2 — (/¿+1 — /¿)12/й?+1/2-
Если построить интерполяционный многочлен Лагранжа по четырем последовательным сеточным узлам ж7-, ж^+1, ж7-+2, ж7-+3 и трижды его продифференцировать, можно получить
¿'з'(ж) = —'6 Л' Д^+1/2(^+3/2 + ^+1/2)(^+3/2 + ^+1/2 + ^'+5/2)] + +6/7+1 /[^+1/2^+3/2(^+3/2 + ^?+5/2)] — 6/7+2 /[(^+3/2 + ^+3/2^+5/2] +
+6/?-+3/[(^^+5/2 + ^7+3/2 + ^+1/2 )( ^+5/2 + ^7+3/2) ^+5/2] •
С использованием этих выражений можно построить граничные условия третьего порядка аппроксимации:
в'''(0) = 0, в'" (а) = 0, (15)
и условия четвертого порядка аппроксимации:
в'''(0) = ¿''(0), в'''(а) = ¿3'' (а). (16)
3. Численные результаты
Для иллюстрации работоспособности построенных сплайнов применим их к задаче интерполяции ступенчатой функции /0(ж), 0 < ж < 1:
/ 0(ж) = 0, если 0 < ж < 0.15,
/ 0(ж) = 1, если 0.15 < ж < 0.45,
/ 0(ж) = 0, если 0.45ж < 0.77,
/ 0(ж) = 1, если 0.77 < ж < 0.83,
/ 0(ж) = 0, если 0.83 < ж < 1.
На рис. 1, а штриховой линией изображается ЕМО-сплайн, полученный интерполяцией приведенной выше функции по ее значениям в узлах равномерной сетки ж¿ = г/1, 1 = 15, отмеченным кружками. В отличие от схем для решения гиперболических уравнений оказалось, что значение параметра сжатия, большее двух, обеспечивает лучшие результаты. Поэтому здесь и в других задачах принято Ь = 6. Недостатком сплайна является провал вблизи левой границы, объясняемый близостью к этой границе разрыва интерполируемой функции (отметим, вблизи границы ЕМО- и ШЕМО-подходы работают на урезанных шаблонах и минимум разделенных разностей отыскивается среди меньшего числа аргументов, причем в данном случае все они являются "разрывными"). Для преодоления этого недостатка в точках, отстоящих от границ на один интервал, т. е. при г =1 и г = 1 — 1, правая часть уравнения (6) вычисляется не по формулам (12), (13), а по монотонным формулам (4), (5). Результат изображен на рис. 1, а сплошной линией. Он представляется приемлемым. На рис. 1, б сплошной линией изображен ШЕМО-сплайн для функции /(ж) = /0(ж) + 0.3ж. Используется неравномерная
сетка, полученная отображением х(п) = п(0.3 + 0.7п) равномерного разбиения интервала 0 < п < 1- Для получения приемлемых вблизи границы результатов здесь также потребовалось использовать формулы (4), (5). Для сравнения на этом рисунке штриховой линией изображен классический сплайн третьей степени, порождающий обычные для таких данных пульсации.
Выше изучалась работоспособность построенных сплайнов на входных данных, описывающих разрывные функции. Применим их также для непрерывных функций с разрывами первой производной. Рассмотрим на интервале 0 < х < 1 функцию
f (х) = V.1875 + у — 4у2/.55 + 0.3х, если у < 0.375;
/(х) = 0.3х, если у > 0.375, у =| х — 0.5 | .
Она имеет изломы в точках х = 0.125, х = 0.5, х = 0.875. Оказалось, что ЕМО-, WENO- и классический сплайн имеют приблизительно одинаковое качество. На рис. 2 приведены графики WENO-сплайна (рис. 2, а) и монотонного сплайна (рис. 2, б), определяемого формулами (3)-(5). Используется равномерная сетка с 15 узлами. При получении результатов, изображенных на рис. 2, а, вблизи границы (при г = 1 и г = I — 1) используются монотонные формулы (4), (5). Монотонный сплайн удовлетворителен везде, кроме одного из имеющихся экстремумов (левого), где он срезает локальный максимум. Возможность этого эффекта отмечалась выше.
Рис. 1. ЕМО-сплайн с граничными условиями (3) (штриховая линия) и условиями (3)-(5) (сплошная линия) (а); ШЕМО- (сплошная линия) и классический (штриховая линия) сплайны (б)
Рис. 2. ШЕМО- (а) и монотонный (б) [7] сплайны
Исследуем точность полученных и исходного сплайнов на гладкой функции:
/(ж) = [ехр(—2ж) — 2ехр(—4ж) + ехр(—6ж)]27/4, 0 < ж < 1.
Она имеет экстремумы на границе интервала при ж = 0 и в точке ж = ж0 , определяемой соотношением ехр(—2ж0) = 1/3, где достигается максимум /(ж0) = 1. Эта функция может служить тестом для оценки влияния различных граничных условий на точность сплайна, так как экстремум на границе интервала может вызвать понижение точности при включении локальной монотонизации (4), (5).
В табл. 1 приведена погрешность 5 = тахх(| / (ж) — в (ж) |) для ЕМО, ШЕМО, классического и монотонного сплайнов на равномерной сетке из 21-го узла.
Можно видеть, что повышение порядка аппроксимации граничных условий приводит, как правило, к повышению точности сплайнов. Следует, однако, заметить, что граничные условия (16) могут вызвать нарушения монотонности в случае входных данных с "разрывами" вблизи границ. Эту немонотонность не удается устранить даже при использовании возле границ формул (4), (5). То есть соотношения (15) представляются предпочтительными по сравнению с соотношением (3) или (16), так как они заметно точнее условий (3) и сравнимы по точности с условиями (16), причем, в отличие от последних, допускают монотонизацию.
Именно с условиями (15), дополненными формулами (4), (5), получаем данные табл. 2.
На рис. 3, а приведены графики ЕМО- (сплошная линия) и классического (штриховая линия) сплайнов, на рис. 3,б — ШЕМО- (сплошная линия) и монотонного ([7], штриховая линия) сплайнов. Наилучшие результаты дает монотонный сплайн. Графики ЕМО- и ШЕМО-сплайнов сходны, они имеют пульсации на серии рядом расположенных изломов данных, когда все фигурирующие в построении этих сплайнов шаблоны являются "негладкими". Тем не менее на изломах, примыкающих к "гладким" участкам, эти сплайны ведут себя лучше, чем классический сплайн (рис. 3, а), который также содержит пульсации на серии изломов.
Резюмируя, можно констатировать, что ЕМО- и ШЕМО-сплайны нигде не уступают классическому. Для входных данных, содержащих серии изломов, эти сплайны хуже монотонного сплайна. Последний, однако, уступает им по точности при интерполяции
Таблица 1.
Граничные условия Сплайны
ЕШ WENO классический монотонный
Формула (3) 5.672 ■ 10-3 5.608 ■ 10-3 4.862 ■ 10-3 8.871 ■ 10-3
Формулы (3) + (4), (5) 9.572 ■ 10-3 9.505 ■ 10-3 8.871 ■ 10-3 8.871 ■ 10-3
Формула (15) 4.808 ■ 10-3 4.468 ■ 10-3 3.979 ■ 10-3 5.871 ■ 10-3
Формулы (15) + (4), (5) 8.752 ■ 10-3 8.435 ■ 10-3 5.871 ■ 10-3 5.871 ■ 10-3
Формула (16) 4.494 ■ 10-3 4.147 ■ 10-3 5.495 ■ 10-3 5.337 ■ 10-3
Формулы (16) + (4), (5) 8.426 ■ 10-3 8.109 ■ 10-3 5.337 ■ 10-3 5.337 ■ 10-3
Таблица 2.
X¿ 0 2 3 5 6 8 9 11 12 14 15
/¿ 10 10 10 10 10 10 10.5 15 56 60 85
Рис. 3. ENO- (сплошная линия), классический (штриховая линия) сплайны (a), WENO-(сплошная линия), монотонный (штриховая линия) [7] сплайны (б)
гладких функций. В целом ситуация сходна с той, которая имеет место при использовании ENO- и WENO-алгоритмов в численных исследованиях течений сжимаемого газа: схемы, построенные на их основе, имеют более высокие аппроксимационные свойства, чем схемы, построенные на лимитировании аппроксимационных поправок, например, TVD-схема Хартена (см. [8]), но уступают им по монотонным качествам. При одинаковом порядке ENO- и WENO-схемы в меньшей степени генерируют осцилляции, чем классические схемы с фиксированным шаблоном (напомним, ENO- и WENO-алгоритмы основаны на локальном выборе оптимального шаблона из набора шаблонов).
Список литературы
[1] Завьялов Ю.С., Квасов Б.И., Мирошниченко В.Л. Методы сплайн-функций. М.: Наука, 1980. 352 с.
[2] ÄKIMA H. A new method of interpolation and smooth curve fitting based on local procedures // J. Assoc. Comput. Machinery. 1970. Vol. 17. P. 589-602.
[3] Квасов Б.И. Методы изогеометрической аппроксимации сплайнами. М.: Физматлит, 2006.
[4] Sapidis N.S., Kaklis P.D. An algorithm for constructing convexity and monotonicity-preserving splines in tensions // Comput. Aided Geometric Design. 1988. Vol. 5. P. 127-137.
[5] Паасонен В.И. Параллельный алгоритм построения гиперболических сплайнов // Вы-числ. технологии. 2006. T. 11, № 6. С. 87-95.
[6] Пинчуков В.И. ENO-модификация нелокального кубического сплайна на равномерной сетке // Вычисл. технологии. 2000. T. 5, № 6. С. 62-69.
[7] Пинчуков Б.И. Монотонный нелокальный кубический сплайн // Журн. вычисл. математики и матем. физики. 2001. T. 41, № 2. C. 200-206.
[8] Пинчуков В.И., Шу Ч.-Ш. Численные методы высоких порядков для задач аэрогидродинамики. Новосибирск: Изд-во СО РАН, 2000. 232 с.
Поступила в редакцию 6 октября 2008 г., в переработанном виде —17 февраля 2009 г.
