home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 5.ddi / SPLINES.DI$ / BSPLIDEM.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  3.4 KB  |  106 lines
  1. %BSPLIDEM    Demonstrate B-splines.
  2. % Copyright (c) 1990-92 by Carl de Boor and The MathWorks, Inc.
  3.                                           % latest change: 24 January 1992
  4. clc;clg;home;echo on;hold off
  5.  
  6. %      The B-spline
  7.  
  8. %                     B(. | t(i), ..., t(i+k))
  9.  
  10. % with knots  t(i) <= ... <= t(i+k)  is positive on the interval 
  11. %  (t(i) .. t(i+k))  and is zero outside the interval.  It is pp of order  k
  12. % with breaks at the points  t(i), ..., t(i+k) . These knots may coincide, 
  13. % and the precise  m u l t i p l i c i t y  of a knot governs the smoothness
  14. % with which the two polynomial pieces join there.
  15.  
  16. %     Here are some pictures of B-splines, together with the
  17. % polynomials whose pieces make up the B-spline. The pictures also show the
  18. % knot locations.
  19.  
  20. pause;                % Touch any key to continue
  21.  
  22. for j=1:7,
  23.    bspline([0:j]); pause
  24. end
  25.  
  26. %     Here is a sequence of pictures of a cubic B-spline as more and more
  27. % knots become coincident.
  28.  
  29. %    First, a double knot develops. 
  30.  
  31. pause                % Touch any key to continue
  32.  
  33. t=[0:4];k=4;bspline(t);
  34. window=221;i=1;propor=[1/2,1/8,1/32,0];delta=t(i+1)-t(i+2);
  35. for j=3:-1:0
  36.    t(i+1)=t(i+2)+propor(4-j)*delta;pp=bspline(t,window+j);pause
  37. end
  38. t(i+1)=t(i+2)+delta;
  39.  
  40. %    In the last picture, t(i+1)=t(i+2) and, correspondingly, the continuity
  41. %  of the second derivative is lost. This can be seen quite nicely by
  42. %  looking at the two polynomials which join there: One is convex, the other
  43. %  concave at that point.
  44.  
  45. pause                % Touch any key to continue
  46.  
  47. xx=[t(i)*10:t(i+k)*10]/10;pl=pppce(pp,1);pr=pppce(pp,2);
  48. clf
  49. plot(xx,fnval(pl,xx),xx,fnval(pr,xx),t(i+2)*[1 1],[-20,10]);pause
  50.  
  51. %    Next, we let a triple knot develop and expect, correspondingly, to see
  52. %  a jump in the first derivative at that break.
  53.  
  54. pause                % Touch any key to continue
  55.  
  56. clg
  57. delta=t(i+1:i+2)-t(i+3);
  58. for j=3:-1:0
  59.    t(i+1:i+2)=t(i+3)+propor(4-j)*delta;pp=bspline(t,window+j);pause
  60. end
  61. t(i+1:i+2)=t(i+3)+delta;
  62.  
  63. %    We were not disappointed.
  64.  
  65. %    Finally, a quadruple knot, at which, for a fourth order spline, we
  66. %  expect a discontinuity in the function values.
  67.  
  68. pause                % Touch any key to continue
  69.  
  70. clg
  71. delta=t(i+1:i+3)-t(i+4);
  72. for j=3:-1:0
  73.    t(i+1:i+3)=t(i+4)+propor(4-j)*delta;pp=bspline(t,window+j);pause
  74. end
  75.  
  76. t=[0 1 1 3 4 6 6 6];
  77. clc;clg;home
  78. %     The rule connecting smoothness across a knot with the multiplicity of
  79. %  that knot is easy to remember if we designate the number of smoothness
  80. %  conditions across a knot as the condition multiplicity there.
  81. %  The rule is:
  82.  
  83. %           knot multiplicity + condition multiplicity = order.
  84.  
  85.  
  86. %  For example, for a B-spline of order  3 , a simple knot would mean  2
  87. %  smoothness conditions, i.e., continuity of function and first derivative,
  88. %  while a double knot would only leave one smoothness condition, i.e., just
  89. %  continuity, and a triple knot would leave no smoothness condition, i.e.,
  90. %  even the function would be discontinuous.
  91.  
  92. pause                %  Touch any key to continue
  93.  
  94. %     Here is a picture of all the third-order B-splines for a certain
  95. %  knot sequence. For each breakpoint, try to determine its multiplicity in
  96. %  the knot sequence, as well as its multiplicity as a knot in each of the
  97. %  B-splines.
  98.  
  99. pause                %  Touch any key to continue
  100.  
  101. x=[-10:70]/10;
  102. c=spcol(t,3,x);[l,m]=size(c);c=c+ones(l,1)*[0:m-1];
  103. axis([-1 7 0 m]);hold on;for tt=t;plot([tt tt],[0 m],'-');end;plot(x,c);pause
  104. hold off;
  105.  
  106.