% InterpErr Script displays error for polynomials
%           interpolating sin(x) on [0 pi/2].

close all
clc
n = input('Please enter number of interpolation points: ');
if n<1, error('Who ya tryin'' to kid?'), end
Cheby = input('Uniform [0] or Chebyshev [1] points? ');

if Cheby
%   x = pi/4*(1-cos(linspace(0,pi,n)'));
  x = pi/4*(1-cos([1:2:2*n-1]*pi/2/n)');
else
  x = linspace(0,pi/2,n)';
end
y = sin(x);

% Construct the interpolant, and 
% compute the error and maximum error.

% Compute divided difference table
% [ y[1] y[12] y[123] ... y[12..n-1] ]'
% Cf. algorithm Kincaid-Cheney 3e p.332.

ddy = y;
for k = 1:n-1
% Order k divided differences in last n-k positions
  ddy(k+1:n) = (ddy(k+1:n)-ddy(k:n-1)) ./ (x(k+1:n)-x(1:n-k));
end

npt = 193;		    % # of points for plot
z = linspace(0,pi/2,npt)';

% Evaluate the interpolant by nested multiplication
% Cf. Eq(6) & algorithm p.310, Kincaid-Cheney 3e.

sI = ddy(n)*ones(size(z));
for k = n-1:-1:1
  sI = (z-x(k)) .* sI + ddy(k);
end

ierr = sin(z) - sI;         % Actual error
ermx = max(abs(ierr));

% Evaluate the error formula:

%   sin(z)-p(z) = sin^(n)(xi) * prod(z-x(i))/n!

% [Replace n-th derivative of sin with +1 or -1.]

berr = (-1)^floor(n/2) * ... 
       [prod([ones(n,1)*z'-x*ones(1,npt)])/prod(1:n)]';
plot(z,ierr,'m',z,berr,'b',x,0,'r*')
xlabel('x'), ylabel('Error (magenta), Bound (blue)'),grid
title(sprintf('Maximum interpolation error with %d points: %4.2e',n,ermx))