home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 5.ddi / SPLINES.DI$ / CSAPS.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  3.7 KB  |  100 lines
  1. function output = csaps(x,y,p,xx)
  2. % CSAPS    Cubic smoothing spline.
  3. %
  4. %      values = csaps(x,y,p,xx)
  5. %
  6. % returns the values at  xx  of the cubic smoothing spline for the given data
  7. %  (x,y)  and depending on the smoothing parameter  p  from  [0,1] . For  p=0,
  8. % this is the least-squares straight line fit to the data, while, on the other
  9. % extreme, i.e., for  p=1 , this is the `natural' cubic spline interpolant.
  10. %
  11. %      The alternative call
  12. %
  13. %      pp = csaps(x,y,p)
  14. %
  15. % returns the pp-form of the cubic smoothing spline instead, for later use 
  16. % with  fnval,  fnder, etc.
  17.  
  18. % C. de Boor / latest change: Sep.2, 1989
  19. % Copyright (c) 1990-92 by Carl de Boor and The MathWorks, Inc.
  20.  
  21. %     Generate the cubic spline interpolant in pp form, depending on
  22. % the number of data points.
  23.  
  24. n=length(x);[xi,ind]=sort(x);xi=xi(:);
  25. output=[];
  26. if n<2,
  27.    error('There should be at least two data points!')
  28. elseif all(diff(xi))==0,
  29.    error('The data abscissae should be distinct!')
  30. elseif n~=length(y),
  31.    error('Abscissa and ordinate vector should be of the same length!')
  32. else   
  33.    yi=y(ind);yi=yi(:);
  34.    if (n==2), % the smoothing spline is the straight line interpolant
  35.       pp=ppmak(xi',[diff(yi)./diff(xi) yi(1)],1);
  36.    else % set up the linear system for solving for the 2nd derivatives at  xi .
  37.         % this is taken from (XIV.6)ff of the `Practical Guide to Splines'
  38.         % with the diagonal matrix  D = eye(n,n)
  39.       dx=diff(xi);divdif=diff(yi)./dx;
  40.       odx=ones(n-1,1)./dx;
  41.      if (n<50),
  42.         diag([dx(1:n-2);0],-1) + 2*diag([0;dx(2:n-1)+dx(1:n-2);0]) ...
  43.                       + diag([0;dx(2:n-1)],1);
  44.         R=ans(2:n-1,2:n-1);
  45.         diag([odx(1:n-2);0],-1) - diag([0;odx(2:n-1)+odx(1:n-2);0])...
  46.                       + diag([0;odx(2:n-1)],1);
  47.         Qt=ans(2:n-1,:);
  48.         % solve for the 2nd derivatives 
  49.         u=(6*(1-p)*Qt*Qt'+p*R)\diff(divdif);
  50.      else, % use explicit Gauss elimination to handle the larger banded
  51.            % matrices. The number 50 was picked without much thought.
  52.            % For very large  n  , iteration might be a good alternative.
  53.         % generate the five diagonals of Qt*Qt'
  54.         P=zeros(n-2,6);
  55.         oddx=odx(1:n-2)+odx(2:n-1);
  56.         odx.*odx;
  57.         P(:,3)=ans(1:n-2)+oddx.*oddx+ans(2:n-1);
  58.         P(:,4)=-[(oddx(1:n-3)+oddx(2:n-2)).*odx(2:n-2);0];
  59.         P(:,5)=[odx(2:n-3).*odx(3:n-2);0;0];
  60.         % combine with the three diagonals of R to form (6*(1-p)*Qt*Qt'+p*R)
  61.         P=(6*(1-p))*P;
  62.         P(:,3)=P(:,3)+(2*p)*(dx(2:n-1)+dx(1:n-2));
  63.         P(:,4)=P(:,4)+p*dx(2:n-1);
  64.         P(:,2)=[0;P(1:n-3,4)];
  65.         P(:,1)=[0;0;P(1:n-4,5)];
  66.         % solve for the 2nd derivatives 
  67.         % u=(6*(1-p)*Qt*Qt'+p*R)\diff(divdif);
  68.         % by Gauss elimination without pivoting
  69.         P(:,6)=diff(divdif);
  70.         for j=1:n-4;
  71.            P(j,4:6)=(1/P(j,3))*P(j,4:6);
  72.            P(j+1,[3 4 6])=P(j+1,[3 4 6])-P(j+1,2)*P(j,[4:6]);
  73.            P(j+2,[2 3 6])=P(j+2,[2 3 6])-P(j+2,1)*P(j,[4:6]);
  74.         end
  75.         j=n-3;
  76.            P(j,4:6)=(1/P(j,3))*P(j,4:6);
  77.            P(j+1,[3 4 6])=P(j+1,[3 4 6])-P(j+1,2)*P(j,[4:6]);
  78.         j=n-2;
  79.            P(j,6)=P(j,6)/P(j,3);
  80.         j=n-3;
  81.            P(j,6)=P(j,6)-P(j,4)*P(j+1,6);
  82.         for j=n-4:-1:1;
  83.            P(j,6)=P(j,6)-P(j,[4 5])*P(j+[1 2],6);
  84.         end
  85.         u=P(:,6);
  86.      end
  87.         % ... and convert to pp form
  88.       % Qt'*u=diff([0;diff([0;u;0])./dx;0])
  89.       yi = yi - 6*(1-p)*diff([0;diff([0;u;0])./dx;0]);c3 = [0;p*u;0];
  90.       c2=diff(yi)./dx-dx.*(2*c3(1:n-1)+c3(2:n));
  91.       pp=ppmak(xi',[diff(c3)./dx,3*c3(1:n-1),c2,yi(1:n-1)],1);
  92.    end
  93.    if nargin==3,
  94.       output=pp;
  95.    else
  96.       output=ppual(pp,xx);
  97.    end
  98. end
  99.  
  100.