## Recherche de minimum de fonction à deux arguments par la méthode de Newton

On souhaite fabriquer une fonction qui recherche le minimum d'une fonction par la méthode de Newton.

Pour cela, prenons déjà une fonction F : R^2 -> R 
                                         (x,y) -> F(x,y)

Les conditions nécessaires d'extréma (F est supposée différentiable 2x) sont, 
pour (xinf,yinf) l'argument minimisant de F,

dF/dx(xinf,yinf) = 0 ; dF/dy(xinf,yinf) = 0

Il s'agit donc de résoudre ces deux équations et on veut que ce soit en utilisant la méthode de 
Newton. 

### Analyse

Si F est connue, trouver ses dérivées partielles

Choisir une expression simple comme x^4+y^4+x*y par exemple.

F:x^4+y^4+x*y;
Fx:diff(F,x);Fy:diff(F,y);Fxx:diff(Fx,x);Fyy:diff(Fy,y);Fxy:diff(Fx,y);

Calculer l'itération de Newton : le passage de (x0,y0) à (x1,y1)

pour les équations dF/dx(x,y)=0 et dF/dy(x,y)=0

[Fx0,Fy0,Fxx0,Fyy0,Fxy0]:ev([Fx,Fy,Fxx,Fyy,Fxy],x=x0,y=y0);
eq:[Fxx0*(x1-x0)+Fxy0*(y1-y0)+Fx0=0,Fxy0*(x1-x0)+Fyy0*(y1-y0)+Fy0=0];
sol:solve(eq,[x1,y1]);

Fabriquer une fonction réalisant cette itération.

On peut utiliser 'define' plutôt que ':=' pour placer des résultats de calculs dans le corps de la fonction

Voir la différence avec les exemples ci-dessous

res:x^2/(1+x);
f(x):=res;
fundef(f);

res:x^2/(1+x);
define(f(x),res);
fundef(f);

define(it(x0,y0),ev([x1,y1],sol));

Essayer la fonction : d'abord sans s'occuper du critère d'arrêt et ensuite en l'introduisant

block([x0:1,y0:1,n:0,N:10],print(n,[x0,y0]),
 while (n<N) do block([],
  [x0,y0]:float(it(x0,y0)),n:n+1,
  print(n,[x0,y0])));

block([x0:1,y0:1,n:0,N:100,x1,y1,eps:1.0e-2],print(n,[x0,y0]),
[x1,y1]:[x0,y0]+2*eps*[1,1],
 while (n<N) and sqrt((x1-x0)^2+(y1-y0)^2)>eps do block([],
  [x0,y0]:[x1,y1],
  [x1,y1]:float(it(x0,y0)),n:n+1,
  print(n,[x0,y0])));

### Synthèse 

Il s'agit de réunir les étapes précédentes

F:x^4+y*sin(x)+y^4-x*cos(y)+x*y;
Fx:diff(F,x)$Fy:diff(F,y)$Fxx:diff(Fx,x)$Fyy:diff(Fy,y)$Fxy:diff(Fx,y)$
[Fx0,Fy0,Fxx0,Fyy0,Fxy0]:ev([Fx,Fy,Fxx,Fyy,Fxy],x=x0,y=y0)$
eq:[Fxx0*(x1-x0)+Fxy0*(y1-y0)+Fx0=0,Fxy0*(x1-x0)+Fyy0*(y1-y0)+Fy0=0]$
sol:solve(eq,[x1,y1])$
define(it(x0,y0),ev([x1,y1],sol))$
block([x0:1,y0:1,n:0,N:100,x1,y1,eps:1.0e-2],
[x1,y1]:[x0,y0]+2*eps*[1,1],
 while (n<N) and sqrt((x1-x0)^2+(y1-y0)^2)>eps do block([],
  [x0,y0]:[x1,y1],
  [x1,y1]:float(it(x0,y0)),n:n+1),
if n<N then print("argument minimisant : ",[x1,y1]) 
else print("non convergé : "))$

Et de les réunir dans une fonction globale 

Nwt(F,x,y,x00,y00,eps,N):=block([Fx,Fy,Fxx,Fxy,Fx0,Fy0,Fxx0,Fyy0,Fxy0,eq,sol,it,x0,y0,x1,y1],
 Fx:diff(F,x),Fy:diff(F,y),Fxx:diff(Fx,x),Fyy:diff(Fy,y),Fxy:diff(Fx,y),
 [Fx0,Fy0,Fxx0,Fyy0,Fxy0]:ev([Fx,Fy,Fxx,Fyy,Fxy],x=x0,y=y0),
 eq:[Fxx0*(x1-x0)+Fxy0*(y1-y0)+Fx0=0,Fxy0*(x1-x0)+Fyy0*(y1-y0)+Fy0=0],
 sol:solve(eq,[x1,y1]),
 define(it(x0,y0),ev([x1,y1],sol)),
 block([x0:x00,y0:y00,n:0,x1,y1],
  [x1,y1]:[x0,y0]+2*eps*[1,1], 
  while (n<N) and sqrt((x1-x0)^2+(y1-y0)^2)>eps do block([],
   [x0,y0]:[x1,y1],
   [x1,y1]:float(it(x0,y0)),n:n+1), 
  if n<N then return([x1,y1]) else return("echec")))$

Nwt(x^4+y*sin(x)+y^4-x*cos(y)+x*y,x,y,1,1,1.0e-2,100);

### Fonction d'un nombre d'arguments quelconque

Voici une fonction qui donne l'itération de Newton pour un nombre quelconque d'arguments

ItNwt(F,X):=block([G,H,X1,sol],
 X1:map(lambda([%x],concat(%x,1)),X),
 G:map(lambda([%x],diff(F,%x)),X),
 H:map(lambda([%g],map(lambda([%x],diff(%g,%x)),X)),G),
 sol:solve(map(lambda([%h,%g],%h.(X1-X)+%g),H,G),X1),
 ev(X1,sol))$

ItNwt(x^4+y^4+z^4,[x,y,z]);

Fabriquer une fonction globale comme dans 1.2

Nwt(F,X,X0,eps,N):=block([it:ItNwt(F,X),X2X0,X1,n:0,d10],  
  X1:X0+2*eps*makelist(1,i,1,length(X)), d10:sqrt((X1-X0).(X1-X0)),
  while (n<N) and d10>eps do block([],
   X0:X1,X2X0:map(lambda([%x,%x0],%x=%x0),X,X0),
   X1:float(ev(it,X2X0),n:n+1),
   d10:sqrt((X1-X0).(X1-X0))),
  if n<N then return(X1) else return("echec"))$

Nwt(x^4+y^2+z^2,[x,y,z],[1,1,1],1.0e-5,100);