spapidem.m

演示matlab曲线拟和与插直的基本方法

%% Splines and real world data
% This a demonstration of a spline function being fit to real-world data.

% Copyright 1987-2005 C. de Boor and The MathWorks, Inc.
% $Revision: 5.24.4.2 $

%% Manual knot choice for spline interpolation
% Here are some data which record a property of
% titanium measured as a function of temperature. We'll use it to
% illustrate some issues with spline interpolation.
%
% The plot of the data shows a rather sharp peak.
hold off
[xx,yy]=titanium;
frame=[-1 1 -.1 .1]+[min(xx),max(xx),min(yy),max(yy)];
plot(xx,yy,'x');
axis(frame);

%%
% We pick a few data points from these somewhat
% rough data, since we want to interpolate. Here is a
% picture of the data, with the selected data points marked.

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');

%%
% 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 site lies in the support
% of the i-th B-spline. We achieve this by using the data
% sites as knots, but adding two simple knots at either end.
%
%   dl = tau(2)-tau(1); dr = tau(end)-tau(end-1);
%   t = [tau(1)-dl*[2 1] tau tau(end)+dr*[1 2]];
%
% We use this knot sequence to construct an interpolating cubic spline:
%
%   sp = spapi(t, tau, y);

dl=tau(2)-tau(1);
dr=tau(end)-tau(end-1);
t=[tau(1)-dl*[2 1] tau tau(end)+dr*[1 2]];   % construct the knot sequence
hold on
axis(frame+[-2*dl 2*dr 0 0])
plot(t,repmat(frame(3)+.03,size(t)),'kx')
hold off
sp=spapi(t,tau,y);                         % This constructs the spline.

%%
% Now, for the plot.
% Since we do not care about the part of the spline
% outside the data interval, we restrict the plot to that interval:
%
%  hold on,  fnplt(sp,[tau(1) tau(end)], 'k'),  hold off
%

plot(xx,yy,'x',tau,y,'ro'), axis(frame), hold on
% Now, for the plot:
fnplt(sp,[tau(1) tau(end)], 'k')
hold off

%%
% A closer look at the left part of the spline fit shows some
% undulations.

xxx = linspace(tau(1),tau(5),41);
plot(xxx, fnval(sp, xxx), 'k', tau, y, 'ro');
axis([tau(1) tau(5) 0.6 1.2]);

%%
% 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, along with its knot sequence.

axis([tau(1) tau(5) 0.6 1.2]);
fnplt(sp,'k');
hold on, plot(t,repmat(.1,size(t)),'kx'), hold off

%%
% Here are again the data points as well.

hold on
plot(tau, y, 'ro');
hold off

%%
% 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(end)+[1 2 3 4]*dr];
xx=[tau(1)-[2 1]*dl tau tau(end)+[1 2]*dr];
yy=[y(1)-[2 1]*(y(2)-y(1)) y y(end)+[1 2]*(y(end)-y(end-1))];
sp2=spapi(tt,xx,yy); fnplt(sp2,'b',tau([1 end]))
hold on
plot(tau,y,'or', xx([1 2 end-1 end]),yy([1 2 end-1 end]),'ko');
axis(frame+[-2*dl 2*dr 0 0]);
hold off

%%
% Here is also the original spline fit, shown in black, to show the
% reduction of the undulation in the first and last interval.

hold on
fnplt(sp,'k',tau([1 end]))
hold off

%%
% 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.

plot(xxx,fnval(sp2,xxx),'b',tau,y,'ro',xxx,fnval(sp,xxx),'k');
axis([tau(1) tau(5) .6 2.2]);
hold off

%% Automatic knot choice for interpolation
% If all this detail turns you off, let the spline toolbox choose the knots
% for you, by using the spline interpolation command SPAPI in the form
%        spapi( order, data_sites, data_values )

autosp = spapi( 4, tau, y);
fnplt(autosp,'g')
hold on, plot(tau, y, 'or'), hold off

%%
% Here is the result of a much better knot choice (plotted above as red x's),
% obtained by shifting the knot at 842 slightly to the right
% and the knot at 985 slightly to the left.

knots = fnbrk(autosp,'knots');
hold on, plot(knots, repmat(.5,size(knots)),'xg')
knots([7 12]) = [851, 971];
plot(knots, repmat(.54,size(knots)),'xr')
adjsp = spapi(knots, tau, y);
fnplt(adjsp,'r',2), hold off

%%
% Else, simply try the standard cubic spline interpolant,
% supplied by CSAPI (which amounts to a slightly different choice of knots):
%
% The next and last slide shows all five interpolants, for comparison.

autocs = csapi(tau, y);
fnplt(autocs,'c')
hold on, plot(tau, y, 'or'), hold off

%%
% With such fast varying data, it is hard to get agreement among all reasonable
% interpolants, even if each of them is a cubic spline.

hold on
fnplt(sp,'k',tau([1 end]))  % black: original
fnplt(sp2,'b',tau([1 end])) % blue:  with special end conditions
fnplt(autosp,'g')           % green: automatic knot choice by SPAPI
fnplt(autocs,'c')           % cyan:  automatic knot choice by CSAPI
fnplt(adjsp,'r',2)          % red:   knot choice by SPAPI slightly changed
hold off


displayEndOfDemoMessage(mfilename)