%% Interpolation HW Chapter 5
% Life Expectancy - Predicting value at 2007.
% First view the data.

clear
figure
year = 1980:5:2010;
t = year-1980;
y = [73.7, 74.7, 75.4, 75.8, 77.0, 77.8, 78.7];
n = length(y)-1;
plot(year,y,'.','MarkerSize',15)
title('Life Expectancy 1980-2010')
xlabel('Year')
ylabel('Life Expectancy')

%
% Do polynomial interpolation first.

c = polyfit(t',y',n);
hold on, 
x=linspace(1980,2010,100)
plot(x,polyval(c,x-1980))
hold off
polyval(c,2007-1980)

%
% Next perform cubic interpolation.

plot(year,y,'.','MarkerSize',15)
title('Life Expectancy 1980-2010')
xlabel('Year')
ylabel('Life Expectancy')
hold on
plot(x,interp1(t,y,x-1980,'spline'))
hold off
interp1(t,y,2007-1980,'spline')

%
% Do piecewise linear interpolation for f(x)=\cos\left(\pi^2x^2\right) for 1600 points in x\in[0,1]for n equispaced nodes, 𝑛=10,20,40,80,160. Record the error. See Example 5.2.3.

clear
f = @(x) cos(pi^2*x.^2) ;
x = linspace(0,1,1600)';  % sample the difference at many points
n_ = [10,20,40,80,160]';
err_ = 0*n_;
for i = 1:length(n_)
    n = n_(i);
    t = linspace(0,1,n+1)';   % interpolation nodes
    p = plinterp(t,f(t));
    err = max(abs( f(x) - p(x) ));
    err_(i) = err;
end
format longE
err_

%
% Show quadratic convergence using a loglog plot.

figure
loglog( n_, err_, '.-' ,'MarkerSize',15)
hold on
C=.473;
loglog( n_, C*(n_/n_(1)).^(-2), '--')  % From example
%C = 47;
%loglog( n_, C*n_.^(-2), '--' )        % Simpler version
xlabel('n'), ylabel('||f-p||_\infty')    
title('Convergence of Interpolation')   
legend('error','2nd order')     
axis tight    
hold off

%

function H = hatfun(x,t,k)
% HATFUN   Hat function/piecewise linear basis function.
% Input: 
%   x      evaluation points (vector)
%   t      interpolation nodes (vector, length n+1)
%   k      node index (integer, in 0,...,n)
% Output:
%   H      values of the kth hat function

n = length(t)-1;
k = k+1;  % adjust for starting with index=1

% Fictitious nodes to deal with first, last funcs.
t = [ 2*t(1)-t(2); t(:); 2*t(n+1)-t(n) ];
k = k+1;  % adjust index for the fictitious first node

H1 = (x-t(k-1))/(t(k)-t(k-1));   % upward slope
H2 = (t(k+1)-x)/(t(k+1)-t(k));   % downward slope

H = min(H1,H2);
H = max(0,H);
end

function p = plinterp(t,y)
% PLINTERP   Piecewise linear interpolation.
% Input:
%   t     interpolation nodes (vector, length n+1)
%   y     interpolation values (vector, length n+1)
% Output:
%   p     piecewise linear interpolant (function)

n = length(t)-1;
p = @evaluate;

    % This function evaluates p when called.
    function f = evaluate(x)
        f = 0;
        for k = 0:n
            f = f + y(k+1)*hatfun(x,t,k);
        end
    end

end