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

Пусть некоторая функция f(x) задана на отрезке $[a,b]$, разбитом на части $[x_{i-1},x_{i}]$, $a=x_{0}<x_{1}<...<x_{N}=b$. Кубическим сплайном называется функция S(x), которая:

• на каждом отрезке $[x_{i-1},x_{i}]$ является полиномом третьей степени;
• имеет непрерывную вторую производную на всём отрезке $[a,b]$;
• в точках $x_{i}$ выполняется равенство $S(x_{i})=f(x_{i})$;
• $S''(a)=S''(b)=0$.

По построению сплайн S(x) интерполирует функцию f в точках $x_{i}$.

Теорема: Существует ровно одна функция S(x), удовлетворяющая этим условиям.

Построение

Обозначим: $h_{i}=x_{i}-x_{i-1}$

На каждом отрезке $[x_{i},x_{i+1}]$ функция $S(x)$ есть полином третьей степени $S_{i}(x)$, коэффициенты которого надо определить. Запишем для удобства $S_{i}(x)$ в виде:

$S_{i}(x)=a_{i}+b_{i}(x-x_{i})+{c_{i} \over 2}(x-x_{i})^{2}+{d_{i} \over 6}(x-x_{i})^{3}$

тогда

$S_{i}\left(x_{i}\right)=a_{i},\quad S'_{i}(x_{i})=b_{i},\quad S''_{i}(x_{i})=c_{i}$

Условия непрерывности всех производных до второго порядка включительно записываются в виде
$S_{i}\left(x_{i-1}\right)=S_{i-1}(x_{i-1})$
$S'_{i}\left(x_{i-1}\right)=S'_{i-1}(x_{i-1})$
$S''_{i}\left(x_{i-1}\right)=S''_{i-1}(x_{i-1})$

а условия инетрполяции в виде 

$S_{i}\left(x_{i-1}\right)=f(x_{i-1})$

Отсюда получаем формулы для вычисления коэффициентов сплайна:

$a_{i}=f\left(x_{i}\right)$

$h_{i}c_{i-1}+2(h_{i}+h_{i+1})c_{i}+h_{i+1}c_{i+1}=6\left({{f_{i+1}-f_{i}} \over {h_{i+1}}}-{{f_{i}-f_{i-1}} \over {h_{i}}}\right)$

$d_{i}={{c_{i}-c_{i-1}} \over {h_{i}}}$

$b_{i}=-{1 \over 2}h_{i}c_{i}-{1 \over 6}h_{i}^{2}d_{i}+{{f_{i}-f_{i-1}} \over {h_{i}}}$

Если учесть, что $c_{0}=c_{n}=0$, то вычисление с можно провести с помощью метода прогонки для трехдиагональной матрицы.

Реализация на языке C++

// Интерполирование функций естественными кубическими сплайнами
// Разработчик: Назар Андриенко Email: nuzikprogrammer@gmail.com

#include <cstdlib>
#include <cmath>
#include <limits>

class cubic_spline
{
private:
	// Структура, описывающая сплайн на каждом сегменте сетки
	struct spline_tuple
	{
		double a, b, c, d, x;
	};

	spline_tuple *splines; // Сплайн
	std::size_t n; // Количество узлов сетки

	void free_mem(); // Освобождение памяти

public:
	cubic_spline();
	~cubic_spline();

	// Построение сплайна
	// x - узлы сетки, должны быть упорядочены по возрастанию, кратные узлы запрещены
	// y - значения функции в узлах сетки
	// n - количество узлов сетки
	void build_spline(const double *x, const double *y, std::size_t n);

	// Вычисление значения интерполированной функции в произвольной точке
	double f(double x) const;
};

cubic_spline::cubic_spline() : splines(NULL)
{

}

cubic_spline::~cubic_spline()
{
	free_mem();
}

void cubic_spline::build_spline(const double *x, const double *y, std::size_t n)
{
	free_mem();

	this->n = n;

	// Инициализация массива сплайнов
	splines = new spline_tuple[n];
	for (std::size_t i = 0; i < n; ++i)
	{
		splines[i].x = x[i];
		splines[i].a = y[i];
	}
	splines[0].c = splines[n - 1].c = 0.;

	// Решение СЛАУ относительно коэфициентов сплайнов c[i] методом прогонки для трехдиагональных матриц
	// Вычисление прогоночных коэфициентов - прямой ход метода прогонки
	double *alpha = new double[n - 1];
	double *beta = new double[n - 1];
	alpha[0] = beta[0] = 0.;
	for (std::size_t i = 1; i < n - 1; ++i)
	{
		double h_i = x[i] - x[i - 1], h_i1 = x[i + 1] - x[i];
		double A = h_i;
		double C = 2. * (h_i + h_i1);
		double B = h_i1;
		double F = 6. * ((y[i + 1] - y[i]) / h_i1 - (y[i] - y[i - 1]) / h_i);
		double z = (A * alpha[i - 1] + C);
		alpha[i] = -B / z;
		beta[i] = (F - A * beta[i - 1]) / z;
	}

	// Нахождение решения - обратный ход метода прогонки
	for (std::size_t i = n - 2; i > 0; --i)
		splines[i].c = alpha[i] * splines[i + 1].c + beta[i];

	// Освобождение памяти, занимаемой прогоночными коэфициентами
	delete[] beta;
	delete[] alpha;

	// По известным коэфициентам c[i] находим значения b[i] и d[i]
	for (std::size_t i = n - 1; i > 0; --i)
	{
		double h_i = x[i] - x[i - 1];
		splines[i].d = (splines[i].c - splines[i - 1].c) / h_i;
		splines[i].b = h_i * (2. * splines[i].c + splines[i - 1].c) / 6. + (y[i] - y[i - 1]) / h_i;
	}
}

double cubic_spline::f(double x) const
{
	if (!splines)
		return std::numeric_limits<double>::quiet_NaN(); // Если сплайны ещё не построены - возвращаем NaN

	spline_tuple *s;
	if (x <= splines[0].x) // Если x меньше точки сетки x[0] - пользуемся первым эл-тов массива
		s = splines + 1;
	else if (x >= splines[n - 1].x) // Если x больше точки сетки x[n - 1] - пользуемся последним эл-том массива
		s = splines + n - 1;
	else // Иначе x лежит между граничными точками сетки - производим бинарный поиск нужного эл-та массива
	{
		std::size_t i = 0, j = n - 1;
		while (i + 1 < j)
		{
			std::size_t k = i + (j - i) / 2;
			if (x <= splines[k].x)
				j = k;
			else
				i = k;
		}
		s = splines + j;
	}

	double dx = (x - s->x);
	return s->a + (s->b + (s->c / 2. + s->d * dx / 6.) * dx) * dx; // Вычисляем значение сплайна в заданной точке по схеме Горнера (в принципе, "умный" компилятор применил бы схему Горнера сам, но ведь не все так умны, как кажутся)
}

void cubic_spline::free_mem()
{
	if (splines)
	{
		delete[] splines;
		splines = NULL;
	}
}

Литература

• Роджерс Д.,Адамс Дж. Математические основы машинной графики. — М.: Мир, 2001. — ISBN 5-03-002143-4.
• Костомаров Д.П., Фаворский А.П. Вводные лекции по численным методам.