Кубический сплайн | это... Что такое Кубический сплайн? (original) (raw)

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

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

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

S''(a) = S''(b) = 0.

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

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

Построение

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

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

$S_i(x) = a_i + b_i(x - x_i) + {c_i\over2}(x-x_i)^2 + {d_i\over6}(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_ic_{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\over2}h_ic_i - {1\over6}h_i^2d_i + {{f_i - f_{i-1}}\over{h_i}} \,\!$

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

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

#include #include #include

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 = 0.;

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

    splines[n - 1].c = (F - A * beta[n - 2]) / (C + A * alpha[n - 2]);

    // Нахождение решения - обратный ход метода прогонки
    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::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() { 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;
    }
}

// Вычисление значения интерполированной функции в произвольной точке
public 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
Костомаров Д. П., Фаворский А. П. Вводные лекции по численным методам.
Волков Е. А. Глава 1. Приближение функций многочленами. § 11. Сплайны // Численные методы. — Учеб. пособие для вузов. — 2-е изд., испр.. — М.: Наука, 1987. — С. 63-68. — 248 с.

Ссылки

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

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