home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 5.ddi / SPLINES.DI$ / SPCOL.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  6.9 KB  |  190 lines
  1. function colloc=spcol(knots,k,tau,slvblk)
  2. % SPCOL    Generate the B-spline collocation matrix. 
  3. %
  4. %       colloc=spcol(knots,k,tau[,slvblk])
  5. %
  6. %  returns the matrix  colloc = (D^m(i)B_j(tau_i)) , with  
  7. %   B_j  the j-th B-spline of order  k  for the knot sequence  knots ,
  8. %   tau  a sequence of points, assumed to be  i n c r e a s i n g , and  
  9. %   m(i) := #{j<i: tau_j = tau_i}.
  10. %
  11. %     If the fourth argument is present, the matrix is returned in the almost
  12. %  block-diagonal format (specialised for splines) required by  slvblk.m .
  13. %
  14. %     The recurrence relations are used to generate the entries of the matrix.
  15.  
  16. % Carl de Boor / latest change: May 13, 1989
  17. % Copyright (c) 1990-92 by Carl de Boor and The MathWorks, Inc.
  18.  
  19. if (~isempty(find(diff(knots)<0))),
  20.    error('The knot sequence should be nondecreasing!')
  21. end
  22.  
  23. %  compute the number  n  of B-splines of order  k  supported by the given
  24. %  knot sequence and return an empty answer in case there aren't any.
  25.  
  26. npk=length(knots);n=npk-k;
  27. if (n<1);b=[];return;end
  28.  
  29. % remove all multiplicities from  tau .
  30. nrows=length(tau);
  31. index=[1 find(diff(tau)>0)+1];m=diff([index nrows+1]);pts=tau(index);
  32.  
  33. %set some abbreviations
  34. nd=min([k,max(m)]);   km1=k-1;
  35.  
  36. %  augment knot sequence to provide a k-fold knot at each end, in order to avoid
  37. % struggles near ends of interval  [knots(1),knots(npk)] . 
  38.  
  39. [augknot,addl]=augknt(knots,k);naug=length(augknot)-k;
  40. pts=pts(:);augknot=augknot(:);
  41.  
  42. %  For each  i , determine  savl(i)  so that  k <= savl(i) < naug+1 , and,
  43. % within that restriction,
  44. %        augknot(savl(i)) <= pts(i) < augknot(savl(i)+1) .
  45.  
  46. savl=max([sorted(augknot(1:naug),pts);k*ones(1,length(pts))]);
  47.  
  48. b=zeros(nrows,k);
  49.  
  50. % first do those without derivatives
  51. index1=find(m==1);
  52. if(~isempty(index1)),
  53.    pt1s=pts(index1);savls=savl(index1);lpt1=length(index1);
  54.    % initialize the  b  array.
  55.    lpt1s=index(index1); b(lpt1s,1)=ones(lpt1,1);
  56.    
  57.    % run the recurrence simultaneously for all  pt1(i) .
  58.    for j=1:km1;
  59.       saved=zeros(lpt1,1);
  60.       for r=1:j;
  61.          tr=augknot(savls+r)-pt1s;
  62.          tl=pt1s-augknot(savls+r-j);
  63.          term=b(lpt1s,r)./(tr+tl);
  64.          b(lpt1s,r)=saved+tr.*term;
  65.          saved=tl.*term;
  66.       end
  67.       b(lpt1s,j+1)=saved;
  68.    end
  69. end
  70.  
  71. % then do those with derivatives, if any:
  72. if (nd>1),
  73.    indexm=find(m>1);ptss=pts(indexm);savls=savl(indexm);lpts=length(indexm);
  74.    % initialize the  bb  array.
  75.    bb=ones(nd*lpts,1)*[1 zeros(1,km1)];
  76.    lptss=nd*[1:lpts];
  77.    
  78.    % run the recurrence simultaneously for all  pts(i) .
  79.    % First, bring it up to the intended level:
  80.    for j=1:k-nd;
  81.       saved=zeros(lpts,1);
  82.       for r=1:j;
  83.          tr=augknot(savls+r)-ptss;
  84.          tl=ptss-augknot(savls+r-j);
  85.          term=bb(lptss,r)./(tr+tl);
  86.          bb(lptss,r)=saved+tr.*term;
  87.          saved=tl.*term;
  88.       end
  89.       bb(lptss,j+1)=saved;
  90.    end
  91.  
  92.    % save the B-spline values in successive blocks in  bb .
  93.    
  94.    for jj=1:nd-1;
  95.       j=k-nd+jj; saved=zeros(lpts,1); lptsn=lptss-1;
  96.       for r=1:j;
  97.          tr=augknot(savls+r)-ptss;
  98.          tl=ptss-augknot(savls+r-j);
  99.          term=bb(lptss,r)./(tr+tl);
  100.          bb(lptsn,r)=saved+tr.*term;
  101.          saved=tl.*term;
  102.       end
  103.       bb(lptsn,j+1)=saved; lptss=lptsn;
  104.    end
  105.  
  106.    % now use the fact that derivative values can be obtained by differencing:
  107.  
  108.    for jj=nd-1:-1:1;
  109.       j=k-jj;
  110.       [jj:nd-1]'*ones(1,lpts)+ones(nd-jj,1)*lptsn;lptss=ans(:); 
  111.       for r=j:-1:1,
  112.          ones(nd-jj,1)*(augknot(savls+r)-augknot(savls+r-j))'/j;
  113.          bb(lptss,r)=-bb(lptss,r)./ans(:);
  114.          bb(lptss,r+1)=bb(lptss,r+1) - bb(lptss,r);
  115.       end
  116.    end
  117.  
  118.    % finally, combine appropriately with  b  by interspersing the multiple 
  119.    % point conditions appropriately: 
  120.    dtau=diff([tau(1)-1 tau(:)' tau(nrows)+1]);
  121.    index=find(min([dtau(2:nrows+1);dtau(1:nrows)])==0); % Determines all rows
  122.                                                     % involving multiple tau.
  123.    dtau=diff(tau(index));index2=find(dtau>0)+1;     % We need to make sure to
  124.    index3=[1 (dtau==0)];                            % skip unwanted derivs:
  125.    if (~isempty(index2)),
  126.                       index3(index2)=1+nd-m(indexm(1:length(indexm)-1));end
  127.    b(index,:)=bb(cumsum(index3),:);
  128.  
  129.    % and appropriately enlarge  savl
  130.    index=cumsum([1 (diff(tau)>0)]);  
  131.    savl=savl(index);
  132. end
  133.  
  134. % Finally, zero out all rows of  b  corresponding to  tau  outside 
  135. %  [knots(1),knots(npk)] .
  136.  
  137. index=find(tau<knots(1)|tau>knots(npk));
  138. if (~isempty(index)),
  139.    [1-nd:0]'*ones(1,length(index))+nd*ones(nd,1)*index(:)';   
  140.    b(ans(:),:)=zeros(nd*length(index),k);
  141. end
  142.  
  143. % The first B-spline of interest begins at knots(1), i.e., at  augknot(1+addl)
  144. % (since  augknot's  first knot has exact multiplicity  k ). If  addl<0 ,
  145. % this refers to a nonexistent index and means that the first  -addl  columns
  146. % of the collocation matrix should be trivial.  This we manage by setting
  147. savl=savl+max([0,-addl]);
  148.  
  149. if (nargin<4), % return the collocation matrix in standard matrix form
  150.    width=max([n,naug])+km1+km1;
  151.    cc=zeros(nrows*width,1);
  152.    cc([1-nrows:0]'*ones(1,k)+nrows*(savl'*ones(1,k)+ones(nrows,1)*[-km1:0]))=b;
  153.    % (This uses the fact that, for a column vector  v  and a matrix  A ,
  154.    %  v(A)(i,j)=v(A(i,j)), all i,j.)
  155.    colloc=zeros(nrows,n);
  156.    colloc(:)=...
  157.           cc([1-nrows:0]'*ones(1,n)+nrows*ones(nrows,1)*(max([0,addl])+[1:n]));
  158.  
  159. else,         % return the collocation matrix in almost block diagonal form.
  160.               % For this, make the blocks out of the entries with the same
  161.               %  savl(i) , with  last  computed from the differences.
  162.    % There are two issues, the change in the B-splines considered because of
  163.    % the use of  augknot  instead of  knots , and the possible drop of B-splines
  164.    % because the extreme  tau  fail to involve the extreme knot intervals.
  165.  
  166.    % savl(j) is the index in  augknot  of the left knot for  tau(j) , hence the
  167.    % corresponding row involves  B-splines to index  savl(j) wrto augknot, i.e.,
  168.    % B-splines to index  savl(j)-addl  wrto  knots. 
  169.    % Those with negative index are removed by cutting out their columns (i.e.,
  170.    % shifting the blocks in which they lie appropriately). Those with index 
  171.    % greater than  n  will be ignored because of  last .
  172.  
  173.    if (addl>0), % if B-splines were added on the left, remove them now:
  174.       width=km1+k;cc=zeros(nrows*width,1);
  175.       index=min([k*ones(1,length(savl));savl-addl]);
  176.     cc([1-nrows:0]'*ones(1,k)+nrows*(index'*ones(1,k)+ones(nrows,1)*[0:km1]))=b;
  177.       b(:)=cc([1-nrows:0]'*ones(1,k)+nrows*(ones(nrows,1)*(k+[0:km1])));
  178.       savl=savl+k-index;
  179.    end
  180.    ds=(diff(savl));
  181.    index=[0 find(ds>0) nrows];
  182.    rows=diff(index);
  183.    nb=length(index)-1;
  184.    last=ds(index(2:nb));
  185.    if(addl<0), nb=nb+1; rows=[0 rows]; last=[-addl last]; end
  186.    addr=naug-n-addl;
  187.    if(addr<0), nb=nb+1; rows=[rows 0]; last=[last -addr]; end
  188.    colloc=[41 nb rows k last n-sum(last) b(:)'];
  189. end
  190.