home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 5.ddi / SPLINES.DI$ / SPALLDEM.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  6.7 KB  |  178 lines
  1. % SPALLDEM    Demonstrate form and usage of 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
  5.  
  6. %  A  s p l i n e  is specified by its (nondecreasing)  k n o t  sequence
  7.  
  8. %                      t
  9.  
  10. % and by its B-spline coefficient sequence
  11.  
  12. %                      a
  13.  
  14. % With  n := length(a), and  n+k := length(t), the spline is
  15. % of  o r d e r   k . This means that its polynomial pieces have
  16. % degree < k .
  17.  
  18. pause;                % Touch any key to continue
  19.  
  20. %     This means, explicitly, that
  21.  
  22. %             n
  23. %  spline := sum  B_(i,k)*a(i) ,
  24. %            i=1
  25.  
  26. % with  B_(i,k)  the i-th B-spline of order  k  for the given knot
  27. % sequence  t  , i.e., the B-spline with knots  t(i), ..., t(i+k) .
  28.  
  29.  
  30. %      The B-spline with knots  t(i) <= ... <= t(i+k)  is
  31. % positive on the interval  (t(i) .. t(i+k))  and is zero outside the interval.
  32. % It is pp of order  k  with breaks at the points  t(i), ..., t(i+k) . These
  33. % knots may coincide, and the precise  m u l t i p l i c i t y  of a knot
  34. % governs the smoothness with which the two polynomial pieces join there.
  35.  
  36. %     Here is the picture of a cubic B-spline, which is typical of B-splines
  37. % with simple knots.
  38.  
  39. pause;                % Touch any key to continue
  40.  
  41. fnplt(spmak([0 1.3 1.9 3.1 4],1));pause
  42.  
  43. % See the demo  bspldem  for various examples, particularly for the interplay
  44. % between knot multiplicity and smoothness.
  45.  
  46. %    To repeat: A spline is a linear combination of B-splines. Usually, a
  47. %  spline is constructed from some information, like function values and/or
  48. %  derivative values, or as the approximate solution of some ordinary
  49. %  differential equation.  But it is also possible to make up a spline from
  50. %  scratch, by providing its knot sequence and its coefficient sequence to
  51. %  the M-file  spmak .
  52.  
  53. pause;                %  Touch any key to continue
  54.  
  55. %    For example, we might supply the uniform knot sequence [1:10] and the
  56. %  coefficient sequence [3:8] . Since there are 10 knots and 6 coefficients,
  57. %  the order must be 4, i.e., we get a cubic spline.
  58.  
  59. sp=spmak([1:10],[3:8])
  60. echo off;
  61. ttype=11;td=1;tn=6;tt=[1:10];tk=4;tcoefs=[3:8];
  62. zcheck(sp-[ttype td tn tcoefs tk tt]);echo on
  63. pause;                %  Touch any key to continue
  64.  
  65. %  The output may be a bit bewildering but, as with pp's, you don't really
  66. %  have to know just how the spline is stored in the vector  sp . If you
  67. %  do ever want to know the constituents, you can always unmake the spline 
  68. %  using  spbrk  as follows.
  69. pause                    % Touch any key to continue
  70.  
  71. [knots,coefs,n,k]=spbrk(sp)
  72.  
  73. pause                                  % Touch any key to continue
  74. echo off;
  75. zcheck(knots-tt);zcheck(coefs-tcoefs);zcheck([n k]-[tn tk]);echo on;
  76.  
  77. %  But the point of the spline toolbox is that there shouldn't be any need
  78. %  for you to look up these details. You simply use  sp  as an argument to
  79. %  M-files that evaluate, differentiate, integrate, convert, or plot the 
  80. %  spline whose description is contained in  sp .
  81.  
  82. %      For example, here is a plot of the spline we just made up. Notice
  83. % how this spline vanishes at the endpoints of its basic interval, due to
  84. % the fact that all its knots are simple.
  85.  
  86. pause                                  % Touch any key to continue
  87. fnplt(sp);pause
  88.  
  89. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  90.  
  91.  
  92. %     A spline is usually constructed as an approximation to some function.
  93. %  Information about that function is typically available as function values.
  94. %  Here, for example, is a data set which has become a standard test case, the
  95. %  t i t a n i u m  h e a t  data. It is challenging since it is relatively
  96. %  flat except for one sharp peak.
  97. pause                                  % Touch any key to continue
  98. [x,y]=titanium;
  99. plot(x,y,'o');
  100. frame=[x(1),x(length(x)),min(y)-1,max(y)+1]; axis(frame); pause
  101. hold on
  102. %    As a first attempt at fitting this data, we try a C^1 piecewise parabolic
  103. %  with six pieces and uniformly placed breaks:
  104.  
  105. k=3; l=6; breaks=x(1)+[0:l]*(x(length(x))-x(1))/l;
  106.  
  107. %    We use the M-file  augknt  to convert the given breakpoint sequence into
  108. %  the knot sequence whose corresponding B-splines provide a basis for the pp
  109. %  space chosen:
  110.  
  111. knots=augknt(breaks,k)
  112. echo off;
  113. zcheck(knots-[595,595,595,675,755,835,915,995,1075,1075,1075]);echo on
  114.  
  115. pause                                  % Touch any key to continue
  116.  
  117. %     The M-file  spap2  provides a least-squares spline fit to given data:
  118.  
  119. splin0=spap2(knots,k,x,y);
  120.  
  121. %  We plot this approximation on top of the data:
  122. pause                                  % Touch any key to continue
  123. fnplt(splin0);pause
  124.  
  125. %   The approximation is not very good, but it is at least correct in the 
  126. % sense that the spline crosses the given function eight times, i.e., at
  127. % least as many times as there are degrees of freedom (i.e., B-splines) in
  128. % our approximating family. Recall that that number is
  129. length(knots)-k
  130. % and here is the plot for another look.
  131. pause                                  % Touch any key to continue
  132. hold off
  133.  
  134. %    The best way to judge the quality of an approximation is to look at the
  135. % error.  Specifically, one looks for systematic behavior and, when using
  136. % pp functions, tries to reduce the error by placing additional breakpoints
  137. % in the area of largest error. Here is the error plot:
  138. pause                                  % Touch any key to continue
  139. error0=y-fnval(splin0,x);
  140. plot(x,error0);hold on;plot(frame(1:2),[0 0]);pause
  141.  
  142. %    Now we add breaks at the point or points of greatest error:
  143. break1=sort([breaks,860,900,950]);
  144. splin1=spap2(augknt(break1,k),k,x,y);
  145. erro1=y-fnval(splin1,x);
  146. pause                                  % Touch any key to continue
  147. plot(x,erro1,'-.');pause
  148.  
  149. %    There is only one systematic part still left. 
  150. % To fit it, we make  900  a double knot, thus allowing the approximating
  151. % spline a discontinuous first derivative (since it is quadratic):
  152. break2=sort([break1,900]);
  153. splin2=spap2(augknt(break2,k),k,x,y);
  154. erro2=y-fnval(splin2,x);
  155. pause                                  % Touch any key to continue
  156.  
  157. plot(x,erro2,':');pause
  158.  
  159. %    There is still some systematic oscillation in the error. We place an 
  160. % additional breakpoint at its maximum:
  161. break3=sort([break2,880]);
  162. splin3=spap2(augknt(break3,k),k,x,y);
  163. erro3=y-fnval(splin3,x);
  164. pause                                  % Touch any key to continue
  165. plot(x,erro3,'-');pause
  166.  
  167. %   This could be improved further, but you get the idea. 
  168.  
  169. %   Here are the data and the last spline fit constructed:
  170. hold off
  171. plot(x,fnval(splin3,x))
  172. axis(frame);
  173. hold on
  174. plot(x,y,'o');pause
  175.  
  176. hold off
  177. splst
  178.