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

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

Одной из задач, возникающих в археологии при исследовании керамических предметов быта (горшков, кувшинов, корчаг и т.д.), является расчет их характеристик (объема, кривизны тулова, кривизны предплечья, максимального диаметра, общей высоты) и описание формы по численным значениям определенных параметров. Эти параметры зависит от типа изделия, например, для кувшинов это 7 поперечных сечений и 6 соответствующих расстояний между ними.

(Археологи считают, что эти параметры позволяют удовлетворительно описать форму сосуда. Второе и пятое сечения соответствуют минимальному и максимальному диаметру сосуда, четвертое и шестое сечения –  окрестности точек максимальной выпуклости предплечья и тулова).

Нетрудно видеть, что контур кувшина представляет собой график полинома третьей степени. Задача состоит в том. Чтобы получить параметры этого полинома, пользуясь значениями параметров, т.е. провести кривую третьего порядка через множество точек с координатами . 

Для решения этой задачи применяется  метод интерполяции кубическими  сплайнами.

Ранее на кафедре математической статистики КГУ был реализован алгоритм этого метода на языке Паскаль.

Теоретическое описание метода

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

В нашем случае

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

В нашем случае n=6.

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

Для вычисления коэффициентов c1_i , c2_i , c3_i , i = 1, 2, ..., n-1 решается система линейных уравнений, построенная из условия непрерывности производной S'(x) в узлах сетки и дополнительных краевых условий на вторую производную, которые имеют вид:

2S''₁  + b1 *S''₂ = b2  
b3 *S''_N-1  + 2S''_N  = b4 

Здесь возможны два случая. Случай первый, когда известны значения первой производной в краевых точках (y'₁  = y'(x₁ ), y'_n  = y'(x_n )), тогда следует положить:

b1  = 1,    b2  = (6/(x₂ -x₁ )) * ((y₂ -y₁ ) / (x₂ -x₁ )-y'₁ ), 
b3  = 1,    b4  = (6/(x_n -x_n-1 )) * (y'_n  - (y_n -y_n-1 )/(x_n -x_n-1 ))

Случай второй, когда  известны значения второй производной (y''₁  = y''(x₁ ), y''_n  = y''(x_n )), тогда полагаем:

b1  = 0,    b2  = 2*y''₁  
b3  = 0,    b4  = 2*y''_n 

Отметим, что точки x_i  должны быть отсортированы по возрастанию.

Таким образом, чтобы получить значение интерполируемой функции f в точке x необходимо в начале определить какому из отрезков разбиения она принадлежит (т.е. найти такое i, что x_i  < x < x_i+1 ) и затем воспользоваться формулой S(x) с соответствующими c1_i , c2_i , c3_i , y_i .