% FINDZERO    - Solve for x in f(x)=0
%
%    FINDZERO( F, INTERVAL ) takes an interval INTERVAL containing a zero
%    of the continuous function F, and returns a subinterval also containing
%    a zero of F. Here F is a function handle associated
%    with a correctly vectorized scalar-valued function.

% Philip D Loewen, 2014-10-19

function subint = findzero( f, interval );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% User input zone:
%
method  = 2;        % Method, chosen from 1 to 6
reltolx = 1.E-12;   % Relative tolerance on final bracket length
reltolf = 1.E-12;   % Relative tolerance on final values f(x)
MaxEval = 100;      % Maximum number of evaluations allowed

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Automatic operation starts here.

% Counter for function evaluations
Nevals = 0;

% Save original interval for later reference
a = min(interval(:));  
b = max(interval(:));
fa= f(a);   Nevals = 1 + Nevals;
fb= f(b);   Nevals = 1 + Nevals;

% Sanity checks on given interval
if (fa*fb>0)
  disp('FINDZERO: Given function has the same sign at both ends of the given interval.');
  subint = [];
  return;
  end;
% Deal correctly with some trivial cases.
if (fa==0), subint=[a,a]; return; end;
if (fb==0), subint=[b,b]; return; end;

% Switch to x0, x1 notation. Use the name x1 for the better endpoint.
if (abs(fa)<abs(fb))
  x0 = a; f0 = fa;
  x1 = b; f1 = fb;
else
  x0 = b; f0 = fb;
  x1 = a; f1 = fa;
end;

% Take note of a typical function value (the better initial guess)
ftyp = abs(f1);

disp('**************************************');
disp(['Starting findzero with method ',num2str(method),' ...']);
disp('evals    x0            x1            f0       f1    ;   c            fc');

close all;
figure('Name','Convergence');
hold on;
xlabel('Number of Function Evaluations');
ylabel('log_1_0(bracket length multiplier)');

while (Nevals < MaxEval)

  % STEP 1: PICK INTERMEDIATE POINT c
  switch method
  case 1
    % Option (i): Random style
    desc = 'Random';
    r  = rand(1);
    c  = x0 + r*(x1-x0);    
  case 2
    % Option (ii): Bisection
    desc = 'Random';
    r  = rand(1);
    c  = x0 + r*(x1-x0);    
  case 3
    % Option (iii): Regula Falsi
    desc = 'Random';
    r  = rand(1);
    c  = x0 + r*(x1-x0);    
  case 4
    % Option (iv): Regula Falsi with bisection
    desc = 'Random';
    r  = rand(1);
    c  = x0 + r*(x1-x0);    
  case 5
    % Option (v): Regula Falsi with biased bisection
    desc = 'Random';
    r  = rand(1);
    c  = x0 + r*(x1-x0);    
  case 6
    % Option (vi): Illinois
    desc = 'Random';
    r  = rand(1);
    c  = x0 + r*(x1-x0);    
  otherwise
    % random
    disp(['Warning: method ',method,' not implemented.']);
    c = x0 + rand(1)*(x1-x0);
  end;

  % STEP 1 (CONTINUED): EVALUATE f(c) AND TREAT CASES
  fc = f(c);  Nevals = 1 + Nevals;

  disp(sprintf('%4d:  %12.5e  %12.5e  %8.1e %8.1e; %12.5e %8.1e',...
               [Nevals,x0,x1,f0,f1,c,fc]));

  if (fc==0)
    x0 = c; f0 = 0;
    x1 = c; f1 = 0;
  elseif (fc*f1<0)
    x0 = x1; f0 = f1;
    x1 =  c; f1 = fc;
  elseif (fc*f0<0)
    x1 =  c; f1 = fc;
  else
    disp('Logic error: this is mathematically impossible.');
  end;

  % STEP 2: ASSESS PROGRESS SO FAR, REPORT, PERHAPS QUIT

  totalxshrink = abs(x1-x0)/(b-a);
  if (totalxshrink < reltolx), break; end;

  fshrink = max([abs(f0),abs(f1)])/ftyp;
  if (fshrink < reltolf), break; end;

  % STEP 3: LOOP BACK AND CONTINUE (nothing to do here)
end;

% Set return value for function findzero
subint=[ min([x0,x1]), max([x0,x1]) ];

lastword = ['Method ',num2str(method),' (',desc,') '...
            'stopped after ',num2str(Nevals),' ',...
            'evaluations.'];
disp(lastword);
title(lastword);
grid on;

