nspline - elseif n~=length(y) error('Abscissa and ordinate...

Info iconThis preview shows pages 1–2. Sign up to view the full content.

View Full Document Right Arrow Icon
function output = nspline(x,y,xx) %NSPLINE Natural cubic spline interpolation. % YY = NSPLINE(X,Y,XX) interpolates the data in vector Y at the % abscissas in vector X using a cubic spline with natural % end conditions, and evaluates it at XX. % PP = NSPLINE(X,Y) returns the pp-form of the cubic spline % interpolant for later use with PPVAL or PPDER. % See also SPLINE, PSPLINE, LSPLINE. % Author: Robert Piche, Tampere University of Technology, Finland % (with code from MatLab toolbox routine SPLINE) % Last updated: Jan. 31, 1994 (requires MatLab 4) if nargin==0, help nspline, return, end x = x(:); y = y(:); n = length(x); [x,ind] = sort(x); if n<2 error('There should be at least two data points.') elseif all(diff(x))==0 error('The data abscissas should be distinct.')
Background image of page 1

Info iconThis preview has intentionally blurred sections. Sign up to view the full version.

View Full DocumentRight Arrow Icon
Background image of page 2
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: elseif n~=length(y) error('Abscissa and ordinate vectors should be the same length.') else y = y(ind); if (n==2), % the interpolant is a straight line pp=mkpp(x',[diff(y)./diff(x) y(1)]); else % % Set up the linear system T*z=v for curvatures (K&amp;C p.319) % h = x(2:n) - x(1:n-1); u = 2*( h(1:n-2) + h(2:n-1) ); b = 6*( y(2:n) - y(1:n-1) )./h; if n&gt;3 T = spdiags([ h(2:n-1) u h(1:n-2)],-1:1,n-2,n-2); else T = u; end v = b(2:n-1) - b(1:n-2); % % Solve for curvatures at the internal knots % z = T\v; % % Add the curvatures at the ends % z = [0; z; 0]; % % Compute coefficients of pp form (K&amp;C Eqn 11) % A = ( z(2:n) - z(1:n-1) )./(6*h); B = z(1:n-1)/2; C = b/6 - ( z(2:n)/6 + z(1:n-1)/3 ).*h; pp = mkpp(x',[ A B C y(1:n-1) ]); end if nargin==2, output = pp; else output = ppval(pp,xx); end end...
View Full Document

Page1 / 2

nspline - elseif n~=length(y) error('Abscissa and ordinate...

This preview shows document pages 1 - 2. Sign up to view the full document.

View Full Document Right Arrow Icon
Ask a homework question - tutors are online