Example: Maple and Gradient Descent. This is doing the process by hand without a loop.

The function is x^2+4y^2-4x-8y, and we start the process at (0,0).

First, define the function f and plot some level curves. Further, define the partial derivatives f_x, f_y as functions:

with(plots):

f:=(x,y)->x^2+4*y^2-4*x-8*y;

AA:=contourplot(f(x,y),x=-1..4,y=-1..4, contours=[-7.96,-7.85,-7.43,-5.9,0,1]):

f(2,1); #This is the minimum, just to compare later:

fx:=unapply(diff(f(x,y),x),x); fy:=unapply(diff(f(x,y),y),y);

Now create a line in the desired direction, and find the point t that minimizes f along that line:

x0:=0-t*fx(0); y0:=0-t*fy(0);

F1:=f(x0,y0);

diff(F1,t);

solve(%=0,t);

x1:=evalf(subs(t=5/34,x0)); y1:=evalf(subs(t=5/34,y0));

The point (x1,y1) is our new point. Now start the steps over again:

(i) Find the gradient at the new point, (ii) move in the direction of negative gradient, (iii) minimize f along that line:

fx1:=fx(x1); fy1:=fy(y1);

x1t:=x1-t*fx1; y1t:=y1-t*fy1;

F2:=simplify(f(x1t,y1t));

dF2:=diff(F2,t); solve(dF2=0,t);

x2:=subs(t=0.3125,x1t); y2:=subs(t=0.3125,y1t);

We did some extraneous computation- Like showing what f looks like when constrained to the line. This time we'll go straight through the three steps.

fx2:=fx(x2); fy2:=fy(y2);

x2t:=x2-t*fx2; y2t:=y2-t*fy2;

F3:=simplify(f(x2t,y2t)): dF3:=diff(F3,t): solve(dF3=0,t);

x3:=subs(t=0.147,x2t); y3:=subs(t=0.147,y2t);

One more time:

x4t:=x3-t*fx(x3); y4t:=y3-t*fy(y3);

F4:=simplify(f(x4t,y4t)): dF4:=diff(F4,t): solve(dF4=0,t);

x4:=subs(t=0.313,x4t); y4:=subs(t=0.313,y4t);

Plot of the contours together with the path we've taken thus far:

AB:=plot([[0,0],[x1,y1],[x2,y2],[x3,y3],[x4,y4]]):

display(AA,AB,scaling=constrained); #Shows the gradient is orthog to level curves: