clear
reset(0)
echo on
%In this example we demonstrate two approaches to
%finding the minimum of a function of two variables.
%We shall use the method of steepest descent, and 
%the simplex method as implemented by the matlab 
%routine fmins.  The function we shall minimize is 
%given by:
%
% f(x)=exp((x(1)-1)^2+(x(2)-1)^2/4)
%
%The contour lines of this function are ellipses
%centered about (1,1).  We must store this function
%under the name fex.m in the appropriate directory.
pause
%Let's plot up the contour lines of this function:
xx1=[.5:.1:1.5];
xx2=[0:.1:2];
for i=1:length(xx2)
 for j=1:length(xx1)
  z(i,j)=fex([xx1(j),xx2(i)]);
  echo off
 end
end
echo on
contour(xx1,xx2,z,10)
hold on
pause
%The first method we shall try is the method of
%steepest descent.  We pick the descent direction
%as the negative of the gradient of the function.
%We then do a line search for the minimum in this
%direction using some algorithm.  In this case we
%shall use successive parabolic interpolation.  We
%will start our iterations at the point (0.5,0.5):
pause
x=[0.5,0.5];
plot(x(1),x(2),'o')
set(gca,'FontSize',18)
xlabel('x(1)')
ylabel('x(2)')
pause(1)
delta=1;
eps=1e-5;
n=length(x);
sdit=0;
flag=1; %This is a flag for the echo command.
while delta>1e-4,
 %We compute the search direction using a center
 %difference formula:
 for i=1:n
  ep=zeros(size(x));
  ep(i)=eps;
  gradf(i)=(fex(x+ep)-fex(x-ep))/2/eps;
  echo off
 end
 if flag==1;echo on;end
 direc=-gradf/norm(gradf);
 %and now for the root finding in this direction:
 clear alpha f
 alpha(1)=0.1;
 alpha(2)=0.2;
 alpha(3)=0.3;
 f(1)=fex(x+alpha(1)*direc);
 f(2)=fex(x+alpha(2)*direc);
 f(3)=fex(x+alpha(3)*direc);
 i=3;
 while abs(alpha(i)-alpha(i-1))>1e-6
  ft=f(i-2:i)';
  at=alpha(i-2:i)';
  a=[ones(3,1),at,at.^2];
  coeff=a\ft;
  i=i+1;
  alpha(i)=-coeff(2)/2/coeff(3);
  f(i)=fex(x+alpha(i)*direc);
  echo off
 end
 if flag==1;echo on;end
 x=x+alpha(i)*direc;
 plot(x(1),x(2),'o')
 pause(1)
 sdit=sdit+1;
 delta=abs(alpha(i));
 flag=0;
 echo off
end
echo on
hold off
%The minimum is given by:
x
%which was found in this number of iterations:
sdit
pause
%Now let's compare this to the modified simplex method
%implemented in the matlab routine fmins.  This is
%alot easier to use, since mathworks has already
%done all of the work for you.  We will also use the
%modified function fexp.m so that we can plot the 
%intermediate steps.  Thus:
x=[0.5,0.5];
contour(xx1,xx2,z,10)
set(gca,'FontSize',18)
xlabel('x(1)')
ylabel('x(2)')
hold on
x=fmins('fexp',x)
%So the minimum is found, but it took nearly 100
%function evaluations
hold off
echo off


_______________________________________


function y=fex(x)
y=exp((x(1)-1)^2+0.25*(x(2)-1)^2);


_______________________________________


function y=fexp(x)
y=exp((x(1)-1)^2+0.25*(x(2)-1)^2);
%This next bit graphs the intermediate steps
plot(x(1),x(2),'o')
drawnow