%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% function gd(alpha,x0)- minimizes f(x)=.5(alpha*x_1^2+x_2^2)
%%%
%%% INPUT: alpha = scalar (alpha > 0 gives elliptic paraboloid)
%%%                       (alpha < 0 gives hyperbolic paraboloid=saddle)
%%%           x0 = 2 by 1 column vector for starting point such that
%%%           x0(1),x0(2) in [-11,11]
%%% EXAMPLE; gd(4,[10;5])
%%% Run this for different alpha values and different starting points
%%%
%%% SAMSI NDHS Undegraduate Workshop May 2006
%%% created: Dr. Amy Langville, North Carolina State University
%%% updated: Anjela Govan, North Carolina State University
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%

function gd(alpha,x0)
hold off
clf

%% points for the surface
X0= [-11: .4: 11];
Y0= [-11: .4: 11];

%% Create mesh plot with contour plot of f(x)=.5(alpha*x_1^2+x_2^2)
figure(1)
[X,Y] = meshgrid(X0,Y0);
Z = .5*(alpha*X.^2 + Y.^2);
surfc(X,Y,Z)    % shows 3D plot of surface f(x)=.5(alpha*x_1^2+x_2^2)
shading('faceted')

%%Contour graph of the surface of f(x) 
figure(2)
v=  alpha*X0.^2 + Y0.^2;
w=[v(1) v(round(length(v)/4) : end-round(length(v)/4) )   v(end)];
[C,H]=contour(X,Y,Z,w);
hold on

xnew=x0;              % user-defined starting point
plot(xnew(1),xnew(2),'xk') %plot the initial point

Q=[alpha 0;0 1];      % Hessian
gradnew=Q*xnew;       % gradient of f at the initial point


disp(['initial point: ',num2str(xnew')]);
disp(['iteration ',num2str(0),', the current point is:   (',num2str([xnew(1)]),',',num2str([xnew(2)]),')']);



i=0;
xs=[];
while (norm(gradnew,2)>1e-10)
  i=i+1;  %% count the iterations
  xs=[xs xnew];   
  %% compute the step length adjustment
  tstep=norm(gradnew)^2/(gradnew'*Q*gradnew);
  %% compute the new point using the multiple of the gradient at the
  %% current point
  xnew=xnew-tstep*gradnew;
  %% compute new gradient
  gradnew=Q*xnew;
  %% plot the new point and the line path between the new and previous
  %% points
  plot(xs(1,:),xs(2,:),xnew(1),xnew(2),'xk');
  disp(['iteration ',num2str(i),', the current point is:   (',num2str([xnew(1)]),',',num2str([xnew(2)]),')' ]);
end
hold off