// G. Chiavassa juin 2012
//methodes de recherche de minimum
// J(x,y)= 2x^3-3 x^2 + 2 y^3 + 3 y^2
// pts critiques: (0,0) (0,-1) (1,-1) (1,0)
// (0,-1) est un max et (1,0) est un min
// les autres pts selles

deff('z=J(x,y)','z=2*x.^3- 3*x.^2 + 2*y.^3 + 3* y.^2');

deff('[z]=grad(x,y)',['z(1)=6*x.^2-6*x';'z(2)=6*y.^2+6*y']);

deff('[z]=grad_approx(x,y,delta)',['z(1)=(J(x+delta,y)-J(x,y))/delta';'z(2)=(J(x,y+delta)-J(x,y))/delta']);
deff('[z]=grad_approxO2(x,y,delta)',['z(1)=(J(x+delta,y)-J(x-delta,y))/2/delta';'z(2)=(J(x,y+delta)-J(x,y-delta))/2/delta']);
xset("window",0)
clf();
//fplot3d(-0.5:0.1:1.5,-1.5:0.1:0.5,J);
xset("window",1)
clf();
contour(-0.5:0.1:1.5,-1.5:0.1:0.5,J,10);

pas=0.03;
delta=0.0001;
B=[1;0];
N=200;
M0=[0.;-0.5];
M=M0;
M1=M0;
M2=M0;

for i=1:N
 dM=grad(M(1),M(2));
 dM1=grad_approx(M(1),M(2),delta);
 dM2=grad_approxO2(M(1),M(2),delta);
 M=M-pas*dM;
 M1=M1-pas*dM1;
 M2=M2-pas*dM2;
 err(i)=norm(M-B,2);
 err_approx(i)=norm(M1-B,2);
 err_approxO2(i)=norm(M2-B,2);
end;

xset("window",2)
clf();
plot2d([1:N],err,3,logflag="nl");
plot2d([1:N],err_approx,logflag="nl");
plot2d([1:N],err_approxO2,2,logflag="nl");
M
M1
M2

function [f,g,ind] = cost(x,ind)
  xref=[1;2;3];
  f=2*x(1)^3- 3*x(1)^2 + 2*x(2)^3 + 3*x(2)^2;
  g=grad(x(1),x(2));
endfunction

[f,xopt]=optim(cost,M0,'gc',imp=2);
norm(xopt-B,2)

