% goldsect.m
%
% Script executes a few steps of golden section search using
% workspace variables named below:
%   fname  - string containing name of function to operate on
%   a,b,c  - real variables a<b<c satisfying f(b)<=min{f(a),f(c)}
%   nstep  - number of steps to iterate

% Philip D. Loewen
% 25 Jan 98:  Added printing for final bracket and fcn vals
% 31 Dec 97:  Basic setup and debugging


gamma = (3.0-sqrt(5.0))/2.0;

disp([' ']);
prtmtx([a b c],'[a b c]');

ya = eval([fname,'(a)']);
yb = eval([fname,'(b)']);
yc = eval([fname,'(c)']);

xx = linspace(a,c,100);
yy = eval([fname,'(xx)']);
plot(xx,yy,'r');
hold on;
zz = axis;
axis(zz);
ymin = zz(3);
ymax = zz(4);
disp([' ']);
disp(['   a          ya         b          yb         c          yc ',...
      '        (c-a)    max{y-yb}']);

for k=1:nstep,
  prtmtx([a,ya,b,yb,c,yc,c-a,max(yc-yb,ya-yb)]);
  plot([a,a],[ymin,ymax],'b');
  plot([c,c],[ymin,ymax],'b');
  plot([b,b],[ymin,ymax],'r');
%  pause
  if ( (b-a)<(c-b) ),  
    % Interval (b,c)  is longer
    z = b + gamma*(c-b);
    yz = eval([fname,'(z)']);
    if (yz < yb),
      % New point is lowest so far, shrink to interval (b,c)
      a = b; ya = yb;
      b = z; yb = yz;
    else
      % Lowest point remains at point b, shrink to interval (a,z)
      c = z; yc = yz;
    end;
  else
    % Interval (a,b) is longer
    z = b + gamma*(a-b);
    yz = eval([fname,'(z)']);
    if (yz < yb),
      % New point is lowest so far, shrink to interval (a,b)
      c = b; yc = yb;
      b = z; yb = yz;
    else
      % Lowest point remains at point b, shrink to interval (z,c)
      a = z; ya = yz;
    end;
  end;
end;

disp(' ');
disp('Golden Section Minimization:  Final Bracket and Function Values:');
disp(['a = ',sprintf('%22.16e',a),'; ',fname,'(a)=',sprintf('%22.16e',ya)]);
disp(['b = ',sprintf('%22.16e',b),'; ',fname,'(b)=',sprintf('%22.16e',yb)]);
disp(['c = ',sprintf('%22.16e',c),'; ',fname,'(c)=',sprintf('%22.16e',yc)]);

disp([' ']);
title(['Golden Section Minimization of Function ',fname],'FontSize',16);
xlabel('x','FontSize',14);
ylabel(['y = ',fname,'(x)'],'FontSize',14);
