home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 5.ddi / SPLINES.DI$ / CHEBLOOP.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.6 KB  |  77 lines
  1. % CHEBLOOP    Used by CHEBDEM, Chebyshev demonstration.
  2. % Copyright (c) 1990-92 by Carl de Boor and The MathWorks, Inc.
  3. %                                          % latest change: February 13, '90
  4. echo on
  5. % For the complete levelling, we use the  Remez algorithm. This means that we 
  6. % construct new tau's as the extrema of our current approximation to  c  
  7. % and try again.  We find the new interior  tau  as the zeros of  c' .
  8.  
  9.      pause;                      % touch any key to continue
  10.  
  11. cp = fnder(c);
  12.  
  13. % We take the zeros of the corresponding control polygon as our first
  14. % guess. For this, we must take apart the spline  cp .
  15.  
  16. [knots,coefs,np,kp] = spbrk(cp);
  17.  
  18. % The control polygon has the vertices  (tstar(i),coefs(i)) , with
  19. %  tstar  the knot averages for the spline as supplied by
  20.  
  21. tstar=aveknt(knots,kp);
  22.  
  23. % Here are the zeros of the resulting control polygon of  cp :
  24. npp=[1:np-1];
  25. guess=tstar(npp) - coefs(npp).*(diff(tstar)./diff(coefs));
  26. plot(guess,zeros(1,np-1),'o');pause
  27.  
  28. % This provides already a very good first guess for the actual zeros
  29.  
  30.      pause;                      % touch any key to continue
  31.  
  32.  
  33. % Now we compute the first derivative at both these points
  34. points=ones(4,1)*tau(2:n-1);
  35. points(1,:)=guess;
  36. values=0*points;
  37. values(1,:)= fnval(cp,points(1,:));
  38. values(2,:)= fnval(cp,points(2,:));
  39.  
  40. % use two steps of the secant method, but guard against division by zero:
  41. for j=2:3
  42.    rows=[j,j-1];cpd=diff(values(rows,:));
  43.    index=find(cpd==0);cpd(index)=ones(1,length(index));
  44.    points(j+1,:)=points(j,:)-values(j,:).*(diff(points(rows,:))./cpd);
  45.    values(j+1,:)=fnval(cp,points(j+1,:));
  46. end
  47.    max(abs(values'))
  48.  
  49.      pause;                      % touch any key to continue
  50.  
  51. % This is a good improvement. Now we take these  points  as our  tau :
  52. tau = [tau(1) points(4,:) tau(n)]
  53. % .. and check the extreme values of our current approximation there
  54. extremes=abs(fnval(c,tau));
  55. % The difference
  56. max(extremes)-min(extremes)
  57. % is an estimate of how far we are from total levelling. 
  58. echo off       % for a quick check
  59. if ~exist('answer'),
  60. % format long
  61. % ans
  62. % format short
  63. zcheck(ans-0.69047027491413);end;echo on
  64. % Here is the picture of the new spline corresponding to that
  65. c=spapi(t,tau,b);
  66. points=sort([tau [0:100]*(t(n+1)-t(k))/100]);
  67. values=fnval(c,points);plot(points,values);pause
  68.  
  69. % If this is not close enough, simply try again, starting from these new  tau .
  70. %
  71. echo off
  72.   fprintf(' If you want to continue, type  1 , then RETURN.')
  73.   answer = input(' Otherwise, just hit RETURN:  ');
  74. if ~isempty(answer)&(answer==1), chebloop
  75. else,hold off
  76. end
  77.