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

Интерполяция функции кубическими сплайнами
Интерполяция кубическими сплайнами - это быстрый, эффективный и устойчивый способ интерполяции функций, который является основным конкурентом полиномиальной интерполяции. В его основе лежит следующая идея - интервал интерполяции разбивается на небольшие отрезки, на каждом из которых функция задается полиномом третьей степени. Коэффициенты полинома подбираются так, что на границах интервалов обеспечивается непрерывность функции, её первой и второй производных. Также есть возможность задать граничные условия - значения первой или второй производной на границах интервала. Если значения одной из производных на границе известны, то задав их, мы получаем крайне точную интерполяционную схему. Если значения неизвестны, то можно положить вторую производную на границе равно нолю и получить достаточно хорошие результаты.
Теперь о математической части. Пусть заданы точки x1 , x2 , ..., xn  и соответствующие им значения y1 , y2 , ..., yn  функции f(x). На каждом из отрезков [xi , xi+1 ], i=1, 2, ..., n-1 функцию приближаем при помощи полинома третьей степени:
S(x) = yi  + c1i (x-xi ) + c2i (x-xi ) 2 + c3i (x-xi ) 3, xi  < x < xi+1
Для вычисления коэффициентов c1i , c2i , c3i , i = 1, 2, ..., n-1 решается система линейных уравнений, построенная из условия непрерывности производной S'(x) в узлах сетки и дополнительных краевых условий на вторую производную, которые имеют вид:
2*S''1  + b1 *S''2  = b2
b3 *S''N-1  + 2*S''N  = b4
здесь возможны два случая. Случай первый, когда известны значения первой производной в краевых точках (y'1  = y'(x1 ), y'n  = y'(xn )), тогда следует положить:
b1  = 1, b2  = (6/(x2 -x1 )) * ((y2 -y1 ) / (x2 -x1 )-y'1 ),
b3  = 1, b4  = (6/(xn -xn-1 )) * (y'n  - (yn -yn-1 )/(xn -xn-1 ))
Случай второй, когда известны значения второй производной (y''1  = y''(x1 ), y''n  = y''(xn )), тогда полагаем:
b1  = 0, b2  = 2*y''1
b3  = 0, b4  = 2*y''N
Описание подпрограммы Spline3BuildTable
Процедура Spline3BuildTable служит для постройки таблицы коэффициентов кубического сплайна по заданным точкам и граничным условиям, накладываемым на производные. В дальнейшем постренная таблица используется функцией Spline3Interpolate.
Входные параметры:
• N - число точек
• DiffN - тип граничного условия. 1 соответствует граничным условиям накладываемым на первые производные, 2 - на вторые.
• xs - массив абсцисс опорных точек с номерами от 0 до N-1. Точки могут быть неупорядочены по возрастанию, алгоритм отсортирует их перед расчетами.
• ys - массив ординат опорных точек с номерами от 0 до N-1.
• BoundL - левое граничное условие. Если DiffN равно 1, то первая производная на левой границе равна BoundL, иначе BoundL равна вторая производная.
• BoundR - аналогично BoundL
Выходные значения:
• ctbl - в этот массив помещается таблица коэффициентов сплайна. Эта таблица используется функцией Spline3Interpolate.
Описание подпрограммы Spline3Interpolate
Подпрограмма вычисляет значение сплайна в данной точке на основе таблицы коэффициентов, расчитанной в подпрограмме Spline3BuildTable.
Входные параметры:
• N - число точек
• C - таблица коэффициентов, построенная процедурой Spline3BuildTable
• X - точка, в которой ведем расчет
Результат - значение кубического сплайна, заданного таблицей C, в точке X.
Алгоритм взят с сайта БЧА НИВЦ МГУ.
Если нашли ошибку в алгоритме - сообщите!
Реализация алгоритма
Исходный код данного алгоритма доступен в версиях на C++, Visual Basic 6 и Delphi. Все версии идентичны по своей функциональности.
Для каждого из языков программирования выводится список файлов, содержащий ссылку на код алгоритма (в начале списка) и вспомогательные алгоритмы, если такие есть (выделены курсивом).
Если вы в первый раз посетили этот сайт, то:
• Скачайте и подключите библиотеку AP (12 КБ) - библиотеку классов и функций, которая необходима для работы программ с сайта. Архив содержит версии библиотеки для каждого из языков, представленных на сайте.
• Обязательно прочитайте FAQ. Обратите внимание на раздел, посвященный выбранному вами языку.
• И ещё - читайте комментарии. Сэкономите много времени, поскольку НУМЕРАЦИЯ ЭЛЕМЕНТОВ МАССИВОВ НЕ ВСЕГДА НАЧИНАЕТСЯ С НОЛЯ, а не все это замечают.
C++
• Интерполяция функции кубическими сплайнами (скачать, открыть в браузере)
Delphi
• Интерполяция функции кубическими сплайнами (скачать, открыть в браузере)
Visual Basic 6
• Интерполяция функции кубическими сплайнами (скачать, открыть в браузере)
Реализация алгоритма на AlgoPascal
Реализация алгоритма на AlgoPascal доступна для лучшего понимания сути алгоритма, если вы захотите разобраться в нем. Автоматический перевод позволяет получить работоспособную программу, но оригинал на AlgoPascal гораздо легче читать, поскольку он набирается человеком.
 открыть AP-файл в браузере
Аппроксимация и интерполяция данных
SPLINE, PPVAL,
MKPP, UNMKPP
 Интерполяция функции одной переменной кубическим сплайном
Синтаксис:
 yi = spline(x, y, xi) v = ppval(pp, xx)
 pp = spline(x, y) [breaks, coefs, l, k] = unmkpp(pp)
  pp = mkpp(breaks, coefs)
Описание:
Функция yi = spline(x, y, xi) интерполирует значения функции y в точках xi внутри области определения функции, используя кубические сплайны [1].
Функция pp = spline(x, y) возвращает pp-форму сплайна, используемую в М-файлах ppval, mkpp, unmkpp.
Функция v = ppval(pp, xx) вычисляет значение кусочно-гладкого полинома pp для значений аргумента xx.
Функция [breaks, coefs, l, k] = unmkpp(pp) возвращает характеристики кусочно гладкого полинома pp:
            breaks - вектор разбиения аргумента;
            coefs - коэффициенты кубических сплайнов;
            l = length(breaks) - 1;
            k = length(coefs)/l.
Функция pp = mkpp(breaks, coefs) формирует кусочно-гладкий полином pp по его характеристикам.
Пример:
Зададим синусоиду всего 10 точками и проведем интерполяцию кубическими сплайнами, используя мелкую сетку.
            x = 0:10; y = sin(x);
            xi = 0:.25:10;
            yi = spline(x, y, xi);
            plot(x, y, 'o', xi, yi, ‘g’), grid
 
Определим pp-форму сплайна.
           pp = spline(x, y);
           [breaks, coeffs, l, k] = unmkpp(pp)
           breaks = 0   1    2   3   4   5   6   7    8   9   10
           coeffs =
 -0.0419 -0.2612 1.1446 0
 -0.0419 -0.3868 0.4965 0.8415
 0.1469 -0.5124 -0.4027 0.9093
 0.1603 -0.0716 -0.9867 0.1411
 0.0372 0.4095 -0.6488 -0.7568
 -0.1234 0.5211 0.2818 -0.9589
 -0.1684 0.1509 0.9538 -0.2794
 -0.0640 -0.3542 0.7506 0.6570
 0.1190 -0.5463 -0.1499 0.9894
 0.1190 -0.1894 -0.8856 0.4121
          l =    10
          k =    4
Вычислим pp-форму в узловых точках сетки.
          v = ppval(pp,x)
          v =
0 0.8415 0.9093 0.1411 -0.7568 -0.9589 -0.2794 0.6570 0.9894 0.4121 -0.5440
Алгоритм:
Интерполяция сплайнами использует вспомогательные функции ppval, mkpp, unmkpp, которые образуют небольшой пакет для работы с кусочно-гладкими полиномами.
Существенно большие возможности для работы со сплайнами предо-ставляет пользователю специализированный пакет Spline Toolbox [2].
Сопутствующие функции: ICUBIC, INTERP1, POLYFIT, Spline Toolbox.
Интерполяция - приближение функции кривой, проходящей через все N точек. Основной недостаток интерполяционных алгоритмов в том, что при изменении значения функции в одной точке необходимо полностью пересчитать интерполяционные формулы.
Аппроксимация - приближение кривой, не обязательно проходящей через все точки. Основные методы аппроксимации обладают (и это очень ценно) свойством 'local control': изменение значения функции в одной точке влечет за собой перевычисление лишь 1-3 формул (это гораздо лучше, чем N формул, особенно в реальных приложениях компьютерной графики).