% dilate.m
% Script to bracket a minimum for a scalar function of one real variable.
%  fname (in workspace) - string for function name
%  aa (returned) - left end of bracketing interval
%  bb (returned) - interior point of bracketing interval
%  cc (returned) - right end of bracketing interval

% Philip D. Loewen
% 30 December 1997

gold = (3.0+sqrt(5.0))/2.0;   % Golden ratio expansion factor
maxswaps = 10;                % Maximum number of dilations allowed

x0 = 0.0; y0 = eval([fname,'(x0)']);   % Starting interval has endpoints
x1 = 1.0; y1 = eval([fname,'(x1)']);   % x0=0 and x1=1.

% Relabel the endpoints according to which has the greater y-value.
if y0<y1,
  xlo = x0; xhi = x1; 
else
  xlo = x1; xhi = x0;
end;
ylo = min([y0,y1]); yhi = max([y0,y1]);

% Do a sample dilation to get the ball rolling.
xnew = xhi + gold*(xlo-xhi);
ynew = eval([fname,'(xnew)']);

disp([' ']);
disp(['dilate:  Table of intermediate results follows:']);
disp([' ']);
disp(['   xlo        ylo        xhi        yhi        xnew       ynew']);
prtmtx([xlo ylo xhi yhi xnew ynew]);

% Loop through these steps until the minimum is bracketed.
nswaps = 0;
while ynew <= min(ylo,yhi),
  nswaps = nswaps+1;
  if nswaps > maxswaps, break, end;
  xlo = xnew; ylo = ynew;
  xnew = xhi + gold*(xlo-xhi);
  ynew = eval([fname,'(xnew)']);
  prtmtx([xlo ylo xhi yhi xnew ynew]);
end;

% Report the results in increasing order.
aa = min([xhi,xnew]);
bb = xlo;
cc = max([xhi,xnew]);

if nswaps > maxswaps,
  aa = min([xlo,xhi,xnew]);  cc = max([xlo,xhi,xnew]);
  disp([' ']);
  disp(['*** dilate:  Function ',fname,' appears to be monotonic on the interval']);
  disp(['*** [',num2str(aa),',',num2str(cc),']:  try ', ...
                 'rescaling the independent variable.']);
end;

disp([' ']);
disp(['dilate:  Finished; returning values ', ...
      'aa=',num2str(aa),', bb=',num2str(bb),', cc=',num2str(cc),'.']);
disp([' ']);
return;
