Главная       Скачать       Коммерческая поддержка       FAQ       Forum       О нас       Английская версия

Интерполяция сплайнами

Интерполяция сплайнами третьего порядка - это быстрый, эффективный и устойчивый способ интерполяции функций. Наравне с рациональной интерполяцией, сплайн-интерполяция является одной из альтернатив полиномиальной интерполяции.

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

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

Содержание

    1 Типы сплайнов
           Линейный сплайн
           Сплайн Эрмита
           Сплайн Катмулла-Рома
           Кубический сплайн
           Сплайн Акимы
    2 Сплайны в ALGLIB
           Основные операции
           Дифференцирование на сетке, переход от сетки к сетке
           Примеры
    3 Аппроксимация кубическими сплайнами
    4 Численные результаты
           Интерполяция различными видами сплайнов
    5 Manual entries
    6 Ссылки по теме

Типы сплайнов

Линейный сплайн

Линейный сплайн - это сплайн, составленный из полиномов первой степени, т.е. из отрезков прямых линий. Точность интерполяции линейными сплайнами невысока, также следует отметить, что они не обеспечивают непрерывности даже первых производных. Однако в некоторых случаях кусочно-линейная аппроксимация функции может оказаться предпочтительнее, чем аппроксимация более высокого порядка. Например, линейный сплайн сохраняет монотонность переданного в него набора точек.

Сплайн Эрмита

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

Сплайн Катмулла-Рома

Сплайн Катмулла-Рома - это сплайн Эрмита, производные которого определяются по формуле:

Как и сплайн Эрмита, сплайн Катмулла-Рома имеет непрерывную первую производную и разрывную вторую. Сплайн Катмулла-Рома локален - значения сплайна зависят только от значений функции в четырех соседних точках (двух слева, двух справа). Можно использовать два типа граничных условий:

Кубический сплайн

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

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

Сплайн Акимы

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

Важным свойством сплайна Акимы является его локальность - значения функции на отрезке [x, xi+1 ] зависят только от fi-2 , fi-1 , f, fi+1 , fi+2 , fi+3 . Вторым свойством, которое следует принимать во внимание, является нелинейность интерполяции сплайнами Акимы - результат интерполяции суммы двух функций не равен сумме интерполяционных схем, построенных на основе отдельных функций. Для построения сплайна Акимы требуется не менее 5 точек. Во внутренней области (т.е. между x и xN-3  при нумерации точек от 0 до N-1) погрешность интерполяции имеет порядок O(h 2).

Сплайны в ALGLIB

Основные операции

Работа с кубическими сплайнами в ALGLIB осуществляется в два этапа:

Для построения сплайна можно использовать следующие подпрограммы:

После того, как сплайн построен, вы можете использовать следующие подпрограммы:

Дифференцирование на сетке, переход от сетки к сетке

Одной из часто встречающихся задач является использование кубических сплайнов для перехода от одной сетки к другой. Мы имеем сетку x, значения y (на сетке x), а также новую сетку x. Нам надо получить значения кубического сплайна, построенного на основе старой сетки, в узлах новой сетки x. Дополнительно могут требоваться значения производных.

Разумеется, можно просто построить кубический сплайн при помощи функциии spline1dbuildcubic, после чего использовать функцию spline1ddiff, вызвав её для каждого из узлов новой сетки (см. пример ниже). Приведенный ниже код имеет трудоемкость O(n·log(n)).

C++ example


... initialization of x, y, x2 ...

alglib::spline1dbuildcubic(x, y, spline);
for(i=0; i<n2; i++)
    y2[i] = alglib::spline1dcalc(spline, x2[i]);
    
... now y2 contains new values ...

Однако есть более удобное и быстрое решение - использовать специальные функции для операций на сетках:

Результат работы этих функций эквивалентен результату, полученному при помощи приведенного выше кода, но принцип работы другой - используется специализированное представление кубического сплайна, оптимизированное для единственной операции - перехода от одной сетки к другой. Это позволяет значительно повысить быстродействие на предварительно отсортированных сетках - до O(max(n,n)).

Замечание #1
В случае, если сетки не отсортированы, трудоемкость увеличивается до O(n·log(n)+n·log(n)) в связи с проведением предварительной сортировки. Однако и это часто оказывается быстрее, чем первый вариант.

С использованием этих функций мы можем написать более компактный код:

C++ example


... initialization of x, y, x2 ...

alglib::spline1dconvcubic(x, y, x2, y2);
    
... now y2 contains new values ...

Примеры

ALGLIB Reference Manual содержит следующие примеры использования кубических сплайнов:

Аппроксимация кубическими сплайнами

Пакет ALGLIB также может использоваться для аппроксимации сплайнами. Более подробно этот вопрос рассмотрен в главе ALGLIB User Guide, посвященной аппроксимации.

Численные результаты

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

Интерполяция различными видами сплайнов

Рассмотрим следующую проблему: интерполяцию функции f(x)=sin(x) на интервале [-π/2,+π/2]. Мы используем равномерную сетку с N узлами и четыре различных типа сплайнов:

Исследуем зависимость погрешности от числа узлов, варьируя N в диапазоне 10...100. В результате мы получим следующий график:

Мы можем видеть, что в данном случае лучшим выбором является кубический сплайн с точными граничными условиями. Если значение производной в узловых точках неизвестно, то лучшим вариантом будет сплайн, завершающийся параболой. Два прочих варианта - сплайн Катмулла-Рома и кубический сплайн с естественными граничными условиями - дают значительно большую ошибку.

Этот вывод верен и для большинства других задач. Если вы знаете производную на границе, лучше использовать её. Если нет - используйте сплайн, завершающийся параболой.

Замечание #2
На графике не показан результат, полученный с помощью сплайна Эрмита. Обычно он очень близок по качеству к кубическому сплайну с точными граничными условиями, но требует больше информации. В частности, в рассмотренном выше случае графики для этих двух сплайнов полностью совпали бы.

Ссылки по теме

  1. 'Spline_interpolation', Wikipedia

Manual entries

C++ spline1d subpackage   
C# spline1d subpackage   

This article is intended for personal use only.

Скачать ALGLIB

C#

Исходный код на C#

Downloads page

 

C++

Исходный код на C++

Downloads page

 

C++, арифметика высокой точности

Исходный код на C++, использующий библиотеки MPFR/GMP.

Исходный код GMP доступен на сайте gmplib.org. Исходный код MPFR доступен на сайте www.mpfr.org.

Downloads page

 

FreePascal

Исходный код на Free Pascal.

Downloads page

 

Delphi

Исходный код на Delphi.

Downloads page

 

VB.NET

Исходный код на VB.NET.

Downloads page

 

VBA

Исходный код на VBA.

Downloads page

 

Python

Исходный код на Python (CPython и IronPython).

Downloads page

 

 

ALGLIB® - numerical analysis library, 1999-2016.
ALGLIB is registered trademark of the ALGLIB Project.