% Analyzing the accuracy of various finite-difference formulae
% Written by Bartosz Protas for the course MATH745 (January 2018)
% This version uses arbitrary-precision computations with Advanpix
% ( available from https://www.advanpix.com/ )

function diff_01mp

echo off, clc
clear all
 
ap = 134;  
mp.Digits(ap)
  
 

x0 = mp('1.0');
hmax = mp('0.09999999');
%hmin = 1.0e-4;
hmin = mp('1.0e-20');


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 - FD2( x0, h )) / df0 );
  Err(i,4) = abs( (df0 - CD( x0, h )) / df0 );
  Err(i,5) = abs( (df0 - CD4( x0, h )) / df0 );
  Err(i,6) = abs( (df0 - CSD( x0, h )) / df0 );
  i = i + 1;
  h = h / mp('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',Err(:,1),Err(:,6),'-s');
legend('1st ord Forward Difference','2nd 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 ) * (mp('2.0')*tan(x)*(cos(x))^2+mp('1.0'))/(mp('2.0')*(tan(x))^(mp('1.5'))*(cos(x))^2);

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

  
% Second-order forward difference formula  
function [fd] = FD2( x0, h )  
  
  fd = ( - f(x0 + mp('2.0')*h) + mp('4.0')*f(x0 + h) - mp('3.0')*f(x0) ) / (mp('2.0')*h);

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

  
