home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 5.ddi / SPLINES.DI$ / SPAPIDEM.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  8.3 KB  |  247 lines
  1. % SPAPIDEM    Demonstrate spline interpolation.
  2. % Copyright (c) 1990-92 by Carl de Boor and The MathWorks, Inc.
  3.                                           %latest change: November 28, 1990
  4. % latest change: May 16, 1991 change x to tau, ll  to  n , and add material 
  5. % on optimal spline interpolation.
  6.  
  7. echo on;clc;clg;home
  8. %    Here are some sample data, much used for testing spline approximation
  9. %  with variable knots, the so-called Titanium Heat Data, which record some
  10. %  property of titanium measured as a function of temperature. We'll use it
  11. %  merely to illustrate the use of spline interpolation.
  12.  
  13. [xx,yy]=titanium;
  14.  
  15. %    The plot of the data shows a rather sharp peak.
  16.  
  17. pause                % Touch any key to continue
  18.  
  19. frame=[-1 1 -.1 .1]+[min(xx),max(xx),min(yy),max(yy)];
  20. plot(xx,yy,'x');
  21. axis(frame);pause
  22. hold on
  23.  
  24. %     We'll pick a few data points from these somewhat rough data, since we
  25. %  want to interpolate it. 
  26.  
  27.  
  28. pick=[1 5 11 21  27 29 31 33 35 40 45 49];
  29. tau=xx(pick);y=yy(pick);
  30.  
  31. % Here is a picture of the data, with the selected data points marked.
  32. pause                % Touch any key to continue
  33. plot(tau,y,'o');pause
  34.  
  35. %    Since a spline of order  k  with  n+k  knots has  n  degrees of
  36. %  freedom, and we have 12 data points, a fit with a cubic spline, i.e., a
  37. %  fourth order spline, requires  12+4  knots.
  38.  
  39. %    Moreover, this knot sequence  t  must satisfy the Schoenberg-Whitney
  40. %  conditions, i.e., must be such that the i-th data
  41. %  abscissa lies in the support of the i-th B-spline, i.e.,
  42.  
  43. %                t(i) < tau(i) < t(i+k) , all  i ,
  44.  
  45. %  (with equality allowed in case of a knot of multiplicity  k ).
  46. %    We can achieve this condition for our example by using the data
  47. %  abscissa as knots, but adding two knots at either end.
  48.  
  49. n=length(tau);dl=tau(2)-tau(1);dr=tau(n)-tau(n-1);
  50. t=[tau(1)-dl*[2 1] tau tau(n)+dr*[1 2]];  % construct the knot sequence
  51. sp=spapi(t,tau,y);     % This constructs the spline.
  52.  
  53. pause                % Touch any key to continue
  54.  
  55. %    Now, for the plot:
  56.  
  57. pp=sp2pp(sp);        % converts to pp form
  58.  
  59. %    We don't care about the part of the curve outside the data interval,
  60. %  so we cut the  pp-rep. down to that interval....
  61.  
  62. pc=ppcut(pp,[tau(1) tau(n)]);
  63.  
  64. %           ....  and only then plot it,
  65.  
  66. values=fnplt(pc,'-.');pause  % saving the plotted points for later use.
  67.  
  68. % (The statement  values=fnplt(pp,'-.',[tau(1) tau(n)])  would have had the
  69. % same effect, without the necessity of introducing  pc .)
  70. hold off
  71.  
  72. %    A closer look at the left part of the spline fit shows some
  73. %  undulations.
  74.  
  75. pause                % Touch any key to continue
  76.  
  77. xxx=[0:40]/40*(tau(5)-tau(1))+tau(1);
  78. clf;
  79. plot(xxx,fnval(pc,xxx),'-.',tau,y,'o');
  80. axis([tau(1) tau(5) .6 1.2]);
  81. pause;
  82.  
  83. %     The unreasonable bump in the first interval stems from the fact that
  84. %  our spline goes smoothly to zero at its first knot. Here is a picture
  85. %  of the entire spline.
  86.  
  87. values2=fnplt(pp,'-.');pause
  88. hold on
  89.  
  90. %     Here are again the data points as well.
  91.  
  92. pause                 % Touch any key to continue
  93.  
  94. plot(tau,y,'o');pause
  95. hold off
  96.  
  97. %     There are many ways to enforce a more reasonable boundary behavior,
  98. %  e.g., by prescribing the slope (or higher derivatives) at the boundary,
  99. %  but they all come down to making sure that the interval of real interest
  100. %  to us does not extend all the way to the endpoints of the data interval.
  101.  
  102. %    Here is a simple idea: We add two more data points outside the given
  103. %  data interval and choose as our data there the values of the straight
  104. %  line through the first two data points.
  105.  
  106. pause                % Touch any key to continue
  107.  
  108. tt=[tau(1)-[4 3 2 1]*dl tau tau(n)+[1 2 3 4]*dr];
  109. xx=[tau(1)-[2 1]*dl tau tau(n)+[1 2]*dr];
  110. yy=[y(1)-[2 1]*(y(2)-y(1)) y y(n)+[1 2]*(y(n)-y(n-1))];
  111. sp2=spapi(tt,xx,yy);
  112. pp2=sp2pp(sp2);
  113. plot(values(1,:),fnval(pp2,values(1,:)),'-',tau,y,'o');
  114. axis(frame);pause
  115.  
  116. %     Here is also the original spline fit, dash-dotted, to show the reduction 
  117. %  of the undulation in the first and last interval.
  118.  
  119. pause                % Touch any key to continue
  120.  
  121. hold on;plot(values(1,:),values(2,:),'-.');pause;hold off;
  122. %    Finally, here is a closer look at the first four data intervals which
  123. %  shows more clearly the reduction of the undulation near the left end.
  124.  
  125. pause                % Touch any key to continue
  126.  
  127. plot(xxx,fnval(pp2,xxx),'-',tau,y,'o',xxx,fnval(pp,xxx),'-.');
  128. axis([tau(1) tau(5) .6 2.2]); pause;
  129. hold off
  130.  
  131. pause                % Touch any key to finish
  132.  
  133. %     In  o p t i m a l  spline interpolation,  at points  
  134. %                         tau(1), ..., tau(n)
  135. % say, one chooses the knots so as to minimize the constant in a standard 
  136. % error formula. Specifically, one chooses the first and the last data point 
  137. % as a k-fold knot. The remaining  n-k  knots are supplied by  optknt
  138.  
  139. help optknt
  140.  
  141. pause                % Touch any key to continue
  142.  
  143. % We try this alternative way for choosing knots for interpolation on our
  144. % example, in which we were interpolating by cubic splines to data 
  145. %  (tau(i),y(i)),  i=1,...,n :
  146.  
  147.  k = 4;   xi = optknt(tau,k);
  148. osp = spapi(augknt([tau(1) xi tau(n)],k),tau,y);
  149.  
  150. % Here is the plot
  151.  
  152. pause                % Touch any key to continue
  153.  
  154. fnplt(osp); pause
  155.  
  156. %  ... a bit disconcerting!
  157.  
  158. % And here is also the earlier interpolant, for comparison:
  159.  
  160. pause                % Touch any key to continue
  161.  
  162. hold on;
  163. plot(values(1,:), values(2,:),'-.'); pause
  164.  
  165. %  One can clearly see the given data, as the points were these two
  166. % interpolants coincide, but here are the data as well, to be sure:
  167.  
  168. pause                % Touch any key to continue
  169.  
  170. plot(tau,y,'o');pause
  171.  
  172. % Finally, here are also the (interior) optimal knots:
  173.  
  174. pause                % Touch any key to continue
  175.  
  176. plot(xi,0*xi,'x');pause
  177. hold off
  178.  
  179. %     The knot choice for optimal interpolation is designed to make the 
  180. % maximum over  a l l  functions  f  of the ratio   norm{f - If}/norm{D^k f}  
  181. % of the norm of the interpolation error  f - If  to the norm of the  k-th 
  182. % derivative  D^k f  of the interpoland as small as  possible.  Since our data 
  183. % imply that  D^k f  is rather large, the interpolation error near the flat 
  184. % part of the data is of acceptable size for such an `optimal' scheme.
  185.  
  186. %  But, to be sure that the optimal knots were calculated correctly, here is
  187. % a test. The (interior) optimal knots  xi  are determined in such a way 
  188. % that the spline function of order  k+1  with interior knots  xi  which
  189. % vanishes at all the data points  tau  is  p e r f e c t , i.e., has an 
  190. % absolutely constant  k-th derivative. 
  191.  
  192. %   Here is the verification: 
  193. pause                % Touch any key to continue
  194.  
  195.  
  196. data = [tau, tau(n)+1]; % add that data point to the right, to pin down a
  197.                         % unique spline of order  k+1  and vanishing at all
  198.                         % the tau(i),
  199. datay=[0*tau,1];   % ... and make the spline nonzero at that point.
  200.  
  201. perfect=spapi(augknt(sort([xi,data([1 n+1])]),k+1),data,datay);
  202.  
  203. %  Now plot the supposedly perfect spline, along with its zeros  tau :
  204. % Areas where it is relatively large are places where, all else being equal,
  205. % the error in optimal spline interpolation can be expected to be relatively 
  206. % large. Ideally, one would want to choose the data points  tau  so as to make
  207. % this perfect spline equi-oscillate.
  208.  
  209. pause                % Touch any key to continue
  210.  
  211. fnplt(perfect,'-',tau([1 n]));
  212.    title('the perfect spline for the given tau')
  213. hold on
  214. plot(tau,zeros(1,n),'o');pause
  215. hold off
  216. % now plot also the derivatives, including the piecewise constant one
  217. % that is supposed to be absolutely constant, always showing  tau  as well:
  218. dp = perfect;
  219. for j=1:k
  220.    dp=fnder(dp);
  221.    fnplt(dp,'-',tau([1 n]));
  222.    title([int2str(j),'. derivative'])
  223.    hold on
  224.    plot(tau,zeros(1,n),'o');pause
  225.    hold off
  226. end
  227.  
  228. %   This looks indeed absolutely constant. Note that its breakpoints  xi 
  229. % satisfy the Schoenberg-Whitney conditions with respect to the  tau  and for
  230. %  k = 4 , i.e.
  231. %
  232. %                 tau(i) < xi(i) < tau(i+4),   i = 1,...,n-4 
  233. %
  234. % as is made explicit in the next overlay, in which the intervals  
  235. %  [tau(i)..tau(i+4)]  are plotted, with  xi(i)  marked by an 'x'.
  236. % As you will see, each  xi(i)  lies even very close to the center of the
  237. % corresponding interval [tau(i)..tau(i+4)] :
  238.  
  239. pause                % Touch any key to continue
  240.  
  241. hold on
  242. for i=1:(n-k)
  243.    plot([tau(i),tau(i+k)],i*.5e-5*[1 1],'-')
  244.    plot(xi(i),i*.5e-5,'x')
  245. end
  246. hold off
  247.