Кубический сплайн: различия между версиями

Интерактивная навигация по истории

(Показать все непатрулированные изменения)

[непроверенная версия]

← Предыдущая правка Следующая правка →

Содержимое удалено Содержимое добавлено

ВизуальныйВики-текст

Линейный

Версия от 12:17, 18 февраля 2018

Некоторая функция f(x) задана на отрезке $[a,b]$ , разбитом на части $[x_{i-1},x_{i}]$ , $a=x_{0}<x_{1}<...<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}$ .

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

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

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

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

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

Построение

На каждом отрезке $[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} \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}\right)=f(x_{i})

Обозначим $:\quad h_{i}=x_{i}-x_{i-1},\quad f_{i}=f(x_{i})$

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

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}}}={{f_{i}-f_{i-1}} \over {h_{i}}}+{{h_{i}(2c_{i}+c_{i-1})} \over 6}

Если учесть, что $c_{0}=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 = 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 (int 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 (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. * 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()
{
	delete[] splines;
	splines = NULL;
}

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

// Интерполирование функций естественными кубическими сплайнами

using System;

class CubicSpline
{
    SplineTuple[] splines; // Сплайн

    // Структура, описывающая сплайн на каждом сегменте сетки
    private 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 hi  = x[i] - x[i - 1];
            double hi1 = x[i + 1] - x[i];
            double A = hi;
            double C = 2.0 * (hi + hi1);
            double B = hi1;
            double F = 6.0 * ((y[i + 1] - y[i]) / hi1 - (y[i] - y[i - 1]) / hi);
            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];
        }

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

    // Вычисление значения интерполированной функции в произвольной точке
    public double Interpolate(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[0];
        }
        else if (x >= splines[n - 1].x) // Если x больше точки сетки x[n - 1] - пользуемся последним эл-том массива
        {
            s = splines[n - 1];
        }
        else // Иначе x лежит между граничными точками сетки - производим бинарный поиск нужного эл-та массива
        {
            int i = 0;
            int 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; 
    }
}

Реализация на языке PHP

class CubicSpline
{
	private static $_splines = []; // Сплайн

	// Структура, описывающая сплайн на каждом сегменте сетки
	private static $_spline_struct = [
				'a' => 0,
				'b' => 0,
				'c' => 0,
				'd' => 0,
				'x' => 0,
			];


	/**
	*   Построение сплайна
    *
	* @param array $x - узлы сетки, должны быть упорядочены по возрастанию, кратные узлы запрещены
	* @param array $y - значения функции в узлах сетки
	* @param integer $n - количество узлов сетки
	*
	* @return array возвращаем страктуру сплайна
	*/
	public static function buildSpline($x, $y, $n)
	{
      	// Инициализация массива сплайнов
        for ($i = 0; $i < $n; ++$i) {
            self::$_splines[$i]['x'] = $x[$i];
            self::$_splines[$i]['a'] = $y[$i];
        }

        self::$_splines[0]['c'] = self::$_splines[$n - 1]['c'] = 0;

        // Решение СЛАУ относительно коэффициентов сплайнов c[i] методом прогонки для трехдиагональных матриц
        // Вычисление прогоночных коэффициентов - прямой ход метода прогонки
        $alpha = [];
        $beta  = [];

        $alpha[0] = $beta[0] = 0;

        for ($i = 1; $i < ($n - 1); ++$i) {
            $hi  = $x[$i] - $x[$i - 1];
            $hi1 = $x[$i + 1] - $x[$i];
            $A = $hi;
            $C = 2.0 * ($hi + $hi1);
            $B = $hi1;
            $F = 6.0 * (($y[$i + 1] - $y[$i]) / $hi1 - ($y[$i] - $y[$i - 1]) / $hi);
            $z = ($A * $alpha[$i - 1] + $C);
            $alpha[$i] = -$B / $z;
            $beta[$i] = ($F - $A * $beta[$i - 1]) / $z;
        }

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

        // По известным коэффициентам c[i] находим значения b[i] и d[i]
        for ($i = ($n - 1); $i > 0; --$i) {
            $hi = $x[$i] - $x[$i - 1];
            self::$_splines[$i]['d'] = (self::$_splines[$i]['c']- self::$_splines[$i - 1]['c']) / $hi;
            self::$_splines[$i]['b'] = $hi * (2.0 * self::$_splines[$i]['c'] + self::$_splines[$i - 1]['c']) / 6.0 + ($y[$i] - $y[$i - 1]) / $hi;
        }

        return self::$_splines;
	}

	// Вычисление значения интерполированной функции в произвольной точке
    public static function interpolate($x)
    {
        if (self::$_splines == null) {
            return null; // Если сплайны ещё не построены - возвращаем NULL
        }

        $n = count(self::$_splines);
        $s = self::$_spline_struct;

        if ($x <= self::$_splines[0]['x']) { // Если x меньше точки сетки x[0] - пользуемся первым эл-тов массива

            $s = self::$_splines[0];

        } elseif ($x >= self::$_splines[$n - 1]['x']) { // Если x больше точки сетки x[n - 1] - пользуемся последним эл-том массива

            $s = self::$_splines[$n - 1];

        } else { // Иначе x лежит между граничными точками сетки - производим бинарный поиск нужного эл-та массива

            $i = 0;
            $j = $n - 1;
            while (($i + 1) < $j) {
                $k = $i + ($j - $i) / 2;
                if ($x <= self::$_splines[$k]['x']) {
                    $j = $k;
                } else {
                    $i = $k;
                }
            }
            $s = self::$_splines[$j];
        }

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

Python 2.7, шаг постоянный

from collections import defaultdict
from matplotlib import mlab
import bisect, pylab, math

def drange(start, stop, step):
    while start < stop: yield start; start += step;

class Dot:
    def __init__(self, x, y): self.x, self.y = [x, y]

class Tuple: a, b, c, d, x = [0., 0., 0., 0., 0.]

def buildSpline(dots):
    for i in range(len(dots)): splines[i].x, splines[i].a = dots[i].x, dots[i].y

    alpha, beta = [defaultdict(lambda: 0.), defaultdict(lambda: 0.)]

    for i in range(1, len(dots)-1):
        C = 4. * in_step
        F = 6. * ((dots[i + 1].y - dots[i].y) / in_step - (dots[i].y - dots[i - 1].y) / in_step)
        z = (in_step* alpha[i - 1] + C)
        alpha[i] = -in_step / z
        beta[i] = (F - in_step* beta[i - 1]) / z

    for i in reversed(range(1, len(dots) - 1)): splines[i].c = alpha[i] * splines[i+1].c + beta[i]

    for i in reversed(range(1, len(dots))):
        hi = dots[i].x - dots[i-1].x
        splines[i].d = (splines[i].c - splines[i-1].c) / hi
        splines[i].b = hi * (2.0 * splines[i].c + splines[i - 1].c) / 6.0 + (dots[i].y - dots[i-1].y) / hi

def calc(x):
    distribution = sorted([t[1].x for t in splines.items()])
    indx = bisect.bisect_left(distribution, x)
    if indx == len(distribution): return 0
    dx = x - splines[indx].x
    return splines[indx].a + splines[indx].b * dx + splines[indx].c * dx**2 / 2. + splines[indx].d * dx**3 / 6.
#============================================

in_func =  lambda x:  math.sin(x)
in_min_x = 0
in_max_x = 25 
in_step = 2.5

out_min_x = 0
out_max_x = 25
out_step = 0.1

#============================================

#build model
splines = defaultdict(lambda: Tuple())
buildSpline(map(lambda x: Dot(x, in_func(x)), [x for x in drange(in_min_x, in_max_x+1, in_step)]))

#print result
for x in drange(out_min_x, out_max_x, out_step):
    print str(x) + ';' + str(calc(x))

#build graphics
xlist = mlab.frange (out_min_x, out_max_x, out_step)
pylab.plot(xlist, [calc(x) for x in xlist])
pylab.plot(xlist, [in_func(x) for x in xlist])
pylab.show()

Литература

Роджерс Д., Адамс Дж. Математические основы машинной графики. — М.: Мир, 2001. — ISBN 5-03-002143-4.
Костомаров Д. П., Фаворский А. П. Вводные лекции по численным методам.
Волков Е. А. Глава 1. Приближение функций многочленами. § 11. Сплайны // Численные методы. — Учеб. пособие для вузов. — 2-е изд., испр.. — М.: Наука, 1987. — С. 63-68. — 248 с.

Ссылки

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

Кубический сплайн: различия между версиями

Версия от 12:17, 18 февраля 2018

Построение

Литература

Ссылки

Навигация

Поиск