home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 5.ddi / SPLINES.DI$ / CSAPE.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  3.9 KB  |  93 lines
  1. function pp = csape(x,y,conds,valconds)
  2. %CSAPE    Cubic spline interpolation with various end-conditions.
  3. %
  4. %        pp = csape(x,y[,conds[,valconds]])
  5. %
  6. % returns the cubic spline interpolant (in pp-form) to the given data  (x,y)
  7. % using the specified end-conditions  conds(i)  with  values  valconds(i) ,
  8. % with i=1 (i=2) referring to the left (right) endpoint. 
  9. %   conds(i)=j  means that the j-th derivative is being specified to be
  10. %   valconds(i) , j=1,2. 
  11. %   conds(1)=0=conds(2)  means periodic end conditions.
  12. %  If conds(i) is not specified or is different from 0, 1 or 2, then the 
  13. %  default value for  conds(i)  is  1  and the default value of valconds(i) 
  14. %  is taken.
  15. %  If  valconds  is not specified, then the default value for valconds(i) is
  16. %       deriv. of cubic interpolant to nearest four points, if  conds(i)=1;
  17. %       0                                                   if  conds(i)=2.
  18.  
  19. % C. de Boor / latest change: December 3, 1990
  20. % C. de Boor / latest change: February 22, 1991 (added comments)
  21. % C. de Boor / latest change: 12 April 1992 (added periodic case and comments)
  22. % Copyright (c) 1990-92 by Carl de Boor and The MathWorks, Inc.
  23.  
  24.  
  25. %     Generate the cubic spline interpolant in pp form.
  26.  
  27. if (nargin<3), conds=[1 1];end
  28. if (nargin<4), valconds=[0 0]; end
  29.  
  30. n=length(x);[xi,ind]=sort(x);xi=xi(:);
  31. if n<2,
  32.    error('There should be at least two data points!')
  33. elseif all(diff(xi))==0,
  34.    error('The data abscissae should be distinct!')
  35. elseif n~=length(y),
  36.    error('Abscissa and ordinate vector should be of the same length!')
  37. else   
  38.    yi=y(ind);yi=yi(:);
  39.       % set up the linear system for solving for the slopes at  xi .
  40.    dx=diff(xi);divdif=diff(yi)./dx; 
  41.    c = diag([dx(2:n-1);0],-1)+2*diag([0;dx(2:n-1)+dx(1:n-2);0]) ...
  42.                  + diag([0;dx(1:n-2)],1);
  43.    b=zeros(n,1);
  44.    b(2:n-1)=3*(dx(2:n-1).*divdif(1:n-2)+dx(1:n-2).*divdif(2:n-1));
  45.    if (~any(conds)),
  46.       c(1,1)=1; c(1,n)=-1; b(1) = 0;
  47.    elseif (conds(1)==2),
  48.       c(1,1:2)=[2 1]; b(1)=3*divdif(1)-dx(1)*valconds(1)/2;
  49.    else, 
  50.       c(1,1:2) = [1 0]; b(1) = valconds(1);
  51.       if (nargin<4|conds(1)~=1), % if endslope was not supplied, 
  52.                                  % get it by local interpolation
  53.          b(1)=divdif(1);
  54.      if (n>2), ddf=(divdif(2)-divdif(1))/(xi(3)-xi(1));
  55.            b(1) = b(1)-ddf*dx(1); end
  56.      if (n>3), ddf2=(divdif(3)-divdif(2))/(xi(4)-xi(2));
  57.            b(1)=b(1)+(ddf2-ddf)*(dx(1)*(xi(3)-xi(1)))/(xi(4)-xi(1));end
  58.       end
  59.    end
  60.    if (~any(conds)),
  61.       c(n,1:2)=dx(n-1)*[2 1]; c(n,n-1:n)= c(n,n-1:n)+dx(1)*[1 2];
  62.       b(n) = 3*(dx(n-1)*divdif(1) + dx(1)*divdif(n-1));
  63.    elseif (conds(2)==2),
  64.       c(n,n-1:n)=[1 2]; b(n)=3*divdif(n-1)+dx(n-1)*valconds(2)/2;
  65.    else, 
  66.       c(n,n-1:n) = [0 1]; b(n) = valconds(2);
  67.       if (nargin<4|conds(2)~=1), % if endslope was not supplied,
  68.                  % get it by local interpolation
  69.          b(n)=divdif(n-1);
  70.      if (n>2), ddf=(divdif(n-1)-divdif(n-2))/(xi(n)-xi(n-2));
  71.        b(n) = b(n)+ddf*dx(n-1); end
  72.      if (n>3), ddf2=(divdif(n-2)-divdif(n-3))/(xi(n-1)-xi(n-3));
  73.        b(n)=b(n)+(ddf-ddf2)*(dx(n-1)*(xi(n)-xi(n-2)))/(xi(n)-xi(n-3));end
  74.       end
  75.    end
  76.  
  77.      % solve for the slopes and convert to pp form
  78.      % The final version should make use of MATLAB's eventual banded matrix
  79.      % capability, by setting up a routine for complete cubic spline inter-
  80.      % polation to vector data, which can then be used, without pivoting,
  81.      % to solve for three splines, one matching function values and zero
  82.      % end slopes, while the other two match zero values and slopes except
  83.      % except for a slope of  1  at one endpoint or the other. Any side 
  84.      % condition can then be obtained as an appropriate linear combination
  85.      % of these three, with coefficients determined from a 2x2 linear 
  86.      % system.
  87.  
  88.    s=c\b;
  89.    c4=(s(1:n-1)+s(2:n)-2*divdif(1:n-1))./dx;
  90.    c3=(divdif(1:n-1)-s(1:n-1))./dx - c4;
  91.    pp=ppmak(xi',[c4./dx c3 s(1:n-1) yi(1:n-1)],1);
  92. end
  93.