% Approximate a Derivative by a Difference Quotient
% The example is for f(x) = sin(x) at x = 1
% Adapted from "Numerical Computing with IEEE Floating Point Arithemtic"
% by Michael Overton 
% Chapter 11 on Cancellation, page 70
% Kartik's MATLAB code (19/2/03)

format long;
n = 1;
x = 1.0;
h = 1.0;
deriv = cos(x);

% The results are stored in a file 'cancellation.res'

fid=fopen('cancellation.res','w');

fprintf(fid,'Derivative = %13.6e\n',deriv);

% Let h range from 10^{-1} down to 10^{-20}

fprintf(fid,'h    diffquo    error\n');

h_vector = [];
err_vector = [];

while ( n <= 20)
   % h = 10^{-n}
   h = h/10;
   h_vector = [h_vector; h];
   % diffquo is the difference quotient  
   diffquo = (sin(x+h)-sin(x))/h; 
   error = abs(deriv - diffquo);
   err_vector = [err_vector; error];
   fprintf(fid,'%5.1e %13.6e % 13.6e \n',h, diffquo, error);
   n = n+1;
end

fclose(fid);

% Plot error as a function of h (logarithmic plot)

figure;

loglog(h_vector,err_vector,'*');
xlabel('h');
ylabel('error');
title('log(error) vs log(h)');