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

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

Одной из задач, возникающих в археологии, является РАСЧЕТ ХАРАКТЕРИСТИК И  ОПИСАНИЕ ФОРМЫ керамических сосудов (объема, показания кривизны тулова, показания кривизны предплечья, максимального  диаметра, общей высоты), найденных при археологических раскопках. Эта задача решается на основе использования статистических методов, которые позволяют получить параметры функций (полиномы 3-ей степени), связывающих вышеперечисленные характеристики с размерами кувшинов, сделанными в определенных точках.

Учитывались тринадцать параметров сосуда: семь поперечных сечений и  шесть соответствующих расстояний между ними. Считается, что это минимальное число размеров, для удовлетворительного описания контура сосуда. Причем второе и пятое сечения берутся в месте минимального и максимального диаметра сосуда соответственно. Четвертое и шестое сечение следует помещать в окрестности максимальной выпуклости предплечья и тулова.

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

Из всех способов одномерной интерполяции можно выделить три основных - полиномиальная, рациональная и сплайн-интерполяция. При полиномиальной интерполяции функция представляется в виде полинома степени N-1, где N - число точек с известными значениями интерполируемой функции. При рациональной интерполяции функция представляется в виде отношения двух полиномов. При сплайн-интерполяции используется кусочно-кубическое представление функции. У каждого способа есть свои преимущества и недостатки.

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

Теперь о математической части. Пусть заданы точки 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 

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

На вход алгоритма (блок-схема см Приложение 1) подаются массивы x, y, на выходе получаем матрицу c:array[1..3,1..n]. Чтобы получить значение интерполируемой функции f в точке x необходимо в начале определить какому из отрезков разбиения она принадлежит (т.е. найти такое i, что xi  < x < xi+1 ) и затем воспользоваться формулой S(x) с сответствующими c1i , c2i , c3i , yi .