Finding relative extreme points and saddle point (Maple)

Question

Searching the site, I have found some problems related to my question, for example this one. But, I want to know how to make codes in Maple to find these points. I could make the following codes as they seems to be (of course it should be completed to be a real program) :

> f:=(x,y)->f(x,y): 
> D1:=diff(f(x,y),x):
> D2:=diff(f(x,y),x):
> S:=solve({D1=0,D2=0},{x,y}):
> ....

But; I feel maybe there is another approach doing the problem in Maple environment. I will be very glad if I know my approach above is enough? Thanks for your time.

Answer 1

This might also give you some ideas, which might be turned into a procedure for re-use on other examples. (...and check it all, and correct my mistakes).

And you might also wish to cover the case where some eigenvalue of the Hessian is zero at a critical point -- and so to differentiate further, say.

You might find that a plot can add some insight. You can overlay the surface, the surface where convex/concave, critical points, solved partial derivatives as spacecurve, etc.

restart:
interface(warnlevel=0):

f:=(x,y)->x^3+3*x*y^2-15*x+y^3-15*y:

critpts:=[solve({D[1](f)(x,y)=0,D[2](f)(x,y)=0},{x,y})]:
critpts:=evalf(map(allvalues,critpts));

     critpts := [{x = 2., y = 1.}, {x = -2., y = -1.}, {x = 0., y = 2.236067977}, 
                 {x = 0., y = -2.236067977}]

H:=VectorCalculus:-Hessian( f(x,y), [x,y] );

                             [6 x     6 y   ]
                        H := [              ]
                             [6 y  6 x + 6 y]

rel_min:=[seq(`if`(LinearAlgebra:-IsDefinite(eval(H,p),
                                         query=positive_definite),
               [(eval([x,y],p)[]),f(eval([x,y],p)[])], NULL), p in critpts)];

                     rel_min := [[2., 1., -30.]]

rel_max:=[seq(`if`(LinearAlgebra:-IsDefinite(eval(H,p),
                                         query=negative_definite),
               [(eval([x,y],p)[]),f(eval([x,y],p)[])], NULL), p in critpts)];

                    rel_max := [[-2., -1., 30.]]

others:=[seq(`if`(LinearAlgebra:-IsDefinite(eval(H,p),
                                         query=indefinite),
               [(eval([x,y],p)[]),f(eval([x,y],p)[])], NULL), p in critpts)];

        others := [[0., 2.236067977, -22.36067978], [0., -2.236067977, 22.36067978]]

posdefconds := (LinearAlgebra:-IsDefinite(H,query=positive_definite)):
posdefconds := simplify(posdefconds) assuming real:
posdefconds := unapply(posdefconds, [x,y]);

                                                        2          2
                posdefconds := (x, y) -> 0 < x and 0 < x  + x y - y 

negdefconds := (LinearAlgebra:-IsDefinite(H,query=negative_definite)):
negdefconds := simplify(negdefconds) assuming real:
negdefconds := unapply(negdefconds, [x,y]);

                                                     2          2    
                negdefconds := (x, y) -> x < 0 and -x  - x y + y  < 0

plots:-display(
    plot3d(f,
         -10..10, -10..10, color=red, transparency=0.5)
  , plot3d(proc(x,y) `if`(negdefconds(x,y),f(x,y),undefined); end proc,
         -10..10, -10..10, color=blue)
  , plot3d(proc(x,y) `if`(posdefconds(x,y),f(x,y),undefined); end proc,
         -10..10, -10..10, color=green)
  , plots:-pointplot3d(rel_min, symbolsize=20, color=red)
  , plots:-pointplot3d(rel_max, symbolsize=20, color=yellow)
  , plots:-pointplot3d(others, symbolsize=20, color=cyan)
  , seq(plots:-spacecurve([rhs(Z[]),y,f(rhs(Z[]),y)],y=-10..10,
                        view=[-10..10,-10..10,-100..100],
                        grid=[40,40,40],color=white),
        Z in solve([D[1](f)(x,y)=0],[x]))
  , seq(plots:-spacecurve([rhs(Z[]),y,f(rhs(Z[]),y)],y=-10..10,
                        view=[-10..10,-10..10,-100..100],
                        grid=[40,40,40],color=black),
        Z in solve([D[2](f)(x,y)=0],[x]))
  , view=-1000..1000, labels=[x,y,f]
  );

Answer 2

f:=(x,y)->x^3+3*x*y^2-15*x+y^3-15*y:
critpts:=[solve({D1(x,y)=0,D2(x,y)=0},{x,y})]:
critpts:=evalf(map(allvalues,critpts)):
dfxx:=(x,y)->[x,y,D1$2](f)(x,y)*D[2$2(x,y)-D1,2(x,y)^2,D1$2(x,y)]:
map(subs,critpts,dfxx(x,y));

[[2., 1., 180., 12.], [-2., -1., 180., -12.], [0., 2.236067977, -179.9999999, 0.], [0., -2.236067977, -179.9999999, 0.]]