<Text-field layout="Heading 1" style="Heading 1"><Font family="Times New Roman"> Gradient</Font></Text-field>

Rendons d'abord disponibles les macro-commandes des extensions linalg, plottools et plots.

with ( linalg ): with ( plottools ): with ( plots ):

Consid\351rons la fonction g d\351finie par NiMvLSUiZ0c2JCUieEclInlHLCgqJkYnIiIiLSUkZXhwRzYjLCYqJCwmRidGK0YrISIiIiIjRjIqJCwmRihGKyQiI0QhIiNGK0YzRjJGK0YrKiZGKEYrLUYtNiMsJiokRidGM0YyKiRGKEYzRjJGK0YrIiIkRis=.

g:=(x,y)->x*exp(-(x-1)^2-(y+0.25)^2)+(y)*exp(-x^2-y^2)+3;

Tra\347ons la surface d'\351quation NiMvJSJ6Ry0lImdHNiQlInhHJSJ5Rw== sur le pav\351 [-2; 3,5] x [-2,5; 2].

Surface:=plot3d([x,y,g(x,y)],x=-2..3.5,y=-2.5..2,grid=[30,30]): display(Surface,labels=[x,y,z], lightmodel=light2,style=hidden,shading=zhue,color=navy, axes=framed,font=[TIMES,ROMAN,8], orientation=[-160,60]);

Cette surface poss\350de des extremums relatifs. Tra\347ons alors des courbes de niveau pour mieux estimer les valeurs x et y qui permettent d'atteindre ces extremums.

Courbes_niveau:=contourplot(g(x,y),x=-2..3,y=-2.5..2,contours=15,numpoints=1300,color=orange): Courbes_niveau;

Les courbes de niveau nous montrent que le maximum est atteint avec (x,y) au voisinage du point (1,4; -0,3) et que le minimum est atteint avec (x,y) dans le voisinage du point (-0,2; -0,7).

Visualisons le fait que le gradient en un point est toujours perpendiculaire \340 la courbe de niveau passant par ce point. Obtenons d'abord un champ de vecteurs gradient. La macro-commande gradplot de l'extension plots sera utilis\351e.

gp:=gradplot(g(x,y),x=-2..3.5,y=-2.5..2,grid=[25,25],arrows=slim,color=navy): gp;

Superposons ensuite ce champ de vecteurs au trac\351 des courbes de niveau.

display([Courbes_niveau,gp]);

On observe effectivement que les fl\350ches sont toujours perpendiculaire aux courbes de niveau. S\351lectionnez le graphique pr\351c\351dent et cliquez sur la plus grosse loupe (ou utilisez le raccourcis-clavier ctrl-n) pour agrandir le graphique. Cette illustration nous indique la direction avec laquelle la fonction poss\350de le plus grand taux de variation. Le sens de ces vecteurs confirme qu'il se pr\351sente effectivement un maximum relatif (fl\350ches rentrantes) dans le voisinage du point (1,4; -0,3) et un minimum relatif (fl\350ches sortantes) dans le voisinage du point (-0,2; -0,7).

Pour rechercher les couples (x,y) avec lesquels la fonction g pr\351sentera un extremum relatif, appliquons d'abord le test d'extremum local. En fait, nous nous contenterons de trouver les points critiques et sur la base du trac\351 de la surface, nous identifierons la nature de chacun de ces points.

gx:=D[1](g): gy:=D[2](g): Systeme:={gx(x,y),gy(x,y)};

Sol_symboliques:=solve({gx(x,y)=0,gy(x,y)=0,x<infinity},{x,y});

Il fallait s'y attendre, l'\351valuateur n'a trouv\351 aucune solution \340 ce syst\350me. Plus souvent qu'autrement, il faut plut\364t r\351soudre num\351riquement de tels syst\350mes d'\351quations. R\351solvons d'abord num\351riquement ce syst\350me au voisinage de (1,4; -0,3) avec la macro-commande fsolve.

xy_max:=fsolve({gx(x,y)=0,gy(x,y)=0},{x,y},{x=0..1.4,y=-0.3..0});

R\351solvons ensuite num\351riquement le syst\350me au voisinage de (-0,2; -0,7).

xy_min:=fsolve({gx(x,y)=0,gy(x,y)=0},{x,y},{x=-0.2..0,y=-0.7..-0.5});

Contr\364lons si chaque r\351sultat annule effectivement le gradient.

Gradient:=[gx(x,y),gy(x,y)];

evalf(subs(xy_max,Gradient)); evalf(subs(xy_min,Gradient));

Les valeurs des composantes sont presque nulles. Obtenons les valeurs de ces extremums relatifs.

z_max:=evalf(subs(xy_max,g(x,y))); z_min:=evalf(subs(xy_min,g(x,y)));

<Text-field layout="Heading 1" style="Heading 1"><Font family="Times New Roman"> Gradient et extremum relatif</Font></Text-field>

Commen\347ons par rechercher un maximum relatif de la fonction (en supposant qu'il existe, bien s\373r). Rappelons que le gradient en un point indique la direction du plus grand taux de variation de la fonction. Alors, avec un point NiMvJiUiUEc2IyIiITckJiUieEdGJiYlInlHRiY= de d\351part, nous allons nous d\351placer l\351g\350rement, dans la direction donn\351e par le gradient en ce point, vers un deuxi\350me point NiMvJiUiUEc2IyIiIjckJiUieEdGJiYlInlHRiY=. Si, avec le point NiMmJSJQRzYjIiIh de d\351part, on n'atteint pas un maximum de la fonction, ce deuxi\350me point nous permettra par contre de se rapprocher du point qui permet \340 la fonction d'atteindre le maximum recherch\351.

Soit donc l'\351quation vectorielle NiMvJSNPUEcsJiZGJDYjIiIhIiIiKiYlJ2xhbWJkYUdGKS0tJSVncmFkRzYjJSJnRzYkJiUieEdGJyYlInlHRidGKUYp de la droite passant par le point NiMmJSJQRzYjIiIh et ayant NiMtLSUlZ3JhZEc2IyUiZ0c2JCYlInhHNiMiIiEmJSJ5R0Yr pour vecteur directeur. La valeur de NiMlJ2xhbWJkYUc=>0 d\351terminera l'importance du d\351placement dans le direction de NiMtLSUlZ3JhZEc2IyUiZ0c2JCYlInhHNiMiIiEmJSJ5R0Yr. Avec une valeur NiMlJ2xhbWJkYUc=>0 donn\351e, on d\351duira un point particulier de cette droite. Ce point particulier, qu'on nommera NiMmJSJQRzYjIiIi, nous permettra de se rapprocher du point (x,y) permettant d'atteindre le maximum cherch\351.

Voyons maintenant comment appliquer ce raisonnement. Comme les calculs seront num\351riques, nous consid\350rerons que le point P(x,y) permet d'atteindre un maximum relatif de la fonction si le gradient en ce point est presque \351gal au vecteur [0,0]. Un fa\347on \351l\351gante de tester ce fait, c'est de v\351rifier si la longueur (la norme) de ce vecteur gradient est presque nulle.

Pour les besoins de la cause, le point P(-2,2) sera d\351sign\351 point de d\351part.

P[0]:=[-2,2];

Contr\364lons (ici, par principe) si ce point permet d'atteindre le maximum de la fonction. Utilisons une tol\351rance de NiMpIiM1LCQiIikhIiI=.

tolerance:=0.00000001; Grad_P:=evalf(subs({x=P[0][1],y=P[0][2]},Gradient)); Grad_longueur:=norm(Grad_P,2); evalb(Grad_longueur<tolerance);

On a utilis\351 la macro-commande norm de l'extension linalg pour calculer la longueur du gradient au point P(-2,2). De plus, rappelons que evalb commande une \351valuation bool\351enne. Le r\351sultat vrai pour l'in\351galit\351 Grad_longueur<tolerance signifiera que le point P(x,y) est, au niveau de tol\351rance admis, suffisamment pr\350s du point permettant d'avoir le maximum cherch\351.

Obtenons maintenant un deuxi\350me point avec un pas NiQvJSdsYW1iZGFHIiIhIiIi.

lambda:=0.1;

Obtenons le point NiMmJSJQRzYjIiIi.

Grad_P:=evalf(subs({x=P[0][1],y=P[0][2]},Gradient)); P[1]:=P[0]+lambda*Grad_P;

Contr\364lons si le point NiMmJSJQRzYjIiIi permet d'atteindre le maximum recherch\351 ( toujours avec la tol\351rance de NiMpIiM1LCQiIikhIiI=).

Grad_P:=evalf(subs({x=P[1][1],y=P[1][2]},Gradient)); Grad_longueur:=norm(Grad_P,2); evalb(Grad_longueur<tolerance);

Il sera beaucoup plus commode d'utiliser une boucle de calcul pour continuer de se d\351placer de point en point. Prenons soins de limiter le nombre d'it\351rations \340 1000 et rappelons les param\350tres qui seront utilis\351s dans cette boucle.

k:=1000; lambda:=0.1; tolerance:=0.00000001;

L'algorithme de cette boucle est le suivant.

1) \311valuer le gradient au point NiMmJSJQRzYjJSJuRw==.

2) Calculer la longueur (la norme) de ce vecteur gradient.

3) Contr\364ler si cette longueur est inf\351rieur \340 la longueur tol\351r\351e.

4) Si tel est le cas, terminer les it\351rations

5) Sinon, calculer un prochain point NiMmJSJQRzYjLCYlIm5HIiIiRihGKA==.

P[0]:=[-2,2]; # Point de d\351part

for n from 0 to k do Grad_P:=evalf(subs({x=P[n][1],y=P[n][2]},Gradient)); Grad_longueur:=norm(Grad_P,2); if Grad_longueur<tolerance then break else P[n+1]:=P[n]+lambda*Grad_P; fi; od: Grad_P:=evalf(subs({x=P[n][1],y=P[n][2]},Gradient)); norm(Grad_P,2); Nombre_iterations:=n; 'P[n]'=P[n]; n:='n':

Comparons ce point avec celui xy_max obtenu \340 la section pr\351c\351dente.

xy_max;

Tra\347ons, dans le plan, sans les relier, le chemin de points ( le point de d\351part plus les 812 points obtenus par les it\351rations pr\351c\351dentes).

Courbe:=pointplot({seq(P[n],n=0..Nombre_iterations)},color=navy): Courbe;

On observe tr\350s bien sur le trac\351 pr\351c\351dent, que plus le taux de variation de la fonction g est important, plus la distance s\351parant les points augmente. Cela d\351coule du fait que le pas de d\351placement NiMlJ2xhbWJkYUc= est constant et que le taux de variation de la fonction g ne l'est pas.

Illustrons ce fait en r\351alisant une animation du trac\351 de ce chemin de points. Pour animer la partie significative du chemin de point, cette animation sera construite \340 partir du 600 i\350me point jusqu'au 750 i\350me. Nous constaterons qu'avec un pas NiMlJ2xhbWJkYUc= constant, plus la distance entre les points augmente, plus l'animation s'acc\351l\350re.

Animation:=seq(pointplot([seq(P[n],n=0..k)],color=navy),k=600..750):

display(Animation,insequence=true);

Superposons maintenant le trac\351 de ces points \340 ceux du champ de vecteurs gradient et des courbes de niveau. On aura alors une id\351e claire du chemin \340 parcourir dans le domaine, \340 partir du point (-2, 2), pour se rapprocher du point (x,y) qui permet d'atteindre le maximum relatif de la fonction. Notez que le chemin obtenu est toujours, \351videmment, perpendiculaire \340 toutes les courbes de niveau passant par ces points.

display([Courbe,Courbes_niveau,gp]);

Obtenons le trac\351 de ce chemin dans l'espace en ajoutant, \340 chacun des points calcul\351s, une troisi\350me composante (cote) dont la valeur sera celle de l'image par la fonction g des deux premi\350res composantes. Relions cette fois-ci les points en pr\351cisant l'option connect=true.

Courbe3d:=pointplot3d([seq([P[n][1],P[n][2],evalf(g(P[n][1],P[n][2]))],n=1..Nombre_iterations)],connect=true,thickness=3,axes=boxed,color=navy,orientation=[-160,60]): Courbe3d;

Superposons finalement le trac\351 de cette courbe \340 celui de la surface. L'effet est spectaculaire.

display([Surface,Courbe3d],labels=[x,y,z], lightmodel=light2,style=hidden,shading=zhue,color=navy, axes=framed,font=[TIMES,ROMAN,8], orientation=[-190,50]);

Pour un effet davantage spectaculaire, animons cette superposition. (Il vous faudra patienter environ 1 min 38 sec avec un PII-400)

Animation:=seq(display([Surface,pointplot3d([seq([P[n][1],P[n][2],evalf(g(P[n][1],P[n][2]))],n=0..k)],color=navy,thickness=3,connect=true)]),k=600..750):

display(Animation,insequence=true, labels=[x,y,z], lightmodel=light2,style=hidden,shading=zhue,color=navy, axes=framed,font=[TIMES,ROMAN,8], orientation=[-190,50]);