This preview has intentionally blurred sections. Sign up to view the full version.
View Full 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&C p.319) % h = x(2:n)  x(1:n1); u = 2*( h(1:n2) + h(2:n1) ); b = 6*( y(2:n)  y(1:n1) )./h; if n>3 T = spdiags([ h(2:n1) u h(1:n2)],1:1,n2,n2); else T = u; end v = b(2:n1)  b(1:n2); % % 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&C Eqn 11) % A = ( z(2:n)  z(1:n1) )./(6*h); B = z(1:n1)/2; C = b/6  ( z(2:n)/6 + z(1:n1)/3 ).*h; pp = mkpp(x',[ A B C y(1:n1) ]); end if nargin==2, output = pp; else output = ppval(pp,xx); end end...
View
Full Document
 Spring '09
 Spline interpolation, cubic spline interpolation, cubic spline, Tampere University of Technology

Click to edit the document details