CC BY
44
УДК 519.67
DOI 10.14529/cmse210101
НАПРАВЛЕННЫЕ СПЛАЙНЫ И ИХ ИСПОЛЬЗОВАНИЕ ДЛЯ СГЛАЖИВАНИЯ ВЫБРОСОВ И ИЗЛОМОВ
ИНТЕРПОЛЯНТА
© 2021 В.А. Коднянко
Сибирский федеральный университет (660074 Красноярск, ул. Академика Киренского, д. 26А) E-mail: kowlad@rambler.ru Поступила в редакцию: 18.07.2020
Сформулирован и предложен метод построения направленного кубического сплайна для набора точек на плоскости. Проведено сравнение сплайна с Б-сплайном Шёнберга, сплайнами Акимы и Катмулла—Рома. Показано, что для неравноотстоящих точек в сравнении с Б-сплайном он дает значительно меньшие выбросы и практически лишен сильных изломов, которые свойственны сплайнам Акимы. Сплайн не дает петель и осцилляций, которые являются характерным недостатком параметрических сплайнов, в частности, эрмитовых, к числу которых относится сплайн Катмулла—Рома. Предложен быстрый метод оптимизации направляющего коэффициента сплайна, цель которой состоит в минимизации разрывов второй производной функции в ее промежуточных точках. Приведен пример оптимизации направленного сплайна третьего порядка. Также предложен направленный сплайн четвертого порядка, который лишен изломов. Сформулирован метод оптимизации направленного сплайна четвертого порядка, изложен алгоритм его оптимизации. Критериями оптимизации являются длина сплайна и наименьшее расстояние между его глобальными максимумом и минимумом. Показано, что в сравнении с сплайна Шёнберга направленный сплайн четвертого порядка имеет меньшие выбросы. Предложен метод автоматического притупления острых пиков кривых, который можно применять ко всем типам сплайнов.
Ключевые слова: сплайн, сплайн Шёнберга, сплайн Акимы, направленный сплайн.
ОБРАЗЕЦ ЦИТИРОВАНИЯ
Коднянко В.А. Направленные сплайны и их использование для сглаживания выбросов и изломов интерполянта // Вестник ЮУрГУ. Серия: Вычислительная математика и информатика. 2021. Т. 10, № 1. С. 5-19. DOI: 10.14529/cmse210101.
Введение
При проведении вычислительных или натурных экспериментов однофакторную связь между входной и выходной величинами обычно определяют посредством вычислений либо измерений. Результатом являются данные в виде набора точек на плоскости
(х$> Уо)< (xh yi)i ■■■ , (xn-h Уп-l)' (хт Уп)> xi-l < xi,i = 1,2, ... ,n. (1)
На его основе зачастую необходимо получить плавную кривую, которая должна проходить через экспериментальные точки.
Задачу можно решить посредством использования интерполяционных методов [1-3]. Для наборов с малым числом точек удовлетворительные результаты дают интерполяционные методы Лагранжа, Ньютона, Стирлинга [2-4]. Однако, с увеличением количества точек набора (1) интерполяционные полиномы в области крайних точек дают осцилляции недопустимо большой амплитуды [4].
Этого недостатка лишены сплайны [5-10, 20-29], которые находят широкое использование в практических приложениях, как-то: при численном решении нелинейных уравнений, сглаживании функций, сжатии и восстановлении графических изображений, ряде других применений. Наиболее известным среди них является Б-сплайн Шёнберга [5, 17]. Сплайн обеспечивает безупречную гладкость интерполянта для равноотстоящих точек х, (i = 0,1, ...,п), для которых
h, = х,— хг_( = h = const. (2)
Если (2) не выполняется, то плавность кривых, как правило, нарушается. В таких случаях 5-сплайн может давать значительные осцилляции кривой (так называемые «выбросы»), которые имеют место в области отрезков с малыми h,.
Избежать выбросов 5-сплайна позволяют сплайн Акимы [11, 16] и эрмитовы сплайны, частным случаем которых является применяемый в графической геометрии параметрический сплайн Катмулла—Рома [12, 13]. Однако и эти сплайны имеют свои недостатки. Так, сплайн Акимы часто дает неприемлемые изломы графика функции в узловых точках, которые видны на рис. 1 (кривая 2).
0,0 0,2 0,4 0,6 0,8 х 0,0 0,2 0,4 0,6 0,8 х
Рис 1. Графики 5-сплайна (1) Рис 2. Кривые 5-сплайна (1)
и сплайна Акимы (2) и сплайна Катмулла—Рома (2)
Более гладкие переходы между соседними полиномами дают сплайны Катмулла— Рома (рис. 2). Однако при малых Л, они могут образовывать петли в узловых точках. Кроме того, на параметрической зависимости х = х(Ь) которая в силу (1), очевидно, должна быть монотонно-возрастающей, могут появляться локальные экстремумы, что свидетельствует о явных недостатках алгоритма интерполяции при обработке наборов данных, содержащих точки, находящиеся на относительно малом расстоянии друг от друга.
Таким образом, актуальной является задача разработки сплайна, который характеризуется меньшими выбросами 5-сплайнов и существенно менее выраженными изломами сплайнов Акимы. Наряду с этим поставлена задача создания такой модели сплайна, которая не уступала бы по скорости вычислений алгоритму сплайна Акимы при изменении одной или нескольких точек, что в сравнении с 5-сплайном является характерным преимуществом данного сплайна.
Статья организована следующим образом. В разделе 1 рассмотрена методика построения направленного кубического сплайна (НСЗ-сплайна). В разделе 2 предложен метод и приведен пример оптимизации такого сплайна. Раздел 3 посвящен построению направленного сплайна четвертого порядка (НС4-сплайна). В разделе 4 изложена методика оптимизации направленного сплайна четвертого порядка по критерии минимума длины кривой и минимума расстояния между ее глобальными максимумом и минимумом. В раз-
деле 5 предложен метод притупления острых пиков сплайнов. Методика может быть применена ко всем типам рассмотренных сплайнов. В заключении приводится краткая сводка результатов, полученных в рамках данного исследования, и указаны направления дальнейших исследований.
Ниже рассмотрена методика построения и алгоритм оптимизации такого сплайна.
1. Построение направленного кубического сплайна (НС3-сплайна)
На каждом отрезке [Х1-(,Х1] будем представлять сплайн Б(х) полиномом третьей степени
Б,(х) = а, + Ь,(х — х,) + с,(х — х,)2 + di(x — х,)3, г = 1,2, ...,п. (3)
На стыках соседних отрезков [х-^х^ для полиномов (3) потребуем выполнения условий непрерывности сплайна и ее первой производной в точках (1)
Б,{х,_() = (4)
Б'1(х,_1) = Б1_1(х,.1). (5)
Используя (3)—(5), получим
а, = У,,
23
'а, — + ^с, — = у,+1, (6)
=
Если считать Ь, известными, то (6) позволяет получить систему уравнений относительно неизвестных коэффициентов с,,
Ь, — и, с, г ,
к (7)
И,
2с, — =
V
где
У,—У,-1
Решив (6), найдем
2Ь, + Ь,-1 — 3и, ^ Ь, + Ь,-1 — 2и,
с,= к, ,й,= и,2 . (8)
В простейшем случае для крайних отрезков [х0,х1] и [хп-1,хп] можно воспользоваться производной линейного интерполяционного сплайна, положив
Ь$ = и^ Ь* = и*.
Эти же коэффициенты можно определить при помощи построенного по трем точкам (х,-1,у,-1), (х,,у,), (х,+1,у,+1) интерполяционного полинома Лагранжа [2]
\и,к,+1 + и,+1к, + (и,+1 — и,)(х — х,)](х — х,) .
К(х) = У,+-г-ГГ-= 1,2, ...,п — 1,
к,+1 +
л к1(и2 — и1) 1 кп(ип — ип-1)
= А (Х$) = и1 Т ГГ ,Ьп = /п-1(хп) = ип + Т ГГ .
к1 + к2 кп-1 + кп
Для промежуточных отрезков [х,--1,х1] будем использовать формулу
Ь, = аи, + (1 — а)и,+1, (9)
где коэффициент а £ [0,1].
Из (4)—(9) следует, что на отрезке [х^+^х^] функция Б,(х) полностью определяется только тремя точками, в то время как локальность сплайна Акимы определятся пятью точками [11], а 5-сплайн всеми точками набора (1). Это свойство характеризует лучшее сравнительное быстродействие рассматриваемого метода при коррекции коэффициентов (3) на случай изменения отдельных точек набора (1).
Разделенная разность и, представляет собой угловой коэффициент соответствующего сегмента линейного интерполяционного сплайна, который показан на рис. 3.
у
2,0 1,0 0,0 -1,0 -2,0
(
I \ л
- -Л 1 л а
I V 1 —о ■ м
\ N \ \ V
\ 1 / п \ \ \
1 2
\ / V
V 1
-1
0,0 0,2 0,4 0,6 0,8 х
Рис 3. 5-сплайн (1) и линейный интерполяционный сплайн (2)
0,0 0,2 0,4 0,6 0,8 X
Рис 4. 5-сплайн (1), сплайн Акимы (2) и направленный сплайн (3)
Очевидно, (9) является угловым коэффициентом касательной в точке (х,,у,) сопряжения сегментов сплайна, находящихся по обе стороны от нее. Варьируя коэффициент а, можно отрегулировать положение касательных в промежуточных точках сплайна, направляя его так, чтобы выбросы были бы минимальными, а изломы не слишком заметными. Назовем такой сплайн направленным, а коэффициент а направляющим.
В простейшем случае можно положить а = 0,5, то есть считать, что направляющая касательная сплайна в промежуточных узловых точках должна занимать среднее положение относительно соседних отрезков линейного интерполяционного сплайна.
На рис. 4 показан пример интерполяции набора данных при помощи 5-сплайна (1), сплайна Акимы (2) и направленного сплайна (3).
Видно, что, в отличие от сплайна Акимы, направленный сплайн не имеет видимых изломов, чего для многих практических применений оказывается достаточным основанием для принятия решения об удовлетворительной аппроксимации исходной табличной функции подобными сплайнами. В сравнении с 5-сплайном менее выражены и выбросы сплайна.
2. Оптимизация НСЗ-сплайна
Оптимизация преследует цель сделать менее заметными изломы направленного сплайна, которые определяются абсолютной величиной разности значений вторых производных в точке сопряжения соседних полиномов. Величину такого разрыва в промежуточных точках можно определить формулой
Э,(а) = ^ Щ'^+г) - Б'+^+Л = \с,- с1_1
3d.il, (I = 1,2, ...,п — 1).
(10)
Вычислительный эксперимент показал, что в большинстве случаев при крайних значениях направляющего коэффициента а = 0 и а = 1, когда углы касательной к кривой в промежуточных точках и одного из отрезков линейного сплайна совпадают, направленный сплайн обычно имеет не только выраженные изломы, но и большие выбросы. При промежуточных значениях а эти недостатки менее заметны. Следовательно, существует такие значения а, при котором эти недостатки будут наименее выражены.
Суть оптимизации заключается в отыскании такого а = аори, при котором наибольший разрыв
Э(а) = тахЭ,(а),1 = 1,2,... ,п —1 (11)
кажется минимальным.
Установлено, что в большинстве случаев И (а) является кусочно-линейной функцией с единственной точкой излома. Однако нередки случаи функций с несколькими изломами.
Найти значение аори можно одним из методов минимизации унимодальных функций [14, 15, 17]. Однако можно воспользоваться свойством кусочной линейности И(а) и на этом основании предложить более быстрый алгоритм.
Методика решения задачи и описание алгоритма нахождения аори приведены ниже. Использованы величины типа Т = (Т.х,Т.у,Т.г), где Т.х, Т.у — абсцисса и ордината точки, Т. г — значение производной функции И(а) в этой точке.
Шаг 1. Зададим достаточно малое число £ — точность определения аори, а также границы интервала А.х = 0; В.х = 1.
Шаг 2. Найдем значение функции А.у = О (А. х) и А.г = [И (А. х + е) — А.у\/£ на левом конце интервала поиска. Определим В.у = И (В.х) и В.г = [В.у — И(В.х-£)]/£ — аналогичные им величины на правом конце интервала поиска.
Шаг 3. По этим точкам и производным построим прямые, первая из которых проходит через точку (А.х, А.у), вторая через (В.х, В. у). Нетрудно показать, что абсциссу точки пересечения этих прямых можно найти по формуле
В.у — А.у — В.гВ.х + А.гА.х , ,
х = —-—А-Б-. (12)
А.г — В.г
Вычислим абсциссу С.х = х + £/3, сдвинутую вправо от х на величину меньшую £. Это необходимо для того чтоб точка х попала в диапазон [А.х, С.х].
Шаг 4. Если \С.х — А.х\ < £ или \С.х — В.х\ < £, то точка минимума аори = х найдена и алгоритм заканчивает работу, иначе находим аналогичные значения С. у и С. г. Если С.г и А.г числа одного знака, то полагаем В = С, иначе А = С и переходим к шагу 2 для выполнения новой итерации.
Продемонстрируем работу алгоритма на примере оптимизации функции c двумя изломами.
Для этого на шаге 1 зададим точность поиска £ = 10+3. На шаге 2 получим А.х = 0;
A.у = 1,2; А.г = —1,7; В.х = 0,999; В.у = 0,6; В.г = 1,7. Разные знаки производных А.г и
B.г указывают на то, что функция И(а) строго унимодальна [18] и, следовательно, ее минимум находится внутри отрезка.
На шаге 3 по формуле (12) найдем х = 0,644; С.х = 0,645. Поскольку на шаге 4 ни одно из его условий не выполнилось, то вычислим С. у и С.г = 0,8 > 0. Это означает, что справа от точки х функция возрастает, следовательно, минимум находится слева от х. Положим В = С и поиск минимума функции продолжим на отрезке [0; 0, 645] меньшей длины.
На новой итерации получим х = 0,425; С.х = 0,426; С.г = 0,8 >0 и новый отрезок [0; 426].
На следующей итерации найдем х = 0,425; С.х = 0,426. Теперь условие 1С.Х — В.х1 < £ выполнилось. Это значит, что минимум функции находится в точке аори = х = 0,425.
Очевидно, число итераций не превышает к + 1, где к — число изломов функции И (а). В частности, данная задача решена за три итерации при двух изломах минимизируемой функции.
Приведем таблицу значений функции у(х) = smx, использованной в работе [19] для оценки погрешности Б-сплайна и сплайна Акимы. Сплайны построены на множестве из 14 случайных точек. Вычисления проведены для равноотстоящих узлов. Для сравнения в таблицу добавлены данные по направленному сплайну. В таблице в колонках Д1, Д2, Дз приведены разности между точным значением функции у(х) и значениями, которые получены с помощью Б-сплайна (Д1), сплайна Акимы (Д2) и направленного сплайна (Дз), соответственно.
Таблица
Сравнительные данные по вычислительной погрешности сплайнов
X У(х) Д1 Д2 Д3
0,00 0,0000000 0,000Е+00 0,000Е+00 0,000Е+00
0,05 0,0499792 7,819Е-08 -6,100Е-05 3,529Е-05
0,10 0,0998334 3,189Е-08 -1,010Е-05 -5,012Е-06
0,15 0,1494381 1,433Е-08 7,158Е-06 9,359Е-06
0,20 0,1986693 8,864Е-09 2,348Е-05 -5,406Е-06
0,25 0,2474040 2,986Е-08 -1,177Е-05 3,179Е-06
0,30 0,2955202 8,974Е-11 -1,624Е-06 6,155Е-07
0,35 0,3428978 2,180Е-08 6,243Е-06 -3,199Е-06
0,40 0,3894183 2,905Е-08 -8,113Е-06 4,885Е-06
0,45 0,4349655 -2,711Е-20 -2,711Е-20 -5,421Е-20
0,50 0,4794255 4,830Е-10 5,141Е-06 -4,438Е-06
0,55 0,5226872 8,824Е-08 -2,120Е-06 2,857Е-06
0,60 0,5646425 -6,664Е-08 -6,810Е-06 2,393Е-06
0,65 0,6051864 -1,044Е-07 4,300Е-06 -1,165Е-06
0,70 0,6442177 8,651Е-07 3,549Е-06 -2,821Е-06
0,75 0,6816388 -2,716Е-06 -1,707Е-06 1,296Е-06
Более высокую точность показал Б-сплайн. Среди двух последних лучшие результаты дал направленный сплайн. Это неочевидный результат, поскольку ожидалось, что якобы имеющий преимущества при интерполяции монотонных функций, а именно таковой в данном примере является функция у(х), сплайн Акимы должен был бы демонстрировать лучшие показатели не только в сравнении с направленным сплайном, но и с Б-сплайном. Однако в данном случае эти ожидания не оправдались.
Изложенная идея позволяет расширить рамки подхода к построению направленного нового сплайна, который будет лишен изломов, и допускает возможность его оптимизации с целью подавления выбросов. Примером может служить направленный сплайн четвертой степени, методика построения которого изложена ниже.
3. Построение направленного сплайна четвертой степени (НС4-сплайна)
На каждом отрезке [х,--1,х1] будем представлять сплайн Б(х) полиномом четвертой степени
Б,(х) = а, + Ь,(х — х,) + с,(х — х,)2 + — х,-1)3 + е,(х — х,)4,1 = 1,2, ...,п. (13)
На стыках соседних отрезков [х-^х^ для полиномов (3) потребуем выполнения условий непрерывности сплайна и ее первой и второй производных в точках (1)
Б,(х,-1) = Б,-1(х,-1), (14)
Б[(х,-1) = Би(х,-1), (15)
Б1'(х1-г>=5Ц1(х1-г>. (16) Полагая по-прежнему Ь, известными и, используя (14), получим
с, = (1,-1Ь., — е,1г2 + ж,, (17) — с, +
е, =-~2-, (18)
12
где
Ь, — и,
wi = ■
Используя (15), найдем
где
1,
3к2
2с1 + 4е112 = у1--I,-1, (19)
, , , , 1, ,
ь, — ь,-1
Подставив (18) в (19), запишем
с1 = 2(11111+3-21(11-1 + 2™1—Ц> (20)
21, 2
2 312-1 /о1 \
= 2" — — Щ. (21)
Условие (16) дает зависимость
с, + 6е1 = с,-1 + Зй-^, (22)
Подставив (20), (21) в (22), найдем
8(11 + (1,-1 ¡15 1011,-1)(1,-2 = + р,-1 — 4(2ы, + ш,-1). (23)
Сдвинув в (23) индекс, получим рекуррентные формулы
у.,(1,+1 + щй, + Л,й,-1 = Ш,Л = 1,2,... ,п 1, (24)
где
12 , Л , 312-1
щ = 8к,+1>= 5к, р3к— + 2),А, = —к—,ш, = 5У,+1 + V, — 4^,+1 + Wi).
Формулы (24) представляют трехдиагональную систему линейных алгебраических уравнений относительно неизвестных коэффициентов , которая с учетом очевидных краевых условий d0 = 0, dп = 0 может быть решена методом прогонки [12]. Далее коэффициенты с,, е, можно найти по формулам (17), (18).
4. Оптимизация НС4-сплайна
При оптимизации НС3-сплайна использована однопараметрическая процедура. Это объясняется тем, что многопараметрическая оптимизация в пределе дает Б-сплайн Шёнберга, что влечет утрату достоинств НС3. Сплайн НС4 лишен изломов, поэтому его оп-
тимизация сводится к только максимальному подавлению выбросов. В этом процессе могут быть задействованы все направляющие коэффициенты сплайна а,, при помощи которых вычисляются коэффициенты НС4-сплайна
Ь, = а,и, + (1 — а,)и,+1,1 = 1,2, ...,п — 1. (25)
В качестве критериев оптимальности данного сплайна использовали
- длину Ь сплайна,
- разность Я между его глобальным максимумом и глобальным минимумом. Критерий Ь определяется суммой длины сегментов сплайна и может быть вычислен
при помощи известной формулы [23]
п __
Ь = Х У ^1 + &(х)]2Лх,
где
Б'(х) = Ь, + 2с,(х — х,) + 3й,(х — х,+()2 + 4е,(х — х,)3. (26)
Для вычисления критерия Я использована формула (26), а также формула
Б''(х) = 2с, + 6а,(х — х,) + 12е,(х — х,)2. (27)
Используя (26), находили корни уравнения
Б'(х) = 0,
их тип контролировали с помощью (27).
Таким образом, как критерий Ь, так и критерий Я являются функциями многих переменных
К = К (а), (28)
где
а = (аъ a.2,..., а*+г).
Очевидно, 5-сплайн Шёнберга является частным случаем НС4-сплайна. Следовательно, с позиций минимума указанных критериев оптимальный сплайн не может быть хуже сплайна Шёнберга.
В процессе минимизации критериев придерживались требования поддержания такой формы сплайнов, изогеомертия которых соответствовала бы сплайну Шёнберга [25].
Расчеты НС4-сплайна показали, что без принятия мер сплайн может терять указанную форму.
Среди причин потери формы выделены следующие:
- сплайн может иметь «выпячивания» отдельных фрагментов кривой,
- на кривой могут появляться новые локальные экстремумы,
- могут также появляться новые точки смены кривизны сплайна.
Данные недостатки связаны с появлением новых локальных экстремумов, как самого сплайна, так и его первой и второй производных, а также новых смен знака его третей производной. При оптимизации варианты таких сплайнов отбраковывались.
Численные эксперименты позволили установить тот факт, что функция (28) является многоэкстремальной, то есть имеет множество локальных минимумов, среди которых следует найти глобальный минимум, который соответствует сплайну оптимальной формы.
Так, при п = 12, для которого проводилось большинство экспериментов, пришлось бы отыскивать глобальный минимум функции 11 переменных, что представляется задачей чрезвычайной сложности. В общем случае трудности получения решения такой задачи ставят под сомнение ценность практического использования обсуждаемого сплайна. Выход из положения был найден в использовании упомянутых свойств сплайна, которые продиктованы жесткими условиями сохранения его формы.
Методика минимизации критериев состоит в следующем. Изначально устанавливается стартовое состояние, для которого выбирается вектор а, все компоненты которого полагаются равными 0,5, и вычисляется стартовый НС-сплайн. Далее назначается шаг т малой длины и проводится переборная однопараметрическая минимизация по каждой компоненте вектора а, для которой процесс начинается со стартового состояния. Полученные таким образом лучшие покоординатные сплайны дают множество а-векторов, число которых т < п. Наблюдения показали, что через фильтр жестких требований сохранения формы сплайна обычно проходит не более половины начальных стартов. Так, например, для п = 12 обычно т <7. Следующим шагом является аналогичная однопараметрическая оптимизация, где в качестве стартов последовательно используются векторы, прошедшие фильтр первого шага. При этом количество новых векторов, которые прошли фильтр сохранения формы также невелико и оно обычно меньше аналогичного количества предыдущего шага. Рекурсия ведется до тех пор, пока отфильтрованные векторы не дадут новые векторы для продолжения процесса. Результатом оптимизации будет сплайн с наименьшим значением критерия К. Расчеты показали, что, например, для п = 12 обычно требуется сформировать сплайн и вычислить его характеристики 1000-2000 раз.
На рис. 5 показаны кривые сплайнов, которые дают типичную картину формы оптимального по критерию длины НС4-сплайна. Сплайн обычно занимает среднее положение между сплайном Шёнберга и НС3 ближе к первому.
Приведем характерный пример процесса оптимизации НС4. Длина сплайна Шёнберга 19,70. Исходный НС4 имеет несколько большую длину — 19,90. Оптимизированный сплайн имеет длину 16,82, что на 14,6% меньше длины сплайна Шёнберга. При проведении оптимизации алгоритм улучшал результат 181 раз. Для решения задачи потребовалось провести вычисление сплайна 1312 раз. Наименьшую длину имел НС3. При длине 11,58 он почти в два раза короче сплайна Шёнберга. Тысячи вычислительных экспериментов, поставленных для п = 12, показали, что оптимизированный НС4 короче сплайна Шёнберга на 5-50%, а НС3 — в 1,5-3 раза.
Примеры процессов оптимизации можно наблюдать в динамике на видео по ссылкам [30-32].
Исследования показывают, НС4 позволяет ослабить проявление выбросов сплайна Шёнберга. Однако лучшим в этом отношении следует признать СП3, на котором выбросов нет, а изломы проявляются весьма слабо.
В процессе изучения свойств предложенных сплайнов обнаружено, что кроме выбросов и изломов сплайны могут иметь острые пики локальных максимумов и минимумов, которые в ряде случаев также следует рассматривать как недостатки интерполяции. Ниже предложен и описан способ, который позволять притупить пики экстремумов.
5. Методика притупления пиков сплайнов
Пусть необходимо притупить локальный минимум какого-либо из рассмотренных сплайнов и пусть (хс,ус) точка такого минимума. Рассмотрим также точки перегиба (ха,Уа) и (%ь,Уь), расположенные по обе стороны от экстремума.
Для сглаживаемого экстремума сплайна введем добавочную функцию
д(х) = в(х — ха)к(хь — х)т, в >0,к>2,т>2. (29)
Для нее имеют место очевидные краевые условия
д'(Ха) = д'(хь) = 0,д"(Ха) = д"(хь) = 0.
Следовательно, на краях интервала функция не привносит изломов и способствует притуплению сплайна в области его минимума.
2,0
1,0
0,0
-1,0
-2,0
2 N
1
/ \
3 \ ч4
0,0
0,2
0,4
0,6
Рис 5. Б-сплайн (1), НСЗ-сплайн (2), НС4-сплайн (3)
Максимум этой функции имеет место при
Рис 6. Линейный интерполяционный сплайн (1), Б-сплайн (2), НС3 (3) и НС3 с притупленными экстремумами (4)
д'{хс) — G(xc - ха)к 1(хь - хс)т 1[к{хь - хс) - т{хс - ха)] — 0,
(30)
откуда следует
т = tk, t —
(xb-xc) (хс—ха)
Если t < 1, то примем т = 2 + £, где £ — малое число. Тогда к — m/t > 2. Иначе при t >1 примем к — 2 + £. Тогда т = kt > 2.
Обратимся к линейному интерполяционному сплайну
L(x) — LaX + Ьь,Ь(Ха) — Уа, Шь) = Уь,
где
Уь-Уа ХьУа - ХаУь
La — , ^ь — .
Хь Ха
Хь Ха
(31)
Для точки рассматриваемого минимума сплайна имеем
в(хс - ха)к{хъ - хс)т = а[Ь(хс) - ус],
где а — коэффициент притупления пика.
При а = 0 добавка д(х) отсутствует, при о = 1 суммарный сплайн ц(х) = з(х) + д(х) коснется отрезка прямой Ь(х), что является предельным случаем, поскольку при а > 1 суммарный сплайн получит новый экстремум-максимум, по обе стороны которого будет образовано два новых минимума. Такие режимы не обеспечивают сохранения формы сплайна, потому 0 < а < 1.
Для сохранения формы суммарный сплайн кроме того должен иметь кривизну одного знака до о = отах < 1. Данная величина может быть определена из условия
ц"(хс) = 0.
Выполнив дифференцирование функции д(х), найдем
5"(хс) + в(к + т)в(хс - ха)к~1(хь - хс)т~1 = 0.
Отсюда
8"(Хс)
G —
(к + т)(хс - ха)к—1(хь - хс)т—1'
(32)
Из (31), (32) вытекает, что
[L(xc)-yc] s"(xc)
(xc-Xa)k~1(xb-Xc)<-(~ (к + т)(хс-Ха)к-1(хь-Хс)<-(- ( )
Отношение (33) позволяет определить предельно допустимую величину а
S c)
°тах = (k + m)[L(Xc)-yc] '
Аналогично проводится притупление локальных максимумов.
Решение о подлежащих притуплению экстремумах, может быть принято, как в автоматизированном, так и автоматическом режимах. При автоматическом принятии решения в качестве критерия оценки экстремума можно использовать величину ке = s"(xc).
Приведем описание методики оценивания экстремумов, подлежащих автоматическому притуплению. В качестве ориентира будем использовать НС3, график которого дан на рис. 4. При визуальной оценке сплайна можно выделить два экстремума, которые необходимо притупить. Один из них является максимумом с абсциссой х = 0,395, второй соседним минимумом с абсциссой х = 0,440. Для первого ке = -2920, для второго ке = 1871. Следующий за ними по убыванию абсолютной величины критерия ке является крайний правый минимум с абсциссой х = 0,946, для которого ке = 1123. Оценим его как экстремум, не требующий притупления. Данный простой анализ позволяет в первом приближении сформулировать следующую методику автоматической экспертной оценки: экстремумы, которые должны быть притуплены автоматически, должны удовлетворять условию 1ке1 > 1200, если это не противоречит приведенным выше условиям сохранения формы сплайна.
На рис. 6 в качестве примера притупления экстремумов показаны линейный интерполяционный сплайн, 5-сплайн, НС3-сплайн и два фрагмента НС3-сплайна, на котором экстремумы притуплены автоматически в соответствии с описанной выше методикой. При вычислениях принята величина коэффициента притупления а = 0,75атах. Методика притупления острых пиков кривых может быть применена к сплайнам любого типа.
Заключение
В статье сформулирован и предложен простой метод построения кубического сплайна для набора точек на плоскости. Проведено сравнение сплайна с 5-сплайном Шёнберга и сплайнами Акимы и Катммула—Рома. Показано, что для неравноотстоящих точек, при которых обычно проявляются недостатки названных сплайнов, в сравнении с 5-сплайном он дает значительно меньшие осцилляции. Сплайн с таким набором точек практически лишен сильных изломов, которые свойственны сплайнам Акимы. Он не дает петель и осцилляций, которые являются характерным недостатком параметрических сплайнов, в частности, эрмитовых, к числу которых относится сплайн Катмулла—Рома. Предложен метод оптимизации направляющего коэффициента сплайна, цель которой состоит в минимизации разрывов второй производной функции в ее промежуточных точках. Также предложен направленный сплайн четвертого порядка, который лишен изломов и в сравнении со сплайном Шёнберга имеет меньшие выбросы. Предложен метод притупления острых пиков кривых, который можно применять ко всем типам сплайнов. Приведенные в статье численные и графические результаты получены на основе предложенных алгоритмов, реализованных при помощи разработанного автором программного обеспечения в среде визуального программирования Delphi с использованием приложения MS Office Excel.
В ходе дальнейшей работы планируется улучшить реализацию алгоритма оптимизации направленного сплайна четвертого порядка с целью сокращения времени его работы и повышения точности процедуры минимизации критериев оптимальности сплайна.
Литература
1. Powell M.J.D. Approximation Theory and Methods. Cambridge University Press, 1981. 352 p. DOI: 10.1017/СВ09781139171502.
2. Atkinson K.A. An Introduction to Numerical Analysis (2nd ed.). John Wiley and Sons, 1988. 615 p.
3. Watson G.A. Approximation Theory and Numerical Methods. John Wiley, 1980. 229 p.
4. Schatzman M. Numerical Analysis: A Mathematical Introduction. Clarendon Press, 2002. 496 p.
5. Schoenberg I.J. Contributions to the problem of approximation of equidistant data by analytic functions // Quart. Appl. Math. 1946. Vol. 4. P. 45-99, 112-141.
6. Ahlberg J.H., Nilson E.N., Walsh J.L. The Theory of Splines and Their Application. Academic Press, 1967. 296 p.
7. David F., Rogers J., Adams A. Mathematical Elements for Computer Graphics. McGraw-Hill Science / Engineering / Math, 2 edition, 1989. 611 p.
8. Cohen D. Incremental Methods for Computer Graphics. PhD Thesis. Harvard University, 1969.
9. Warnock J.E. A Hidden Surface Algorithm for Computer-Generated Halftone Pictures. Computer Science Department, University of Utah, TR 1-15. 1969.
10. Watkins G.S. A Real-Time Visible Surface Algorithm. Computer Science Department, University of Utah, UTECH-CSC-70-101, 1970.
11. Akima H. A New Method of Interpolation and Smooth Curve Fitting Based on Local Procedures // Journal of the ACM. 1970. Vol. 17, no. 4. P. 589-602. DOI: 10.1145/321607.321609.
12. Catmull E., Rom R. A class of local interpolating splines // Computer Aided Geometric Design. 1974. P. 317-326.
13. Barry P.J., Goldman R.N. Recursive evaluation algorithm for a class of Catmull-Rom splines // Computer Graphics. 1988. Vol. 22, no. 4. P. 199-204. DOI: 10.1145/378456.378511.
14. Kiefer J.K. Sequential minimax search for a maximum // P. Am. Math. Soc. 1953. P. 502506.
15. Brent R.P. Algorithms for Minimization without Derivatives. Dover. 2002. 195 p.
16. Круковец А.С., Горелкин Г.А. Разработка метода интерполяции значений номограммы // Современные научные исследования и инновации. 2015. № 5(2). URL: http://web.snauka.ru/issues/2015/05/53846 (дата обращения: 14.06.2020).
17. Dobson A.J. An Introduction to Statistical Modelling. Chapman and Hall, London, 1983.
18. Коднянко В.А. О вычислительной избыточности метода дихотомии и условной минимизации унимодальных функций методом экономной дихотомии // Системы и средства информатики. 2019. Т. 29, № 1. С. 164-173. DOI: 10.14357/08696527190113.
19. Ruckdeschel F.R. Basic scientific subroutines. Vol. 2. BYTE/McGRAW-HILL, 1981.
20. Яненко Н.Н., Квасов Б.И. Итерационный метод построения поликубических сплайн-функций // Докл. АН СССР. 1970. Т. 195, № 5. С. 1055-1057.
21. Constantini P., Morandi R. An algorithm for computing shape-preserving cubic spline interpolation to data // Calcolo. 1984. Vol. 21. P. 295-305.
22. Рябенький В.С. Локальные формулы гладкого восполнения и гладкой интерполяции функций по их значениям в узлах неравномерной прямоугольной сетки // Препринт. Ин-т прикл. математики АН СССР. 1974. № 21. 35 с.
23. Dietze S., Schmidt J.W. Determination of shape preserving spline interpolants with minimal curvature via dual programs // J. Approxim. Theory. 1988. Vol. 52, no. 1. P. 43-57.
24. Завьялов Ю.С., Квасов Б.И., Мирошниченко В.Л. Методы сплайн-функций. М.: Наука, 1980. 352 с.
25. Мирошниченко В.Л. Изогеометрические свойства и погрешность аппроксимации взвешенных кубических сплайнов // Вычислительные системы. Новосибирск: ИМ СО РАН, Вып. 154: Сплайны и их приложения. 1995. С. 127-154.
26. Корнейчук Н.П., Бабенко В.Ф., Лигун A.A. Экстремальные свойства полиномов и сплайнов. Киев: Наукова думка. 1992. 304 с.
27. Dzyubenko G.A., Gilewicz J., Shevchuk I.A. New phenomena in coconvex approximation // Analysis Mathematica. 2006. Vol. 32, no. 2. P. 113-121. DOI: 10.1007/s10476-006-0005-x.
28. Стечкин С.Б., Субботин Ю.Н. Сплайны в вычислительной математике. М.: Наука, 1976. 248 с.
29. Волков Ю.С. Новый способ построения интерполяционных кубических сплайнов // Журн. вычисл. матем. и матем. физ. 2004. Т. 44, № 2. C. 231-241.
30. Коднянко В.А. Video DS4-spline optimization 1. URL: http://smiuk.sfu-kras.ru/ kodnyanko/site/science/Video1.mp4 (дата обращения: 12.07.2020).
31. Коднянко В.А. Video DS4-spline optimization 2. URL: http://smiuk.sfu-kras.ru/ kodnyanko/site/science/Video2.mp4 (дата обращения: 12.07.2020).
32. Коднянко В.А. Video DS4-spline optimization 3. URL: http://smiuk.sfu-kras.ru/ kodnyanko/site/science/Video3.mp4 (дата обращения: 12.07.2020).
Коднянко Владимир Александрович, д.т.н, профессор, кафедра стандартизации, метрологии и управления качеством, Сибирский федеральный университет (Красноярск, Российская Федерация)
DOI: 10.14529/cmse210101
DIRECTIONAL SPLINES AND THEIR USE FOR SMOOTHING EJECTIONS AND FRACTURES
OF INTERPOLANT
© 2021 V.A. Kodnyanko
Siberian Federal University (Kirensky 26A, Krasnoyarsk, 660074 Russia)
E-mail: kowlad@rambler.ru Received: 18.07.2020
A method for constructing a directional cubic spline for a set of points on a plane is formulated and proposed. The spline is compared with the Schoenberg B-spline, Akima and Catmull-Rom splines. It is shown that for unequally spaced points, in comparison with the B-spline, it gives significantly lower overshoots and is practically free of strong kinks, which are characteristic of Akima splines. The spline does not give loops and oscillations, which are a characteristic drawback of parametric splines, in particular, Hermitian ones, which include the Catmull-Rom spline. A fast method for optimizing the spline guiding coefficient is proposed, the purpose of which is to minimize the discontinuities of the second derivative of the function at its intermediate points. An example of optimization of a directional third-order spline is given. A fourth-order directional spline, which is free of kinks, is also proposed. The method of optimization of the directional spline of the fourth order is formulated, the algorithm of its optimization is stated. The optimization criteria are the spline length and the smallest distance between its global maximum and minimum. It is shown that, in comparison with the Schoenberg spline, the fourth-order directional spline has smaller outliers. A method for automatic blunting of sharp peaks of curves is proposed, which can be applied to all types of splines.
Keywords: spline, Schoenberg spline, Akima spline, directional spline.
Направленные сплайны и их использование для сглаживания выбросов и изломов... FOR CITATION
Kodnyanko V.A. Directional Splines and Their Use for Smoothing Ejections and Fractures of Interpolant. Bulletin of the South Ural State University. Series: Computational Mathematics and Software Engineering. 2021. Vol. 10, no. 1. P. 5-19. (in Russian) DOI: 10.14529/cmse210101.
This paper is distributed under the terms of the Creative Commons Attribution-Non Commercial 3.0 License which permits non-commercial use, reproduction and distribution of the work without further permission provided the original work is properly cites.
References
1. Powell M.J.D. Approximation Theory and Methods. Cambridge University Press, 1981. 352 p. DOI: 10.1017/CB09781139171502.
2. Atkinson K.A. An Introduction to Numerical Analysis (2nd ed.). John Wiley and Sons, 1988. 615 p.
3. Watson G.A. Approximation Theory and Numerical Methods. John Wiley, 1980. 229 p.
4. Schatzman M. Numerical Analysis: A Mathematical Introduction. Clarendon Press, 2002. 496 p.
5. Schoenberg I.J. Contributions to the problem of approximation of equidistant data by analytic functions. Quart. Appl. Math. 1946. Vol. 4. P. 45-99, 112-141.
6. Ahlberg J.H., Nilson E.N., Walsh J.L. The Theory of Splines and Their Application. Academic Press, 1967. 296 p.
7. David F., Rogers J., Adams A. Mathematical Elements for Computer Graphics. McGraw-Hill Science / Engineering / Math, 2 edition, 1989. 611 p.
8. Cohen D. Incremental Methods for Computer Graphics. PhD Thesis. Harvard University, 1969.
9. Warnock J.E. A Hidden Surface Algorithm for Computer-Generated Halftone Pictures. Computer Science Department, University of Utah, TR 1-15. 1969.
10. Watkins G.S. A Real-Time Visible Surface Algorithm. Computer Science Department, University of Utah, UTECH-CSC-70-101, 1970.
11. Akima H. A New Method of Interpolation and Smooth Curve Fitting Based on Local Procedures. Journal of the ACM. 1970. Vol. 17, no. 4. P. 589-602. DOI: 10.1145/321607.321609.
12. Catmull E., Rom R. A class of local interpolating splines. Computer Aided Geometric Design. 1974. P. 317-326.
13. Barry P.J., Goldman R.N. Recursive evaluation algorithm for a class of Catmull-Rom splines. Computer Graphics. 1988. Vol. 22, no. 4. P. 199-204. DOI: 10.1145/378456.378511.
14. Kiefer J.K. Sequential minimax search for a maximum. P. Am. Math. Soc. 1953. P. 502506.
15. Brent R.P. Algorithms for Minimization without Derivatives. Dover, 2002. 195 p.
16. Krukovets A.S., Gorelkin G.A. Development of a method for interpolating the values of the nomogram. Modern scientific research and innovations. 2015. No. 5(2). URL: http://web.snauka.ru/issues/2015/05/53846 (accessed: 14.06.2020) (in Russian)
17. Dobson A.J. An Introduction to Statistical Modelling. Chapman and Hall, London, 1983.
18. Kodnyanko V.A. On computational redundancy of the dichotomous search and conditional minimization of unimodal functions by the economical dichotomous search. Systems and
means of informatics. 2019. Vol. 29, no. 1. P. 164-173. DOI: 10.14357/08696527190113 (in Russian)
19. Ruckdeschel F.R. Basic scientific subroutines. Vol. 2. BYTE/McGRAW-HILL, 1981.
20. Yanenko N.N., Kvasov B.I. An Iterative Method for Constructing Polycubic Spline Functions. Dokl. USSR Academy of Sciences. 1970. Vol. 195, no. 5. P. 1055-1057. (in Russian)
21. Constantini P., Morandi R. An algorithm for computing shape-preserving cubic spline interpolation to data. Calcolo. 1984. Vol. 21. P. 295-305.
22. Ryabenky V.S. Local formulas for smooth completion and smooth interpolation of functions by their values at nodes of an uneven rectangular grid. Preprint, Inst. mathematics of the Academy of Sciences of the USSR. IPM. 1974. No. 21. (in Russian)
23. Dietze S., Schmidt J.W. Determination of shape preserving spline interpolants with minimal curvature via dual programs. J. Approxim. Theory. 1988. Vol. 52, no. 1. P. 43-57.
24. Zavyalov Yu.S., Kvasov B.I., Miroshnichenko V.L. Methods of spline functions. Moscow, Nauka, 1980. (in Russian)
25. Miroshnichenko V.L. Isogeometric properties and approximation error for weighted cubic splines. Computing Systems. Novosibirsk: IM SB RAS, 1995. Vol. 154. P. 127-154. (in Russian)
26. Korneychuk N.P., Babenko V.F., Ligun A.A. Extreme properties of polynomials and splines. Kiev, Naukova Dumka, 1992. 304 p. (in Russian)
27. Dzyubenko G.A., Gilewicz J., Shevchuk I. A. New phenomena in coconvex approximation. Analysis Mathematica. 2006. Vol. 32, no. 2. P. 113-121. DOI: 10.1007/s10476-006-0005-x.
28. Stechkin S.B., Subbotin Yu.N. Splines in computational mathematics. Moscow, Nauka, 1976. 248 p. (in Russian)
29. Volkov Yu.S. A new method for constructing interpolation cubic splines. Computational Mathematics and Mathematical Physics. 2004. Vol. 44, no. 2. P. 231-241. (in Russian)
30. Kodnyanko V.A. Video DS4-spline optimization 1. URL: http://smiuk.sfu-kras.ru/ kodnyanko/site/science/Video1.mp4 (accessed: 12.07.2020).
31. Kodnyanko V.A. Video DS4-spline optimization 2. URL: http://smiuk.sfu-kras.ru/ kodnyanko/site/science/Video2.mp4 (accessed: 12.07.2020).
32. Kodnyanko V.A. Video DS4-spline optimization 3. URL: http://smiuk.sfu-kras.ru/ kodnyanko/site/science/Video3.mp4 (accessed: 12.07.2020).
