% minimum beräknar minimum med gradientmetoden till
% funktionen fun 

%fun = @(x,y) (7*x + 15* y).*exp(-(0.01*x.^2 + 0.02* y.^2));
%[X,Y]=meshgrid(linspace(-6,3,100));figure, mesh(X,Y,fun(X,Y)),hold on
% Funktionens gradient beräknas numeriskt.

%funx=@(x,y) (fun(x+10^(-6),y) - fun(x-10^(-6),y))/(2*10^(-6));
%funy = @(x,y) (fun(x,y+10^(-6)) - fun(x,y-10^(-6)))/(2*10^(-6));
%För att programmet skall fungera skall fun, funx och funy definieras först.
%Av okänd anledning kan inte de senare definieras i skriptifilen. 
%Gör man det i kommandofönstret är det inga problem.
% (x0,y0) är startpunkten som anges här:
x0 = 1; y0= 1;

 


z0=fun(x0,y0); % Funktionsvärdet i startpunkten
plot3(x0,y0,z0,'*')

% while-snurrans kontrollvariabel q skall vara 1 om avståndet mellan ,
% den senast beräknade funktionsvärdet och den näst senast beräknade 
% är större än toleransen tol som anges här:

tol = 10^-3;

% När avståndet är mindre än tol skall q sättas till 0.
% För att starta snurran sätts q=1.

q=1;
while q

% beräknar  gradienten i (x0,y0).

   px =funx(x0,y0); py=funy(x0,y0);

   % beräknar koordinater för tre punkter i riktningen -gradf(x0,y0) från (x0,y0). 


   r=[0.5, 1, 2];

   xt=x0-r*px;yt=y0-r*py;zt=zeros(1,3); 
    

   for i=1:3
      zt(i)=fun(xt(i),yt(i)); % Motsvarande z-värden.
   end

   % Välj ut bästa punkten (xt1,yt1,zt1) och motsvarande r-värde
   
   k=find(zt==min(zt));
   xt1=xt(k);yt1=yt(k);zt1=zt(k);r=r(k); 

   % Leta efter ännu bättre punkt  (xt2,yt2,zt2) genom att testa med potenser av detta r.

   f=r;
   r=f*r;
   xt2=x0-r*px;yt2=y0-r*py;zt2=fun(xt2,yt2);
% Om denna punkt ger bättre värde upprepas proceduren. 
   
      while (zt2 <zt1)&abs(zt2-zt1)>tol
       xt1=xt2;yt1=yt2;zt1=zt2; 
       r=f*r;
       xt2=x0-r*px;yt2=y0-r*py;zt2=fun(xt2,yt2);
      end
      x1=xt1;y1=yt1;z1=zt1;
   
% Vi har nu hittat en bra punkt (x1,y1) i riktningen -gradf(x0,y0) 
   % men undersöker först om den är bättre än punkten (x0,y0)
   % 
   if z1 > z0  q=0;
   else
   plot3([x0 x1],[y0 y1],[z0 z1],'r'), plot3(x1,y1,z1,'r*'),pause,hold on

   d=abs(z1-z0); %Normalt är z1 < z0 men om man väljer liten tol
   q=(d>tol);
x0=x1;y0=y1;z0=z1;
 end
end