www.pudn.com > NumericalComputing.rar > splinetx.m, change:2003-12-20,size:1525b


function v = splinetx(x,y,u) 
%SPLINETX  Textbook spline function. 
%  v = splinetx(x,y,u) finds the piecewise cubic interpolatory 
%  spline S(x), with S(x(j)) = y(j), and returns v(k) = S(u(k)). 
% 
%  See SPLINE, PCHIPTX. 
 
%  First derivatives 
 
   h = diff(x); 
   delta = diff(y)./h; 
   d = splineslopes(h,delta); 
 
%  Piecewise polynomial coefficients 
 
   n = length(x); 
   c = (3*delta - 2*d(1:n-1) - d(2:n))./h; 
   b = (d(1:n-1) - 2*delta + d(2:n))./h.^2; 
 
%  Find subinterval indices k so that x(k) <= u < x(k+1) 
 
   k = ones(size(u)); 
   for j = 2:n-1 
      k(x(j) <= u) = j; 
   end 
 
%  Evaluate spline 
 
   s = u - x(k); 
   v = y(k) + s.*(d(k) + s.*(c(k) + s.*b(k))); 
 
 
% ------------------------------------------------------- 
 
function d = splineslopes(h,delta) 
%  SPLINESLOPES  Slopes for cubic spline interpolation. 
%  splineslopes(h,delta) computes d(k) = S'(x(k)). 
%  Uses not-a-knot end conditions. 
 
%  Diagonals of tridiagonal system 
 
   n = length(h)+1; 
   a = zeros(size(h)); b = a; c = a; r = a; 
   a(1:n-2) = h(2:n-1); 
   a(n-1) = h(n-2)+h(n-1); 
   b(1) = h(2); 
   b(2:n-1) = 2*(h(2:n-1)+h(1:n-2)); 
   b(n) = h(n-2); 
   c(1) = h(1)+h(2); 
   c(2:n-1) = h(1:n-2); 
 
%  Right-hand side 
 
   r(1) = ((h(1)+2*c(1))*h(2)*delta(1)+h(1)^2*delta(2))/c(1); 
   r(2:n-1) = 3*(h(2:n-1).*delta(1:n-2)+h(1:n-2).*delta(2:n-1)); 
   r(n) = (h(n-1)^2*delta(n-2)+(2*a(n-1)+h(n-1))*h(n-2)*delta(n-1))/a(n-1); 
 
%  Solve tridiagonal linear system 
 
   d = tridisolve(a,b,c,r);