function output = cspline(x,y,xx) % CSPLINE Cubic spline interpolation % % cubic spline interpolation for given nodes x % Two types of spline interpolation are available: % With given slopes at the endpoints: % y = [left_slope y_values right_slope] : given values % Without given slopes at the endpoints: % y = [y_values] : given values % xx - where to interpolate % % usage: [a, b, c, d] = cspline(x,y,xx) or yy = cspline(x,y,xx) % % See also the matlab-function SPLINE % % Example: % x = 0:10; y = cos(x); % xx = 0:.25:10; % yy = cspline(x,y,xx); % plot(x,y,'o',xx,yy) % % Knowing the slopes at the endpoints, % a better result can be obtained by % % x = 0:10; y = cos(x); % xx = 0:.25:10; % yy = cspline(x,[-sin(xx(1)) y -sin(xx(end))],xx); % plot(x,y,'o',xx,yy) % (c) Rolf Krause, 1998, krause@math.fu-berlin.de if nargin~=2 & nargin~=3, error('wrong number of arguments'); end % we need the stepsize output = []; n = length(x); if n ~= length(y)-2 & n ~= length(y) error(['y has to be of length length(x) + 2 or length(x)']); end if n < 2, error('only one value given, can not interpolate'); end % check for the slopes at the endpoints being given or not [nr, nc] = size(y); if nr == 1, y = reshape(y, nc, 1); nr = nc; end [nr, nc] = size(x); if nr == 1, x = reshape(x, nc, 1); nr = nc; end if(length(y) == length(x)) naturalInterpolation = 1; dy_l = 0; dy_r = 0; else naturalInterpolation = 0; % y consists of the slopes at the endpoints and of the values of y dy_l = y(1); dy_r = y(n+2); y = y(2:n+1); end if size(x) ~= size(y), error('x and y are of different size'); end dx = [0; diff(x); 0]; dxx = dx(1:n) + dx(2:n+1); % d_xx_j = h_j+h_{j+1} % assemble matrix and rhs M = spdiags([[dx(2:n)./dxx(2:n); 0] 2*ones(n,1) [0; dx(2:n)./dxx(1:n-1)]], -1:1, n,n); % compute the rhs using aitken-neville scheme % c : second derivative % a = y: values of y % b : first derivative % d : third derivative b = diff(y) ./ dx(2:n); c = 6 * diff([dy_l; b; dy_r])./ dxx; %% For natural spline interpolation if(naturalInterpolation == 1) c(1) = 0; c(n) = 0; M(1,2) = 0; M(n,n-1) = 0; end c = M\c; d = diff(c)./dx(2:n); b = b - dx(2:n).* (c(1:n-1)/3 + c(2:n)/6); if nargin == 2 output = [y,b,c,d] return; end % now compute the values yy yy = zeros(size(xx)); for i=1:nr-1 I = find(xx <= x(i+1) & xx >= x(i)); yy(I) = y(i) + b(i)*(xx(I)-x(i))+c(i)/2*(xx(I)-x(i)).^2 + ... d(i)/6*(xx(I)-x(i)).^3; end output = yy;