home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 5.ddi / SPLINES.DI$ / CSAPIDEM.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  14.4 KB  |  380 lines
  1. % CSAPIDEM    Demonstrate cubic spline interpolation.
  2. % Copyright (c) 1990-92 by Carl de Boor and The MathWorks, Inc.
  3.                                     % latest change: March 22, 1991
  4.                                     % latest change: 27 mar 92 (periodic)
  5.                        % latest change: 13 apr 92 (other end conditions)
  6. clc;clg;echo on;home
  7. %     The M-file  csapi  constructs the cubic spline interpolant to given
  8. % data, using the "not-a-knot" end condition. The call
  9.  
  10. %     values = csapi(x,y,xx)
  11.  
  12. % returns the values at  xx  of that interpolating cubic spline. 
  13.  
  14. %     With  xi  the data abscissae in strictly increasing order, the cubic
  15. % spline interpolant is a piecewise cubic with breakpoints  xi  whose cubic 
  16. % pieces join together to form a function with two continuous derivatives. 
  17. % The "not-a-knot" end condition means that, at the first and last interior
  18. % breakpoint, even the third derivative is continuous (up to roundoff). 
  19.  
  20. %      Here are some trial runs.
  21.  
  22. pause                                             %Touch any key to continue
  23. %     Just two data points result in a straight line interpolant:
  24.  
  25. pause                                             %Touch any key to continue
  26. xx=[0:60]/10;x=[0 1];y=[2 0];
  27. plot(xx,csapi(x,y,xx),'-',x,y,'o');pause
  28.  
  29. %     Three data points give a parabola:
  30.  
  31. pause                                             %Touch any key to continue
  32. x=[2 3 5];y=[1 0 4];
  33. plot(xx,csapi(x,y,xx),'-',x,y,'o');pause
  34.  
  35. %     Four or more data points give in general a cubic spline.
  36.  
  37. pause                                             %Touch any key to continue
  38.  
  39. x=[1 1.5 2 4.1 5];y=[1 -1 1 -1 1];
  40. plot(xx,csapi(x,y,xx),'-',x,y,'o');pause
  41.  
  42. pause                                             %Touch any key to continue
  43.  
  44. %     These look like nice interpolants, but how do we check it?
  45.  
  46. % We already saw that we interpolate, for we always plotted the data points
  47. % and the interpolant went right through those points.
  48.  
  49. %     But, for checking purposes, it is better to choose specific data taken
  50. % from a cubic spline so that the entire interpolant should coincide with
  51. % that function.
  52. %     One such is the truncated third power, i.e., the function  
  53. %  x ]-> subplus(x - knot)^3 , with  knot  one of the breaks and  subplus
  54. % the  t r u n c a t i o n  f u n c t i o n , i.e., 
  55. help subplus
  56. pause                                             %Touch any key to continue
  57.  
  58.  
  59. % We try it with  knot = 2:
  60. x=[0:6]; y=subplus(x-2).^3;
  61. % This is data taken from a cubic spline with just one break, at 2.
  62.  
  63. values = csapi(x,y,xx);
  64.  
  65. % we plot the interpolant:
  66.  
  67. pause                                             %Touch any key to continue
  68. plot(xx,values);pause
  69. % The function is zero to the left of  2  and rises like (x-2)^3 to the 
  70. % right of 2. To compare it with subplus(.-2)^3, we plot the difference:
  71.  
  72. pause                                             %Touch any key to continue
  73. plot(xx,values-subplus(xx-2).^3);pause
  74. % Since the maximum of the given function values is 
  75. max_y=max(abs(y))
  76. % the size of the absolute error (as indicated by the plot) is acceptably low.
  77.  
  78. pause                                             %Touch any key to continue
  79.  
  80. %   As a further test, we try to interpolate to a truncated power which
  81. % cannot be interpolated exactly. For example, the first interior 
  82. % breakpoint of the interpolating spline is not really a knot since the
  83. % interpolant has three continuous derivatives there. Hence we should not be
  84. % able to reproduce the truncated power centered at that point. We try 
  85.  
  86. values = csapi(x,subplus(x-1).^3,xx);
  87. pause                                             %Touch any key to continue
  88. plot(xx,values-subplus(xx-1).^3);pause
  89. %The difference is as large as .2, but decays rapidly as we move away from
  90. % 1. This ilustrates that cubic spline interpolation is essentially local.
  91.  
  92. pause                                             %Touch any key to continue
  93.  
  94. %     It is possible to retain the interpolating cubic spline in a form
  95. % suitable for later, repeated evaluation, or for calculating its 
  96. % derivatives, or for other manipulations.
  97. %     This is done by calling  csapi  in the form
  98. %
  99. %          pp = csapi(x,y)
  100. %
  101. % which returns, in  pp , the pp-form of the interpolant. This form can be 
  102. % evaluated at some points  xx  by
  103. %
  104. %         values = fnval(pp,xx)
  105. %
  106. % It can be differentiated by
  107. %
  108. %        dpp = fnder(pp)     
  109. %
  110. % or integrated by
  111. %
  112. %        ipp = fnint(pp)     
  113. %
  114. % which return , in  dpp  or  ipp , the pp-form of the derivative or primitive.
  115. pause                                             %Touch any key to continue
  116. %     For example, we try again the interpolant to the truncated power.
  117. pp = csapi(x,subplus(x-2).^3);
  118. % and look at its derivatives:
  119. pause                                             %Touch any key to continue
  120. dpp = fnder(pp);plot(xx,fnval(dpp,xx));pause
  121. % This should coincide with  3*subplus(.-2).^3. We look at the error:
  122. pause                                             %Touch any key to continue
  123. plot(xx,fnval(dpp,xx)-3*subplus(xx-2).^2);pause
  124. % and that looks ok.
  125.  
  126. %     The second derivative of the interpolant should be  6*subplus(.-2).
  127. % We try it.
  128. pause                                             %Touch any key to continue
  129. ddpp = fnder(dpp);plot(xx,fnval(ddpp,xx)-6*subplus(xx-2));pause
  130.  
  131. % There are jumps, but they are within roundoff.
  132.  
  133. pause                                             %Touch any key to continue
  134.  
  135. % Here is also the integral
  136. ipp = fnder(pp);plot(xx,fnval(ipp,xx));pause
  137. % and here is the difference between the constructed and the ideal integral:
  138. ipp=fnint(pp);plot(xx,fnval(ipp,xx)-subplus(xx-2).^4/4);pause
  139.  
  140. pause                                             %Touch any key to continue
  141.  
  142. %  The M-file  csape  also provides a cubic spline interpolant to given data,
  143. % but permits various other end conditions. It does not directly return values 
  144. % of that interpolant, but only its pp-form.  Its simplest version,
  145.  
  146. %       pp = csape(x,y)
  147.  
  148. % uses the Lagrange end conditions, which are better at times than the
  149. % not-a-knot condition used by  csapi . 
  150.  
  151. pause                                             %Touch any key to continue
  152.  
  153. % For example, here is again the interpolant to the truncated power which  
  154. %  csapi  fails to reproduce:
  155.  
  156. values = csapi(x,subplus(x-1).^3,xx);
  157. pause                                             %Touch any key to continue
  158.  
  159. plot(xx,values-subplus(xx-1).^3);pause;
  160. hold on
  161.  
  162. % For comparison, here is the error when we use the Lagrange end conditions:
  163. pause                                             %Touch any key to continue
  164.  
  165. plot(xx, fnval(csape(x,subplus(x-1).^3),xx) - subplus(xx-1).^3,'-.' )
  166. title('not-a-knot (-) vs. Lagrange (-.)');pause;
  167. hold off
  168.  
  169. % Well, in this case, there isn't much difference.
  170. pause                                             %Touch any key to continue
  171.  
  172. %     The M-file  csape  only provides the interpolating cubic spline in 
  173. % pp-form, but for several different end conditions. For example, the call
  174. %
  175. %         pp = csape(x,y,[2 2])
  176. %        
  177. % uses the `natural' end conditions. This means that the second derivative is 
  178. % zero at the two extreme breakpoints. This is indicated here by specifying
  179. % a  2  for both ends, meaning specification of 2nd derivatives, but leaving
  180. % the 2nd derivative value to the default, which is  0 .
  181.  
  182. %   We apply `natural' cubic spline interpolation to the truncated power, and
  183. % plot the error:
  184.  
  185. pp = csape(x,subplus(x-1).^3,[2 2]);
  186. pause                                             %Touch any key to continue
  187.  
  188. plot(xx, fnval(pp,xx) - subplus(xx-1).^3 );pause;
  189.  
  190. % Not surprisingly, the error is not small, for the simple reason that we are
  191. % interpolating to a cubic spline function whose second derivative happens not 
  192. % to be zero at the right endpoint. 
  193.  
  194. pause                                             %Touch any key to continue
  195.  
  196. % When we use explicitly the correct second derivatives, we get a small error:
  197.  
  198. pp=csape(x,subplus(x-1).^3,[2 2],6*subplus(x([1 length(x)])-1));
  199. %                           ^ ^  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  200. %                           | |  value of 2.deriv.at endpoints
  201. plot(xx, fnval(pp,xx) - subplus(xx-1).^3 );pause;
  202.  
  203. pause                                             %Touch any key to continue
  204. %  csape  also permits specification of endpoint  s l o p e s . This is the
  205. %  c l a m p e d  spline  (or,  c o m p l e t e  cubic spline interpolant).
  206. % The statement
  207. %
  208. %         pp = csape(x,y,[1 1],[sl,sr])
  209. %        
  210. % supplies the cubic spline interpolant to the data  x , y  which also
  211. % has slope  sl  at the first data point and slope  sr  at the last data 
  212. % point.
  213. pause                                             %Touch any key to continue
  214. %   It is even possible to mix these conditions. For example, our much
  215. % exercised truncated power function  f(x) = subplus(x-1)^3  has slope  0
  216. % at  x = 0   and second derivative  30  at  x = 6 (our last data abscissa).
  217.  
  218. % We therefore expect no error in the following interpolant:
  219.  
  220. pp=csape(x,subplus(x-1).^3,[1 2],[0,30]);
  221. %                           ^ ^   ^ ^^
  222. %                           | |   | ||
  223.  
  224. plot(xx, fnval(pp,xx) - subplus(xx-1).^3 );pause;
  225.  
  226. %   ... and we were not disappointed.
  227.  
  228. pause                                             %Touch any key to continue
  229.  
  230. %    It is also possible to prescribe  p e r i o d i c  end conditions.
  231. % For example, a periodic interpolant to the ordinates 
  232. y = [0 -1 0 1 0];
  233. % at the abscissae
  234. x = (pi/2)*[-2:2];
  235. % should be a good approximation to the sine function:
  236. pp = csape(x,y,[0 0]);
  237. %               ^ ^ 
  238. %               | |  specification of periodic end conditions
  239.  
  240. %   Here is a plot of the difference between the sine function and this
  241. % periodic piecewise cubic function, ...
  242. xx = (pi/50)*[-50:50];
  243. plot(xx, sin(xx) - fnval(pp,xx), 'x'); pause;
  244.  
  245. %   ... showing a relative error of only 2 percent. Not bad!
  246.  
  247. pause                                             %Touch any key to continue
  248.  
  249. %     Finally, any end condition not covered explicitly by  csapi  or  csape
  250. % can be handled by constructing the interpolant with the default side 
  251. % conditions, and then adding to it an appropriate scalar multiple of
  252. % an interpolant to zero values and some side conditions. If there are two
  253. % `nonstandard' side conditions to be satisfied, one may have to solve a 
  254. % 2x2 linear system first.
  255.  
  256. pause                                             %Touch any key to continue
  257.  
  258.  
  259. %     For example, suppose that you want to enforce the condition
  260. %  
  261. %         lambda(s)  :=  a Ds(e)  + b D^2 s(e) = c
  262. %
  263. % on the cubic spline interpolant to data  
  264.  
  265. x = 0:.25:3;
  266. y = x.*(-1 + x.*(-1+x.*x/5));
  267.  
  268. % with
  269.  
  270. e = x(1); a = 2; b = -3; c = 4;
  271.  
  272. % For simplicity, the data are taken from a quartic polynomial which happens to 
  273. % satisfy the specific side condition.
  274.  
  275. pause                                             %Touch any key to continue
  276.  
  277. % Then, in addition to the interpolant with the default end conditions ...
  278.  
  279. pp1 = csape(x,y);
  280. p1 = pppce(pp1,1); dp1 = fnder(p1);
  281.  
  282. % ... and its first derivative, we also construct the cubic
  283. % spline interpolant to zero data, and with a slope of  1  at  e .
  284.  
  285. pp0 = csape(x,0*y,[1,0],[1,0]);
  286. p0 = pppce(pp0,1); dp0 = fnder(p0);
  287.  
  288. % Then we compute  lambda(s)  for both  s = pp1  and  s = pp0  ...
  289.  
  290. lam1 = a*fnval(dp1,e) + b*fnval(fnder(dp1),e);
  291. lam0 = a*fnval(dp0,e) + b*fnval(fnder(dp0),e);
  292.  
  293. pause                                             %Touch any key to continue
  294.  
  295. % ... and construct the right linear combination of  pp1  and
  296. %  pp0  to get a cubic spline 
  297.  
  298. %  pp  :=  pp1 + (c - lambda(pp1))/lambda(pp0) *pp0
  299.  
  300. % which does satisfy the desired condition (as well as the default end 
  301. % condition at the right endpoint). 
  302. %   We form the linear combination with the help of  fncmb :
  303.  
  304. pp = fncmb(pp0,(c-lam1)/lam0,pp1);
  305.  
  306. pause                                             %Touch any key to continue
  307.  
  308. % We expect  pp  to fit our quartic polynomial slightly better near  e  
  309. % than does the interpolant  pp1  with the default conditions ...
  310.  
  311. xx = (-.3):.05:.7;
  312. yy = xx.*(-1 + xx.*(-1+xx.*xx/5));
  313. plot(xx, fnval(pp1,xx) - yy, 'x')
  314. hold on
  315. plot(xx, fnval(pp,xx) - yy, 'o')
  316. title('default (x) vs. nonstardard (o)');pause
  317. hold off
  318.  
  319. %  ... and we are not disappointed:
  320. pause                                             %Touch any key to continue
  321.  
  322.  
  323. % If we also want to enforce the condition
  324.  
  325. %    mu(s) := D^3(3) = 14.6
  326.  
  327. % (which our quartic also satisfies), then we construct an additional cubic
  328. % spline interpolating to zero values, and with zero first derivative at
  329. % the left endpoint, hence certain to be independent from  pp0, ...
  330.  
  331. pp2 = csape(x,0*y,[0,1],[0,1]);
  332.  
  333. % ... and solve the linear system for the coefficients in the linear
  334. % combination    pp :=  pp1 + d0 * pp0  + d2 * pp2
  335. % for which  lambda(pp) = c, mu(pp) = 14.6 . (Note that both  pp0   and
  336. %  pp2  vanish at all interpolation points, hence   pp  will match the 
  337. % given data for any choice of  d0  and  d2 ).
  338.  
  339. pause                                             %Touch any key to continue
  340.  
  341. %     For amusement, we use MATLAB's encoding facility to write a loop
  342. % for computing  lambda(ppj), and mu(ppj), for j=0:2
  343.  
  344. dd = zeros(2,3);
  345. for j=0:2
  346.    J = num2str(j);
  347.    eval(['dpp',J,'=fnder(pp',J,');']);
  348.    eval(['ddpp',J,'=fnder(dpp',J,');']);
  349.    eval(['dd(1,1+',J,')=a*fnval(dpp',J,',e)+b*fnval(ddpp',J,',e);']);
  350.    eval(['dd(2,1+',J,')=fnval(fnder(ddpp',J,'),3);']);
  351. end
  352.  
  353. d = dd(:,[1,3])\([c;14.6]-dd(:,2));
  354. pp = fncmb(fncmb(pp0,d(1),pp2,d(2)),pp1);
  355.  
  356. xxx = 0:.05:3;
  357. yyy = xxx.*(-1 + xxx.*(-1+xxx.*xxx/5));
  358. plot(xxx, yyy - fnval(pp,xxx),'x');pause
  359. pause                                             %Touch any key to continue
  360.  
  361.  
  362. %    For reassurance, we compare this error with the one obtained in complete
  363. % cubic spline interpolation to this function:
  364.  
  365. hold on
  366. plot(xxx, yyy - fnval(csape(x,y,[1 1],[-1,-7+(4/5)*27]),xxx),'o')
  367. title('       nonstandard (x) vs endslopes (o)');pause
  368. hold off
  369.  
  370. %    The errors differ (and not by much) only near the end points, testifying
  371. % to the fact that both  pp0  and  pp2  are sizable only near their respective
  372. % end points.
  373.  
  374. %    As a final check, here is the third derivative of  pp  at  3 (which
  375. % should be  14.6 ):
  376.  
  377. fnval(fnder(fnder(fnder(pp))),3)
  378. pause                                             %Touch any key to continue
  379.  
  380.