% FCDIFF Script demonstrates forward- and 
% centered-difference approximations to derivative
% on the function fexmp(x) = x/(x+1), at x = 4.
% 
close all
% 
f4   = fexmp(4);               % Function value at x = 4
df4  = fexmp(4,1);             % Exact derivative
d2f4 = fexmp(4,2);             % Second derivative
d3f4 = fexmp(4,3);             % Third derivative
% 
% Forward difference approximation
% 
h   = logspace(-12,-4);	       
dq  = ( fexmp(4+h) - f4 ) ./ h;    % Difference quotient
te  = abs(h*d2f4/2);               % Truncation error
% 
% log-log plot with
%   Magenta x's: actual error
%   Blue:    truncation error
%   
loglog(h,abs(df4-dq),'mx', h,te,'b-')  
xlabel('h'), ylabel('Error'), grid on
title('Error in difference quotient, truncation error')
disp('press key to continue')
pause
% 
tot = eps*abs(f4)./h + te;     % Error including roundoff
% 
% log-log plot with
%   Magenta x's: actual error
%   Blue:    truncation error
%   Red:      error estimate including roundoff, assuming 
%             1 ulp in each evaluation of f, no cancellation.
% 
loglog(h,abs(df4-dq),'mx', h,te,'b-', h, tot,'r-')  
xlabel('h'), ylabel('Error'), grid on
title('Error in difference quotient, truncation error + roundoff estimate')
disp('press key to continue')
pause
% 
% Centered difference approximation
% 
figure
h = logspace(-8,-2);
cd = ( fexmp(4+h) - fexmp(4-h) ) ./ (2*h); % Centered Difference
te = abs(h.^2 * d3f4/6);                   % Truncation error
% 
% log-log plot
%   Magenta x's: actual error
%   Blue:    truncation error
loglog(h,abs(df4-cd),'mx', h,te,'b')   
xlabel('h'), ylabel('Error'), grid on
title('Error in centered difference, truncation error')
disp('press key to continue')
pause
% 
tot = eps*abs(f4) ./ (2*h) + te;       % Error including roundoff
% 
% log-log plot
%   Magenta x's: actual error
%   Blue:    truncation error
%   Red:      error estimate including roundoff, assuming 
%             1 ulp in each evaluation of f, no cancellation.
% 
loglog(h,abs(df4-cd),'mx', h,te,'b', h, tot, 'r-')   
xlabel('h'), ylabel('Error'), grid on
title('Error in centered difference, truncation + roundoff')