home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / POLYFUN.DI$ / SPLINE.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  2.5 KB  |  71 lines
  1. function output = spline(x,y,xx)
  2. %SPLINE    Cubic spline data interpolation.
  3. %    Given data vectors X and Y, and a new abscissa vector XI, the
  4. %    function YI = SPLINE(X,Y,XI) uses cubic spline interpolation
  5. %    to find a vector YI corresponding to XI.
  6. %
  7. %    Here's an example that generates a coarse sine curve, then
  8. %    interpolates over a finer abscissa:
  9. %
  10. %        x = 0:10;  y = sin(x);
  11. %        xi = 0:.25:10;
  12. %        yi = spline(x,y,xi);
  13. %        plot(x,y,'o',xi,yi)
  14. %
  15. %    PP = spline(x,y) returns the pp-form of the cubic spline interpolant
  16. %    instead, for later use with  ppval, etc.
  17. %
  18. %    See also INTERP1, INTERP2, PPVAL, MKPP, UNMKPP, the Spline Toolbox.
  19.  
  20. %    Carl de Boor 7-2-86
  21. %    Revised 11-24-87 JNL, 6-16-92 CBM.
  22. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  23.  
  24. % Generate the cubic spline interpolant in pp form, depending on
  25. % the number of data points.  (using the not-a-knot end condition).
  26.  
  27. n=length(x);[xi,ind]=sort(x);xi=xi(:);
  28. output=[];
  29. if n<2,
  30.    fprintf('There should be at least two data points!\n')
  31. elseif all(diff(xi))==0,
  32.    fprintf('The data abscissae should be distinct!\n')
  33. elseif n~=length(y),
  34.    fprintf('Abscissa and ordinate vector should be of the same length!\n')
  35. else   
  36.    yi=y(ind);yi=yi(:);
  37.    if (n==2), % the interpolant is a straight line
  38.       pp=mkpp(xi',[diff(yi)./diff(xi) yi(1)]);
  39.    elseif (n==3), % the interpolant is a parabola
  40.       yi(2:3)=diff(yi)./diff(xi);
  41.       yi(3)=(yi(3)-yi(2))/(xi(3)-xi(1));
  42.       yi(2)=yi(2)-yi(3)*(xi(2)-xi(1));
  43.       pp = mkpp([xi(1),xi(3)],yi([3 2 1])');
  44.    else % set up the sparse, tridiagonal, linear system for the slopes at  xi .
  45.       dx=diff(xi);divdif=diff(yi)./dx;xi31=xi(3)-xi(1);xin=xi(n)-xi(n-2);
  46.       c = spdiags([ [dx(2:n-1);xin;0] ...
  47.            [dx(2);2*[dx(2:n-1)+dx(1:n-2)];dx(n-2)] ...
  48.            [0;xi31;dx(1:n-2)] ],[-1 0 1],n,n);
  49.       b=zeros(n,1);
  50.       b(1)=((dx(1)+2*xi31)*dx(2)*divdif(1)+dx(1)^2*divdif(2))/xi31;
  51.       b(2:n-1)=3*(dx(2:n-1).*divdif(1:n-2)+dx(1:n-2).*divdif(2:n-1));
  52.       b(n)=...
  53.       (dx(n-1)^2*divdif(n-2)+(2*xin+dx(n-1))*dx(n-2)*divdif(n-1))/xin;
  54.  
  55.         % sparse linear equation solution for the slopes
  56.       mmdflag = spparms('autommd');
  57.       spparms('autommd',0);
  58.       s=c\b;
  59.       spparms('autommd',mmdflag);
  60.         % convert to pp form
  61.       c4=(s(1:n-1)+s(2:n)-2*divdif(1:n-1))./dx;
  62.       c3=(divdif(1:n-1)-s(1:n-1))./dx - c4;
  63.       pp=mkpp(xi',[c4./dx c3 s(1:n-1) yi(1:n-1)]);
  64.    end
  65.    if nargin==2,
  66.       output=pp;
  67.    else
  68.       output=ppval(pp,xx);
  69.    end
  70. end
  71.