home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 5.ddi / SPLINES.DI$ / OPTKNT.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  1.4 KB  |  37 lines
  1. function xi = optknt(tau,k)
  2. %OPTKNT    Optimal knots.
  3. %
  4. %        xi = optknt(tau,k)
  5. %
  6. % returns the (interior) optimal knots (according to the optimal recovery 
  7. % theory of Micchelli/Rivlin/Winograd and Gaffney/Powell) for interpolation 
  8. % at data points  tau(1), ..., tau(n)  by splines of order  k .
  9. %
  10. % These are the  n-k  sign-changes in the absolutely constant function  h  
  11. % which satisfies  integral{ f(x)h(x) : tau(1) < x < tau(n) } = 0  for all
  12. % splines  f  of order  k  with knot sequence  tau .
  13.  
  14. % C. de Boor: March 23, 1991
  15. % Copyright (c) 1990-92 by Carl de Boor and The MathWorks, Inc.
  16.  
  17. n = length(tau);nmk=n-k;
  18. signs=(-ones(1,nmk)).^[nmk-1:-1:0];
  19. tauext=[tau, tau(n)*ones(1,k)];
  20. % initial guess:
  21. xi = aveknt(tau,k);  
  22. tol = 1.e-7*(tau(n)-tau(1))/(nmk);maxiter=10;
  23.  
  24. % Newton iteration, a la `Computational aspects of optimal recovery'
  25. for iter=1:maxiter,
  26.    ais = [signs*spcol(tauext,k+1,xi) -1/2];     % a_1,...,a_n  from (18)
  27.    cumsum(ais(n:-1:1)')'; cumais=ans(n:-1:1);      % sum{a_j: j=i,...,n}
  28.    % now solve for the changes in  xi  (see eqns (20) and (21))
  29.    dx = ...
  30.   signs.*((-cumais(1:nmk).*(tau(k+[1:nmk])-tau([1:nmk]))/k)/spcol(tau,k,xi));
  31.    % make sure that change does not destroy order
  32.    xi = sort(xi + dx);
  33.    if (max(abs(dx))<tol), return, end
  34. end
  35. error('Newton iteration failed to converge in %f steps\n',maxiter)
  36.  
  37.