clear

%%
% For illustration, here is a spline interpolant using just a few nodes. 
f = @(x) exp(sin(7*x));         
fplot(f, [0,1])
t = [0, 0.075, 0.25, 0.55, 0.7, 1]';  % nodes
y = f(t);                            % values at nodes
hold on 
plot(t, y, 'r.')
tt = linspace(0., 1., 100);
plot(tt, spline(t, y, tt))
xlabel('x') 
ylabel('f(x)')
title('Function and cubic spline interpolant')
legend('function', 'data points', 'spline interpolation')
grid on

%%
% Now we look at the convergence rate as the number of nodes increases.
x = linspace(0, 1, 10001)';  % sample the difference at many points
nn = 2.^(4:0.5:9)';
err = 0*nn;

for i = 1:length(nn)
    n = nn(i);
    t = linspace(0, 1, n+1)';   % interpolation nodes
    err(i) = norm(f(x) - spline(t, f(t), x), inf);
end

%%
% Since we expect convergence that is $O(h^4)=O(n^{-4})$, we use a log-log
% graph of error and expect a straight line of slope $-4$.
clf
loglog(nn, err, '.-' )
hold on
loglog(nn, 0.1*(nn/nn(1)).^(-4), '--')
xlabel('n')
ylabel('||f-S||_\infty')
title('Convergence of spline interpolation')
legend('error','4th order')
grid on