echo off
%CSAPIDEM Cubic spline interpolation demo.

%   Copyright 1987-2000 by C. de Boor and The MathWorks, Inc.
%   cb: March 22, 1991; 27 mar 92 (periodic); 13 apr 92 (other end conditions)
%   21 jun 92 (replace zeros(y)); 28 aug 94 (replace pppce by ppbrk); 28sep96
%   11apr97 (correct text), 8aug97
%   $Revision: 1.10 $

home, close, echo on
%     The M-file  csapi  constructs the cubic spline interpolant to given
% data, using the "not-a-knot" end condition. The call

%     values = csapi(x,y,xx)

% returns the values at  xx  of that interpolating cubic spline.

%     With  xi  the data sites in strictly increasing order, the cubic
% spline interpolant is a piecewise cubic with break sequence  xi  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).

%      Here are some trial runs.

pause                                             %Touch any key to continue
%     Just two data points result in a straight line interpolant:

pause                                             %Touch any key to continue
xx=[0:60]/10;x=[0 1];y=[2 0];
plot(xx,csapi(x,y,xx),'k-',x,y,'ro'), grid off
title('interpolant to two points'), pause

%     Three data points give a parabola:

pause                                             %Touch any key to continue
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'), pause

%     Four or more data points give in general a cubic spline.

pause                                             %Touch any key to continue

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'), pause

pause                                             %Touch any key to continue

%     These look like nice interpolants, but how do we check it?

% We already saw that we interpolate, for we always plotted the data points
% and the interpolant went right through those points.

%     But, for checking purposes, it is better to choose specific data taken
% from a cubic spline so that the entire interpolant should coincide with
% that function.

pause                                             %Touch any key to continue
%     One such is the truncated third power, i.e., the function
%  x ]-> subplus(x - knot)^3 , with  knot  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 , i.e.,
help subplus
pause                                             %Touch any key to continue


% We try it with  knot = 2:
x=[0:6]; y=subplus(x-2).^3;
% This is data taken from a cubic spline with just one break, at 2.

values = csapi(x,y,xx);

% we plot the interpolant:

pause                                             %Touch any key to continue
plot(xx,values,'k',x,y,'ro'), grid off
title('interpolant to subplus(x-2)^3'), pause
% The function is zero to the left of  2  and rises like (x-2)^3 to the
% right of 2. To compare it with subplus(.-2)^3, we plot the difference:

pause                                             %Touch any key to continue
plot(xx,values-subplus(xx-2).^3)
title('error in cubic spline interpolation to subplus(x-2)^3'),pause
% Since the maximum of the given function values is
max_y=max(abs(y))
% the size of the absolute error (as indicated by the plot) is acceptably low.

pause                                             %Touch any key to continue

%   As a further test, we try to interpolate to a truncated power which
% cannot be interpolated exactly. 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 point. We try

values = csapi(x,subplus(x-1).^3,xx);
pause                                             %Touch any key to continue
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'])
pause
% The difference is as large as .2, but decays rapidly as we move away from
% 1. This illustrates that cubic spline interpolation is essentially local.

pause                                             %Touch any key to continue

%     It is possible to retain the interpolating cubic spline in a form
% suitable for later, repeated 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 primitive.
pause                                             %Touch any key to continue
%     For example, we try again the interpolant to the truncated power.
pp = csapi(x,subplus(x-2).^3);
% and look at its derivatives:
pause                                             %Touch any key to continue
dpp = fnder(pp);plot(xx,fnval(dpp,xx)), grid off
title('derivative of interpolant to subplus(x-2)^3'), pause
% This should coincide with  3*subplus(.-2).^3. We look at the error:
pause                                             %Touch any key to continue
plot(xx,fnval(dpp,xx)-3*subplus(xx-2).^2)
title('(derivative of interpolant to subplus(x-2)^3 )  - 3*subplus(x-2)^2')
pause
% and that looks ok.

%     The second derivative of the interpolant should be  6*subplus(.-2).
% We try it.
pause                                             %Touch any key to continue
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)')
pause

% There are jumps, but they are within roundoff.

pause                                             %Touch any key to continue

% Here is also the integral
ipp = fnint(pp);plot(xx,fnval(ipp,xx)), grid off
title('integral of interpolant to subplus(x-2)^3'), pause
% and here is the difference between the constructed and the ideal integral:
plot(xx,fnval(ipp,xx)-subplus(xx-2).^4/4)
title('(integral of interpolant to subplus(x-2)^3 ) - subplus(x-2)^4/4'), pause


%  The M-file  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 .

pause                                             %Touch any key to continue

% For example, here is again the interpolant to the truncated power which
%  csapi  fails to reproduce:

values = csapi(x,subplus(x-1).^3,xx);
pause                                             %Touch any key to continue

plot(xx,values-subplus(xx-1).^3)
title('error in spline interpolation to subplus(x-1)^3'), pause
hold on

% For comparison, here is the error when we use the Lagrange end conditions:
pause                                             %Touch any key to continue

plot(xx, fnval( csape(x,subplus(x-1).^3), xx ) - subplus(xx-1).^3, '-.' )
title('not-a-knot (-) vs. Lagrange (-.)'),pause
hold off

% Well, in this case, there isn't much difference.
pause                                             %Touch any key to continue

%     The M-file  csape  only provides the interpolating cubic spline in
% ppform, but for several different end conditions. For example, the call
%
%         pp = csape(x,y,'second')
%
% uses the so-called `natural' end conditions. This means that the 
% second derivative is zero at the two extreme breaks. This is specified
% here by asking that second derivatives be matched at the endpoints, but then
% letting the default value, 0, be used.

%   We apply `natural' cubic spline interpolation to the truncated power, and
% plot the error:

pp = csape(x,subplus(x-1).^3,'second');
pause                                             %Touch any key to continue

plot(xx, fnval(pp,xx) - subplus(xx-1).^3 )
title('error in ''natural'' spline interpolation to subplus(x-1)^3')
xlabel(' note the large error near the end'),pause

% Not surprisingly, the error is not small, for the simple reason that we are
% interpolating to a cubic spline function whose second derivative happens not
% to be zero at the right endpoint.

pause                                             %Touch any key to continue

% When we use explicitly the correct second derivatives, we get a small error:

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


pause                                             %Touch any key to continue
%  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,y,'clamped',[sl,sr])
%
% supplies the cubic spline interpolant to the data  x , y  that also
% has slope  sl  at the first data point and slope  sr  at the last data
% point.
pause                                             %Touch any key to continue
%   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,subplus(x-1).^3, [1 2], [0,30] );
%                            ^ ^    ^ ^^
%                            | |    | ||

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.')
pause

%   ... and we were not disappointed.

pause                                             %Touch any key to continue

%    It is also possible to prescribe  p e r i o d i c  end conditions.
% For example, a periodic interpolant to the values
y = [0 -1 0 1 0];
% at the sites
x = (pi/2)*[-2:2];
% should be a good approximation to the sine function:
pp = csape(x,y, 'periodic' );

%   Here is a plot of the difference between the sine function and this
% periodic piecewise cubic function, ...
xx = linspace(-pi,pi,101);
plot(xx, sin(xx) - fnval(pp,xx), 'x')
title('error in periodic cubic spline interpolation to sin(x)'), pause

%   ... showing a relative error of only 2 percent. Not bad!

pause                                             %Touch any key to continue
clf

%     Finally, any end condition not covered explicitly by  csapi  or  csape
% can be handled by constructing the interpolant with the 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.

pause                                             %Touch any key to continue


%     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 to data

x = 0:.25:3;
y = x.*(-1 + x.*(-1+x.*x/5));

% with

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

% For simplicity, the data are taken from a quartic polynomial that happens to
% satisfy the specific side condition.

pause                                             %Touch any key to continue

% Then, in addition to the interpolant with the default end conditions ...

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

% ... and its first derivative, we also construct the cubic
% spline interpolant to zero data, and with a slope of  1  at  e .

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

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

pause                                             %Touch any key to continue

% ... 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 the linear combination with the help of  fncmb :

pp = fncmb(pp0,(c-lam1)/lam0,pp1);

pause                                             %Touch any key to continue

% We expect  pp  to fit our quartic polynomial slightly better near  e
% than does the interpolant  pp1  with the default conditions ...

xx = (-.3):.05:.7;
yy = xx.*(-1 + xx.*(-1+xx.*xx/5));
plot(xx, fnval(pp1,xx) - yy, 'x')
hold on
plot(xx, fnval(pp,xx) - yy, 'o')
title('default (x) vs. nonstardard (o)');pause
hold off

%  ... and we are not disappointed:
pause                                             %Touch any key to continue


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

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

% ... and solve the linear system for the coefficients 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 points, hence   pp  will match the
% given data for any choice of  d0  and  d2 ).

pause                                             %Touch any key to continue

%     For amusement, we use MATLAB's 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 = xxx.*(-1 + xxx.*(-1+xxx.*xxx/5));
plot(xxx, yyy - fnval(pp,xxx),'x')
title('error in spline interpolant to y = x.*(-1 + x.*(-1+x.*x/5))')
pause                                             %Touch any key to continue


%    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,y,'clamped',[-1,-7+(4/5)*27]),xxx),'o')
title('nonstandard (x) vs endslopes (o)');pause
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)
pause                                             %Touch any key to continue
echo off, close
