csapidem.html

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

dp0 = fnder( fnbrk(pp0,1) );
</pre><p>Then we compute  lambda(s)  for both  s = PP1  and  s = PP0  ...</p><pre> lam1 := a*fnval(dp1,e) + b*fnval(fnder(dp1),e);
 lam0 := a*fnval(dp0,e) + b*fnval(fnder(dp0),e);</pre><p>... and construct the right linear combination of PP1 and PP0 to get a cubic spline</p><pre>pp  :=  pp1 + ( (c - lambda(pp1))/lambda(pp0) )* pp0</pre><p>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.
         </p><pre class="codeinput">lam1 = a*fnval(dp1,e) + b*fnval(fnder(dp1),e);
lam0 = a*fnval(dp0,e) + b*fnval(fnder(dp0),e);

pp = fncmb(pp0,(c-lam1)/lam0,pp1);
</pre><p>The plot shows that PP fits our quartic polynomial slightly better near E than does the interpolant PP1 with the default conditions
            ...
         </p><pre class="codeinput">xx = (-.3):.05:.7; yy = q(xx);
plot(xx, fnval(pp1,xx) - yy, <span class="string">'x'</span>)
hold <span class="string">on</span>
plot(xx, fnval(pp,xx) - yy, <span class="string">'o'</span>)
title(<span class="string">'Default (x) vs. nonstardard (o)'</span>)
hold <span class="string">off</span>
</pre><img vspace="5" hspace="5" src="csapidem_22.png"> <p>If we also want to enforce the condition</p><pre>  mu(s) := (D^3 s)(3) = 14.6</pre><p>(which our quartic also satisfies), 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, as
         </p><pre>  pp2 = csape(x,[0,zeros(size(y)),1],[0,1]);</pre><p>and solve the linear system for the coefficients D0 and D1 in the linear combination    pp :=  pp1 + d0 * pp0  + d2 * pp2
            for which  lambda(pp) = c, mu(pp) = 14.6 . (Note that both PP0 and PP2 vanish at all interpolation sites, hence PP will match
            the given data for any choice of  D0  and  D2 ).
         </p><pre class="codeinput">cla, title(<span class="string">''</span>), axis <span class="string">off</span>

pp2 = csape(x,[0,zeros(size(y)),1],[0,1]);
</pre><img vspace="5" hspace="5" src="csapidem_23.png"> <p>For amusement, we use the MATLAB encoding facility to write a loop for computing  lambda(PPj), and mu(PPj), for j=0:2</p><pre class="codeinput">dd = zeros(2,3);
<span class="keyword">for</span> j=0:2
   J = num2str(j);
   eval([<span class="string">'dpp'</span>,J,<span class="string">'=fnder(pp'</span>,J,<span class="string">');'</span>]);
   eval([<span class="string">'ddpp'</span>,J,<span class="string">'=fnder(dpp'</span>,J,<span class="string">');'</span>]);
   eval([<span class="string">'dd(1,1+'</span>,J,<span class="string">')=a*fnval(dpp'</span>,J,<span class="string">',e)+b*fnval(ddpp'</span>,J,<span class="string">',e);'</span>]);
   eval([<span class="string">'dd(2,1+'</span>,J,<span class="string">')=fnval(fnder(ddpp'</span>,J,<span class="string">'),3);'</span>]);
<span class="keyword">end</span>

d = dd(:,[1,3])\([c;14.6]-dd(:,2));
pp = fncmb(fncmb(pp0,d(1),pp2,d(2)),pp1);

xxx = 0:.05:3;
yyy = q(xxx);
plot(xxx, yyy - fnval(pp,xxx),<span class="string">'x'</span>)
title(<span class="string">'Error in spline interpolant to y = x.*(-1 + x.*(-1+x.*x/5))'</span>)
</pre><img vspace="5" hspace="5" src="csapidem_24.png"> <p>For reassurance, we compare this error with the one obtained in complete cubic spline interpolation to this function:</p><pre class="codeinput">hold <span class="string">on</span>
plot(xxx, yyy - fnval(csape(x,[-1,y,-7+(4/5)*27],<span class="string">'clamped'</span>),xxx),<span class="string">'o'</span>)
title(<span class="string">'Nonstandard (x) vs endslopes (o)'</span>)
hold <span class="string">off</span>
</pre><img vspace="5" hspace="5" src="csapidem_25.png"> <p>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.
         </p><pre class="codeinput"><span class="comment">% As a final check, here is the third derivative of PP at 3 (which</span>
<span class="comment">% should be  14.6 ):</span>

fnval(fnder(pp,3),3)
</pre><pre class="codeoutput">
ans =

   14.6000

</pre><p class="footer">Copyright 1987-2005 C. de Boor and The MathWorks, Inc.<br>
            Published with MATLAB&reg; 7.1<br></p>
      </div>
      <!--
##### SOURCE BEGIN #####
%%             Cubic spline interpolation
%
% Illustration of the use of CSAPI and CSAPE.

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

%%  CSAPI
% The spline toolbox 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 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 roundoff).

%% Two data points
%
% Just two data points result in a straight line interpolant:

xx = linspace(0,6,121); x = [0 1]; y = [2 0];
plot(xx,csapi(x,y,xx),'k-',x,y,'ro'), grid off
title('Interpolant to two points')

%% Three data points
%
% Three data points give a parabola:

x = [2 3 5]; y = [1 0 4];
plot(xx,csapi(x,y,xx),'k-',x,y,'ro'), grid off
title('Interpolant to three points')

%% Four or more points
% Four or more data points give in general 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'), grid off
title('Cubic spline interpolant to five points')

%% How to check CSAPI
% These look like nice interpolants, but how do we check that CSAPI
% performs as advertised?
%
% We already saw that we interpolate, for we always 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 gives us
% back that cubic spline.

%% The truncated power
% One such is the truncated third power, i.e., the function
% x |-> subplus(x - xi)^3 , with XI one of the breaks and SUBPLUS
% the  t r u n c a t i o n  f u n c t i o n, provided by the spline toolbox
% command SUBPLUS:
cla, title(''), axis off

help subplus

%% A particular truncated power
% For the particular choice XI = 2 , we get the function plotted above.
% 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), grid off
axis([0,6,-10,70])

%% The particular truncated power interpolated
% 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 subplus(x-2)^3')

%% The error
% When comparing two functions, it is usually much more informative to plot
% their difference. To put the size of their difference into context,
% we also compute the maximum data value.
% This shows the error to be no worse than the inevitable round-off.

plot(xx,values-subplus(xx-2).^3)
title('Error in cubic spline interpolation to subplus(x-2)^3')

max_y=max(abs(y))

%% A truncated power that cannot be reproduced
% As a further test, we interpolate a truncated power whose CSAPI 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 the interpolant has three continuous derivatives there.
% Hence we should not be able to reproduce the truncated power centered at
% that site. We try

values = csapi(x,subplus(x-1).^3,xx);

plot(xx,values-subplus(xx-1).^3)
title('Error in not-a-knot interpolant to subplus(x-1)^3')
xlabel(['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.

%% 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, in PP, the ppform of the interpolant. This form can be
% evaluated at some points XX by
%
%         values = fnval(pp,xx)
%
% It can be differentiated by
%
%        dpp = fnder(pp)
%
% or integrated by
%
%        ipp = fnint(pp)
%
% which return, in DPP or IPP, the ppform of the derivative or the primitive,
% respectively.

%% Differentiating the interpolant
% For example, we plot the derivative of this truncated power, i.e.,
% x |-> 3*subplus(x-2)^2 (again in yellow) and, on top of it,
% the derivative of our interpolant (again in black).

pp = csapi(x,subplus(x-2).^3); dpp = fnder(pp);
plot(xx,fnval(dpp,xx),'y','linewidth',3), hold on
plot(xx,3*subplus(xx-2).^2,'k'), hold off, grid off
title('Derivative of interpolant to subplus(x-2)^3')

%% Error
% Again, the more informative comparison is to plot their difference, and this
% is again no bigger than round-off.

plot(xx,fnval(dpp,xx)-3*subplus(xx-2).^2)
title('(derivative of interpolant to subplus(x-2)^3 )  - 3*subplus(x-2)^2')

%% Second derivative
%
% The second derivative of the interpolant should be  6*subplus(.-2).
% We try it, plotting above the difference between this function and the second
% derivative of the interpolant to subplus(.-2)^3.
% Now there are jumps, but they are still within roundoff.

ddpp = fnder(dpp); plot(xx,fnval(ddpp,xx)-6*subplus(xx-2))
title('(second derivative of interpolant to subplus(x-2)^3 )  - 6*subplus(x-2)')

%% Integral
%
% The integral of subplus(.-2)^3 is subplus(.-2)^4/4.
% We plot the difference between it and the integral of the interpolant to
% subplus(.-2)^3:

ipp = fnint(pp);
plot(xx,fnval(ipp,xx)-subplus(xx-2).^4/4)
title('(integral of interpolant to subplus(x-2)^3 ) - subplus(x-2)^4/4')

%% CSAPE