% CSAPIDEM	Demonstrate cubic spline interpolation.
% Copyright (c) 1990-92 by Carl de Boor and The MathWorks, Inc.
                                    % latest change: March 22, 1991
                                    % latest change: 27 mar 92 (periodic)
                       % latest change: 13 apr 92 (other end conditions)
clc;clg;echo on;home
%     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 abscissae in strictly increasing order, the cubic
% spline interpolant is a piecewise cubic with breakpoints  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
% breakpoint, 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),'-',x,y,'o');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),'-',x,y,'o');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),'-',x,y,'o');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.
%     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);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);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 
% breakpoint 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);pause
%The difference is as large as .2, but decays rapidly as we move away from
% 1. This ilustrates 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 pp-form 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 pp-form 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));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);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));pause

% There are jumps, but they are within roundoff.

pause                                             %Touch any key to continue

% Here is also the integral
ipp = fnder(pp);plot(xx,fnval(ipp,xx));pause
% and here is the difference between the constructed and the ideal integral:
ipp=fnint(pp);plot(xx,fnval(ipp,xx)-subplus(xx-2).^4/4);pause

pause                                             %Touch any key to continue

%  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 pp-form.  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);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 
% pp-form, but for several different end conditions. For example, the call
%
%         pp = csape(x,y,[2 2])
%        
% uses the `natural' end conditions. This means that the second derivative is 
% zero at the two extreme breakpoints. This is indicated here by specifying
% a  2  for both ends, meaning specification of 2nd derivatives, but leaving
% the 2nd derivative value to the default, which is  0 .

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

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

plot(xx, fnval(pp,xx) - subplus(xx-1).^3 );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,[2 2],6*subplus(x([1 length(x)])-1));
%                           ^ ^  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
%                           | |  value of 2.deriv.at endpoints
plot(xx, fnval(pp,xx) - subplus(xx-1).^3 );pause;

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,[1 1],[sl,sr])
%        
% supplies the cubic spline interpolant to the data  x , y  which 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 abscissa).

% 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 );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 ordinates 
y = [0 -1 0 1 0];
% at the abscissae
x = (pi/2)*[-2:2];
% should be a good approximation to the sine function:
pp = csape(x,y,[0 0]);
%               ^ ^ 
%               | |  specification of periodic end conditions

%   Here is a plot of the difference between the sine function and this
% periodic piecewise cubic function, ...
xx = (pi/50)*[-50:50];
plot(xx, sin(xx) - fnval(pp,xx), 'x'); pause;

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

pause                                             %Touch any key to continue

%     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 which 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 = pppce(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 = pppce(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

% which 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(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');pause
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,[1 1],[-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(fnder(fnder(pp))),3)
pause                                             %Touch any key to continue

