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

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

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

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

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

Теорема: Для любой функции f и любого разбиения отрезка [a,b] cуществует ровно один естественный сплайн S(x), удовлетворяющий перечисленным выше условиям.

Эта теорема является следствием более общей теоремы Шёнберга-Уитни об условиях существования интерполяционного сплайна.

Построение

Обозначим: h_i = x_i x_{i 1}

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

тогда

Условия непрерывности всех производных до второго порядка включительно записываются в виде



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

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

Если учесть, что c₀ = c_n = 0, то вычисление c можно провести с помощью метода прогонки для трёхдиагональной матрицы.

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

#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;
}
}

Реализация на языке C# Платформа .NET  

// Интерполирование функций естественными кубическими сплайнами
 
using System;
 
class CubicSpline
{
    SplineTuple[] splines; // Сплайн
 
    // Структура, описывающая сплайн на каждом сегменте сетки
    struct SplineTuple
    {
        public double a, b, c, d, x;
    }
 
    // Построение сплайна
    // x - узлы сетки, должны быть упорядочены по возрастанию, кратные узлы запрещены
    // y - значения функции в узлах сетки
    // n - количество узлов сетки
    public void BuildSpline(double[] x, double[] y, int n)
    {
        // Инициализация массива сплайнов
        splines = new SplineTuple[n];
        for (int i = 0; i < n; ++i)
        {
            splines[i].x = x[i];
            splines[i].a = y[i];
        }
        splines[0].c = splines[n - 1].c = 0.0;
 
        // Решение СЛАУ относительно коэффициентов сплайнов c[i] методом прогонки для трехдиагональных матриц
        // Вычисление прогоночных коэффициентов - прямой ход метода прогонки
        double[] alpha = new double[n - 1];
        double[] beta = new double[n - 1];
        alpha[0] = beta[0] = 0.0;
        for (int 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.0 * (h_i + h_i1);
            double B = h_i1;
            double F = 6.0 * ((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 (int i = n - 2; i > 0; --i)
            splines[i].c = alpha[i] * splines[i + 1].c + beta[i];
 
        // Освобождение памяти, занимаемой прогоночными коэффициентами
        beta = null;
        alpha = null;
 
        // По известным коэффициентам c[i] находим значения b[i] и d[i]
        for (int 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.0 * splines[i].c + splines[i - 1].c) / 6.0 + (y[i] - y[i - 1]) / h_i;
        }
    }
 
    // Вычисление значения интерполированной функции в произвольной точке
    double Func(double x)
    {
        if (splines == null)
            return double.NaN; // Если сплайны ещё не построены - возвращаем NaN
 
        int n = splines.Length;
        SplineTuple 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 лежит между граничными точками сетки - производим бинарный поиск нужного эл-та массива
        {
            int i = 0, j = n - 1;
            while (i + 1 < j)
            {
                int 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.0 + s.d * dx / 6.0) * dx) * dx;
    }
}

Литература

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

п·о·р Кривые
Плоские
Алгебраические
Конические сечения (2-й порядок)	Гипербола • Парабола • Эллипс (Окружность)
3-й порядок	Эллиптические: Эллиптическая кривая • Функции Якоби • Интеграл • Функции Другие: Верзьера Аньези • Декартов лист • Полукубическая парабола • Строфоида • Циссоида Диокла
Лемнискаты (2n порядок)	Бернулли (Овал Кассини) • Бута • Жероно
Аппроксимационные	Сплайны (B-сплайн · Кубический сплайн · Моносплайн · Сплайн Эрмита) • Кривая Безье
Трансцендентные
Спирали	Архимедова (Ферма) • Гиперболическая • «Жезл» • Клотоида • Логарифмическая
Циклоидальные (порождённые катящейся окружностью)	Циклоида • Эпициклоида (Кардиоида · Нефроида) • Гипоциклоида (Дельтоида (кривая Штейнера) · Астроида) • Трохоида (Удлинённая циклоида · Укороченная циклоида) • Эпитрохоида (Удлинённая эпициклоида · Укороченная эпициклоида · Улитка Паскаля · «Роза») • Гипотрохоида • Скорейшего спуска (Брахистохрона, дуга циклоиды)
Фрактальные	Кривая Коха • Кривая Леви • Кривая Минковского • Кривая Пеано • Топологические: салфетка и ковёр Серпинского, губка Менгера
Другие	Квадратриса • Кривая погони (Трактриса) • Трохоида • Цепная линия (перевёрнутая арочная) • Кривая Безье • Кривая постоянной ширины • Синусоида
Неплоские
Винтовая линия • Линия откоса • Локсодрома • Ортодромия • Губка
Связанные понятия
Определения кривой	Аналитическая кривая • Кривая Жордана • Канторова кривая • Кривая Урысона
Преобразованные кривые	Эволюта • Эвольвента • Каустика