home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 5.ddi / SPLINES.DI$ / PPALLDEM.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  4.3 KB  |  112 lines
  1. % PPALLDEM    Demonstrate form and usage of piecewise polynomials.
  2. % Copyright (c) 1990-92 by Carl de Boor and The MathWorks, Inc.
  3.                                             % latest change: September 28, '90
  4. clc;clg;echo on;home                             
  5. %       A piecewise polynomial, or   pp   for short, is characterized by its
  6. %  b r e a k  sequence,   break   say, and the   c o e f f i c i e n t 
  7. %  array,   coefs   say, of the local power form of its polynomial pieces.
  8. %
  9. %     The break sequence is assumed to be strictly increasing,
  10. %               break(1) < break(2) < ... < break(l+1),
  11. %  with  l  the number of polynomial pieces which make up  pp .
  12. %
  13. %      While these polynomials may be of varying degrees, they are
  14. %  all recorded as polynomials of the same  o r d e r   k , i.e., the
  15. %  coefficient array  c  is of size  (l,k) , with  coefs(j,:)  containing
  16. %  the  k  coefficients in the local power form for the j-th polynomial
  17. %  piece.
  18. pause                    % Touch any key to continue
  19. %
  20. %       With this, the precise description of  pp  in terms of the break
  21. %   sequence  break  and the coefficient array  coefs  is
  22. % 
  23. %   pp(t) = polyval(coefs(j,:),t-break(j))  for   break(j) <= t < break(j+1) ,
  24. %
  25. %  Here,  polyval(a,x)  is the MATLAB function; it returns the number 
  26. %     sum{a(j)x^(k-j):j=1,..,k} = a(1)*x^(k-1) + a(2)*x^(k-2) + ... + a(k)*x^0 .
  27. %
  28. %  For  t not in [break(1),break(l+1)),  pp(t)  is defined by extending the 
  29. % first, resp. last, polynomial piece.  
  30. %
  31. pause                    % Touch any key to continue
  32. %    A pp is usually constructed in some M-file, through a process of
  33. %  interpolation or approximation, or conversion from some other form (e.g.,
  34. %  from a spline), and is output as a single (row) vector. But it is also
  35. %  possible to make one up from scratch, using the statement
  36.  
  37. %           pp=ppmak(breaks,coefs).
  38.  
  39. pause                    % Touch any key to continue
  40.  
  41. %    Let's try that.
  42.  
  43. pp=ppmak([-5:-1],[-22:-11])
  44. echo off;
  45. ttype=10;td=1;tl=4;tbreak=[-5:-1];tk=3;tcoef=[-22:-20;-19:-17;-16:-14;-13:-11];
  46. zcheck(pp-[ttype td tl tbreak tk tcoef(:)']);
  47. % Check also integration and differentiation. 
  48. % Evaluation and taking a piece are checked below.
  49. zcheck(pp-fnder(fnint(pp)));echo on
  50. pause                    % Touch any key to continue
  51.  
  52. %    Back comes this row vector in which you could actually find all those
  53. %  numbers input, since they are all negative and distinct. But you thereby
  54. %  also see some other numbers. In fact, it is not worthwhile to wonder just
  55. %  how the pp is described in this vector.
  56.  
  57. pause                    % Touch any key to continue
  58.  
  59. %    If you really want to know the constituent parts of pp, use
  60. %           [breaks,coefs,l,k]=ppbrk(pp)
  61.  
  62. pause(5);
  63.  
  64. [breaks,coefs,l,k]=ppbrk(pp)
  65. pause                    % ... and there they are. Touch any key to continue
  66. echo off
  67. zcheck(breaks-tbreak);zcheck(coefs-tcoef);zcheck([l k]-[tl tk]);echo on
  68. %    Here are some operations you can perform on  pp.
  69.  
  70. %       v=fnval(pp,x)     evaluates,
  71.  
  72. %       dpp=fnder(pp)     differentiates,
  73.  
  74. %       dpp=fnint(pp)     integrates,
  75.  
  76. %       pj=pppce(pp,j)   pulls out the j-th polynomial piece.
  77.  
  78. %       pc=ppcut(pp,[a b])  pulls out the part relevant for the interval
  79. %                           [a,b] .
  80.  
  81. %       (values = ) fnplt(pp)  plots pp between its extreme breaks.
  82.  
  83. %     For example, here is a plot of the particular pp we just made up.
  84.  
  85. pause                     % Touch any key to continue
  86.  
  87. x=[-55:-5]/10; plot(x,fnval(pp,x),'-.');pause
  88.  
  89.  
  90. %   Here, added to the plot, are the breaklines
  91.  
  92. pause                     % Touch any key to continue
  93.  
  94. yy=axis;hold on;
  95. for j=1:l+1;plot([breaks(j) breaks(j)],yy(3:4));end;pause
  96.  
  97. %   And here is, added to that plot, the plot of its third piece
  98.  
  99. pause                     % Touch any key to continue
  100.  
  101. plot(x,fnval(pppce(pp,3),x));pause
  102. hold off; echo off
  103. zcheck(-14+(x+3).*(-15 + (x+3)*(-16))-fnval(pppce(pp,3),x));echo on
  104.  
  105. %     While the pp-representation of a pp is efficient for evaluation, the
  106. %  c o n s t r u c t i o n  of a pp from some data is usually more
  107. %  efficiently handled by determining first its
  108. %   B - r e p r e s e n t a t i o n , i.e., its representation as a linear
  109. %  combination of B-splines.    For this, look at      spalldemo  .
  110.  
  111. pause                     % Touch any key to finish
  112.