home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 5.ddi / SPLINES.DI$ / TSPDEM.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  10.8 KB  |  313 lines
  1. %TSPDEM    Demonstrate tensor product splines.
  2. % Copyright (c) 1990-92 by Carl de Boor and The MathWorks, Inc.
  3.                                     % latest change: May 24, 1991
  4. clc;clg;echo on;home
  5. % Since the toolbox can handle splines with  v e c t o r  coefficients, it is
  6. % easy to implement interpolation or approximation to gridded data by tensor
  7. % product splines.
  8. % 
  9. % Consider, for example, least-squares approximation to given data
  10. %  z(i,j) = f(x(i),y(j)), i=1,...,I, j=1,...,J , e.g.,
  11.  
  12. x = sort([[0:10]/10,.03 .07, .93 .97]);
  13. y = sort([[0:6]/6,.03 .07, .93 .97]);
  14. z = franke(x'*ones(1,length(y)),ones(length(x),1)*y);
  15.  
  16. %  Here is a plot of the data ..
  17.  
  18. pause                                               %Touch any key to continue
  19.  
  20. mesh(z,[50,50]);pause
  21.  
  22. % It doesn't give an accurate image of the data because the MATLAB function
  23. %  mesh  has no facilities for data from a  n o n uniform rectangular grid.
  24. % I chose point density somewhat higher near the boundary, to pin down the
  25. % approximation there.
  26.  
  27. pause                                               %Touch any key to continue
  28.  
  29. % Next, we choose some appropriate knot sequence and spline order for the 
  30. % y-direction, e.g.,
  31.  
  32. ky = 3; knotsy = augknt([0,.25,.5,.75,1],ky);
  33.  
  34. % Then
  35.  
  36. sp=spap2(knotsy,ky,y,z);
  37.  
  38. % provides a spline approximation to all the curves  (y,z(i,:)) .
  39. % In particular, the statement
  40.  
  41. yy=[-.1:.05:1.1];
  42. vals = fnval(sp,yy);
  43.  
  44. % provides the array  vals , whose entry  vals(i,j)  can be taken as an
  45. % approximation to the value  f(x(i),yy(j))  of the underlying function
  46. %  f  at the mesh-point (x(i),yy(j)) .  This is evident when we plot  vals  
  47. % using  mesh :
  48.  
  49. pause                                               %Touch any key to continue
  50.  
  51. mesh(vals,[50,50]);pause
  52.  
  53. % Note that, for each  x(i) , both the first two and the last two values
  54. % are zero since both the first two and the last two points in  yy  are
  55. % outside the basic interval for the spline in  sp .
  56.  
  57. % Note also the ridges. They confirm that we are plotting here smooth curves in
  58. % one direction only.
  59.  
  60. pause
  61.  
  62. % To get an actual surface, you now have to go a step further.  Look at the
  63. % coefficients  coefsy  of the spline in  sp :
  64.  
  65. [ignored,coefsy]=spbrk(sp);
  66.  
  67. % Abstractly, you can think of the spline in  sp  as the function
  68.  
  69. %       y |--> sum_r coefsy(r,:) B_{r,ky}(y) 
  70.  
  71. % with the i-th entry  coefsy(r,i)  of the vector coefficient  coefsy(r,:)
  72. % corresponding to  x(i) , all  i .  This suggests approximating each 
  73. % coefficient vector  coefsy(r,:)  by a spline of the same order  kx  and with 
  74. % the same appropriate knot sequence  knotsx :
  75.  
  76. kx = 4; knotsx=augknt([0:.2:1],kx);
  77. sp2=spap2(knotsx,kx,x,coefsy');
  78.  
  79. pause                                               %Touch any key to continue
  80.  
  81. % Note that  spap2(knots,k,x,fx)  expects  fx(:,j)  to be the datum at  x(j) ,
  82. % i.e., expects each  c o l u m n  of  fx  to be a datum. Since we wanted to 
  83. % fit the data  coefsy(r,:)  at  x(r) , all  r , we have to present  spap2  
  84. % with the  t r a n s p o s e  of  coefsy .
  85.  
  86. % Now consider the transpose   coefs = cxy'  of the coefficients  cxy  of the 
  87. % resulting spline `curve' :
  88.  
  89. [ignored,cxy]=spbrk(sp2);
  90. coefs=cxy';
  91.  
  92. % It provides the  b i v a r i a t e  spline approximation
  93. % 
  94. %     (x,y) |--> sum_q sum_r coefs(q,r) B_{q,kx}(x) B_{r,ky}(y)
  95. % 
  96. % to the original data
  97. % 
  98. %    (x(i),y(j)) |-->  z(x(i),y(j)) .
  99.  
  100. pause                                               %Touch any key to continue
  101.  
  102. % To evaluate this spline surface at the grid (xv(i),yv(j)), e.g., the grid
  103.  
  104. xv=[0:.025:1];yv=[0:.025:1];
  105.  
  106. % you can do the following:
  107.  
  108. values = spcol(knotsx,kx,xv)*coefs*spcol(knotsy,ky,yv)';
  109.  
  110. % Here is the mesh-plot of these values: 
  111.  
  112. pause                                               %Touch any key to continue
  113.  
  114. mesh(values,[50,50]);pause
  115.  
  116. % This makes good sense since  spcol(knotsx,kx,xv)  is the matrix whose
  117. % (i,q)-entry equals the value  B_{q,kx}(xv(i))  at  xv(i)  of the  q-th 
  118. % B-spline of order  kx  for the knot sequence  knotsx .
  119.  
  120. pause                                               %Touch any key to continue
  121.  
  122. % Since the matrices  spcol(knotsx,kx,xv)  and  spcol(knotsy,ky,yv)  are 
  123. % banded, it may be more efficient (though perhaps more memory-consuming)
  124. % for `large'  xv  and  yv  to make use of   fnval , as follows:
  125.  
  126. value2=fnval(spmak(knotsx,fnval(spmak(knotsy,coefs),yv)'),xv)';
  127.  
  128.  
  129. % Here is a check, viz. the  r e l a t i v e  difference between the values 
  130. % computed in these two different ways:
  131.  
  132. disp( max(max(abs(values-value2)))/max(max(abs(values))) )
  133.  
  134. pause                                               %Touch any key to continue
  135.  
  136. % Here is also a plot of the error, i.e., the difference between the given data
  137. % and the value of the approximation at those data points:
  138.  
  139. errors = z -  spcol(knotsx,kx,x)*coefs*spcol(knotsy,ky,y)';
  140.  
  141. pause                                               %Touch any key to continue
  142.  
  143. mesh(errors,[50,50]);pause
  144.  
  145. % The  r e l a t i v e  error is of size
  146. disp( max(max(abs(errors)))/max(max(abs(z))) )
  147.  
  148. % This is perhaps not too impressive. On the other hand, we used only ...
  149. disp(size(coefs))
  150. % ... coefficients to fit a data array of size
  151. disp(size(z))
  152.  
  153.  
  154. pause                                               %Touch any key to continue
  155.  
  156. %    The approach followed here seems  b i a s e d :  We first think of the
  157. % given data  z  as describing a vector-valued function of  y , and then we
  158. % treat the matrix formed by the vector coefficients of the approximating 
  159. % curve as describing a vector-valued function of  x .
  160.  
  161. %    What happens when we take things in the opposite order, i.e., think of
  162. %  z  as describing a vector-valued function of  x , and then treat the matrix
  163. % made up from the vector coefficients of the approximating curve as describing
  164. % a vector-valued function of  y ?
  165.  
  166. %    Perhaps surprisingly, the final approximation is the same (up to
  167. % roundoff). Here is the numerical experiment.
  168.  
  169. %    First, we fit a spline curve to the data, but this time with  x  as the
  170. % independent variable, hence it is the   r o w s  of  z  which now become the
  171. % data points. Correspondingly, we must supply  z' (rather than  z ) to spap2:
  172.  
  173. pause                                               %Touch any key to continue
  174.  
  175. spb=spap2(knotsx,kx,x,z');
  176.  
  177. % ... thus obtaining a spline approximation to all the curves  (x,z(:,j)) .
  178. % In particular, the statement
  179.  
  180. valsb = (fnval(spb,xv))';
  181.  
  182. % provides the array  valsb , whose entry  valsb(i,j)  can be taken as an
  183. % approximation to the value  f(xv(i),y(j))  of the underlying function
  184. %  f  at the mesh-point (xv(i),y(j)) .  This is evident when we plot  valsb  
  185. % using  mesh :
  186.  
  187. pause                                               %Touch any key to continue
  188.  
  189. mesh(valsb,[50,50]);pause
  190.  
  191. % Note the ridges. They confirm that we are plotting here smooth curves in
  192. % one direction only.
  193.  
  194. pause                                               %Touch any key to continue
  195.  
  196. % Here, for a quick comparison is the earlier mesh of curves, with the smooth
  197. % curves running in the other direction:
  198.  
  199. pause                                               %Touch any key to continue
  200.  
  201. mesh(vals,[50,50]);pause
  202.  
  203. %    Now comes the second step, to get the actual surface. First, extract 
  204. % the coefficients:
  205.  
  206. [ignored,coefsx]=spbrk(spb);
  207.  
  208. % Abstractly, you can think of the spline in  spb  as the function
  209.  
  210. %       x |--> sum_r coefsx(r,:) B_{r,kx}(x) 
  211.  
  212. % with the j-th entry  coefsy(r,j)  of the vector coefficient  coefsx(r,:)
  213. % corresponding to  y(j) , all  j .  Thus we now fit each
  214. % coefficient vector  coefsx(r,:)  by a spline of the same order  ky  and with 
  215. % the same (appropriate) knot sequence  knotsy :
  216.  
  217. spb2=spap2(knotsy,ky,y,coefsx');
  218.  
  219. pause                                               %Touch any key to continue
  220.  
  221. % Note that we need again to transpose the coefficient array from  spb , 
  222. % since  spap2  takes the columns of that last input argument as the data 
  223. % points.  
  224.  
  225. % Correspondingly, there is no need now to transpose the coefficient array
  226. %  coefsb  of the resulting  `curve':
  227.  
  228. [ignored,coefsb]=spbrk(spb2);
  229.  
  230. %    The claim is that  coefsb  equals the earlier coefficient array  coefs
  231. % (up to round-off), and here is the test:
  232.  
  233. disp( max(max(abs(coefs - coefsb))) )
  234.  
  235. %     Not bad, eh? - See the User's Guide for an explanation.
  236.  
  237. %  Thus, the  b i v a r i a t e  spline approximation
  238. % 
  239. %     (x,y) |--> sum_q sum_r coefsb(q,r) B_{q,kx}(x) B_{r,ky}(y)
  240. % 
  241. % to the original data
  242. % 
  243. %    (x(i),y(j)) |-->  z(x(i),y(j)) 
  244. %
  245. % obtained coincides with the earlier one (which generated  coefs  rather than
  246. %  coefsb ).
  247.  
  248.  
  249. pause                                               %Touch any key to continue
  250.  
  251. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  252.  
  253. %    Since the data come from a smooth function, we should be interpolating
  254. % it, i.e., use  spapi  instead of  spap2 , or, equivalently, use  spap2
  255. % with the appropriate knot sequences. For illustration, here is the same 
  256. % process done with  spapi .
  257.  
  258. pause                                               %Touch any key to continue
  259.  
  260. % To recall, the data points were
  261. % x = sort([[0:10]/10,.03 .07, .93 .97]);
  262. % y = sort([[0:6]/6,.03 .07, .93 .97]);
  263. % 
  264.  
  265. % We use again quadratic splines in  y , hence use knots midway between data
  266. % points:
  267.  
  268. ly = length(y);
  269. knotsy = augknt(sort([0 1 (y(2:(ly-2))+y(3:(ly-1)))/2]),ky);
  270.  
  271. spi=spapi(knotsy,y,z);
  272. [ignored,coefsy] = spbrk(spi);
  273. pause                                               %Touch any key to continue
  274.  
  275. % We use again cubics in  x , and use the not-a-knot condition, therefore use
  276. % all but the second and the second-last data point as knots:
  277.  
  278. lx = length(x);
  279. knotsx = augknt(x([1,[3:(lx-2)],lx]),kx);
  280. spi2 = spapi(knotsx,x,coefsy');
  281.  
  282. [ignored,cxy] = spbrk(spi2);
  283. icoefs = cxy';
  284.  
  285. pause                                               %Touch any key to continue
  286.  
  287. % We are ready to evaluate the interpolant at a fine mesh:
  288.  
  289. ivalues = spcol(knotsx,kx,xv)*icoefs*spcol(knotsy,ky,yv)';
  290.  
  291. % Here is the mesh-plot of these values: 
  292.  
  293. pause                                               %Touch any key to continue
  294.  
  295. mesh(ivalues,[50,50]);pause
  296.  
  297.  
  298. % Its error, as an approximation to the Franke function, is computed next:
  299.  
  300. fvalues = franke(xv'*ones(1,length(yv)),ones(length(xv),1)*yv);
  301. error =  fvalues - ivalues;
  302.  
  303. % The error is of  r e l a t i v e  size ... 
  304. disp( max(max(abs(error)))/max(max(abs(fvalues))) )
  305.  
  306. %  ... which is better than the error at the data points of the
  307. % earlier least-squares approximation. For comparison, here is also the 
  308. % overall relative error of that least-squares approximation:
  309.  
  310. disp( max(max(abs(fvalues-values)))/max(max(abs(fvalues))) )
  311. pause                                               %Touch any key to continue
  312.  
  313.