clear
echo on
%This program demonstrates the machine round off
%error and algorithm error of a higher order forward
%difference algorithm for derivatives.  The 
%algorithm is:
%      df/dx = (-3f(x)+4f(x+h)-f(x+2h))/2h
%where h is some small quantity.  We apply the formula
%to the function f(x)=exp(x) over the range x=0:.1.
%We average the error for a number of different values
%of x in order to average out the random round off
%error.  We also pick values of h in the range of
%e-4 to e-7.
pause

x=[0:.001:.1];
for i=1:30
  h(i)=10^(-i/10)/10^4;
  deriv=(4*exp(x+h(i))-3*exp(x)-exp(x+2*h(i)))/(2*h(i));
  error(i)=sum(abs(deriv-exp(x)))/101;
echo off
end;
echo on
loglog(h,error)
xlabel('h')
ylabel('error')

[minerror,i]=min(error);

%The minimum error is:
minerror
%The optimum h is:
hopt=h(i)
pause

%As can be seen, the minimum error occurs for a larger
%value of h than was the case for the first order
%forward difference approximation and the error is much
%smaller.  A detailed analysis would show that the
%algorithm error varies as h^2 rather than h^1, thus
%the optimum value of h is O(eps*f/f''')^(1/3).

echo off