home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 5.ddi / SPLINES.DI$ / CSAPI.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.1 KB  |  61 lines
  1. function output = csapi(x,y,xx)
  2. % CSAPI    Cubic spline interpolation.
  3. %
  4. %      values = csapi(x,y,xx)
  5. %
  6. % returns the values at  xx  of the cubic spline interpolant to the
  7. % given data  (x,y)  (using the not-a-knot end condition).
  8. %      The alternative call
  9. %
  10. %      pp = csapi(x,y)
  11. %
  12. % returns the pp-form of the cubic spline interpolant instead, for later
  13. % use with  fnval,  fnder, etc.
  14.  
  15. % C. de Boor / latest change: Feb.25, 1989
  16. % Copyright (c) 1990-92 by Carl de Boor and The MathWorks, Inc.
  17.  
  18. %     Generate the cubic spline interpolant in pp form, depending on
  19. % the number of data points.
  20.  
  21. n=length(x);[xi,ind]=sort(x);xi=xi(:);
  22. if n<2,
  23.    error('There should be at least two data points!')
  24. elseif all(diff(xi))==0,
  25.    error('The data abscissae should be distinct!')
  26. elseif n~=length(y),
  27.    error('Abscissa and ordinate vector should be of the same length!')
  28. else   
  29.    yi=y(ind);yi=yi(:);
  30.    if (n==2), % the interpolant is a straight line
  31.       pp=ppmak(xi',[diff(yi)./diff(xi) yi(1)],1);
  32.    elseif (n==3), % the interpolant is a parabola
  33.       yi(2:3)=diff(yi)./diff(xi);
  34.       yi(3)=(yi(3)-yi(2))/(xi(3)-xi(1));
  35.       yi(2)=yi(2)-yi(3)*(xi(2)-xi(1));
  36.       pp = ppmak([xi(1),xi(3)],yi([3 2 1])',1);
  37.    else % set up the linear system for solving for the slopes at  xi .
  38.       dx=diff(xi);divdif=diff(yi)./dx;xi31=xi(3)-xi(1);xin=xi(n)-xi(n-2);
  39.       c = diag([dx(2:n-1);0],-1)+2*diag([0;dx(2:n-1)+dx(1:n-2);0]) ...
  40.                     + diag([0;dx(1:n-2)],1);
  41.       c(1,1:2)=[dx(2) xi31];
  42.       c(n,n-1:n) = [xin dx(n-2)];
  43.       b=zeros(n,1);
  44.       b(1)=((dx(1)+2*xi31)*dx(2)*divdif(1)+dx(1)^2*divdif(2))/xi31;
  45.       b(2:n-1)=3*(dx(2:n-1).*divdif(1:n-2)+dx(1:n-2).*divdif(2:n-1));
  46.       b(n)=...
  47.       (dx(n-1)^2*divdif(n-2)+(2*xin+dx(n-1))*dx(n-2)*divdif(n-1))/xin;
  48.  
  49.         % solve for the slopes and convert to pp form
  50.       s=c\b;
  51.       c4=(s(1:n-1)+s(2:n)-2*divdif(1:n-1))./dx;
  52.       c3=(divdif(1:n-1)-s(1:n-1))./dx - c4;
  53.       pp=ppmak(xi',[c4./dx c3 s(1:n-1) yi(1:n-1)],1);
  54.    end
  55.    if nargin==2,
  56.       output=pp;
  57.    else
  58.       output=ppual(pp,xx);
  59.    end
  60. end
  61.