CoCalc Public FilesComplementsSolD5.sagewsOpen with one click!
Author: Juan Carlos Bustamante
Views : 55
Compute Environment: Ubuntu 18.04 (Deprecated)

Descente de gradient

Dans les pages 108 à 110 du livre de Stewart on présente une méthode de minimisation, connue sous le nom de "descente de gradient". En gros, on part d'un point puis on fait des pas dans la direction oposée au gradient. La longueur du pas est celle qui permet de minimiser la valeur de la fonction objectif (si elle existe...). Ci après, une implémentation très naïve de la chose, en SAGEMath.

L'implémentation ci bas fait appel à la fonction find_local_minimum(...) qui retourne une liste de valeurs.

L'exemple utilisé est celui de la page 109 - 110.

var('x,y') f(x,y) = x^4+y^2 -2*x^2 *y+2*y+x
(x, y)
var('t') def myDescente(f, prec, x0,y0):# prec = précision, x0, y0 sont les coordonnées du point de départ. Df = f.gradient() print "Point d'évalution :", (x0,y0) print "Valeur de la fonction objectif:", f(x0,y0).n(digits=6), "Norme du gradient = ", norm(Df(x0,y0)).n(digits=6) if norm(Df(x0,y0)) < prec : print "Valeur minimale trouvée. Arrêt." else: F(t) = f(x0-t*Df(x0,y0)[0],y0-t*Df(x0,y0)[1]) t0 = F.find_local_minimum(0,2)[1] print "On va ailleurs. Longueur du pas (en fraction du gradient) : ", t0 myDescente(f,prec,x0-t0*Df(x0,y0)[0], y0-t0*Df(x0,y0)[1])
tt
myDescente(f,0.01, 0,0)
Point d'évalution : (0, 0) Valeur de la fonction objectif: 0.000000 Norme du gradient = 2.23607 On va ailleurs. Longueur du pas (en fraction du gradient) : 0.380408919629 Point d'évalution : (-0.38040891962887613, -0.7608178392577523) Valeur de la fonction objectif: -1.08206 Norme du gradient = 0.422488 On va ailleurs. Longueur du pas (en fraction du gradient) : 0.345102226762 Point d'évalution : (-0.2499999964962533, -0.8260222923108163) Valeur de la fonction objectif: -1.11257 Norme du gradient = 0.249272 On va ailleurs. Longueur du pas (en fraction du gradient) : 0.301870420298 Point d'évalution : (-0.28365182325811644, -0.8933259383163823) Valeur de la fonction objectif: -1.12205 Norme du gradient = 0.117240 On va ailleurs. Longueur du pas (en fraction du gradient) : 0.320913060169 Point d'évalution : (-0.2499999930467315, -0.910151862449617) Valeur de la fonction objectif: -1.12425 Norme du gradient = 0.0611523 On va ailleurs. Longueur du pas (en fraction du gradient) : 0.302753560077 Point d'évalution : (-0.25827974829437783, -0.9267113565644658) Valeur de la fonction objectif: -1.12482 Norme du gradient = 0.0294276 On va ailleurs. Longueur du pas (en fraction du gradient) : 0.314569974107 Point d'évalution : (-0.24999999981318474, -0.9308512327243619) Valeur de la fonction objectif: -1.12496 Norme du gradient = 0.0148671 On va ailleurs. Longueur du pas (en fraction du gradient) : 0.30296365379 Point d'évalution : (-0.25201433489340147, -0.93487990243502) Valeur de la fonction objectif: -1.12499 Norme du gradient = 0.00719511 Valeur minimale trouvée. Arrêt.

Voyons maintenant l'exercice n.31, section 3.1 (celui du devoir).

D'abord avec le point de départ (1/4,1/4)(1/4, 1/4), le critère d'arêt est que la norme du gradient soit inférieure à 10610^{-6}.

g(x,y) = x^2+ y^2 - (x^2+y^2)^(3/2) myDescente(g, 10^-6, 0.25,0.25)
Point d'évalution : (0.250000000000000, 0.250000000000000) Valeur de la fonction objectif: 0.0808058 Norme du gradient = 0.332107 On va ailleurs. Longueur du pas (en fraction du gradient) : 1.06457745783 Point d'évalution : (-1.59089003082258e-9, -1.59089003082258e-9) Valeur de la fonction objectif: 5.06186e-18 Norme du gradient = 4.49972e-9 Valeur minimale trouvée. Arrêt.

À partir du point (1,1)

Les choses sont très différentes : la fonction h(t)=f((1,1)+tf(1,1))h(t) = f((1,1) + t \nabla f(1,1)) n'a pas de minimums locaux, la descente ne s'arête jamais. L'algorithme ne converge pas (on peut essayer de l'exécuter, voir ce que ça donne)

#myDescente(g, 10^-6, 1,1) Dg = g.gradient() h(t) = g(1 - t* Dg(1,1)[0], 1 - t* Dg(1,1)[1]) plot(h,0,1)

Voyons la surface, on comprendra ce qui est arrivé.

Surf = plot3d(g,(x,-1.1,1.1), (y,-1.1,1.1), mesh=0.7) show(Surf)
3D rendering not yet implemented

Multiplicateurs de Lagrange (ex. 38, section 3.3)

On doit maximiser f(x,y)=2x+3yf(x,y) = 2x+3y avec la contrainte g(x,y)=x+y=5g(x,y) = \sqrt{x} + \sqrt{y} = 5. Voyons quelques courbes de niveau de la fonction objectif, ainsi que la courbe de contrainte

var('x,y') C = contour_plot(2*x+3*y, (x,0,15), (y,0,15), cmap = "Blues", colorbar=True, fill=True) Cg = plot((5-sqrt(x))^2, (x,0,15),ymin = 0, ymax = 15, color = "red" ) show(Cg+C)
(xx, yy)

Exercice n.25, page 246

Il s'agit de donner des approximations de la solution de y(x)=6x23x2yy'(x) = 6x^2 -3x^2 y, avec y(0)=3y(0)= 3. On utilise la méthode d'Euler.

V=[3]
var('x,y') def F(x,y) : return 6*x^2 - 3 * x^2 * y #Définition de la fonction h = 0.0001 #Le pas N = 1/h+1 #Le nombre de points x0 = 0 #L'abscisse du point de départ y0 = 3 #L'ordonnée du point de départ X = [x0 + j*h for j in range(N)] # Les abscisses des points Y = [y0] for j in range(N) : Y = Y + [Y[-1] + h * F(X[j],Y[-1]) ] #La boucle pour calculer les ordonnées Points = [(X[j], Y[j]) for j in range(N)] Points[-1] V = V + [Points[-1][1]] # Juste pour avoir la liste des valeurs estimées de y(1). Pas très élégant comme code V
(xx, yy)
(1.000000000000001.00000000000000, 2.367901516308502.36790151630850)
[33, 33, 2.392793913961662.39279391396166, 2.370110818325782.37011081832578, 2.368100406247932.36810040624793, 2.367901516308502.36790151630850]
f(x)= 2+e^(-x^3)

On évalue maintenant l'erreur commise entre les approximations et la "vraie valeur"

[(v-f(1)).n(digits=10) for v in V]
[0.63212055880.6321205588, 0.63212055880.6321205588, 0.024914472790.02491447279, 0.0022313771550.002231377155, 0.00022096507750.0002209650775, 0.000022075135350.00002207513535]

La question n. 26, section 6.3

L'équation différentielle s'écrit yeyy=xx2+1y e^y y' = x\sqrt{x^2 +1}, il s'agit d'une équation à variables séparables. On peut intégrar sans trop de mal et trouver la solution exacte. Faisons le avec SAGE

var('x,y') integrate(y*e^y,y) integrate(x*sqrt(x^2 +1),x)
(xx, yy)
(y1)ey{\left(y - 1\right)} e^{y}
13(x2+1)32\frac{1}{3} \, {\left(x^{2} + 1\right)}^{\frac{3}{2}}

Les courbes solution sont donc les courbes (y1)ey=13(x2+1)3/2+C(y-1)e^y = \frac{1}{3}(x^2 + 1 )^{3/2} + C dessinons-en queleques unes

contour_plot((y-1)*e^y - (1/3)*(x^2+1)^(3/2), (x,-5,5), (y,-5,5), cmap="autumn", fill = False, contours = range(-10,10) + [10 * (k+1) for k in range(5)], labels = True, colorbar = True)

Exercice n. 31, section 6.3.

Il s'agit de trouver les courbes orthogonales à la famille de courbes y=kxy = \frac{k}{x}. Les courbes de la famille cherchée ont pour équation y2=x2+cy^2 = x^2 +c. Voyons quelques courbes de la famille originale, et de la famille des trajectoires orthogonales.

var('x,y') ListeCourbes = [implicit_plot( x*y - c , (x,-3,3), (y,-3,3), color="red") for c in range(-5,5)] FamilleOrthogonale = [implicit_plot(y^2 - x^2 - c, (x,-3,3), (y,-3,3),color="blue") for c in range(-5,5) ] show(sum(ListeCourbes) + sum(FamilleOrthogonale))
(xx, yy)