% HW4 Problem 3 
% 
% Sample driver for interpolation
%
% You need to implement:
% - newtondd.m - Newton divided differences algorithm
% - newtonev.m - divided differences polynomial evaluation
% - lagrange.m - evaluate interpolation polynomial in Lagrange form
%
% Please make sure that your functions newtonev and lagrange can evaluate the
% interpolation polynomial at all elements of a vector (i.e. that you use
% the vectorized version of operators), otherwise the plotting part may not work.



% B&F 3.1.5 c and 3.1.7 c

fprintf(' --- B&F 3.1.5 c and 3.1.7 c --- \n');

% interpolation node absicssas (we reordered st the node at which we evaluate
% the interpolation is always inside the interval)
xd = [0.2, 0.3, 0.4, 0.1];
% ordinates (BOOK HAS A MISSING MINUS SIGN)
yd = [-0.28398668, 0.00660095, 0.24842440, -0.62049958];
% evaluation point
x  = 0.25;                    
% true function (USING VECTORIZED OPERATIONS)
f = @(x) x.*cos(x) - 2*x.^2 + 3*x -1;

% The function and first few derivatives are:
%  f(z)     = -2*z^2 + z*cos(z) + 3*z -1
%  f'(z)    = cos(z) -4*z -z*sin(z) + 3
%  f''(z)   = -2*sin(z) -z*cos(z) - 4
%  f'''(z)  = z*sin(z) -3*cos(z)
%  f''''(z) = 4*sin(z) + z*cos(z)

% here we give the part of the bound depending on f. It would be better
% to find the max of this functions over the interpolation interval [0.1,0.4]

% for interp poly of degree 1 we bound  f''
fmax(1) =  abs(-2*sin(0.3)-0.3*cos(0.3)-4);
% for interp poly of degree 2 we bound  f'''
fmax(2) = abs(0.2*sin(0.2) -3*cos(0.2))
% for interp poly of degree 2 we bound  f''''
fmax(3) = abs(4*sin(0.4) + 0.4*cos(0.4));

% vectors to store interpolation polynomial evaluated at x
% and error
lagr = zeros(3,1); newt = zeros(3,1);  err = zeros(3,1);
for k=1:3, % loop on degree
 % Lagrange interpolation
 % little trick to take only  the first k nodes
 lagr(k) = lagrange(xd(1:k+1),yd(1:k+1),x);

 % Newton interpolation using divided differences
 d = newtondd(xd(1:k+1),yd(1:k+1));
 newt(k) = newtonev(xd(1:k+1),d,x);

 % compute error bound
 err(k) = fmax(k)/factorial(k+1);
 for j=1:k+1;
  err(k) = err(k)*(x-xd(j));
 end;% for j
end;% for k

% display errors
fprintf('%10s %13s %13s %13s\n',...
  'degree','Lagrange','Newton','err bound');
for k=1:3,
 fprintf('%10d %13.4e %13.4e %13.4e\n',...
  k,abs(lagr(k)-f(x)),abs(newt(k)-f(x)),abs(err(k)));
end;% for k

% Show all three interpolation polynomials in a plot


figure(1); % IF DOING MORE THAN ONE FIGURE, CHANGE FIGURE NUMBER
clf;       % clears figure
xs = linspace(0,1);  % equally space points in an interval containing interp. points

% we are going to do several plots, so we tell Matlab to 
% accumulate all the plot commands
hold on; 

% Plots the true function
plot(xs,f(xs),'r-');

% Plots the interpolation nodes
plot(xd,yd,'*');

% Plot the interpolation polynomials
for k=1:3,
 % Newton interpolation using divided differences
 d = newtondd(xd(1:k+1),yd(1:k+1));
 newt = newtonev(xd(1:k+1),d,xs);
 plot(xs,newt,'b--');
end;

% We don't want to add anymore plots to the current figure
hold off;

% put some labels
xlabel('x');
title('polynomial interpolation 3.1.5 c');

% print to a file
% print('-depsc2','prob3_sample.eps');

% Here is some sample output (compare with book's solutions)
%  
%  >> prob3_sample
%   --- B&F 3.1.5 c and 3.1.7 c --- 
%  
%  fmax =
%  
%     4.877641360060361
%     2.900465867364713
%     1.926097766835756
%  
%      degree      Lagrange        Newton     err bound
%           1    5.9210e-03    5.9210e-03    6.0971e-03
%           2    1.7455e-04    1.7455e-04    1.8128e-04
%           3    2.8798e-06    2.8798e-06    4.5143e-06
%  >> 
%  
