Skip to content

Instantly share code, notes, and snippets.

@cmaes
Created October 7, 2011 07:40
Show Gist options
  • Save cmaes/1269709 to your computer and use it in GitHub Desktop.
Save cmaes/1269709 to your computer and use it in GitHub Desktop.
Small cubic spline implementation in Matlab
function cubic_driver(num_points)
% cubic_driver(n) computes a cubic spline
% interpolant of the runge function
%
% f(x) = 1/(1+25*x^2)
%
% based on n linearly spaced
% points in the interval [-1,1]
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);
function [s0,s1,s2,s3]=cubic_spline(x,y)
% [s0,s1,s2,s3]=cubic_spline(x,y)
%
% computes the coefficents of a cubic spline
% interpolant through the data points (x,y)
%
% The spline is defined as the piecewise cubic
% polynomial
%
% S(x) = { Sk(x) x(k) <= x <= x(k+1)
% = { 0 otherwise
%
% The cubic polynomial Sk(x) is given by
%
% Sk(x) = sk0 + sk1*(x-x(k)) + sk2*(x-x(k))^2 + sk3*(x-x(k))^3
%
% The coefficents sk0, sk1, sk2, and sk3 for each of the
% polynomials are returned in the vectors s0,s1,s2, and s3
% respectively.
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)
% plot_cubic_spline(x,s0,s1,s2,s3)
%
% plots a cubic spline with break points x
% and coefficents 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
Display the source blob
Display the rendered blob
Raw
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
@cmaes
Copy link
Author

cmaes commented Oct 7, 2011

A tiny Matlab implementation of cubic spline interpolation, based on work done for the 18.310 class at MIT. This gist was motivated by the answer to this question on stack overflow.

@ouboub
Copy link

ouboub commented Feb 18, 2017

Hi
thanks for the code, however when i run it with an example of mine it does not coincide. I found
https://stackoverflow.com/questions/7642921/cubic-spline-program

lower = h(2:end);
main = 2*(h(1:end-1) + h(2:end));
upper = h(1:end-1);

which gives the correct result. Could you please correct here also?

Your code and my calculations do not agree with matlab. here is a simple example
x=[-3 -2 1 4 6]
y=[2 0 3 1 4]
[s0,s1,s2,s3]=cubic_spline(x',y')
comment the length of the vectors is not consistent
s3 s2 s1 s0
0.5044 0 -2.5044 2.0000
-0.2831 1.5132 -0.9912 0
0.2217 -1.0351 0.4430 3.0000
-0.1601 0.9605 0.2193 1.0000

however matlab gives using

clear all
x=[-3 -2 1 4 6]
y=[2 0 3 1 4]
pspl=spline(x,[0 y 0])
[breaks,coefs,l,k,d] = unmkpp(pspl)

 s3           s2            s1       s0
2.0595   -4.0595         0    2.0000

-0.3796 2.1190 -1.9405 0
0.3003 -1.2976 0.5238 3.0000
-0.5387 1.4048 0.8452 1.0000

any ideas why this is so?

thanks

Uwe Brauer

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment