Created
October 7, 2011 07:40
-
-
Save cmaes/1269709 to your computer and use it in GitHub Desktop.
Small cubic spline implementation in Matlab
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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); | |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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); |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
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)
-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