|
Если учесть, что <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 (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()
{
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>
== Литература ==
| 