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

    %% Cubic Spline Interpolation
%
% This example shows how to use the |csapi| and |csape| commands from Curve
% Fitting Toolbox(TM) to construct cubic spline interpolants.

%   Copyright 1987-2012 The MathWorks, Inc.

%%  The CSAPI Command
% The command
%
%     values = csapi(x,y,xx)
%
% returns the values at |xx| of the cubic spline interpolant to the given data
% (|x,y|), using the not-a-knot end condition. This interpolant is a
% piecewise cubic function, with break sequence |x|, whose cubic pieces join
% together to form a function with two continuous derivatives. The
% "not-a-knot" end condition means that, at the first and last interior break,
% even the third derivative is continuous (up to round-off error).

%%
% Specifying only two data points results in a straight line interpolant.
x = [0 1];
y = [2 0];
xx = linspace(0,6,121);
plot(xx,csapi(x,y,xx),'k-',x,y,'ro')
title('Interpolant to Two Points')

%%
% Specifying three data points gives a parabola.
x = [2 3 5];
y = [1 0 4];
plot(xx,csapi(x,y,xx),'k-',x,y,'ro')
title('Interpolant to Three Points')

%%
% More generally, four or more data points give a cubic spline.
x = [1 1.5 2 4.1 5];
y = [1 -1 1 -1 1];
plot(xx,csapi(x,y,xx),'k-',x,y,'ro')
title('Cubic Spline Interpolant to Five Points')

%% How to Check the Output of CSAPI
% These look like nice interpolants, but how do we check that |csapi|
% performs as advertised?
%
% We already saw that |csapi| interpolates, because we plotted the data points
% and the interpolant went right through those points. But to be sure that we
% get a cubic spline, it is best to start with data from a cubic spline of the
% expected sort and check whether |csapi| _reproduces_ that cubic spline, i.e.,
% gives back that cubic spline from which the data were taken.

%% Example: The Truncated Power Function
% One simple example of a cubic spline function to check against is the
% truncated third power, i.e., the function
%
% $$f(x) = ((x - xi)_{+})^3,$$
%
% where |xi| is one of the breaks and the "+" subscript indicates the
% _truncation function_, provided by the command |subplus|:
help subplus

%%
% The truncated 3rd power is plotted below for the particular choice |xi| = |2|.
% As expected, it is zero to the left of 2, and rises like (x-2)^3 to the right
% of 2.
plot(xx, subplus(xx-2).^3,'y','LineWidth',3)
axis([0,6,-10,70])

%%
% Now we interpolate this particular cubic spline at the data sites 0:6, and
% plot the interpolant on top of the spline, in black.
x = 0:6;
y = subplus(x-2).^3;
values = csapi(x,y,xx);
hold on
plot(xx,values,'k',x,y,'ro')
hold off
title('Interpolant to ((x-2)_+)^3')

%% The Interpolation Error
% When comparing two functions, it is usually much more informative to plot
% their difference.
plot(xx, values - subplus(xx-2).^3)
title('Error in Cubic Spline Interpolation to ((x-2)_+)^3')
%%
% To put the size of their difference into context, you can also compute the
% maximum data value. This shows the error to be no worse than the inevitable
% round-off.
max_y = max(abs(y))

%% A Truncated Power That Cannot be Reproduced
% As a further test, we interpolate a truncated power whose |csapi|-produced
% interpolant at the sites 0:6 cannot coincide with it. For example, the first
% interior break of the interpolating spline is not really a knot since |csapi|
% uses the "not-a-knot" condition, hence the interpolant has three continuous 
% derivatives at that site. This implies that we should not be able to reproduce
% the truncated 3rd power centered at that site since its third derivative
% is discontinuous across that site.
values = csapi(x,subplus(x-1).^3,xx);
plot(xx, values - subplus(xx-1).^3)
title('Error in Not-a-Knot Interpolant to ((x-1)_+)^3')
%%
% Since 1 is a first interior knot, it is not active for this interpolant.
%
% The difference is as large as .18, but decays rapidly as we move away from
% 1. This illustrates that _cubic spline interpolation is essentially local_.

%% Using the ppform Instead of Values
% It is possible to retain the interpolating cubic spline in a form suitable
% for subsequent evaluation, or for calculating its derivatives, or for other
% manipulations. This is done by calling |csapi| in the form
%
%    pp = csapi(x,y)
%
% which returns the ppform of the interpolant. You can evaluate this form at
% some new points |xx| by the command
%
%    values = fnval(pp,xx)
%
% You can differentiate the interpolant by the command
%
%    dpp = fnder(pp)
%
% or integrate it by the command
%
%    ipp = fnint(pp)
%
% which return the ppform of the derivative or the integral, respectively.

%% Example: Differentiating and Integrating the Interpolant
% To show differentiation of an interpolant, we plot the derivative of this
% truncated power
%
% $$f'_2(x) = 3((x - 2)_{+})^2,$$

%%
% (again in yellow) and then, on top of it, the derivative of our interpolant
% to the original truncated third power function (again in black).
plot(xx,3*subplus(xx-2).^2,'y','LineWidth',3)
pp = csapi(x,subplus(x-2).^3);
dpp = fnder(pp);
hold on
plot(xx,fnval(dpp,xx),'k')
hold off
title('Derivative of Interpolant to ((x-2)_+)^3')

%%
% Again, the more informative comparison is to plot their difference, and as
% before this is no bigger than round-off.
plot(xx, fnval(dpp,xx) - 3*subplus(xx-2).^2)
title('Error in Derivative of interpolant to ((x-2)_+)^3')

%%
% The second derivative of the truncated power is
%
% $$f''_2(x) = 6(x-2)_+$$
%
% A plot of the difference between this function and the second derivative of
% the interpolant to the original function shows that there are now jumps, but
% they are still within roundoff.
ddpp = fnder(dpp);
plot(xx, fnval(ddpp,xx) - 6*subplus(xx-2))
title('Error in Second Derivative of Interpolant to ((x-2)_+)^3')

%%
% The integral of the truncated power is
%
% $$F_2(x) = ((x-2)_+)^4/4.$$
%
% A plot of the difference between this function and the integral of the
% interpolant to the original function again shows that the errors are within
% round-off.
ipp = fnint(pp);
plot(xx, fnval(ipp,xx) - subplus(xx-2).^4/4)
title('Error in Integral of Interpolant to ((x-2)_+)^3')

%% The CSAPE Command
%
% Like |csapi|, the |csape| command provides a cubic spline interpolant to
% given data.  However, it permits various additional end conditions. Its
% simplest version,
%
%    pp = csape(x,y)
%
% uses the Lagrange end condition, which is a common alternative to the
% not-a-knot condition used by |csapi|.  |csape| does not directly return
% values of the interpolant, but only its ppform.

%%
% For example, consider again interpolation to the function
%
% $$f_1(x) = ((x-1)_+)^3,$$
%
% which |csapi| fails to reproduce well. We plot the error of the not-a-knot
% interpolant returned by |csapi| (in black), along with the error of the
% interpolant obtained from |csape| (in red).
exact = subplus(xx-1).^3;
plot(xx, fnval(csapi(x,subplus(x-1).^3),xx) - exact,'k')
hold on
plot(xx, fnval(csape(x,subplus(x-1).^3),xx) - exact,'r')
title('Error in Not-a-Knot vs. Lagrange End Conditions')
legend({'Not-a-Knot' 'Lagrange'});
hold off
%%
% There is not much difference between the two interpolants in this case.

%% Other End Conditions: The `Natural' Spline Interpolant
%
% The |csape| command also provides ways to specify several other types of
% end conditions for an interpolating cubic spline. For example, the command
%
%    pp = csape(x,y,'variational')
%
% uses the so-called `natural' end conditions. This means that the
% second derivative is zero at the two extreme breaks.

%%
% This step shows how to apply `natural' cubic spline interpolation to the
% function
%
% $$f_2(x) = ((x-2)_+)^3,$$
%
% and plot the error. The code below computes the `natural' spline interpolant
% with an alternative argument syntax that is equivalent to the 'variational'
% string argument: using the string 'second' specifies that |csape| should set
% the second derivative at the extreme data sites to the default value of 0.
pp = csape(x,subplus(x-2).^3,'second');
plot(xx, fnval(pp,xx) - subplus(xx-2).^3)
title('Error in ''Natural'' Spline Interpolation to ((x-2)_+)^3')
%%
% Note the large error near the right end.  This is due to the fact that the
% 'natural' end conditions implicitly insist on having a zero second
% derivative there.

%% Other End Conditions: Prescribing Second Derivatives
% We can also explicitly use the correct second derivatives to get a small
% error. First, we compute the correct second derivative values of the
% truncated power at the endpoints.
endcond = 6*subplus(x([1 end])-2);

%%
% Then we create the interpolant, specifying that second derivatives at the
% endpoints are to be matched to the second derivative values we just
% computed.  We do this by providing |endcond(1)| for the left endpoint
% condition, and |endcond(2)| for the right, along with the data values.
pp = csape(x,[endcond(1) subplus(x-2).^3 endcond(2)], 'second');
plot(xx, fnval(pp,xx) - subplus(xx-2).^3,'r')
title(['Error in Spline Interpolation to ((x-1)_+)^3'; ...
       '  When Matching the 2nd Derivative at Ends  '])

%% Other End Conditions: Prescribing Slopes
% |csape| also permits specification of endpoint _slopes_. This is the
% _clamped_ (or, _complete_) cubic spline interpolant.
% The statement
%
%         pp = csape(x,[sl,y,sr],'clamped')
%
% creates the cubic spline interpolant to the data (|x|, |y|) that also has
% slope |sl| at the leftmost data site and slope |sr| at the rightmost data
% site.

%% Other End Conditions: Mixed End Conditions
% It is even possible to mix these conditions. For example, our much-exercised
% truncated power function 
%
% $$f_1(x) = ((x-1)_+)^3$$
%
% has slope 0 at |x|=0 and second derivative 30 at |x|=6 (the last data site).
%
% Therefore, by matching the slope at the left end and the curvature at the
% right, we expect no error in the resulting interpolant.
pp = csape(x, [0 subplus(x-1).^3 30], [1 2]);
plot(xx, fnval(pp,xx) - subplus(xx-1).^3)
title(['Error in Spline Interpolation to ((x-1)_+)^3'; ...
       '        with Mixed End Conditions.          '])

%% Other End Conditions: Periodic Conditions
% It is also possible to prescribe _periodic_ end conditions. For example, the
% sine function is 2*pi-periodic and has the values |[0 -1 0 1 0]| at the
% sites |(pi/2)*(-2:2)|. The difference, between the sine function and its
% periodic cubic spline interpolant at these sites, is only 2 percent. Not
% bad.
x = (pi/2)*(-2:2);
y = [0 -1 0 1 0];
pp = csape(x,y, 'periodic' );
xx = linspace(-pi,pi,201);
plot(xx, sin(xx) - fnval(pp,xx), 'x')
title('Error in Periodic Cubic Spline Interpolation to sin(x)')

%% End Conditions Not Explicitly Covered by CSAPI or CSAPE
% Any end condition not covered explicitly by |csapi| or |csape| can be
% handled by constructing the interpolant with the |csape| default side
% conditions, and then adding to it an appropriate scalar multiple of
% an interpolant to zero values and some side conditions. If there are two
% `nonstandard' side conditions to be satisfied, you may have to solve a
% 2-by-2 linear system first.
%
% For example, suppose that you want to compute the cubic spline interpolant
% |s| to the data
x = 0:.25:3;
q = @(x) x.*(-1 + x.*(-1+x.*x/5));
y = q(x);
%%
% and enforce the condition
%
%    lambda(s) := a * (Ds)(e) + b * (D^2 s)(e) = c
%
% on the first and second derivatives of |s| at the point |e|.
%%
% The data were generated from a quartic polynomial that happens to satisfy
% this side condition with specific parameters
e = x(1);
a = 2; b = -3; c = 4;

%%
% To construct the interpolant that satisfies this specific condition, we first
% construct the interpolant with the default end conditions
pp1 = csape(x,y);
%%
% and the first derivative of its first polynomial piece.
dp1 = fnder(fnbrk(pp1,1));
%%
% In addition, we construct the cubic spline interpolant to zero data values,
% specifying that it have a slope of 1 at |e|,
pp0 = csape(x,[1,zeros(size(y)),0], [1,0]);
%%
% as well as constructing the first derivative of its first polynomial piece.
dp0 = fnder(fnbrk(pp0,1));

%%
% Then we compute |lambda| for both |pp1| and |pp0|,
lam1 = a*fnval(dp1,e) + b*fnval(fnder(dp1),e);
lam0 = a*fnval(dp0,e) + b*fnval(fnder(dp0),e);
%%
% and construct the correct linear combination of |pp1| and |pp0| to get a
% cubic spline
%
%    s := pp1 + ((c - lambda(pp1))/lambda(pp0)) * pp0
%
% that does satisfy the desired condition, as well as the default end
% condition at the right endpoint. We form this linear combination with
% the help of |fncmb|.
s = fncmb(pp0,(c-lam1)/lam0,pp1);

%%
% A plot of the interpolation error shows that |s| fits the quartic polynomial
% slightly better near |e| than the interpolant |pp1| with the default
% conditions does.
xx = (-.3):.05:.7; yy = q(xx);
plot(xx, fnval(pp1,xx) - yy, 'x')
hold on
plot(xx, fnval(s,xx) - yy, 'o')
hold off
legend({'Default conditions' 'Nonstandard conditions'},'location','SE')

%%
% If we want to enforce the condition
%
%    mu(s) := (D^3 s)(3) = 14.6
%
% on the third derivative of the interpolant (the quartic satisfies this
% condition), then we construct an additional cubic spline interpolating
% to zero values, and with zero first derivative at the left endpoint, hence
% certain to be independent from |pp0|.
pp2 = csape(x,[0,zeros(size(y)),1],[0,1]);
%%
% Then we find the coefficients |d0| and |d2| in the linear combination
%
%    s := pp1 + d0*pp0 + d2*pp2
%
% that solves the linear system 
%
%    lambda(s) = c
%    mu(s) = 14.6
%
% Note that both |pp0| and |pp2| vanish at all interpolation sites, hence |s|
% will match the given data for any choice of |d0| and |d2|.

%%
% For amusement, we use the MATLAB(R) encoding facility to write a loop to
% compute |lambda(pp_j)| and |mu(pp_j)|, for |j|=0:2.
dd = zeros(2,3);
for j=0:2
   J = num2str(j);
   eval(['dpp',J,'=fnder(pp',J,');']);
   eval(['ddpp',J,'=fnder(dpp',J,');']);
   eval(['dd(1,1+',J,')=a*fnval(dpp',J,',e)+b*fnval(ddpp',J,',e);']);
   eval(['dd(2,1+',J,')=fnval(fnder(ddpp',J,'),3);']);
end

%%
% Given the values of |lambda| and |mu| for |pp0|, |pp1|, and |pp2|, we then
% solve for the coefficients that define the correct linear combination.
d = dd(:,[1,3])\([c;14.6]-dd(:,2));
s = fncmb(fncmb(pp0,d(1),pp2,d(2)),pp1);

xxx = 0:.05:3;
yyy = q(xxx);
plot(xxx, yyy - fnval(s,xxx),'x')
title('Error in Spline Interpolant to y = x*(-1 + x*(-1+x*x/5))')

%%
% For reassurance, we compare this error with the one obtained in complete cubic
% spline interpolation to this function:
hold on
plot(xxx, yyy - fnval(csape(x,[-1,y,-7+(4/5)*27],'clamped'),xxx),'o')
hold off
legend({'Nonstandard conditions' 'Endslope conditions'})

%%
% The errors differ (and not by much) only near the end points, testifying
% to the fact that both |pp0| and |pp2| are sizable only near their respective
% end points.

%%
% As a final check, we verify that |s| satisfies the desired third derivative
% condition at 3.
fnval(fnder(s,3),3)


displayEndOfDemoMessage(mfilename)