%SPAPIDM2 Demonstrate spline interpolation. % This Expo demo adapted from ... % SPAPIDEM Demonstrate spline interpolation. % Carl de Boor 1990 % Ned Gulley, 6-21-93 % Denise Chen, 9-2-93 % % Copyright (c) 1984-94 by The MathWorks, Inc. % Demo initialization ==================== if ~exist('SlideShowGUIFlag'), figNumber=0; end; [xx,yy]=alonso; % The plot of the data shows a rather sharp peak. frame=[-1 1 -.1 .1]+[min(xx),max(xx),min(yy),max(yy)]; plot(xx,yy,'x'); axis(frame); drawnow; if ssinit(figNumber), str= str2mat( ... ' ', ... ' Press the "Start" button to see a demonstration of a spline', ... ' curve being used to fit real world data.'); ssdisp(figNumber,str); if figNumber, return; end end % Beginning of the demo ================== % Here are some sample data, much used for testing spline approximation % with variable knots, the so-called Titanium Heat Data, which record some % property of alonso measured as a function of temperature. We'll use it % merely to illustrate the use of spline interpolation. str= str2mat( ... ' Here are some sample data which record a property of', ... ' alonso measured as a function of temperature. We''ll use it to', ... ' illustrate the use of spline interpolation.', ... ' ', ... ' >> [xx, yy] = alonso;', ... ' >> frame = [-1 1 -.1 .1]+[min(xx), max(xx), min(yy), max(yy)];', ... ' >> plot(xx, yy, ''x'');', ... ' >> axis(frame);', ... ' >> drawnow;'); ssdisp(figNumber,str); if sspause(figNumber), return; end; str= ... [' We will pick a few data points from these somewhat ' ' rough data, since we want to interpolate it. Here is a ' ' picture of the data, with the selected data points marked.']; str1 = str2mat( ... ' >> pick = [1 5 11 21 27 29 31 33 35 40 45 49];', ... ' >> tau = xx(pick);', ... ' >> y = yy(pick);', ... ' >> plot(tau, y, ''ro'');'); str = str2mat(str,'',str1); ssdisp(figNumber,str); hold on pick=[1 5 11 21 27 29 31 33 35 40 45 49]; tau=xx(pick);y=yy(pick); plot(tau,y,'ro'); drawnow if sspause(figNumber), return; end; % Since a spline of order k with n+k knots has n degrees of % freedom, and we have 12 data points, a fit with a cubic spline, i.e., a % fourth order spline, requires 12+4 knots. % Moreover, this knot sequence t must satisfy the Schoenberg-Whitney % conditions, i.e., must be such that the i-th data % abscissa lies in the support of the i-th B-spline, i.e., % t(i) < tau(i) < t(i+k) , all i , % (with equality allowed in case of a knot of multiplicity k ). % We can achieve this condition for our example by using the data % abscissa as knots, but adding two knots at either end. str= str2mat( ... ' Since a spline of order k with n+k knots has n degrees of',... ' freedom, and we have 12 data points, a fit with a fourth order',... ' spline requires 12+4 = 16 knots. Moreover, this knot sequence',... ' t must be such that the i-th data abscissa lies in the support',... ' of the i-th B-spline. We achieve this by using the data',... ' abscissa as knots, but adding two knots at either end.',... '',... ' >> n = length(tau); dl = tau(2)-tau(1); dr = tau(n)-tau(n-1);',... ' >> t = [tau(1)-dl*[2 1] tau tau(n)+dr*[1 2]];',... ' >> sp = spapi(t, tau, y);'); ssdisp(figNumber,str); n=length(tau);dl=tau(2)-tau(1);dr=tau(n)-tau(n-1); t=[tau(1)-dl*[2 1] tau tau(n)+dr*[1 2]]; % construct the knot sequence sp=spapi(t,tau,y); % This constructs the spline. if sspause(figNumber), return; end; str = str2mat( ... ' Now, for the plot. We first convert the spline from B-form to ', ... ' pp-form. Since we do not care about the part of the curve', ... ' outside the data, we cut the pp-form down to that interval and ', ... ' then plot it.', ... '', ... ' >> pp = sp2pp(sp); ', ... ' >> pc = ppcut(pp, [tau(1) tau(n)]);', ... ' ', ... ' Evaluate and plot "pc" ', ... ' >> values = fnplt(pc);'); ssdisp(figNumber,str); % Now, for the plot: pp=sp2pp(sp); % converts to pp form % We don't care about the part of the curve outside the data interval, % so we cut the pp-rep. down to that interval.... pc=ppcut(pp,[tau(1) tau(n)]); % .... and only then plot it, values=fnplt(pc); drawnow hold off % saving the plotted points for later use. % (The statement values=fnplt(pp,'-.',[tau(1) tau(n)]) would have had the % same effect, without the necessity of introducing pc .) if sspause(figNumber), return; end; str= str2mat( ... ' A closer look at the left part of the spline fit shows some', ... ' undulations.', ... '', ... ' >> xxx = [0:40] / 40 * (tau(5)-tau(1)) + tau(1);', ... ' >> plot(xxx, fnval(pc, xxx), tau, y, ''ro'');', ... ' >> axis([tau(1) tau(5) 0.6 1.2]);'); ssdisp(figNumber,str); xxx=[0:40]/40*(tau(5)-tau(1))+tau(1); plot(xxx,fnval(pc,xxx),tau,y,'ro'); axis([tau(1) tau(5) 0.6 1.2]); drawnow; if sspause(figNumber), return; end; str= str2mat( ... ' The unreasonable bump in the first interval stems from the fact', ... ' that our spline goes smoothly to zero at its first knot.', ... ' ', ... ' Here is a picture of the entire spline.', ... ' ', ... ' >> values2 = fnplt(pp);'); ssdisp(figNumber,str); values2=fnplt(pp); drawnow if sspause(figNumber), return; end; str= str2mat( ... ' Here are again the data points as well.' ,... ' ', ... ' >> plot(tau, y, ''ro'');'); ssdisp(figNumber,str); hold on plot(tau,y,'ro');drawnow hold off if sspause(figNumber), return; end; % There are many ways to enforce a more reasonable boundary behavior, % e.g., by prescribing the slope (or higher derivatives) at the boundary, % but they all come down to making sure that the interval of real interest % to us does not extend all the way to the endpoints of the data interval. % Here is a simple idea: We add two more data points outside the given % data interval and choose as our data there the values of the straight % line through the first two data points. str = str2mat( ... ' Here is a simple way to enforce a more reasonable boundary', ... ' behavior. We add two more data points outside the given data', ... ' interval and choose as our data there the values of the straight', ... ' line through the first two data points.', ... ' ', ... ' >> tt = [tau(1)-[4 3 2 1]*dl tau tau(n)+[1 2 3 4]*dr];', ... ' >> xx = [tau(1)-[2 1]*dl tau tau(n)+[1 2]*dr];', ... ' >> yy = [y(1)-[2 1]*(y(2)-y(1)) y y(n)+[1 2]*(y(n)-y(n-1))];', ... ' >> sp2 = spapi(tt, xx, yy); pp2 = sp2pp(sp2);', ... ' >> plot(values(1,:), fnval(pp2, values(1,:)), ''-'', tau, y, ''ro'');'); ssdisp(figNumber,str); tt=[tau(1)-[4 3 2 1]*dl tau tau(n)+[1 2 3 4]*dr]; xx=[tau(1)-[2 1]*dl tau tau(n)+[1 2]*dr]; yy=[y(1)-[2 1]*(y(2)-y(1)) y y(n)+[1 2]*(y(n)-y(n-1))]; sp2=spapi(tt,xx,yy); pp2=sp2pp(sp2); plot(values(1,:),fnval(pp2,values(1,:)),'-',tau,y,'ro'); axis(frame); drawnow if sspause(figNumber), return; end; str = str2mat( ... ' Here is also the original spline fit, shown in blue to show the', ... ' reduction of the undulation in the first and last interval.', ... ' ', ... ' >> plot(values(1,:), values(2,:), ''c'');'); ssdisp(figNumber,str); hold on; plot(values(1,:),values(2,:),'c'); hold off; drawnow; if sspause(figNumber), return; end; str = str2mat( ... ' Finally, here is a closer look at the first four data intervals', ... ' which shows more clearly the reduction of the undulation near', ... ' the left end.', ... '', ... ' >> hold on', ... ' >> plot(xxx, fnval(pp2,xxx), ''-'');', ... ' >> plot(tau, y, ''ro'');', ... ' >> plot( xxx, fnval(pp,xxx), ''c'');', ... ' >> axis([tau(1) tau(5) .6 2.2]);'); ssdisp(figNumber,str); plot(xxx,fnval(pp2,xxx),'-',tau,y,'ro',xxx,fnval(pp,xxx),'c'); axis([tau(1) tau(5) .6 2.2]); hold off % End of the demo ========================