function [P0,y0,h,err,P] = grads(Fn,Gn,P0,max1,delta,epsilon,show)
% [P0,y0,P] = grads(Fn,Gn,P0,max1,delta,epsilon,show)
% Gradient search for a minimum.
% Fn is the function, input.
% Gn is the gradient for Fn, input.
% P0 is the starting point, input.
% max1 is the maximum number of iterations, input
% delta is the tolerance, input;
% epsilon is the tolerance, input;
% show if show==1 the iterations are displayed, input.
% P0 is the point for the minimum, output.
% y0 is the function value Fn(V0), output.
% dV is the size of the final simplex, output.
% P is a vector containing the iterations, output.
if nargin==5, show = 0; end
[mm n] = size(P0);
maxj = 20;
big = 1e8;
h = 1;
len = norm(P0);
y0 = Fn(P0);
if (len>1e4), h = len/1e4; end
err = 1;
cnt = 0;
cond = 0;
while (cnt<max1 & cond~=5 & (h>delta | err>epsilon))
S = feval(Gn,P0); echo off;
P1 = P0 + h*S;
P2 = P0 + 2*h*S;
y1 = feval(Fn,P1);
y2 = feval(Fn,P2);
cond = 0;
j = 0;
while (j<maxj & cond==0)
len = norm(P0);
if (y0<y1),
P2 = P1;
y2 = y1;
h = h/2;
P1 = P0 + h*S;
y1 = feval(Fn,P1);
else
if (y2<y1),
P1 = P2;
y1 = y2;
h = 2*h;
P2 = P0 + 2*h*S;
y2 = feval(Fn,P2);
else
cond = -1;
end
end
j = j+1;
if (h<delta), cond=1; end
if (abs(h)>big | len>big), cond=5; end
end
if (cond==5),
Pmin = P1;
ymin = y1;
else
d = 4*y1 - 2*y0 - 2*y2; % Start of a long block:
if (d<0),
hmin = h*(4*y1-3*y0-y2)/d;
else
cond = 4;
hmin = h/3;
end
Pmin = P0 + hmin*S;
ymin = feval(Fn,Pmin);
h0 = abs(hmin);
h1 = abs(hmin-h);
h2 = abs(hmin-2*h);
if (h0<h), h = h0; end
if (h1<h), h = h1; end
if (h2<h), h = h2; end
if (h==0), h = hmin; end
if (h<delta), cond=1; end
e0 = abs(y0-ymin);
e1 = abs(y1-ymin);
e2 = abs(y2-ymin);
if (e0~=0 & e0<err), err = e0; end
if (e1~=0 & e1<err), err = e1; end
if (e2~=0 & e2<err), err = e2; end
if (e0==0 & e1==0 & e2==0), err = 0; end
if (err<epsilon), cond=2; end
if (cond==2 & h<delta), cond=3; end
end % End of the long block.
cnt = cnt+1;
P0 = Pmin;
y0 = ymin;
home;
if cnt==1,
diary output,...
disp('The results for a Gradient search.'),...
disp(' p q f(p,q)'),...
diary off,clc;
end;
if show==1,
Mx1 = 'Gradient search iteration No. ';
Mx2 = ' p q f(p,q)';
disp([Mx1,int2str(cnt)]),disp(Mx2),...
diary output,disp([P0,y0]),diary off;
pause(0.5);
plot(P0(1),P0(2),'or');
end
end