www.gusucode.com > curvefit 案例源码程序 matlab代码 > curvefit/spapidem.m

    %% Fitting a Spline to Titanium Test Data
%
% This example shows how to use commands from Curve Fitting Toolbox(TM) to fit a
% spline to titanium test data with manual and automatic selection of knots.

% Copyright 1987-2012 The MathWorks, Inc.

%% Manual Knot Choice for Spline Interpolation
% Here are some data that record a certain property of titanium, measured as a
% function of temperature.  We'll use it to illustrate some issues with spline
% interpolation.
[xx,yy] = titanium;
%%
% A plot of the data shows a rather sharp peak.
plot(xx,yy,'bx');
frame = [-10 10 -.1 .3]+[min(xx),max(xx),min(yy),max(yy)];
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.
pick = [1 5 11 21 27 29 31 33 35 40 45 49];
tau = xx(pick);
y = yy(pick);
hold on
plot(tau,y,'ro');
hold off

%%
% 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 add 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]];  % construct the knot sequence
plot(tau,y,'ro');
hold on
axis(frame+[-2*dl 2*dr 0 0])
plot(t,repmat(frame(3)+.03,size(t)),'kx')
hold off
legend({'Data Values' 'Knots'},'location','NW')

%%
% We use this knot sequence to construct an interpolating cubic spline.
sp = spapi(t,tau,y);

%%
% 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.
plot(tau,y,'ro')
axis(frame)
hold on
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.  To see that, here is a
% picture of the entire spline, along with its knot sequence and the data
% points.
fnplt(sp,'k');
hold on
plot(tau,y,'ro', t,repmat(.1,size(t)),'kx');
hold off
legend({'Spline Interpolant' 'Data Values' 'Knots'},'location','NW')

%%
% 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);
plot(tau,y,'ro', xx([1 2 end-1 end]),yy([1 2 end-1 end]),'bo');
axis(frame+[-2*dl 2*dr 0 0]);
hold on
fnplt(sp2,'b',tau([1 end]))
hold off
legend({'Original Data' 'Data Added for End Conditions' ...
        'Fit with Added Data'},'location','NW')

%%
% Here is a comparison of the two spline fits, to show the reduction of the
% undulation in the first and last interval.
hold on
fnplt(sp,'k',tau([1 end]))
hold off
legend({'Original Data' 'Data Added for End Conditions' ...
        'Fit with Added Data' 'Original Fit'},'location','NW')

%%
% Finally, here is a closer look at the first four data intervals
% that shows more clearly the reduction of the undulation near
% the left end.
plot(tau,y,'ro',xxx,fnval(sp2,xxx),'b',xxx,fnval(sp,xxx),'k');
axis([tau(1) tau(5) .6 1.2]);
legend({'Original Data' 'Fit with Added Data' ...
        'Original Fit'},'location','NW')

%% Automatic Knot Choice for Interpolation
% If all this detail turns you off, let Curve Fitting Toolbox choose the
% knots for you. Specify the desired order of the interpolant as the first
% input argument to the spline interpolation command |spapi|, rather than a
% knot sequence.
autosp = spapi(4, tau, y);
knots = fnbrk(autosp,'knots');
plot(tau, y, 'ro')
hold on
fnplt(autosp,'g')
plot(knots, repmat(.5,size(knots)),'gx')
hold off
legend({'Data Values' 'Fit With Knots Chosen by SPAPI' ...
        'Knots Chosen by SPAPI'}, 'location','NW')

%%
% Below is the result of a much better knot choice, obtained by shifting the
% knot at 842 slightly to the right and the knot at 985 slightly to the left.
knots([7 12]) = [851, 971];
adjsp = spapi(knots, tau, y);
hold on
fnplt(adjsp,'r',2)
plot(knots, repmat(.54,size(knots)),'rx')
hold off
legend({'Data Values' 'Fit With Knots Chosen by SPAPI' ...
        'Knots Chosen by SPAPI' 'Fit With Knots Adjusted' ...
        'Adjusted Knots'}, 'location','NW')

%%
% Or else, simply try the standard cubic spline interpolant,
% supplied by |csapi|.  This amounts to a slightly different choice of knots.
autocs = csapi(tau, y);
plot(tau, y, 'ro')
hold on
fnplt(autocs,'c')
hold off

%%
% With such rapidly-varying data, it is hard to get agreement among all
% reasonable interpolants, even if each of them is a cubic spline.  The plot
% below shows all five interpolants, for comparison.
plot(tau, y, 'ro')
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
legend({'Data Values' 'Original Fit' 'Special End Conditions' ...
        'With Knots Chosen by SPAPI' 'With Knots Chosen by CSAPI' ...
        'With Adjusted Knots'},'location','NW')

displayEndOfDemoMessage(mfilename)