CC BY
66
Вычислительные технологии
Том 5, № 6, 2000
eno-модификация нелокального кубического сплайна на равномерной сетке*
В. И. Пинчуков Институт вычислительных технологий СО РАН Новосибирск, Россия e-mail: pinchvi@net.ict.nsc.ru
The modification of unlocal cubic spline is supposed. Monotone increasing or decreasing of spline is garanteed, if values of initial function in the points of interpolation are sutisfied to this condition. Modified spline is identical to classical cubic spline for smooth functions without extremas. Numerical examples of test problems solution are presented.
Как известно, классический нелокальный сплайн третьего порядка [1] позволяет получать очень хорошие результаты для гладких функций. Если же входные данные соответствуют разрывам функции или ее производной, то качество сплайна значительно ухудшается. Поэтому при построении сплайна налагаются различные требования, дополняющие аппрок-симационные [1-5]. Проблема монотонности эффективно решена в локальных сплайнах [3, 4], обладающих свойствами изогеометрии. Предложены также модификации нелокальных сплайнов с сохранением монотонности интерполируемой функции, построение которых основано на двухэтапной процедуре [5]. На первом этапе строится классический нелокальный кубический сплайн и вычисляются его производные. Альтернативный подход включает определение приближенных значений производных функции на основе явных интерполяционных формул Лагранжа без расчета сплайна. На втором этапе численные значения его производных корректируются для обеспечения монотонности и других нужных свойств. Далее на основе известных значений функции и ее производных строится эрмитов сплайн с натяжением, обладающий свойствами изогеометрии.
В настоящей работе предлагается другая модификация нелокального кубического сплайна, в котором обобщено уравнение для расчета его коэффициентов. Надо отметить, что это уравнение видоизменяется локально — в окрестности разрывов интерполируемой функции или ее производных. В остальных областях он представляет собой традиционный нелокальный кубический сплайн. Таким образом, предлагаемая процедура расчета сплайна эквивалентна автоматической сборке кусков классического сплайна третьей степени и видоизмененного возле разрывов интерполируемой функции сплайна. Сборка производится при условии непрерывности первой производной, причем полученный сплайн гарантирует сохранение монотонности входных данных. Такой сплайн обеспечивает интерполяцию
* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, грант №99-01-00560.
© В. И. Пинчуков, 2000.
а б
О 0.5 1.0 0 0.5 1.0
Рис. 1. Классический нелокальный сплайн (а) и его БМС-модификация (б).
достаточно высокого качества в результате простой и удобной численной процедуры, он также может использоваться как база для построения эрмитова сплайна с натяжением.
Рассмотрим стандартную процедуру построения кубического сплайна для однородного распределения узлов интерполяции хг = гН. Пусть /г — значения функции в узлах интерполяции, иг — значения второй производной сплайна в этих узлах. Поскольку его вторая производная должна быть кусочно-линейной функцией, используем формулу
в(2) = иг_1(1 - £) + 0 < £ = (х - гН + Н)/Н < 1.
Проинтегрировав дважды, с учетом условий интерполяции в \х=гн_н= /г_1, в \х=гн= /г получаем
в = /г—1 (1 - £) + /г£ + [иг_1(3£2 - £3 - 2£) + иг(£3 - £)]Н2/6. (1)
Формула (1) гарантирует непрерывность второй производной сплайна, поскольку выведена из непрерывного кусочно-линейного распределения. Продифференцировав (1) и подставив полученную формулу в условие непрерывности первой производной в(1) \х=гн+о= в(1) \ х=гъ,_о, получим
иг_1 + 4щ + иг+1 = 6(/г+1 - 2/г + /_)/Н. (2)
Совокупность этих соотношений и граничных условий, в качестве которых используем равенство нулю второй производной (и0 = и/ = 0), образует замкнутую линейно-независимую систему уравнений для определения величин иг, 0 < г < I. Эта система обладает диагональным преобладанием и может эффективно решаться с помощью стандартной трехточечной прогонки.
Результаты использования описанного сплайна для аппроксимации функции в виде ступеньки /г = 0, если 0 < г < г0, /г = 1, если г0 < г < I, приведены на рис. 1, а (здесь и далее кривая — значения сплайна, значки — узлы интерполяции).
Видно, что качество аппроксимации невысоко из-за наличия осцилляционых структур, типичных для "разрывных" данных. Отметим, что порождение паразитных осцилляций
имеет место также при использовании схем высоких порядков для численного решения гиперболических уравнений, в частности уравнений газовой динамики. В последние годы найдены эффективные методы борьбы с этим негативным эффектом [6-9]. Методы, описанные в работах [6, 7], основаны на переходе от аппроксимаций высоких порядков на гладких участках решения к монотонным аппроксимациям первого порядка возле разрывов. Некоторые из алгоритмов этой группы [7] гарантируют получение монотонных численных решений. Методы другой группы исследований [8, 9] базируются на алгоритмах ENO-аппроксимации, заключающихся в локальном выборе аппроксимационной формулы с наименьшими осцилляционными свойствами из набора формул одного и того же порядка. Теоретически монотонность их не доказана, однако на практике обычно удается получить результаты высокого качества.
В настоящей работе демонстрируется эффективность этих алгоритмов и при использовании их для модификации классического кубического сплайна. Однако для этой цели рассмотренная выше процедура построения сплайна непригодна, так как модификация уравнений (2) приведет к нарушению непрерывности первой производной сплайна и оставит непрерывной вторую, в то время как естественно допустить разрывность прежде всего старшей (второй) производной.
Поэтому используем другую формулировку в узлах интерполяции, которая порождает эквивалентный сплайн, основываясь на значениях его первой производной vi, 0 < i < I. Поскольку первая производная сплайна — кусочно-квадратичная непрерывная функция, рассмотрим формулу
s(i) = Vi-i(1 - О + + c(1 - £)£, i = (x - ih + h)/h,
где c — произвольная константа, которую определим после интегрирования этой формулы с учетом условий интерполяции s \x=íh-h= fi—i, s \x=íh= fi- В результате запишем
s = fi-i(1 - 3i2 + 2i3) + fi(3i2 - 2i3) + [Vi-i(i - 2i2 + i3) + Vi(i3 - i2)]h. (3)
Выражение (3) гарантирует непрерывность первой производной сплайна. Продифференцировав его дважды, подставим полученную формулу в условие непрерывности второй производной s(2) \x=ih+o= s(2) \x=ih—о и получим соотношение
Vi-i + 4vi + vi+i = 3(fi+i - fi—i)/h. (4)
В качестве граничных условий используем равенство нулю второй производной, эквивалентное соотношениям
2vo + vi = 3(fi - fo)/h, 2vj + vi—i = 3(f - fj—i)/h. (5)
Итак, выведена замкнутая линейно-независимая система уравнений с диагональным преобладанием для определения величин vi, 0 < i < I. Ввиду единственности решения задачи о построении сплайна в результате вычислений получаем тот же сплайн, что порождается системой (1), (2).
В соответствии с идеями ENO [8, 9] центрально-разностную аппроксимацию правой части заменяем на минимальное по модулю значение различных аппроксимаций — центральной, левой и правой:
vi—i + 4vi + vi+i = 3mx(a, b, c), mx(a, b, c) = max[-d, min(b, d)],
(6)
d = шт( \ а \, \ с \ ),
а = (3/г - 4/г—1 + /г_2)/Н, Ь = (/г+1 - /г_1 )/Н, С = -(3/ - 4/г+1 + /г+2)/Н.
Таким образом в узлах интерполяции введены разрывы второй производной, причем они невелики в областях гладкого решения и возрастают на скачках аппроксимируемой функции. На рис. 1, б представлены результаты решения задачи об интерполяции ступенчатой функции. Качество полученного сплайна вполне удовлетворительное, однако для более сложных случаев оно оставляет желать лучшего. Например, на рис. 2, б приведены результаты интерполяции функции, описывающей последовательно треугольный (у = шт[(х - х0 + а)/а, (х0 + а - х)/а]), прямоугольный (у = 1, х1 - а < х < х1 + а) и параболический (у = - (х - х2)2/а2) выступы с единичной высотой и одинаковой шириной основания. Как видно, сплайн имеет осцилляционные структуры возле точек разрывов функции или ее производной, однако амплитуда осцилляции у них значительно меньше, чем у исходного сплайна. Для сравнения на рис. 2, а приведены результаты интерполяции этой же функции с помощью сплайна (3)-(5).
Итак, сплайн (5) - (6) для улучшения его монотонных свойств нуждается в дальнейшем развитии. Рассмотрим модификацию исходного сплайна, которая позволяет гарантировать монотонное поведение сплайна при монотонных значениях функции в узлах интерполяции. На наш взгляд, имеется определенная аналогия между этой модификацией и алгоритмами построения конечно-разностных схем [6, 7]. Запишем уравнение для значений производной сплайна в узлах интерполяции:
Рг^г_1 + (6 - 2рг+ ^+1 = Зшх(27^ , 6_ + 27^+), (7)
рг = Ш1П[1, 27 Ш1П( \ \ , \ ¿г_ \ )/( \ \ + \ ¿г_ \ )], (8)
= (/г+1 - /г)/Н, ¿г_ = (/г - /г_х)/Н
(7 — параметр, выбор которого обсуждается далее, здесь же потребуем, чтобы он был больше единицы). Можно показать, что для монотонных входных данных вновь построенный сплайн совпадает с исходным (3)-(5), если в (7) во всех узлах минимальным по
модулю из аргументов функции тх является второй аргумент (в этом случае тх принимает его значение) и = 1. Эти условия выполняются, если во всех узлах
27ш1п(\ \ \ ) > ( \ 5+\ + \ \ ). (9)
Поскольку \ а \ + \ Ь \ = шах( \ а \ , \ Ь \ ) + ш1п( \ а \ , \ Ь \ ), то последнее условие (9) можно преобразовать к виду
(27 - 2) Ш1П(\ ¿г+ \ , \ ¿г_ \ ) > шах(\ \ , \ \ ) - Ш1П(\ ¿г+ \ , \ ¿г_ \ ), (10)
где правая часть неравенства эквивалентна разности модулей величин и . Заменяя ее на модуль разности этих величин, увеличиваем правую часть. Заменив конечные разности производными в некоторых промежуточных точках и приняв для производных в левой части минимальное, а в правой — максимальное по интервалу значения, получим более жесткое условие совпадения модифицированного и исходного сплайнов:
(27 - 2)ш1п( \ /(1) \ ) > Ншах(\ /(2) \ ). (11)
х х
В случае гладких монотонных функций без экстремумов эта оценка выполняется для достаточно малого шага сетки Н.
Докажем, что построенный сплайн является невозрастающим (неубывающим), если этим свойством обладают входные данные. Общая схема рассуждений такова. Установим, что для монотонных входных данных значения производной сплайна в узлах интерполяции, полученные в результате решения системы (7), (8), неотрицательны (неположительны) , если значения функции в узлах не убывают (не возрастают). Для этого построим сходящийся итерационный процесс ее решения и путем индукции докажем нужное свойство. Затем покажем, что нужный знак производной сплайна в узлах интерполяции гарантирует такой же знак производной во всех точках интервала интерполяции.
Пусть, например, /г_ 1 < /г, 1 < г < I. Отметим, что для неотрицательных аргументов функция тх может быть представлена в более простом виде тх(а,Ь,с) = ш1п(а,Ь,с). Построим итерационный процесс вида
(6 - 2^ = 3 ш1п(27^г_, + 5+, 2Т5+) - _1 - р^1, (12)
2^ = 35+ - 1, 2^ = - (13)
В качестве начального приближения используем значения производных: V0 = 0, 0 < г < I.
Покажем, что если на какой-либо итерации выполнены оценки
0 < _1 < вш1п(^_,5+), 0 < «к _1 < 0 < «к _1 < в^"", (14)
то на следующей итерации они также выполнены. Отметим, что из соотношений (14) как частный случай вытекают неравенства, которые будем использовать в дальнейшем анализе:
Заменяя в соотношениях (12) величины ö_ , ö+ их наименьшими значениями min(öj , ö+), а вычитаемые — их максимальными значениями из неравенства (15), получим
(6 - 2pi)vk > 3min(ö_, ö+) min(27, 2, 2y) - pß(ö_ + ö+). (16)
Из определения (8) вытекает соотношение p < 2ymin(ö+, ö_)/(ö+ + ö_). Увеличивая с его помощью вычитаемое в формуле (16) и учитывая что min(27, 2, 27) = 2, запишем
(6 - 2p>k > 6 min(ö_, ö+) - ß27 min(ö_, ö+).
Таким образом, если 6 - ß27 > 0, то левая часть оценок (14), т.е. ограниченность параметров v снизу нулем, выполнена на новой итерации во внутренних узлах.
Покажем положительность v для граничных узлов. Завышая вычитаемое в первой из формул (13) с помощью соотношения (14), получим
2v0fc > 3ö+ - 0öo+.
Аналогичное соотношение справедливо и для правого граничного узла. Следовательно, если 3 - ß > 0, то нужные неравенства имеются во всех узлах.
Докажем правую часть оценок (14), т.е. ограниченность производных сверху. В соотношении (12) отбросим неотрицательное вычитаемое и второй аргумент функции min. Поскольку и то и другое может лишь увеличить значение выражения, то имеем
(6 - 2рг)^к < 3min(27ö", 27ö+). (17)
Полагая в левой части соотношения (17) p = 1 и тем самым уменьшая ее, получим
vk < 1-57 min(ö_, ö+).
Таким образом, если ß > 1.57, то правая часть оценок (14) выполнена. Аналогично для граничных узлов в соотношениях (13), отбросив вычитаемое, получим правую часть оценок (14) на границе, если выполнено условие ß > 1.5.
Итак, итерационный процесс, сохраняющий исходные нужные знаки производных в узлах интерполяции, построен. Уточним допустимые значения параметров, встречающихся при определении сплайна. Сведя воедино все принятые ограничения на параметры сплайна, запишем систему неравенств
6 - ß27 > 0, ß > 1.5, ß > 1.5Y, 3 - ß > 0, Y > 1.
В соответствии с условиями эквивалентности модифицированного и исходного сплайнов (10) ,(11) следует выбирать максимальные допустимые значения параметра y. Это обеспечивает реализацию немодифицированного сплайна третьего порядка для максимальных значений шага сетки. Наибольшее значение этого параметра, допускаемое приведенной выше системой неравенств, дается значением y = л/2. Оно имеет место лишь в случае
ß = 1.5Y = 1.5л/2.
Построенный итерационный процесс можно описать формулой Vk = AVk-1 + B, где V — вектор-столбец с компонентами v; A = {а^ } — матрица, имеющая лишь две отличные от нуля диагонали, расположенные рядом с главной диагональю. Диагональ, расположенная выше главной, содержит числа -1/2, q1, ..., qi_1, где q = -p/(6 - 2p), а расположенная ниже — числа q1, ..., qi_1, -1/2. Прямое вычисление позволяет легко установить, что
норма = тахк^</ I I = тах(1/2, | 2^1 ..., | 2^/-1 |). Поскольку максималь-
но возможное значение | qi |= 1/4, то ||А||те = 1/2. Таким образом, итерационный процесс сходится. Выше показано, что на каждой итерации для неубывающих значений функции в узлах интерполяции производные сплайна в них неотрицательны. Следовательно, это имеет место и для предела итерационной последовательности. Аналогично доказывается неположительность производных для невозрастающих значений функции в узлах интерполяции.
Рис. 3. Монотонный кубический нелокальный сплайн.
Наконец, покажем, что нужный знак производной сплайна в узлах интерполяции обеспечивает тот же знак производной в любой точке рассматриваемого отрезка. Согласно формуле (3), запишем выражение для производной сплайна:
s(1) = 6£ (1 - £ Кг + (1 - £ )(1 - 3£) + v £ (3£ - 2), 0 < £ < 1.
Тождественные преобразования приводят к следующей форме этого соотношения: s(1) = £(1 - £)[65г- 2.5Vi_i£- 2.5Vi(1 - £)] + vt_i(1 - £)(1 - 3£ + 2.5£2) + v£(2.5£2- 2£ + 0.5).
Покажем, что каждое из слагаемых имеет знак производных сплайна в узлах интерполяции. Рассмотрим случай неубывающей интерполируемой функции. Для последних двух слагаемых нужный знак следует из неотрицательности многочленов, являющихся сомножителями этих слагаемых. Чтобы показать неотрицательность первого слагаемого, рассмотрим два неравенства, вытекающие как частный случай из соотношений (14) для
в = 1.5y = 1.5л/2:
Vi_i < 1.5л/25_, v < 1.5л/25_.
Умножая первое неравенство на неотрицательный множитель 2.5£, второе — на неотрицательный множитель 2.5(1 - £) и складывая их, получаем
2.5vi_1£ + 2.5^(1 - £) < 2.5 х 1.5Л/2£_ < 6¿_.
Таким образом, первое слагаемое, а значит, и производная сплайна неотрицательны во всех точках для неубывающей интерполируемой функции. Аналогично этому производная сплайна неположительна для невозрастающей функции, что и требовалось доказать.
На рис. 3 представлены результаты решения задачи об интерполяции функции, описывающей последовательно треугольный, прямоугольный и параболический участки. Как видно, качество интерполяции выше, чем у сплайнов (3)-(5) или (5), (6).
Список литературы
[1] Завьялов Ю.С., Квлсов Б. И., Мирошниченко В. Л. Методы сплайн-функций. М.: Наука, 1980.
[2] KOBKOV V. V. Interpolation Hermite Spline-function in Interval Analysis // Abstr. of Intern. Conf. on Interval Methods and Their Appl. in Global Optimization (INTERVAL'98), Nanjing, China, Apr. 20-23. 1998. P. 55.
[3] KVASOV B. I. Shape preserving C2 spline interpolation // Proc. Second Asian Math. Conf. 1995 / S. Tangmanee and E. Schulz (Eds). Singapore: World Sci. Publ. Co. Pte. Ltd., 1998. P. 430-439.
[4] DELBOURGO R., GREGORY J. A. Shape preserving piecewise rational interpolation // SIAM J. Numer. Sci. and Statist. Comput. 1985. V. 6, No. 4. P. 967-976.
[5] SPATH H. Spline Algorithms for Curves and Surfaces. Winnipeg: Utilitas Math. Publ. 1974.
[6] SOD G. A. Numerical Methods in Fluid Dynamics. Cambridge: Cambridge Univ. Press, 1985.
[7] Пинчуков В. И. О построении монотонных схем типа предиктор — корректор произвольного порядка аппроксимации // Мат. моделирование. 1991. Т. 3, № 9. С. 95-104.
[8] Harten A., ENGQUIST B., Osher S., CHAKRAVARTHY S. Uniformly high-order essentially non-oscillatory schemes. Pt III //J. Comput. Phys. 1987. V. 71. P. 279-309.
[9] JIANG G., SHU C.-W. Efficient implementation of weighted ENO schemes //J. Comput. Phys. 1996. V. 126. P. 202-228.
Поступила в редакцию 28 июня 1999 г., в переработанном виде — 23 марта 2000 г.
