//
// calcul du minimum de la fonction de Rosenbrock
//  par la méthode du gradient et de Newton
//  comparaison des vitesses de convergence
//
function y=R(x)
    y=10*(x(2)-x(1)^2)^2+(x(1)-1)^2
endfunction
//
function y=nablaR(x)
    y=[-40*x(1)*(x(2)-x(1)^2)+2*(x(1)-1);20*(x(2)-x(1)^2)]    
endfunction
function y=HR(x)
    y=[-40*x(2)+120*x(1)^2+2,-40*x(1);-40*x(1),20];    
endfunction
//
x1=linspace(-1,2,100)
x2=linspace(-1,2,100)
z=zeros(100,100)
//
for i=1:100
    for j=1:100
     z(i,j)=R([x1(i);x2(j)])
    end    
end
figure(1)
clf()
contour(x1,x2,z,[0:0.1:1,2:10,20,40])
//
function [Utot,u, nb_iter, residu] = GPF2 (J,nablaJ,u0,epsilon,max_iter,rho)
u=u0
Utot=u;
nb_iter=0;
residu=norm(nablaJ(u))
while (nb_iter<max_iter)&(residu>epsilon)
  nb_iter=nb_iter+1
  u=u-rho*(nablaJ(u))
  Utot=[Utot,u]
  residu=norm(nablaJ(u))
end
endfunction
//////
function [Utot,u,nb_iter, residu] = Newton(J,nablaJ,HJ,u0,epsilon,max_iter)    
u=u0
Utot=u;
nb_iter=0;
residu=norm(nablaJ(u))
while (nb_iter<max_iter)&(residu>epsilon)
  nb_iter=nb_iter+1
  u=u-inv(HJ(u))*nablaJ(u)
  Utot=[Utot,u]
  residu=norm(nablaJ(u))
end
endfunction
//
//
u0=rand(2,1);
epsilon=1E-4;
max_iter=20000;
rho=0.02;
[Utot1,u, nb_iter1, residu] = GPF2 (R,nablaR,u0,epsilon,max_iter,rho)
disp('minimum GPF2');disp(u)
disp('precision GPF2');disp(residu)
disp('nb iterations GPF2');disp(nb_iter1)
//
plot(Utot1(1,:),Utot1(2,:),'b*')
//
[Utot2,u, nb_iter2, residu] = Newton(R,nablaR,HR,u0,epsilon,max_iter)
disp('minimum Newton');disp(u)
disp('precision Newton');disp(residu)
disp('nb iterations Newon');disp(nb_iter2) 
//
plot(Utot2(1,:),Utot2(2,:),'black-')
//
//
// convergence de l'approximation vers la solution exacte
//
figure(2)
clf()
err1=Utot1-ones(Utot1);
err2=Utot2-ones(Utot2);
plot((1:20),log(sqrt(err1(1,1:20).^2+err1(2,1:20).^2)),'b*')
plot(1:(nb_iter2+1),log(sqrt(err2(1,:).^2+err2(2,:).^2)),'black-')
