%% Example for interpolation with piecewise cubic functions
%% Function 1/(1+x^2) on [-5,5] with 7 equidistant nodes
% Consider the function f(x)=(1+x^2)^{-1}. The derivative is f'(x)=-(1+x^2)^{-2}\cdot 2x.
% We compute the function values and the derivatives at the nodes.
% We will use 1000 equidistant points xe on [-5,5] for plotting and finding the maximal errors.

clearvars; clf                  % clear variables and figures
a = -5; b = 5;
x = linspace(a,b,7);            % 7 equidistant nodes 
y = 1./(1+x.^2);                % function values at nodes
s = -2*x./(1+x.^2).^2;          % derivatives at nodes

xe = linspace(a,b,1e3);         % evaluation points
ye = 1./(1+xe.^2);              % function values at xe

plot(x,y,'o',xe,ye,'k:')
title('1/(1+x^2)')

%
%% Piecwise cubic Hermite interpolation
% The piecewise cubic Hermite function uses function values and derivatives at all nodes.
% You need to download the m-file hermite.m from the course web page.

yherm = hermite(x,y,s,xe);      % find pw cubic Hermite interpolation and evaluate at points xe
subplot(2,1,1);
plot(x,y,'o',xe,ye,'k:',xe,yherm,'b');
title('pw cubic Hermite interpolation')
legend('given points','f(x)','Hermite')
subplot(2,1,2)
plot(x,0*y,'o',xe,yherm-ye,'b'); grid on
title('Errors')
maxerr_herm = max(abs(yherm-ye))

%
%% Complete cubic spline
% The complete spline uses function values at the nodes, and the derivative at the first and last node.

sa = s(1); sb = s(end);         % slopes at left and right endpoint
ycomp = spline(x,[sa,y,sb],xe); % find complete cubic spline and evaluate at points xe
subplot(2,1,1);
plot(x,y,'o',xe,ye,'k:',xe,ycomp,'r');
title('Complete cubic spline')
legend('given points','f(x)','compl. spline')
subplot(2,1,2)
plot(x,0*y,'o',xe,ycomp-ye,'r'); grid on
title('Errors')
maxerr_complete = max(abs(ycomp-ye))

%
%% Not-a-knot cubic spline
% The not-a-knot spline uses only the function values at the nodes.

ynak = spline(x,y,xe);          % find not-a-knot cubic spline and evaluate at points xe
subplot(2,1,1);
plot(x,y,'o',xe,ye,'k:',xe,ynak,'g');
title('Not-a-knot cubic spline')
legend('given points','f(x)','not-a-knot spline')
subplot(2,1,2)
plot(x,0*y,'o',xe,ynak-ye,'g'); grid on
title('Errors')
maxerr_nak = max(abs(ynak-ye))

%
%% Comparison
% The Hermite function gives the smallest error (it uses the most information about function f).
% The not-a-knot spline performs worse than the complete spline near the endpoints (the complete spline uses the derivatives at the endpoints).

subplot(2,1,1)
plot(x,y,'o',xe,ye,'k:',xe,yherm,'b',xe,ycomp,'r',xe,ynak,'g')
legend('given points','f(x)','Hermite','complete spline','not-a-knot spline')
subplot(2,1,2)
plot(x,0*y,'o',xe,yherm-ye,'b',xe,ycomp-ye,'r',xe,ynak-ye,'g'); grid on
title('Errors')
legend('nodes','Hermite','complete spline','not-a-knot spline')