Method of Steepest Descent

> Steep := proc(x)
local phi, phiP, t, v;
v := gradf(x);
phi := f( evalm( x-t*v ) );
phiP := -dotprod( gradf( evalm( x -t*v ) ), v );
if phiP<>0 then
t := op(1,[fsolve( phiP=0, t )])
else
t := 0
end if;
convert( evalm(
x - t*gradf(x)
), vector );
end proc:

>

> s[0] := evalm( x0 ):
for k from 0 to N-1 do
s[k+1] := Steep( s[k] );
if norm( s[k+1] - s[k] ) < tol then break end if;
end do:

>

> results[s] := convert( [seq( [ i, evalm(s[i]), f(s[i]) ], i=0..k )] , matrix );

>

> values[s] := plot( [ seq( [i,f(s[i])], i=0..k ) ], color=green,
title="Steepest Descent: Function Values" ):
points[s] := plot( [ seq( convert(s[i],list), i=0..k ) ], color=green,
title="Steepest Descent: Iterates" ):
display( array( [values[s], points[s]] ) );

>