% 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-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 - 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 );
  Err(i,7) = abs( (df0 - Pade4( x0, h, 15 )) / 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',Err(:,1),Err(:,7),'-s');
legend('1st ord Forward Difference','2nd ord Forward Difference', ...
       '2nd ord Central Difference', '4th ord Central Difference',...
      '4th ord Pade');

figure(2); 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 ) * (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 forward difference formula  
function [fd] = FD2( x0, h )  
  
  fd = ( - f(x0 + 2.0*h) + 4.0*f(x0 + h) - 3.0*f(x0) ) / (2.0*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;

  
% Computes the Pade approximation at x0 using N2 nodes at each side
% with spacing h

function [fd] = Pade4( x0, h, N2 )  
  
% total number of equations in the system; 
  N = 2 * N2 + 1;
  
% Using a tridiagonal solver;  
% Assembling the system matrix;
  for i=1:N
    k = - (N2 + 1 - i);
    b(i) = 1.0;
    if ( i > 1 & i < N )
      c(i) = 0.25;
    else
      c(i) = 0.0;
    end;
    if ( i < N - 1 )
      a(i) = 0.25;
    elseif ( i == N - 1 )
	a(i) = 0.0;
    end;
    if ( i > 1 & i < N )
% RHS terms in the interior;      
      rhs(i) = 3.0 * ( f(x0 + (k + 1) * h) - f(x0 + (k - 1) * h) ) / (4.0 * h);
    else
% RHS terms at endpoints;      
      rhs(i) = ( f(x0 + (k + 1) * h) - f(x0 + (k - 1) * h) ) / (2.0 * h);      
    end;
  end;
  
% Solving the system;
  X = tridiag1( a, b, c, rhs );
  fd = X( N2 + 1 );

  
  
  
%      x = tridiag1( a, b, c, y )
%
%  Solve the  N x N  tridiagonal system for x:
%
%  [ b(1)  c(1)                                  ] [  x(1)  ]   [  y(1)  ] 
%  [ a(1)  b(2)  c(2)                            ] [  x(2)  ]   [  y(2)  ] 
%  [       a(2)  b(3)  c(3)                      ] [        ]   [        ] 
%  [            ...   ...   ...                  ] [  ...   ] = [  ...   ]
%  [                    ...    ...    ...        ] [        ]   [        ]
%  [                        a(N-2) b(N-1) c(N-1) ] [ x(N-1) ]   [ y(N-1) ]
%  [                               a(N-1)  b(N)  ] [  x(N)  ]   [  y(N)  ]
%
%  y must be a vector (row or column); N is determined from its length.
%  a, b, c must be vectors of lengths at least N-1, N, and N-1 respectively,
function x = tridiag1( a, b, c, y )

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Check that the input arrays have acceptable sizes

if min(size(y))~=1 | min(size(a))~=1 | min(size(b))~=1 | min(size(c))~=1
   error('a, b, c, y must be vectors');
end
N = length(y);
if length(a)<N-1 | length(b)<N | length(c)<N-1
   error('a, b, c must be vectors of length at least N-1, N, N-1');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Solve the problem by back substitution

gam = zeros(1,N);   %  hold the LU decomposition
x = zeros(size(y));

%  phase 1:  LU decomposition
beta = b(1);
if beta==0
   error('beta = 0 at j=1:  matrix is singular');
end

x(1) = y(1) / beta;
for j=2:N
   gam(j) = c(j-1) / beta;
   beta = b(j) - a(j-1)*gam(j);
   if beta==0
      error( [ 'beta = 0 at j=' num2str(j) ':  matrix is singular' ]);
   end
   x(j) = ( y(j) - a(j-1)*x(j-1) )/beta;
end

%  phase 2:  back-substitution
for j=N-1:-1:1
   x(j) = x(j) - gam(j+1)*x(j+1);
end

