csapidem.html

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

%
% The command CSAPE also provides a cubic spline interpolant to given data,
% but permits various other end conditions. It does not directly return values
% of that interpolant, but only its ppform.  Its simplest version,
%
%       pp = csape(x,y)
%
% uses the Lagrange end conditions, which are better at times than the
% not-a-knot condition used by CSAPI.

cla, title(''), axis off
%% csapi(x,y) vs csape(x,y)
% For example, consider again interpolation to the truncated power which
% CSAPI fails to reproduce. We plot here (in black), the error of the
% not-a-knot interpolant, along with the error (in red) of the interpolant
% obtained from CSAPE; REPLACE_WITH_DASH_DASH not much difference between the two in this case.

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 (black) vs. Lagrange (red)')
hold off

%% Other end conditions: the `natural' spline interpolant
%
% The command CSAPE only provides the interpolating cubic spline in
% ppform, but for several different end conditions. For example, the call
%
%         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.

cla, title(''), axis off
%% The `natural' spline interpolant: an example
% We apply `natural' cubic spline interpolation to the truncated power we
% were able to reproduce, and plot the error. We now get a large error near
% the right end, due to the fact that we insisted on a zero second derivative
% there. For variety, we get the `natural' spline interpolant here by the
% command  csape(x,y,'second') which uses the default value 0 for the second
% derivative at the extreme data sites.

pp = csape(x,subplus(x-2).^3,'second');
plot(xx, fnval(pp,xx) - subplus(xx-2).^3 )
title('Error in ''natural'' spline interpolation to subplus(x-2)^3')
xlabel(' note the large error near the right end')

%% Other end conditions: prescribing second derivatives
% When we use explicitly the correct second derivatives, we get a small error,
% as the above plot shows.
% The command CSAPE(X,[VALS(1),Y VALS(2)],'second') specifies that second
% derivatives are to be matched, with VALS(1) the second derivative at the
% left endpoint, and VALS(2) the second derivative at the right.
% We compute these values directly from the second derivative of the truncated
% power subplus(.-2)^3.

pp = csape(x,[6*subplus(x(1)-2),subplus(x-2).^3,6*subplus(x(end)-2)], 'second');
plot(xx, fnval( pp, xx ) - subplus(xx-2).^3,'r')
title('Error in spline interpolation to subplus(x-1)^3 when matching f'''' at ends ')

%% Other end conditions: prescribing slopes
% CSAPE also permits specification of endpoint  s l o p e s . This is the
%  c l a m p e d  spline  (or,  c o m p l e t e  cubic spline interpolant).
% The statement
%
%         pp = csape(x,[sl,y,sr],'clamped')
%
% supplies 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.

cla, title(''), axis off
%% Other end conditions: mixed end conditions
% It is even possible to mix these conditions. For example, our much
% exercised truncated power function  f(x) = subplus(x-1)^3  has slope  0
% at  x = 0   and second derivative  30  at  x = 6 (our last data site).
%
% We therefore expect no error in the following 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 subplus(x-1)^3 when matching  ...')
xlabel(' ... slope at left end and curvature at right.')

%% Other end conditions: periodic conditions
% It is also possible to prescribe  p e r i o d i c  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, plotted above, 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, one may have to solve a
% 2x2 linear system first.
%
% For example, suppose that you want to enforce the condition
%
%         lambda(s)  :=  a (Ds)(e)  + b (D^2 s)(e) = c
%
% on the cubic spline interpolant  s  to the following data (taken from a
% a quartic polynomial that happens to satisfy this specific side condition):

cla, title(''), axis off
q = inline('x.*(-1 + x.*(-1+x.*x/5))');
x = 0:.25:3; y = q(x);

e = x(1); a = 2; b = -3; c = 4;

%%
% Then, in addition to the interpolant PP1 with the default end conditions ...
%  pp1 = csape(x,y);
% ... and first derivative DP1 of its first polynomial piece ...
%  dp1 = fnder( fnbrk(pp1,1) );
% ... we also construct the cubic spline interpolant PP0 to zero data,
% and with a slope of  1  at E, as well as the first derivative DP0 of its
% first polynomial piece.

pp1 = csape(x,y);
dp1 = fnder(fnbrk(pp1,1));

pp0 = csape(x,[1,zeros(size(y)),0], [1,0] );
dp0 = fnder( fnbrk(pp0,1) );

%%
% Then we compute  lambda(s)  for both  s = PP1  and  s = PP0  ...
%
%   lam1 := a*fnval(dp1,e) + b*fnval(fnder(dp1),e);
%   lam0 := a*fnval(dp0,e) + b*fnval(fnder(dp0),e);
%
% ... and construct the right linear combination of PP1 and
% PP0 to get a cubic spline
%
%  pp  :=  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.

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

%%
% The plot shows that PP fits our quartic polynomial slightly better near E
% than does the interpolant PP1 with the default conditions ...

xx = (-.3):.05:.7; yy = q(xx);
plot(xx, fnval(pp1,xx) - yy, 'x')
hold on
plot(xx, fnval(pp,xx) - yy, 'o')
title('Default (x) vs. nonstardard (o)')
hold off

%%
% If we also want to enforce the condition
%
%    mu(s) := (D^3 s)(3) = 14.6
%
% (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
%
%    pp2 = csape(x,[0,zeros(size(y)),1],[0,1]);
%
% 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 ).
cla, title(''), axis off

pp2 = csape(x,[0,zeros(size(y)),1],[0,1]);

%%
% For amusement, we use the MATLAB encoding facility to write a loop
% for computing  lambda(PPj), and mu(PPj), 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

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),'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')
title('Nonstandard (x) vs endslopes (o)')
hold off

%%
% 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, here is the third derivative of PP at 3 (which
% should be  14.6 ):

fnval(fnder(pp,3),3)


displayEndOfDemoMessage(mfilename)

##### SOURCE END #####
-->
   </body>
</html>