Я написал пакет с кубическим сплайном в Mathematica очень давно.Вот мой перевод этого пакета на Matlab.Заметьте, я не смотрел кубические сплайны около 7 лет, поэтому я основываюсь на своей собственной документации.Вы должны проверить все, что я говорю.
Основная проблема в том, что нам даны n точки данных (x(1), y(1)) , ... , (x(n), y(n)), и мы хотим вычислить кусочно-кубический интерполант.Интерполант определяется как
S(x) = { Sk(x) when x(k) <= x <= x(k+1)
{ 0 otherwise
Здесь Sk (x) - кубический многочлен вида
Sk(x) = sk0 + sk1*(x-x(k)) + sk2*(x-x(k))^2 + sk3*(x-x(k))^3
Свойства сплайна:
- Сплайн проходит через точку данных
Sk(x(k)) = y(k) - Сплайн является непрерывным в конечных точках и, следовательно, непрерывным повсюду в интервале интерполяции
Sk(x(k+1)) = Sk+1(x(k+1)) - Сплайн имеет непрерывную первую производную
Sk'(x(k+1)) = Sk+1'(x(k+1)) - Сплайн имеет непрерывную вторую производную
Sk''(x(k+1)) = Sk+1''(x(k+1))
Чтобы построить кубический сплайн из набора точек данных, нам нужно найти коэффициенты sk0, sk1, sk2 и sk3 для каждого из n-1 кубических полиномов.Это всего 4*(n-1) = 4*n - 4 неизвестных.Свойство 1 содержит ограничения n, а свойства 2,3,4 - дополнительные ограничения n-2.Таким образом, у нас есть n + 3*(n-2) = 4*n - 6 ограничений и 4*n - 4 неизвестных.Это оставляет две степени свободы.Мы фиксируем эти степени свободы, устанавливая в качестве начальной и конечной вершины вторую производную равной нулю.
Let m(k) = Sk''(x(k)), h(k) = x(k+1) - x(k) и d(k) = (y(k+1) - y(k))/h(k).Имеет место следующее трехчленное рекуррентное соотношение
h(k-1)*m(k-1) + 2*(h(k-1) + h(k))*m(k) + h(k)*m(k+1) = 6*(d(k) - d(k-1))
m (k) - неизвестные, для которых мы хотим найти решение.h(k) и d(k) определяются входными данными.Это трехчленное рекуррентное соотношение определяет трехдиагональную линейную систему.Как только m(k) определены, коэффициенты для Sk задаются как
sk0 = y(k)
sk1 = d(k) - h(k)*(2*m(k) + m(k-1))/6
sk2 = m(k)/2
sk3 = m(k+1) - m(k)/(6*h(k))
Хорошо, это все, что вам нужно знать, чтобы полностью определить алгоритм для вычисления кубического сплайна.Вот это в Matlab:
function [s0,s1,s2,s3]=cubic_spline(x,y)
if any(size(x) ~= size(y)) || size(x,2) ~= 1
error('inputs x and y must be column vectors of equal length');
end
n = length(x)
h = x(2:n) - x(1:n-1);
d = (y(2:n) - y(1:n-1))./h;
lower = h(1:end-1);
main = 2*(h(1:end-1) + h(2:end));
upper = h(2:end);
T = spdiags([lower main upper], [-1 0 1], n-2, n-2);
rhs = 6*(d(2:end)-d(1:end-1));
m = T\rhs;
% Use natural boundary conditions where second derivative
% is zero at the endpoints
m = [ 0; m; 0];
s0 = y;
s1 = d - h.*(2*m(1:end-1) + m(2:end))/6;
s2 = m/2;
s3 =(m(2:end)-m(1:end-1))./(6*h);
Вот некоторый код для построения кубического сплайна:
function plot_cubic_spline(x,s0,s1,s2,s3)
n = length(x);
inner_points = 20;
for i=1:n-1
xx = linspace(x(i),x(i+1),inner_points);
xi = repmat(x(i),1,inner_points);
yy = s0(i) + s1(i)*(xx-xi) + ...
s2(i)*(xx-xi).^2 + s3(i)*(xx - xi).^3;
plot(xx,yy,'b')
plot(x(i),0,'r');
end
Вот функция, которая создает кубический сплайн и строит на знаменитом Rungeфункция:
function cubic_driver(num_points)
runge = @(x) 1./(1+ 25*x.^2);
x = linspace(-1,1,num_points);
y = runge(x);
[s0,s1,s2,s3] = cubic_spline(x',y');
plot_points = 1000;
xx = linspace(-1,1,plot_points);
yy = runge(xx);
plot(xx,yy,'g');
hold on;
plot_cubic_spline(x,s0,s1,s2,s3);
Вы можете увидеть это в действии, выполнив следующую команду Matlab
>> cubic_driver(5)
>> clf
>> cubic_driver(10)
>> clf
>> cubic_driver(20)
К тому времени, когда у вас будет двадцать узлов, ваш интерполант визуально неотличим от функции Runge,
Некоторые комментарии к коду Matlab: я не использую циклы for или while.Я могу векторизовать все операции.Я быстро формирую разреженную трехдиагональную матрицу с spdiags.Я решаю это с помощью оператора обратной косой черты.Я рассчитываю на UMFPACK Тима Дэвиса, чтобы справиться с разложением, а также решить задачи вперед и назад.
Надеюсь, это поможет.Код доступен в виде gist на github https://gist.github.com/1269709