% ### EXsplines1.m ### 10.09.14 % Determine cublic spline fit to a set of (specified) data. Does this three % ways: % Method 0 - via built-in Matlab spline.m function % Method 1 - via built-in Matlab interp1.m function % Method 2 - directly coded % Method 3 - polynomial of order n-1 (where n=numel(x)) % source code for doing cubic spline (sans Matlab spline.m) is % http://m2matlabdb.ma.tum.de/cspline.m?MP_ID=14 clear; figure(1); clf; % ----------------------------- % User Inputs N= 70; % # of points for fit (over interval [min(x) max(x)]) % define 'data' (see pg.69 from Kuts, 2013) x= [0 0.5 1.1 1.7 2.1 2.5 2.9 3.3 3.7 4.2 4.9 5.3 6.0 6.7 7.0]; y= [1.1 1.6 2.4 3.8 4.3 4.7 4.8 5.5 6.1 6.3 7.1 7.1 8.2 6.9 5.3]; % ------------------------------------------------- xx= linspace(min(x),max(x),N); % determine fit x-values % ------------ % Method 0: Built-in Matlab function spline.m yspline0= spline(x,y,xx); % ------------ % Method 1: Built-in Matlab function interp1.m yspline1= interp1(x,y,xx,'spline'); % Note: Other options for interp1 exists, such as: % > yspline1= interp1(x,y,xx); % straight line nearest neighbor interp (default for plot.m) % > yspline1= interp1(x,y,xx,'nearest'); % % ------------ % Method 2: Direct computation of splines (see comments at top for source) n = numel(x); % # of points to be fit 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); % 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); % now compute the values yy (i.e., the fit) yspline2 = zeros(size(xx)); for i=1:nr-1 I = find(xx <= x(i+1) & xx >= x(i)); yspline2(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 % ------------ % Method 3: polynomial fit via built-in Matlab function polyfit.m coeff= polyfit(x,y,numel(x)-1); ypoly= polyval(coeff,xx); % ------------ % plot the data figure(1); clf hold on; grid on; xlabel('x');ylabel('y'); plot(xx,yspline0,'ks','MarkerSize',6) plot(xx,yspline1,'g.') plot(xx,yspline2); plot(xx,ypoly,'m--','LineWidth',2) plot(x,y,'ro','MarkerSize',7,'LineWidth',2); axis([-0.2 7.2 0.5 8.5]) legend('spline.m','interp1.m','direct calculation','polynomial','original data','Location','South')