checkDerivative.m

% function [fh,ph,th] = checkDerivative(f0tn,x0);
%
% checks the implementation of a derivative by comparing the function 
% with the Taylor-poly
%
%   \| f(x0 + h ) - TP_p(x0,h) \|   !=   O( h^{p+1} ) 
%
% Input:
%   f0tn    function handle
%  x0      expanding point
%
% Output:
%  fh      figure handle to graphical output
%  ph      plot handle
%  th      text handle
% call checkDerivative for a minimal example.
%==============================================================================

function [fh,ph,th] = checkDerivative(fctn,x0,varargin)

if nargin == 0, % help and minimal example
  help(mfilename); 
  fctn = @xSquare;  x0 = 1;  checkDerivative(fctn,x0);
  return;
end;


fig = [];
for k=1:2:length(varargin),     % overwrites default parameter
  eval([varargin{k},'=varargin{',int2str(k+1),'};']);
end;

fprintf('%s: derivative test\n',mfilename);
fprintf('T0 = |f0 - ft|, T1 = |f0+h*f0'' - ft|\n');
[f0,df] = feval(fctn,x0);  
if isstruct(df),
  try
    [f0,para,df] = feval(fctn,x0);  
  catch
    if isfield(df,'Q'), df = df.Q;  end;
  end;
end;

h = logspace(-1,-10,10);
v = randn(size(x0));
if isnumeric(df),
  dvf = df*v; 
elseif isa(df,'function_handle'),
  dvf = df(v);
else
  keyboard;
end;

for j=1:length(h),
  ft = feval(fctn,x0+h(j)*v);      % function value
  T0(j) = norm(f0-ft);             % TaylorPoly 0
  T1(j) = norm(f0+h(j)*dvf - ft);  % TaylorPoly 1
 
  f0-ft ;
  h(j)*dvf ;

  fprintf('h=%12.4e     T0=%12.4e    T1=%12.4e\n',h(j),T0(j),T1(j));
end;
if isempty(fig),
  fh = figure;
else
  fh = figure(fig);
end;
ph = loglog(h,[T0;T1]); set(ph(2),'linestyle','--')
th  = title(sprintf('%s: |f-f(h)|,|f+h*dvf -f(h)| vs. h',mfilename));