Кубический сплайн: различия между версиями
[непроверенная версия] | [отпатрулированная версия] |
отмена правки 93244503 участника 176.99.198.6 (обс.) Метка: отмена |
Oleg4280 (обсуждение | вклад) 6 правок возвращено к версии 86694403 OneLittleMouse: убрал программную реализацию |
||
Строка 40: | Строка 40: | ||
Если учесть, что <math>c_0 = c_n = 0</math>, то вычисление <math>c</math> можно провести с помощью [[Метод прогонки|метода прогонки]] для [[трёхдиагональная матрица|трёхдиагональной матрицы]]. |
Если учесть, что <math>c_0 = c_n = 0</math>, то вычисление <math>c</math> можно провести с помощью [[Метод прогонки|метода прогонки]] для [[трёхдиагональная матрица|трёхдиагональной матрицы]]. |
||
Реализация на языке [[C++]] |
|||
<source lang="cpp"> |
|||
#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; |
|||
} |
|||
</source> |
|||
Реализация на языке [[C Sharp|C#]] Платформа [[.NET Framework|.NET]] |
|||
<source lang="csharp"> |
|||
// Интерполирование функций естественными кубическими сплайнами |
|||
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; |
|||
} |
|||
} |
|||
</source> |
|||
Реализация на языке [[PHP]]<syntaxhighlight lang="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; |
|||
} |
|||
} |
|||
</syntaxhighlight> |
|||
Python 2.7, шаг постоянный |
|||
<source lang="python"> |
|||
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() |
|||
</source> |
|||
== Литература == |
== Литература == |
Версия от 19:04, 11 июня 2018
Некоторая функция f(x) задана на отрезке , разбитом на части , . Кубическим сплайном дефекта 1 называется функция , которая:
- на каждом отрезке является многочленом степени не выше третьей;
- имеет непрерывные первую и вторую производные на всём отрезке ;
- в точках выполняется равенство , т. е. сплайн интерполирует функцию f в точках .
Для однозначного задания сплайна перечисленных условий недостаточно, для построения сплайна необходимо наложить какие-то дополнительные требования.
Естественным кубическим сплайном называется кубический сплайн, удовлетворяющий также граничным условиям вида:
Теорема: Для любой функции и любого разбиения отрезка существует ровно один естественный сплайн S(x), удовлетворяющий перечисленным выше условиям.
Эта теорема является следствием более общей теоремы Шёнберга-Уитни об условиях существования интерполяционного сплайна.
Построение
На каждом отрезке функция есть полином третьей степени , коэффициенты которого надо определить. Запишем для удобства в виде:
тогда
Условия непрерывности всех производных до второго порядка включительно
записываются в виде
а условия интерполяции в виде
Обозначим
Отсюда получаем формулы для вычисления коэффициентов сплайна:
Если учесть, что , то вычисление можно провести с помощью метода прогонки для трёхдиагональной матрицы.
Литература
- Роджерс Д., Адамс Дж. Математические основы машинной графики. — М.: Мир, 2001. — ISBN 5-03-002143-4.
- Костомаров Д. П., Фаворский А. П. Вводные лекции по численным методам.
- Волков Е. А. Глава 1. Приближение функций многочленами. § 11. Сплайны // Численные методы. — Учеб. пособие для вузов. — 2-е изд., испр.. — М.: Наука, 1987. — С. 63-68. — 248 с.
Ссылки
Для улучшения этой статьи желательно:
|