% Analyzing the accuracy of various finite-difference formulae
% Written by Bartosz Protas for the course MATH745 (January 2004)

function diff_01
  
echo off, clc
clear all
 

x0 = 1.0;
hmax = 0.09999999;
%hmin = 1.0e-4;
hmin = 1.0e-10;


h = hmax;
df0 = df( x0 );
i = 1;
while ( h >= hmin );  
  Err(i,1) = h;
  Err(i,2) = abs( (df0 - FD( x0, h )) / df0 );
  Err(i,3) = abs( (df0 - CD( x0, h )) / df0 );
  Err(i,4) = abs( (df0 - CD4( x0, h )) / df0 );
  Err(i,5) = abs( (df0 - CSD( x0, h )) / df0 );
  i = i + 1;
  h = h / 2.0;
end;


figure(1); loglog(Err(:,1),Err(:,2),'-x',Err(:,1),Err(:,3),'-+',Err(:,1),Err(:,4),'-o',Err(:,1),Err(:,5),'-d');
legend('1st ord Forward Difference','2nd ord Central Difference', ...
       '4th ord Central Difference','2nd ord Complex Step');


% Implementation of the functions

% The function to be differentiated;
function [f] = f( x )

  f = exp( -x ) / sqrt(tan( x ));


% The exact derivative of the function f(x) = exp( -x ) / sqrt( tan(x) )
function [df] = df( x )

  df = - exp( -x ) * (2.0*tan(x)*(cos(x))^2+1.0)/(2.0*(tan(x))^(1.5)*(cos(x))^2);

  
% First-order forward difference formula
function [fd] = FD( x0, h )  
  
  fd = ( f(x0 + h) - f(x0) ) / h;

  
% Second-order central difference formula
function [fd] = CD( x0, h )  
  
  fd = ( f(x0 + h) - f(x0 - h) ) / (2.0 * h);
  
  
% Fourth-order central difference formula  
function [fd] = CD4( x0, h )  
  
  fd = ( - f(x0 + 2.0*h) + 8.0*f(x0 + h) - 8.0*f(x0 - h) + f(x0 - 2.0*h) ) / (12.0*h);
  
% Complex Step Derivative
function [fd] = CSD( x0, h )  
  
  fd = imag( f(x0 + j*h)  ) / h;

  
