Classifying Critical Points and Lagrange Multipliers 

 

Here we explore quadratic approximations, the "second derivative test" (i.e. classifying min's, max's, and saddle points), and the method of Lagrange multipliers. 

Maple worksheet: minmax.mw 

 

> restart; with(plots): with(VectorCalculus): with(LinearAlgebra):
 

 

Example: Let  f(x, y) = exp(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(4, y), VectorCalculus:-`-`(`*`(`^`(x, 2)))), VectorCalculus:-`-`(`*`(`^`(y, 2))))). We wish to find and classify all of the critical points of f. 

 

> f := (x,y) -> exp(4*y-x^2-y^2):
'f(x,y)' = f(x,y);

contourplot3d(f(x,y),x=-4..4,y=2-sqrt(16-x^2)..2+sqrt(16-x^2),filledregions=true,coloring=[black,red],axes=boxed,contours=15,orientation=[140,70],viewpoint=circleleft);
 

 

f(x, y) = exp(`+`(`*`(4, `*`(y)), `-`(`*`(`^`(x, 2))), `-`(`*`(`^`(y, 2)))))
Plot_2d
 

 

First, we should find the critical points. This is done by solving the system of equations:  f[x](x, y) = 0   and   f[y](x, y) = 0 

 

> solve({diff(f(x,y),x)=0,diff(f(x,y),y)=0});
 

{x = 0, y = 2} (1)
 

 

There is one critical point: (0, 2). 

 

To classify this point we can either compute eigenvalues or use the 2nd derivative test (compute the determinant of the Hessian etc). 

 

> H := subs({x=0,y=2},Hessian(f(x,y),[x,y]));
eigenValues = Eigenvalues(H);
det = Determinant(H);
 

 

 

Matrix(%id = 18446744078099562734)
eigenValues = Vector[column](%id = 18446744078099563934)
det = `+`(`*`(4, `*`(`^`(exp(4), 2)))) (2)
 

 

(0,2) is a relative maximum [negative eigenvalues OR positive determinant and negative diagonal entries] 

 

 

Example: Let  f(x, y) = VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(3, x), VectorCalculus:-`-`(`*`(`^`(x, 3)))), VectorCalculus:-`-`(VectorCalculus:-`*`(2, `*`(`^`(y, 2))))), `*`.... Let's find and classify all of the critical points of f. 

 

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

contourplot3d(f(x,y),x=-1.5..1.5,y=-1.5..1.5,filledregions=true,coloring=[blue,green],axes=boxed,contours=15,orientation=[140,70],viewpoint=circleleft);
 

 

f(x, y) = `+`(`*`(3, `*`(x)), `-`(`*`(`^`(x, 3))), `-`(`*`(2, `*`(`^`(y, 2)))), `*`(`^`(y, 4)))
Plot_2d
 

 

First, we should find the critical points. This is done by solving the system of equations:  f[x](x, y) = 0   and   f[y](x, y) = 0 

 

> solve({diff(f(x,y),x)=0,diff(f(x,y),y)=0});
 

{x = 1, y = 0}, {x = -1, y = 0}, {x = 1, y = 1}, {x = 1, y = -1}, {x = -1, y = 1}, {x = -1, y = -1} (3)
 

 

There are 6 critical points: (1,0), (-1,0), (1,1), (1,-1), (-1,1), (-1,-1). To classify these points we can either compute eigenvalues or use the 2nd derivative test (compute the determinant of the Hessian etc). 

 

> # the critical point (x,y)=(1,0)...
H := evalf(subs({x=1,y=0},Hessian(f(x,y),[x,y])));
eigenValues = evalf(Eigenvalues(H));
det = evalf(Determinant(H));
 

 

 

Matrix(%id = 18446744078122148198)
eigenValues = Vector[column](%id = 18446744078122150598)
det = 24. (4)
 

 

(1,0) is a relative maximum [negative eigenvalues OR positive determinant and negative diagonal entries] 

 

> # the critical point (x,y)=(-1,0)...
H := evalf(subs({x=-1,y=0},Hessian(f(x,y),[x,y])));
eigenValues = evalf(Eigenvalues(H));
det = evalf(Determinant(H));
 

 

 

Matrix(%id = 18446744078122150958)
eigenValues = Vector[column](%id = 18446744078122145182)
det = -24. (5)
 

 

(-1,0) is a saddle point [eigenvalues have mixed signs OR determinant is negative] 

 

> # the critical point (x,y)=(1,1)...
H := evalf(subs({x=1,y=1},Hessian(f(x,y),[x,y])));
eigenValues = evalf(Eigenvalues(H));
det = evalf(Determinant(H));
 

 

 

Matrix(%id = 18446744078122145542)
eigenValues = Vector[column](%id = 18446744078122139766)
det = -48. (6)
 

 

(1,1) is a saddle point [eigenvalues have mixed signs OR determinant is negative] 

 

> # the critical point (x,y)=(-1,1)...
H := evalf(subs({x=-1,y=1},Hessian(f(x,y),[x,y])));
eigenValues = evalf(Eigenvalues(H));
det = evalf(Determinant(H));
 

 

 

Matrix(%id = 18446744078122140126)
eigenValues = Vector[column](%id = 18446744078122142526)
det = 48. (7)
 

 

(-1,1) is a relative minimum [positive eigenvalues OR positive determinant and positive diagonal entries] 

 

> # the critical point (x,y)=(1,-1)...
H := evalf(subs({x=1,y=-1},Hessian(f(x,y),[x,y])));
eigenValues = evalf(Eigenvalues(H));
det = evalf(Determinant(H));
 

 

 

Matrix(%id = 18446744078122142886)
eigenValues = Vector[column](%id = 18446744078122137110)
det = -48. (8)
 

 

(1,-1) is a saddle point [eigenvalues have mixed signs OR determinant is negative] 

 

> # the critical point (x,y)=(-1,-1)...
H := evalf(subs({x=-1,y=-1},Hessian(f(x,y),[x,y])));
eigenValues = evalf(Eigenvalues(H));
det = evalf(Determinant(H));
 

 

 

Matrix(%id = 18446744078122137470)
eigenValues = Vector[column](%id = 18446744078122131694)
det = 48. (9)
 

 

(-1,-1) is a relative minimum [positive eigenvalues OR positive determinant and positive diagonal entries] 

 

 

Example: Consider the following function. 

 

> h := (x,y) -> (3/2)*(x^2+y^2)-(x^4+y^4)/2:
'h(x,y)' = h(x,y);

plot3d(h(x,y),x=-2..2,y=-2..2);
 

 

h(x, y) = `+`(`*`(`/`(3, 2), `*`(`^`(x, 2))), `*`(`/`(3, 2), `*`(`^`(y, 2))), `-`(`*`(`/`(1, 2), `*`(`^`(x, 4)))), `-`(`*`(`/`(1, 2), `*`(`^`(y, 4)))))
Plot_2d
 

 

The following command finds where f[x] = 0 and f[y] = 0 (critical points). Maple finds 4 critical points -- Maple lost several solutions. There are visibally 4 peaks, 1 valley, and 4 saddles betwen peaks. So we should have a total of 9 critical points. 

 

> evalf(solve({diff(h(x,y),x), diff(h(x,y),y)}));
 

{x = 0., y = 0.}, {x = 1.224744872, y = 0.}, {x = 0., y = 1.224744872}, {x = 1.224744872, y = 1.224744872} (10)
 

 

Note: 1.224744872 is actually sqrt(VectorCalculus:-`*`(3, `/`(1, 2))) 

 

The following command computes the Hessian matrix of h(x, y) and then evaluates the Hessian at a critical point and finds its eigenvalues. 

 

Note: If all of the eigenvalues are positive, the critical point is a relative minimum. 

         If all of the eigenvalues are negative, the critical point is a relative maximum. 

         If some eigenvalues are positive and others negative, the critical point is a saddle point. 

         If any of the eigenvalues are 0, this "second derivative test" does not apply. 

 

> H := Hessian(h(x,y),[x,y]);

Eigenvalues(subs({x=sqrt(3/2),y=sqrt(3/2)},H));
 

 

Matrix(%id = 18446744078122132654)
Vector[column](%id = 18446744078122133854) (11)
 

 

The critical point (x, y) = (1.224744872, 1.224744872) is a local maximum since both eigenvalues are negative (they're both -6). 

 

Now let's compute the eigenvalues for all of the other critical points and determine whether they are relative mins, maxes, saddle points, or if the test does not apply. 

 

> Eigenvalues(subs({x=-sqrt(3/2),y=sqrt(3/2)},H));
Eigenvalues(subs({x=sqrt(3/2),y=-sqrt(3/2)},H));
Eigenvalues(subs({x=-sqrt(3/2),y=-sqrt(3/2)},H));
Eigenvalues(subs({x=0,y=0},H));
Eigenvalues(subs({x=0,y=sqrt(3/2)},H));
Eigenvalues(subs({x=0,y=-sqrt(3/2)},H));
Eigenvalues(subs({x=sqrt(3/2),y=0},H));
Eigenvalues(subs({x=-sqrt(3/2),y=0},H));
 

 

 

 

 

 

 

 

Vector[column](%id = 18446744078122135054)
Vector[column](%id = 18446744078122128078)
Vector[column](%id = 18446744078122129278)
Vector[column](%id = 18446744078122130478)
Vector[column](%id = 18446744078122123502)
Vector[column](%id = 18446744078122124702)
Vector[column](%id = 18446744078122125902)
Vector[column](%id = 18446744078122127102) (12)
 

 

Both the eigenvalues are negative when (x,y)=(0,0), so this is a relative minimum. All of the points which have 2 positive eigenvalues are relative maximums. The remaining points which have eigenvalues of mixed signs are saddle points. 

 

local minima:  (0,0) 

 

local maxima:         

 

saddle points:         

 

The test does not fail for any of the critical points. 

 

 

The following code plots the surface along with its quadratic approximation based at (1.224744872, 1.224744872). 

 

> a := 1.224744872:
b := 1.224744872:

Q := h(a,b) +
    subs({x=a,y=b},diff(h(x,y),x))*(x-a) +
    subs({x=a,y=b},diff(h(x,y),y))*(y-b) +
    (1/2)*subs({x=a,y=b},diff(h(x,y),x,x))*(x-a)^2 +
    (1/2)*subs({x=a,y=b},diff(h(x,y),x,y))*(x-a)*(y-b) +
    (1/2)*subs({x=a,y=b},diff(h(x,y),y,x))*(y-b)*(x-a) +
    (1/2)*subs({x=a,y=b},diff(h(x,y),y,y))*(y-b)^2;

A := plot3d(h(x,y),x=-2..2,y=-2..2):
B := plot3d(Q,x=a-1.25..a+1.25,y=-sqrt((1.25)^2-(x-a)^2)+b..sqrt((1.25)^2-(x-a)^2)+b,color=red,thickness=3):
display({A,B},viewpoint=circleleft);
 

 

`+`(2.250000009, `-`(`*`(0.4e-8, `*`(x))), `-`(`*`(0.4e-8, `*`(y))), `-`(`*`(3.000000003, `*`(`^`(`+`(x, `-`(1.224744872)), 2)))), `-`(`*`(3.000000003, `*`(`^`(`+`(y, `-`(1.224744872)), 2)))))
Plot_2d
 

 

Now let's see what happens at a saddle point and the relative min at the origin. 

 

> a := 0:
b := 1.224744872:

Q := h(a,b) +
    subs({x=a,y=b},diff(h(x,y),x))*(x-a) +
    subs({x=a,y=b},diff(h(x,y),y))*(y-b) +
    (1/2)*subs({x=a,y=b},diff(h(x,y),x,x))*(x-a)^2 +
    (1/2)*subs({x=a,y=b},diff(h(x,y),x,y))*(x-a)*(y-b) +
    (1/2)*subs({x=a,y=b},diff(h(x,y),y,x))*(y-b)*(x-a) +
    (1/2)*subs({x=a,y=b},diff(h(x,y),y,y))*(y-b)^2;

A := plot3d(h(x,y),x=-2..2,y=-2..2):
B := plot3d(Q,x=a-1.25..a+1.25,y=-sqrt((1.25)^2-(x-a)^2)+b..sqrt((1.25)^2-(x-a)^2)+b,color=red,thickness=3):
display({A,B},viewpoint=circleleft);

a := 0:
b := 0:

Q := h(a,b) +
    subs({x=a,y=b},diff(h(x,y),x))*(x-a) +
    subs({x=a,y=b},diff(h(x,y),y))*(y-b) +
    (1/2)*subs({x=a,y=b},diff(h(x,y),x,x))*(x-a)^2 +
    (1/2)*subs({x=a,y=b},diff(h(x,y),x,y))*(x-a)*(y-b) +
    (1/2)*subs({x=a,y=b},diff(h(x,y),y,x))*(y-b)*(x-a) +
    (1/2)*subs({x=a,y=b},diff(h(x,y),y,y))*(y-b)^2;

A := plot3d(h(x,y),x=-2..2,y=-2..2):
B := plot3d(Q,x=a-1.25..a+1.25,y=-sqrt((1.25)^2-(x-a)^2)+b..sqrt((1.25)^2-(x-a)^2)+b,color=red,thickness=3):
display({A,B},viewpoint=circleleft);
 

 

 

 

`+`(1.125000005, `-`(`*`(0.4e-8, `*`(y))), `*`(`/`(3, 2), `*`(`^`(x, 2))), `-`(`*`(3.000000003, `*`(`^`(`+`(y, `-`(1.224744872)), 2)))))
Plot_2d
`+`(`*`(`/`(3, 2), `*`(`^`(x, 2))), `*`(`/`(3, 2), `*`(`^`(y, 2))))
Plot_2d
 

 

Example: Consider the following function... 

 

> h := (x,y) -> 3*x^2*y+y^3-3*x^2-3*y^2+2+exp(-x^2-y^2):
'h(x,y)' = h(x,y);

contourplot3d(h(x,y),x=-2..2,y=1-sqrt(4-x^2)..1+sqrt(4-x^2),filledregions=true,contours=15,axes=boxed,viewpoint=circleleft);
 

 

h(x, y) = `+`(`*`(3, `*`(`^`(x, 2), `*`(y))), `*`(`^`(y, 3)), `-`(`*`(3, `*`(`^`(x, 2)))), `-`(`*`(3, `*`(`^`(y, 2)))), 2, exp(`+`(`-`(`*`(`^`(x, 2))), `-`(`*`(`^`(y, 2))))))
Plot_2d
 

 

First, we will plot the plane tangent to z = h(x, y) at the point  

 

> a :=  0.5:
b := -0.5:

L := h(a,b) +
    subs({x=a,y=b},diff(h(x,y),x))*(x-a) +
    subs({x=a,y=b},diff(h(x,y),y))*(y-b);

A := plot3d(h(x,y),x=-2..2,y=1-sqrt(4-x^2)..1+sqrt(4-x^2)):
B := plot3d(L,x=a-1.5..a+1.5,y=-sqrt((1.5)^2-(x-a)^2)+b..sqrt((1.5)^2-(x-a)^2)+b,color=red,thickness=3):
display({A,B},viewpoint=circleleft);
 

 

`+`(.6065306597, `*`(`+`(`-`(4.50), `-`(`*`(1.0, `*`(exp(-.50))))), `*`(`+`(x, `-`(.5)))), `*`(`+`(4.50, `*`(1.0, `*`(exp(-.50)))), `*`(`+`(y, .5))))
Plot_2d
 

 

h(x, y) has 4 critical points. Let's determine what type of critical point each one is.
 

The graph helps us estimate where the critical points are. It look like they are near: (1,1), (-1,1), (0,0), and (0,2). Notice the use of "fsolve" with initial guess of (x, y) = (1, 1). Maple tracks down the (nearest) solution (x, y) = (1.0385 .. (), 1.0385 .. ()). We'll to use our other estimates of point locations to find the other 3 points. In addition, let's classify each point using 2 different methods: (1) We'll compute eigenvalues for the Hessian matrix. (2) We'll compute the Hessian and its determinant and check out h[xx] if necessary. 

 

> fsolve({diff(h(x,y),x)=0,diff(h(x,y),y)=0},{x=1,y=1});
"Eigenvalues" = evalf(Eigenvalues(subs({x= 1.038550488,y=1.038550488},Hessian(h(x,y),[x,y]))));

"Hessian" = evalf(subs({x= 1.038550488,y=1.038550488},Hessian(h(x,y),[x,y])));
"D" = evalf(Determinant(subs({x= 1.038550488,y=1.038550488},Hessian(h(x,y),[x,y]))));
 

 

 

 

{x = 1.038550488, y = 1.038550488}
Eigenvalues
Hessian
D (13)
 

 

Point #1 is (. At this point the Hessian's eigenvalues have mixed signs. So it must be a saddle point. Alternativley, the determinant of the Hessian matrix is negative at this point, thus this is a saddle point. 

 

It looks like there is another saddle point near -1, 1 as well. Let's use this as our second guess. 

 

> fsolve({diff(h(x,y),x)=0,diff(h(x,y),y)=0},{x=-1,y=1});
"Eigenvalues" = evalf(Eigenvalues(subs({x=-1.038550488,y=1.038550488},Hessian(h(x,y),[x,y]))));

"Hessian" = evalf(subs({x=-1.038550488,y=1.038550488},Hessian(h(x,y),[x,y])));
"D" = evalf(Determinant(subs({x=-1.038550488,y=1.038550488},Hessian(h(x,y),[x,y]))));
 

 

 

 

{x = -1.038550488, y = 1.038550488}
Eigenvalues
Hessian
D (14)
 

 

Point #2 is (-. At this point the Hessian's eigenvalues have mixed signs. So it is a saddle point as well. (Also, again, the determinant of the Hessian matrix is negative at this point.) 

 

It looks like there is a relative max near 0, 0. Let's use this as our third guess. 

 

> fsolve({diff(h(x,y),x)=0,diff(h(x,y),y)=0},{x=0,y=0});
"Eigenvalues" = evalf(Eigenvalues(subs({x=0,y=0},Hessian(h(x,y),[x,y]))));

"Hessian" = evalf(subs({x=0,y=0},Hessian(h(x,y),[x,y])));
"D" = evalf(Determinant(subs({x=0,y=0},Hessian(h(x,y),[x,y]))));
 

 

 

 

{x = 0., y = 0.}
Eigenvalues
Hessian
D (15)
 

 

Point #3 is (0,0). At the origin the Hessian's eigenvalues are both negative. Thus this is a relative maximum. Alternatively, we could note that (at this point) the determinant of the Hessian is postive while the diagonal entries are negiative, so (again) it's a relative maximum. 

 

For our final guess. It looks like there is a relative minimum near (0,2). 

 

> fsolve({diff(h(x,y),x)=0,diff(h(x,y),y)=0},{x=0,y=2});
"Eigenvalues" = evalf(Eigenvalues(subs({x=0,y=2.011652763},Hessian(h(x,y),[x,y]))));

"Hessian" = evalf(subs({x=0,y=2.011652763},Hessian(h(x,y),[x,y])));
"D" = evalf(Determinant(subs({x=0,y=2.011652763},Hessian(h(x,y),[x,y]))));
 

 

 

 

{x = 0., y = 2.011652763}
Eigenvalues
Hessian
D (16)
 

 

Point #4 is (0, 2.011652763). At this point the Hessian's eigenvalues are both positive. Thus this is a relative minimum. Alternatively, (at this point) we have a positive determinant and positive values on the diagonal of the Hessian, so (again) it's a relative minimum. 

 

 

Summary: (- is a saddle point. ( is a saddle point. (0,0) is a relative maximum. (0, 2.011652763) is a relative minimum. 

 

 

Let's plot z = h(x, y) together with its quadratic approximation at each of the 4 critical points.  

 

> # Point #1: Saddle Point
a := -1.038550488:
b :=  1.038550488:

Q := h(a,b) +
    subs({x=a,y=b},diff(h(x,y),x))*(x-a) +
    subs({x=a,y=b},diff(h(x,y),y))*(y-b) +
    (1/2)*subs({x=a,y=b},diff(h(x,y),x,x))*(x-a)^2 +
    (1/2)*subs({x=a,y=b},diff(h(x,y),x,y))*(x-a)*(y-b) +
    (1/2)*subs({x=a,y=b},diff(h(x,y),y,x))*(x-a)*(y-b) +
    (1/2)*subs({x=a,y=b},diff(h(x,y),y,y))*(y-b)^2;


A := plot3d(h(x,y),x=-2..2,y=1-sqrt(4-x^2)..1+sqrt(4-x^2)):
B := plot3d(Q,x=a-1.5..a+1.5,y=-sqrt((1.5)^2-(x-a)^2)+b..sqrt((1.5)^2-(x-a)^2)+b,color=red,thickness=3):
display({A,B},viewpoint=circleleft);

# Point #2: Saddle Point
a :=  1.038550488:
b :=  1.038550488:

Q := h(a,b) +
    subs({x=a,y=b},diff(h(x,y),x))*(x-a) +
    subs({x=a,y=b},diff(h(x,y),y))*(y-b) +
    (1/2)*subs({x=a,y=b},diff(h(x,y),x,x))*(x-a)^2 +
    (1/2)*subs({x=a,y=b},diff(h(x,y),x,y))*(x-a)*(y-b) +
    (1/2)*subs({x=a,y=b},diff(h(x,y),y,x))*(x-a)*(y-b) +
    (1/2)*subs({x=a,y=b},diff(h(x,y),y,y))*(y-b)^2;


A := plot3d(h(x,y),x=-2..2,y=1-sqrt(4-x^2)..1+sqrt(4-x^2)):
B := plot3d(Q,x=a-1.5..a+1.5,y=-sqrt((1.5)^2-(x-a)^2)+b..sqrt((1.5)^2-(x-a)^2)+b,color=red,thickness=3):
display({A,B},viewpoint=circleleft);


# Point #3: Relative Maximum
a :=  0:
b :=  0:

Q := h(a,b) +
    subs({x=a,y=b},diff(h(x,y),x))*(x-a) +
    subs({x=a,y=b},diff(h(x,y),y))*(y-b) +
    (1/2)*subs({x=a,y=b},diff(h(x,y),x,x))*(x-a)^2 +
    (1/2)*subs({x=a,y=b},diff(h(x,y),x,y))*(x-a)*(y-b) +
    (1/2)*subs({x=a,y=b},diff(h(x,y),y,x))*(x-a)*(y-b) +
    (1/2)*subs({x=a,y=b},diff(h(x,y),y,y))*(y-b)^2;


A := plot3d(h(x,y),x=-2..2,y=1-sqrt(4-x^2)..1+sqrt(4-x^2)):
B := plot3d(Q,x=a-1.5..a+1.5,y=-sqrt((1.5)^2-(x-a)^2)+b..sqrt((1.5)^2-(x-a)^2)+b,color=red,thickness=3):
display({A,B},viewpoint=circleleft);


# Point #4: Relative Minimum
a :=  0:
b :=  2.011652763:

Q := h(a,b) +
    subs({x=a,y=b},diff(h(x,y),x))*(x-a) +
    subs({x=a,y=b},diff(h(x,y),y))*(y-b) +
    (1/2)*subs({x=a,y=b},diff(h(x,y),x,x))*(x-a)^2 +
    (1/2)*subs({x=a,y=b},diff(h(x,y),x,y))*(x-a)*(y-b) +
    (1/2)*subs({x=a,y=b},diff(h(x,y),y,x))*(x-a)*(y-b) +
    (1/2)*subs({x=a,y=b},diff(h(x,y),y,y))*(y-b)^2;


A := plot3d(h(x,y),x=-2..2,y=1-sqrt(4-x^2)..1+sqrt(4-x^2)):
B := plot3d(Q,x=a-1.5..a+1.5,y=-sqrt((1.5)^2-(x-a)^2)+b..sqrt((1.5)^2-(x-a)^2)+b,color=red,thickness=3):
display({A,B},viewpoint=circleleft);
 

 

 

 

 

 

 

 

`+`(.1247974709, `*`(`+`(`-`(.240219768), `*`(2.077100976, `*`(exp(-2.157174232)))), `*`(`+`(x, 1.038550488))), `*`(`+`(.240219768, `-`(`*`(2.077100976, `*`(exp(-2.157174232))))), `*`(`+`(y, `-`(1.038...
`+`(.1247974709, `*`(`+`(`-`(.240219768), `*`(2.077100976, `*`(exp(-2.157174232)))), `*`(`+`(x, 1.038550488))), `*`(`+`(.240219768, `-`(`*`(2.077100976, `*`(exp(-2.157174232))))), `*`(`+`(y, `-`(1.038...
`+`(.1247974709, `*`(`+`(`-`(.240219768), `*`(2.077100976, `*`(exp(-2.157174232)))), `*`(`+`(x, 1.038550488))), `*`(`+`(.240219768, `-`(`*`(2.077100976, `*`(exp(-2.157174232))))), `*`(`+`(y, `-`(1.038...
`+`(.1247974709, `*`(`+`(`-`(.240219768), `*`(2.077100976, `*`(exp(-2.157174232)))), `*`(`+`(x, 1.038550488))), `*`(`+`(.240219768, `-`(`*`(2.077100976, `*`(exp(-2.157174232))))), `*`(`+`(y, `-`(1.038...
Plot_2d
`+`(.1247974709, `*`(`+`(.240219768, `-`(`*`(2.077100976, `*`(exp(-2.157174232))))), `*`(`+`(x, `-`(1.038550488)))), `*`(`+`(.240219768, `-`(`*`(2.077100976, `*`(exp(-2.157174232))))), `*`(`+`(y, `-`(...
`+`(.1247974709, `*`(`+`(.240219768, `-`(`*`(2.077100976, `*`(exp(-2.157174232))))), `*`(`+`(x, `-`(1.038550488)))), `*`(`+`(.240219768, `-`(`*`(2.077100976, `*`(exp(-2.157174232))))), `*`(`+`(y, `-`(...
`+`(.1247974709, `*`(`+`(.240219768, `-`(`*`(2.077100976, `*`(exp(-2.157174232))))), `*`(`+`(x, `-`(1.038550488)))), `*`(`+`(.240219768, `-`(`*`(2.077100976, `*`(exp(-2.157174232))))), `*`(`+`(y, `-`(...
`+`(.1247974709, `*`(`+`(.240219768, `-`(`*`(2.077100976, `*`(exp(-2.157174232))))), `*`(`+`(x, `-`(1.038550488)))), `*`(`+`(.240219768, `-`(`*`(2.077100976, `*`(exp(-2.157174232))))), `*`(`+`(y, `-`(...
Plot_2d
`+`(3, `*`(`+`(`-`(3), `-`(exp(0))), `*`(`^`(x, 2))), `*`(`+`(`-`(3), `-`(exp(0))), `*`(`^`(y, 2))))
Plot_2d
`+`(`-`(1.982111915), `*`(`+`(0.7032394e-1, `-`(`*`(4.023305526, `*`(exp(-4.046746839))))), `*`(`+`(y, `-`(2.011652763)))), `*`(`+`(3.034958290, `-`(exp(-4.046746839))), `*`(`^`(x, 2))), `*`(`+`(3.034...
`+`(`-`(1.982111915), `*`(`+`(0.7032394e-1, `-`(`*`(4.023305526, `*`(exp(-4.046746839))))), `*`(`+`(y, `-`(2.011652763)))), `*`(`+`(3.034958290, `-`(exp(-4.046746839))), `*`(`^`(x, 2))), `*`(`+`(3.034...
Plot_2d
 

 

Example: Let's find the points on the level surface VectorCalculus:-`*`(VectorCalculus:-`*`(`*`(`^`(x, 2)), `*`(`^`(y, 2))), z) = 1 which are closest to the origin. It's helpful to begin with a graph of this surface. 

 

> F := (x,y,z) -> x^2*y^2*z: 'F(x,y,z)' = F(x,y,z); F(x,y,z)=1; implicitplot3d(F(x,y,z)=1,x=-2..2,y=-2..2,z=-1..2,scaling=constrained,axes=normal,numpoints=5000,orientation=[-81,77],viewpoint=circleleft);
 

 

 

F(x, y, z) = `*`(`^`(x, 2), `*`(`^`(y, 2), `*`(z)))
`*`(`^`(x, 2), `*`(`^`(y, 2), `*`(z))) = 1
Plot_2d
 

 

From the graph we can see that the problem should have FOUR solutions. We need to find a function describing the distance squared (to make things look nicer) from the surface to the origin. Then we mustcompute some partial derivatives, find critical points, and test them. 

 

Note: Minimizing the distanced squared is the same as minimizing the distance. 

 

> distSq := (x,y,z) -> (x-0)^2 + (y-0)^2 + (z-0)^2: 'distSq(x,y,z)' = distSq(x,y,z); f := (x,y) -> distSq(x,y,solve(F(x,y,z)=1,z)): 'f(x,y)' = f(x,y);
 

 

distSq(x, y, z) = `+`(`*`(`^`(x, 2)), `*`(`^`(y, 2)), `*`(`^`(z, 2)))
f(x, y) = `+`(`*`(`^`(x, 2)), `*`(`^`(y, 2)), `/`(1, `*`(`^`(x, 4), `*`(`^`(y, 4))))) (17)
 

 

So f(x, y) is the distance from (x,y,z) to (0,0,0) where F(x,y,z)=1 (we solved F(x,y,z)=1 for z and plugged that into the distance function and thus eliminated the third variable). Our next task is to find the critical points of f. Four of these points should be our solutions. Recall that critical points occur when either a partial fails to exist or if both partials are equal to zero. Let's compute the (first and second) partials of f  and analyze. 

 

> fx := diff(f(x,y),x); fy := diff(f(x,y),y); fxx := diff(f(x,y),x,x); fxy := diff(f(x,y),x,y); fyx := diff(f(x,y),y,x);
fyy := diff(f(x,y),y,y);
 

 

 

 

 

 

`+`(`*`(2, `*`(x)), `-`(`/`(`*`(4), `*`(`^`(x, 5), `*`(`^`(y, 4))))))
`+`(`*`(2, `*`(y)), `-`(`/`(`*`(4), `*`(`^`(x, 4), `*`(`^`(y, 5))))))
`+`(2, `/`(`*`(20), `*`(`^`(x, 6), `*`(`^`(y, 4)))))
`+`(`/`(`*`(16), `*`(`^`(x, 5), `*`(`^`(y, 5)))))
`+`(`/`(`*`(16), `*`(`^`(x, 5), `*`(`^`(y, 5)))))
`+`(2, `/`(`*`(20), `*`(`^`(x, 4), `*`(`^`(y, 6))))) (18)
 

 

To find the critical points we need to know where the first partials are equal to zero. 

 

Note: First, notice that the mixed partials f[xy] and f[yx] are equal (not surprising given Clairaut's theorem). Second, notice that these partials are defined (and continuous) everywhere except for the origin (which isn't relevant since `<>`(VectorCalculus:-`*`(VectorCalculus:-`*`(`^`(0, 2), `^`(0, 2)), z), 1) -- such points aren't on the surface F(x,y,z)=1). So we aren't going to get any (relevant) critical points from the non-existence of partials. 

 

> critPts := solve({fx=0,fy=0},[x,y]);
 

[[x = `+`(`-`(RootOf(`+`(`*`(`^`(_Z, 10)), `-`(2))))), y = RootOf(`+`(`*`(`^`(_Z, 10)), `-`(2)))], [x = RootOf(`+`(`*`(`^`(_Z, 10)), `-`(2))), y = RootOf(`+`(`*`(`^`(_Z, 10)), `-`(2)))]] (19)
 

 

Note: "RootOf(VectorCalculus:-`+`(`*`(`^`(Z, 10)), -2))" stands for a number Z which is a solution of VectorCalculus:-`+`(`*`(`^`(Z, 10)), -2) = 0 this means Z = `*`(`^`(2, `/`(1, 10))) or Z = VectorCalculus:-`-`(`*`(`^`(2, `/`(1, 10)))). So we actually have FOUR solutions! 

 

> plot(Z^10-2,Z=-1.1..1.1);
 

Plot_2d
 

 

The above graph shows that VectorCalculus:-`+`(`*`(`^`(Z, 10)), -2) has TWO x-intercepts.  

 

At this point the problem is essentially solved. Since we have four critical points and are expecting four solutions, these must be our solutions. We could plug them into the function f and see that they all give the same (minimal) value and thus are minimums. But let's apply our critical point test anyway. Next, we need to construct the "D" operator which appears in our second derivative test. (We can't use the name "D" itself -- Maple's all ready using that for something.) 

 

> testD := fxx*fyy - (fxy)^2;
 

`+`(`*`(`+`(2, `/`(`*`(20), `*`(`^`(x, 6), `*`(`^`(y, 4))))), `*`(`+`(2, `/`(`*`(20), `*`(`^`(x, 4), `*`(`^`(y, 6))))))), `-`(`/`(`*`(256), `*`(`^`(x, 10), `*`(`^`(y, 10)))))) (20)
 

 

Let's plug-in our critical points and see what we have. 

 

> 0 < simplify(subs(critPts[1],testD)); 0 < simplify(subs(critPts[1],fxx)); 0 < simplify(subs(critPts[2],testD));
0 < simplify(subs(critPts[2],fxx));
 

 

 

 

`<`(0, 80)
`<`(0, 12)
`<`(0, 80)
`<`(0, 12) (21)
 

 

Therefore, at the critical points, we have D=80 > 0 and `and`(f[xx] = 12, `>`(12, 0)). So these are local minima. 

 

We get that the surface is closest to the origin when... 

 

> op(critPts[1]),z=subs(critPts[1],solve(F(x,y,z)=2,z)); op(critPts[2]),z=subs(critPts[2],solve(F(x,y,z)=2,z));
 

 

x = `+`(`-`(RootOf(`+`(`*`(`^`(_Z, 10)), `-`(2))))), y = RootOf(`+`(`*`(`^`(_Z, 10)), `-`(2))), z = `+`(`/`(`*`(2), `*`(`^`(RootOf(`+`(`*`(`^`(_Z, 10)), `-`(2))), 4))))
x = RootOf(`+`(`*`(`^`(_Z, 10)), `-`(2))), y = RootOf(`+`(`*`(`^`(_Z, 10)), `-`(2))), z = `+`(`/`(`*`(2), `*`(`^`(RootOf(`+`(`*`(`^`(_Z, 10)), `-`(2))), 4)))) (22)
 

 

Note... 

 

> z=2/(2^(1/10))^4;
 

z = `*`(`^`(2, `/`(3, 5))) (23)
 

 

So the four solutions are (x, y, z) = (`&+-`(`*`(`^`(2, `/`(1, 10)))), `*`(`^`(2, `/`(1, 10))), `^`(2, VectorCalculus:-`*`(3, `/`(1, 5))))and `&+-`(`*`(`^`(2, `/`(1, 10)))), VectorCalculus:-`-`(`*`(`^`(2, `/`(1, 10)))), `^`(2, VectorCalculus:-`*`(3, `/`(1, 5))). 

 

Next, let's solve the problem again. But this time we'll use Lagrange multipliers. First: our objective function f(x, y, z)is the distance from the origin (squared) and our constraint g(x, y, z) = K is the level surface defined by VectorCalculus:-`*`(VectorCalculus:-`*`(`*`(`^`(x, 2)), `*`(`^`(y, 2))), z) = 1. 

 

> f := (x,y,z) -> x^2+y^2+z^2: 'f(x,y,z)'=f(x,y,z); g := (x,y,z) -> x^2*y^2*z: g(x,y,z)=1;
 

 

f(x, y, z) = `+`(`*`(`^`(x, 2)), `*`(`^`(y, 2)), `*`(`^`(z, 2)))
`*`(`^`(x, 2), `*`(`^`(y, 2), `*`(z))) = 1 (24)
 

 

Next, we need to compute the gradient of f  and the gradient of g to find our system of equations. 

 

> gradf := Gradient(f(x,y,z),[x,y,z]); gradg := Gradient(g(x,y,z),[x,y,z]); gradf = lambda * gradg;
 

 

 

Vector[column](%id = 18446744078120402566)
Vector[column](%id = 18446744078120402686)
Vector[column](%id = 18446744078120402566) = Vector[column](%id = 18446744078120402806) (25)
 

 

We can rip off the components of each side of the vector equation by using something like this (dotting with the i, j, and k vectors)... 

 

> Eqn1 := gradf.<1,0,0> = lambda * gradg.<1,0,0>; Eqn2 := gradf.<0,1,0> = lambda * gradg.<0,1,0>; Eqn3 := gradf.<0,0,1> = lambda * gradg.<0,0,1>;
 

 

 

`+`(`*`(2, `*`(x))) = `+`(`*`(2, `*`(lambda, `*`(x, `*`(`^`(y, 2), `*`(z))))))
`+`(`*`(2, `*`(y))) = `+`(`*`(2, `*`(lambda, `*`(`^`(x, 2), `*`(y, `*`(z))))))
`+`(`*`(2, `*`(z))) = `*`(lambda, `*`(`^`(x, 2), `*`(`^`(y, 2)))) (26)
 

 

Now we solve (don't forget the constraint: g(x,y,z)=1) and then plug out solution(s) into the objective function to see which critical points give us the minimum. We see that all critical points are minima. Evaluating our objective function at these points tells us that the minimum distance to the origin is about 2.87. 

 

> soln := solve({Eqn1,Eqn2,Eqn3,g(x,y,z)=1}, [x,y,z,lambda]); evalf(subs(soln[1],f(x,y,z))); evalf(subs(soln[2],f(x,y,z)));
 

 

 

[[x = `+`(`-`(RootOf(`+`(`*`(`^`(_Z, 2)), `-`(RootOf(`+`(`*`(`^`(_Z, 5)), `-`(2)))))))), y = RootOf(`+`(`*`(`^`(_Z, 2)), `-`(RootOf(`+`(`*`(`^`(_Z, 5)), `-`(2)))))), z = `+`(`*`(`/`(1, 2), `*`(`^`(Roo...
[[x = `+`(`-`(RootOf(`+`(`*`(`^`(_Z, 2)), `-`(RootOf(`+`(`*`(`^`(_Z, 5)), `-`(2)))))))), y = RootOf(`+`(`*`(`^`(_Z, 2)), `-`(RootOf(`+`(`*`(`^`(_Z, 5)), `-`(2)))))), z = `+`(`*`(`/`(1, 2), `*`(`^`(Roo...
[[x = `+`(`-`(RootOf(`+`(`*`(`^`(_Z, 2)), `-`(RootOf(`+`(`*`(`^`(_Z, 5)), `-`(2)))))))), y = RootOf(`+`(`*`(`^`(_Z, 2)), `-`(RootOf(`+`(`*`(`^`(_Z, 5)), `-`(2)))))), z = `+`(`*`(`/`(1, 2), `*`(`^`(Roo...
2.871745890
2.871745890 (27)
 

 

Although they look different from the our previous answers (finding critical points and using the second derivative test), they are in fact the same. For example: The root of VectorCalculus:-`+`(`*`(`^`(Z, 5)), -2) is Z = `*`(`^`(2, `/`(1, 5))). Then the roots of VectorCalculus:-`+`(`*`(`^`(Z, 2)), VectorCalculus:-`-`(`*`(`^`(2, `/`(1, 5))))) = 0 are So x = `&+-`(`*`(`^`(2, `/`(1, 10)))).  

 

 

Example: The temperature at a point x, y on a certain metal plate is T(x, y) = VectorCalculus:-`+`(VectorCalculus:-`+`(`*`(`^`(x, 4)), `*`(`^`(y, 4))), VectorCalculus:-`-`(VectorCalculus:-`*`(VectorCalculus:-`*`(4, x), y))) degrees Fahrenheit. Let begin by defining the function T(x, y) and graphing the surface z = T(x, y). 

 

> T := (x,y) -> x^4+y^4-4*x*y: 'T(x,y)' = T(x,y); plot3d(T(x,y),x=-2..2,y=-2..2,orientation=[100,75],axes=boxed,viewpoint=circleleft);
 

 

T(x, y) = `+`(`*`(`^`(x, 4)), `*`(`^`(y, 4)), `-`(`*`(4, `*`(x, `*`(y)))))
Plot_2d
 

 

Now let's find all of the function T(x, y)'s critical points. Define "critPts" to be the list of all REAL critical points (notice that the roots of VectorCalculus:-`+`(`*`(`^`(Z, 2)), 1) = 0 are imaginary). 

 

> Tx := diff(T(x,y),x); Ty := diff(T(x,y),y); soln := solve({Tx=0,Ty=0},[x,y]); critPts := [soln[1],soln[2],soln[3]];
 

 

 

 

`+`(`*`(4, `*`(`^`(x, 3))), `-`(`*`(4, `*`(y))))
`+`(`*`(4, `*`(`^`(y, 3))), `-`(`*`(4, `*`(x))))
[[x = 0, y = 0], [x = 1, y = 1], [x = -1, y = -1], [x = RootOf(`+`(`*`(`^`(_Z, 2)), 1)), y = `+`(`-`(RootOf(`+`(`*`(`^`(_Z, 2)), 1))))], [x = RootOf(`+`(`-`(RootOf(`+`(`*`(`^`(_Z, 2)), 1))), `*`(`^`(_...
[[x = 0, y = 0], [x = 1, y = 1], [x = -1, y = -1], [x = RootOf(`+`(`*`(`^`(_Z, 2)), 1)), y = `+`(`-`(RootOf(`+`(`*`(`^`(_Z, 2)), 1))))], [x = RootOf(`+`(`-`(RootOf(`+`(`*`(`^`(_Z, 2)), 1))), `*`(`^`(_...
[[x = 0, y = 0], [x = 1, y = 1], [x = -1, y = -1]] (28)
 

 

Next, let's define the discriminant "D = VectorCalculus:-`+`(VectorCalculus:-`*`(T[xx], T[yy]), VectorCalculus:-`-`(`*`(`^`(T[xy], 2))))" and test all of the REAL critical points. 

 

Note: We will assign the discriminant the name "DD" since "D" is reserved by Maple for another operator. 

 

> Txx := diff(T(x,y),x,x); Txy := diff(T(x,y),x,y); Tyy := diff(T(x,y),y,y); DD := Txx*Tyy-(Txy)^2; 'DD' = simplify(subs(critPts[1],DD)); 'Txx' = simplify(subs(critPts[1],Txx)); 'DD' = simplify(subs(critPts[2],DD));
'Txx' = simplify(subs(critPts[2],Txx)); 'DD' = simplify(subs(critPts[3],DD));
'Txx' = simplify(subs(critPts[3],Txx));
 

 

 

 

 

 

 

 

 

 

`+`(`*`(12, `*`(`^`(x, 2))))
-4
`+`(`*`(12, `*`(`^`(y, 2))))
`+`(`*`(144, `*`(`^`(x, 2), `*`(`^`(y, 2)))), `-`(16))
DD = -16
Txx = 0
DD = 128
Txx = 12
DD = 128
Txx = 12 (29)
 

 

The discriminant is negative for our first critial point 0, 0, so it is a saddle point. We can see that the discriminant is positive and the second partial of T(x, y) is positive for our last two critical points 1, 1 and -1, -1. Therefore, these points are local minima. 

 

 

Suppose that our metal plate is circular with radius 2 and is centered at 0, 0. If an ant is walking around on the plate, what point will it find to be the hottest and what point will it find to be the coldest? 

 

To solve this problem we'll need to use Lagrange multipliers to find the coldest and hottest points on the edge of the plate (our constraint is the circle VectorCalculus:-`+`(`*`(`^`(x, 2)), `*`(`^`(y, 2))) = `^`(2, 2)). And then we'll need to combine these results with our previous critical points results. 

 

> g := (x,y) -> x^2+y^2: g(x,y) = 2^2; EqX := diff(T(x,y),x) = lambda*diff(g(x,y),x); EqY := diff(T(x,y),y) = lambda*diff(g(x,y),y); EqConstraint := g(x,y) = 2^2; lagrangeCrit := solve({EqX,EqY,EqConstraint},[x,y,lambda]); evalf(subs(lagrangeCrit[1],T(x,y))); evalf(subs(lagrangeCrit[2],T(x,y)));
evalf(subs(lagrangeCrit[3],T(x,y))); evalf(subs(critPts[1],T(x,y))); evalf(subs(critPts[2],T(x,y))); evalf(subs(critPts[3],T(x,y)));
 

 

 

 

 

 

 

 

 

 

 

`+`(`*`(`^`(x, 2)), `*`(`^`(y, 2))) = 4
`+`(`*`(4, `*`(`^`(x, 3))), `-`(`*`(4, `*`(y)))) = `+`(`*`(2, `*`(lambda, `*`(x))))
`+`(`*`(4, `*`(`^`(y, 3))), `-`(`*`(4, `*`(x)))) = `+`(`*`(2, `*`(lambda, `*`(y))))
`+`(`*`(`^`(x, 2)), `*`(`^`(y, 2))) = 4
[[x = RootOf(`+`(`*`(`^`(_Z, 2)), `-`(2))), y = `+`(`-`(RootOf(`+`(`*`(`^`(_Z, 2)), `-`(2))))), lambda = 6], [x = RootOf(`+`(`*`(`^`(_Z, 2)), `-`(2))), y = RootOf(`+`(`*`(`^`(_Z, 2)), `-`(2))), lambda...
[[x = RootOf(`+`(`*`(`^`(_Z, 2)), `-`(2))), y = `+`(`-`(RootOf(`+`(`*`(`^`(_Z, 2)), `-`(2))))), lambda = 6], [x = RootOf(`+`(`*`(`^`(_Z, 2)), `-`(2))), y = RootOf(`+`(`*`(`^`(_Z, 2)), `-`(2))), lambda...
[[x = RootOf(`+`(`*`(`^`(_Z, 2)), `-`(2))), y = `+`(`-`(RootOf(`+`(`*`(`^`(_Z, 2)), `-`(2))))), lambda = 6], [x = RootOf(`+`(`*`(`^`(_Z, 2)), `-`(2))), y = RootOf(`+`(`*`(`^`(_Z, 2)), `-`(2))), lambda...
15.99999999
-0.4e-8
HFloat(18.00000000000001)
0.
-2.
-2. (30)
 

 

First, notice that all three critical points ("critPts") from before are inside the circle of radius 2, so we should consider all of those points. Now comparing the temperature at all of the internal critical points and all of the points found by the method of Lagrange multipliers, we find that the maximum temperature is approximately 18 degrees and occurs at "lagrangeCrit[3]" which is located approximately at .5176, -1.9319. The lowest temperature is -2 degrees and this occurs at "critPts[2]" and "critPts[3]" which are the points 1, 1and -1, -1. 

 

 

Example: Let's solve the following optimization problem using Lagrange multipliers: We want to find the min/max values of f(x, y) = xy subject to the constraint `and`(g(x, y) = VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`+`(`*`(`^`(y, 4)), VectorCalculus:-`-`(VectorCalculus:-`*`(2, xy))), VectorCalculus:-`-`(VectorCalculus:-`*`(6, `*`(`^`(x, 2))).... Moreover, we want to find where the min/max values occur and create a plot showing the relevant level curves of f and gas well as a few gradient vectors. 

      

 

The objective function is  f(x, y) and the constraint equation is g(x, y) = 5.  

 

Note: If we solve the Lagrange multiplier equations exactly, Maple's (exact symbolic) solution will be ugly and unuseful. 

 

> f := (x,y) -> x*y:
'f(x,y)'=f(x,y);

g := (x,y) -> y^4-2*x*y-6*x^2+x^4:
g(x,y)=5; fGrad := Gradient(f(x,y),[x,y]);
gGrad := Gradient(g(x,y),[x,y]); gGrad = lambda*fGrad; solve({-2*y-12*x+4*x^3=lambda*y,4*y^3-2*x=lambda*x,g(x,y)=5});
 

 

 

 

 

 

f(x, y) = `*`(x, `*`(y))
`+`(`*`(`^`(y, 4)), `-`(`*`(2, `*`(x, `*`(y)))), `-`(`*`(6, `*`(`^`(x, 2)))), `*`(`^`(x, 4))) = 5
Vector[column](%id = 18446744078120428238)
Vector[column](%id = 18446744078120428358)
Vector[column](%id = 18446744078120428358) = Vector[column](%id = 18446744078120428478)
{lambda = `+`(`-`(`/`(459, 2)), `-`(`*`(`/`(11787, 10), `*`(RootOf(`+`(625, `*`(4500, `*`(_Z)), `*`(11150, `*`(`^`(_Z, 2))), `*`(9228, `*`(`^`(_Z, 3))), `-`(`*`(2575, `*`(`^`(_Z, 4)))), `-`(`*`(3672, ...
{lambda = `+`(`-`(`/`(459, 2)), `-`(`*`(`/`(11787, 10), `*`(RootOf(`+`(625, `*`(4500, `*`(_Z)), `*`(11150, `*`(`^`(_Z, 2))), `*`(9228, `*`(`^`(_Z, 3))), `-`(`*`(2575, `*`(`^`(_Z, 4)))), `-`(`*`(3672, ...
{lambda = `+`(`-`(`/`(459, 2)), `-`(`*`(`/`(11787, 10), `*`(RootOf(`+`(625, `*`(4500, `*`(_Z)), `*`(11150, `*`(`^`(_Z, 2))), `*`(9228, `*`(`^`(_Z, 3))), `-`(`*`(2575, `*`(`^`(_Z, 4)))), `-`(`*`(3672, ...
{lambda = `+`(`-`(`/`(459, 2)), `-`(`*`(`/`(11787, 10), `*`(RootOf(`+`(625, `*`(4500, `*`(_Z)), `*`(11150, `*`(`^`(_Z, 2))), `*`(9228, `*`(`^`(_Z, 3))), `-`(`*`(2575, `*`(`^`(_Z, 4)))), `-`(`*`(3672, ...
{lambda = `+`(`-`(`/`(459, 2)), `-`(`*`(`/`(11787, 10), `*`(RootOf(`+`(625, `*`(4500, `*`(_Z)), `*`(11150, `*`(`^`(_Z, 2))), `*`(9228, `*`(`^`(_Z, 3))), `-`(`*`(2575, `*`(`^`(_Z, 4)))), `-`(`*`(3672, ...
{lambda = `+`(`-`(`/`(459, 2)), `-`(`*`(`/`(11787, 10), `*`(RootOf(`+`(625, `*`(4500, `*`(_Z)), `*`(11150, `*`(`^`(_Z, 2))), `*`(9228, `*`(`^`(_Z, 3))), `-`(`*`(2575, `*`(`^`(_Z, 4)))), `-`(`*`(3672, ...
{lambda = `+`(`-`(`/`(459, 2)), `-`(`*`(`/`(11787, 10), `*`(RootOf(`+`(625, `*`(4500, `*`(_Z)), `*`(11150, `*`(`^`(_Z, 2))), `*`(9228, `*`(`^`(_Z, 3))), `-`(`*`(2575, `*`(`^`(_Z, 4)))), `-`(`*`(3672, ...
{lambda = `+`(`-`(`/`(459, 2)), `-`(`*`(`/`(11787, 10), `*`(RootOf(`+`(625, `*`(4500, `*`(_Z)), `*`(11150, `*`(`^`(_Z, 2))), `*`(9228, `*`(`^`(_Z, 3))), `-`(`*`(2575, `*`(`^`(_Z, 4)))), `-`(`*`(3672, ...
{lambda = `+`(`-`(`/`(459, 2)), `-`(`*`(`/`(11787, 10), `*`(RootOf(`+`(625, `*`(4500, `*`(_Z)), `*`(11150, `*`(`^`(_Z, 2))), `*`(9228, `*`(`^`(_Z, 3))), `-`(`*`(2575, `*`(`^`(_Z, 4)))), `-`(`*`(3672, ...
{lambda = `+`(`-`(`/`(459, 2)), `-`(`*`(`/`(11787, 10), `*`(RootOf(`+`(625, `*`(4500, `*`(_Z)), `*`(11150, `*`(`^`(_Z, 2))), `*`(9228, `*`(`^`(_Z, 3))), `-`(`*`(2575, `*`(`^`(_Z, 4)))), `-`(`*`(3672, ...
{lambda = `+`(`-`(`/`(459, 2)), `-`(`*`(`/`(11787, 10), `*`(RootOf(`+`(625, `*`(4500, `*`(_Z)), `*`(11150, `*`(`^`(_Z, 2))), `*`(9228, `*`(`^`(_Z, 3))), `-`(`*`(2575, `*`(`^`(_Z, 4)))), `-`(`*`(3672, ...
{lambda = `+`(`-`(`/`(459, 2)), `-`(`*`(`/`(11787, 10), `*`(RootOf(`+`(625, `*`(4500, `*`(_Z)), `*`(11150, `*`(`^`(_Z, 2))), `*`(9228, `*`(`^`(_Z, 3))), `-`(`*`(2575, `*`(`^`(_Z, 4)))), `-`(`*`(3672, ...
{lambda = `+`(`-`(`/`(459, 2)), `-`(`*`(`/`(11787, 10), `*`(RootOf(`+`(625, `*`(4500, `*`(_Z)), `*`(11150, `*`(`^`(_Z, 2))), `*`(9228, `*`(`^`(_Z, 3))), `-`(`*`(2575, `*`(`^`(_Z, 4)))), `-`(`*`(3672, ...
{lambda = `+`(`-`(`/`(459, 2)), `-`(`*`(`/`(11787, 10), `*`(RootOf(`+`(625, `*`(4500, `*`(_Z)), `*`(11150, `*`(`^`(_Z, 2))), `*`(9228, `*`(`^`(_Z, 3))), `-`(`*`(2575, `*`(`^`(_Z, 4)))), `-`(`*`(3672, ...
{lambda = `+`(`-`(`/`(459, 2)), `-`(`*`(`/`(11787, 10), `*`(RootOf(`+`(625, `*`(4500, `*`(_Z)), `*`(11150, `*`(`^`(_Z, 2))), `*`(9228, `*`(`^`(_Z, 3))), `-`(`*`(2575, `*`(`^`(_Z, 4)))), `-`(`*`(3672, ...
{lambda = `+`(`-`(`/`(459, 2)), `-`(`*`(`/`(11787, 10), `*`(RootOf(`+`(625, `*`(4500, `*`(_Z)), `*`(11150, `*`(`^`(_Z, 2))), `*`(9228, `*`(`^`(_Z, 3))), `-`(`*`(2575, `*`(`^`(_Z, 4)))), `-`(`*`(3672, ...
{lambda = `+`(`-`(`/`(459, 2)), `-`(`*`(`/`(11787, 10), `*`(RootOf(`+`(625, `*`(4500, `*`(_Z)), `*`(11150, `*`(`^`(_Z, 2))), `*`(9228, `*`(`^`(_Z, 3))), `-`(`*`(2575, `*`(`^`(_Z, 4)))), `-`(`*`(3672, ...
{lambda = `+`(`-`(`/`(459, 2)), `-`(`*`(`/`(11787, 10), `*`(RootOf(`+`(625, `*`(4500, `*`(_Z)), `*`(11150, `*`(`^`(_Z, 2))), `*`(9228, `*`(`^`(_Z, 3))), `-`(`*`(2575, `*`(`^`(_Z, 4)))), `-`(`*`(3672, ...
{lambda = `+`(`-`(`/`(459, 2)), `-`(`*`(`/`(11787, 10), `*`(RootOf(`+`(625, `*`(4500, `*`(_Z)), `*`(11150, `*`(`^`(_Z, 2))), `*`(9228, `*`(`^`(_Z, 3))), `-`(`*`(2575, `*`(`^`(_Z, 4)))), `-`(`*`(3672, ...
(31)
 

 

Let's actually solve the equations using fsolve (the numerical solver). To help Maple out we should provide it with some guesses. We'll our guesses by looking at the contour plot.  

 

Notes: Since we need Typesetting:-delayGradient(f, cartesian) = VectorCalculus:-`*`(lambda, Typesetting:-delayGradient(g, cartesian)), as a guess for lambda we should use λ=VectorCalculus:-`*`(abs(Typesetting:-delayGradient(f, cartesian)), `/`(1, `*`(abs(Typesetting:-delayGradient(g, cartesian))))). Also, We will plot normalized gradient vectors (so their lengths don't overwhelm the picture). 

 

We will find all the points at which min/max values occur and then add the corresponding gradient vectors of f and g. Also, we'll only show contours for level curves corresponding to min/max values. 

 

> fsolve({-2*y-12*x+4*x^3=lambda*y,4*y^3-2*x=lambda*x,g(x,y)=5},{x= 2,y=-2,lambda=Norm(evalVF(gGrad,< 2,-2>))/Norm(evalVF(fGrad,< 2,-2>))});
fsolve({-2*y-12*x+4*x^3=lambda*y,4*y^3-2*x=lambda*x,g(x,y)=5},{x= 2,y= 2,lambda=Norm(evalVF(gGrad,< 2, 2>))/Norm(evalVF(fGrad,< 2, 2>))});
fsolve({-2*y-12*x+4*x^3=lambda*y,4*y^3-2*x=lambda*x,g(x,y)=5},{x=-2,y= 2,lambda=Norm(evalVF(gGrad,<-2, 2>))/Norm(evalVF(fGrad,<-2, 2>))});
fsolve({-2*y-12*x+4*x^3=lambda*y,4*y^3-2*x=lambda*x,g(x,y)=5},{x=-2,y=-2,lambda=Norm(evalVF(gGrad,<-2,-2>))/Norm(evalVF(fGrad,<-2,-2>))});
x1 :=  2.081433147; y1 := -1.550018748; 'f(x1,y1)'=f(x1,y1);
x2 :=  2.401211104; y2 :=  1.998347190; 'f(x2,y2)'=f(x2,y2);
x3 := -2.081433147; y3 :=  1.550018748; 'f(x3,y3)'=f(x3,y3);
x4 := -2.401211104; y4 := -1.998347190; 'f(x4,y4)'=f(x4,y4); constraint := contourplot(g(x,y),x=-3..3,y=-3..3,contours=[5]):
objective := contourplot(f(x,y),x=-3..3,y=-3..3,contours=[f(x1,y1),f(x2,y2),f(x3,y3),f(x4,y4)],color=black):
gGradplot1 := plots[arrow](<x1,y1>,Normalize(evalVF(gGrad,<x1,y1>)),shape=arrow,thickness=5,color=green):
fGradplot1 := plots[arrow](<x1,y1>,Normalize(evalVF(fGrad,<x1,y1>)),shape=arrow,thickness=1,color=blue):
gGradplot2 := plots[arrow](<x2,y2>,Normalize(evalVF(gGrad,<x2,y2>)),shape=arrow,thickness=5,color=green):
fGradplot2 := plots[arrow](<x2,y2>,Normalize(evalVF(fGrad,<x2,y2>)),shape=arrow,thickness=1,color=blue):
gGradplot3 := plots[arrow](<x3,y3>,Normalize(evalVF(gGrad,<x3,y3>)),shape=arrow,thickness=5,color=green):
fGradplot3 := plots[arrow](<x3,y3>,Normalize(evalVF(fGrad,<x3,y3>)),shape=arrow,thickness=1,color=blue):
gGradplot4 := plots[arrow](<x4,y4>,Normalize(evalVF(gGrad,<x4,y4>)),shape=arrow,thickness=5,color=green):
fGradplot4 := plots[arrow](<x4,y4>,Normalize(evalVF(fGrad,<x4,y4>)),shape=arrow,thickness=1,color=blue):
display(constraint,objective,gGradplot1,fGradplot1,gGradplot2,fGradplot2,gGradplot3,fGradplot3,gGradplot4,fGradplot4);
 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

{lambda = -9.156626931, x = 2.081433147, y = -1.550018748}
{lambda = 11.29359614, x = 2.401211104, y = 1.998347190}
{lambda = -9.156626931, x = -2.081433147, y = 1.550018748}
{lambda = 11.29359614, x = -2.401211104, y = -1.998347190}
2.081433147
-1.550018748
f(x1, y1) = -3.226260401
2.401211104
1.998347190
f(x2, y2) = 4.798453462
-2.081433147
1.550018748
f(x3, y3) = -3.226260401
-2.401211104
-1.998347190
f(x4, y4) = 4.798453462
Plot_2d
 

 

The maximum value of f(x, y) constrained to g(x, y) = 5 is 4.798453462. This occurs when ) and ). 

 

The minimum value of f(x, y) constrained to g(x, y) = 5 is -3.226260401. This occurs when (x1, y1) = (2.081433147, -1.550018748) and (x3, y3) = (-2.081433147, 1.550018748). 

 

 

Example: Here we classify critical points of a function of 3 variables. We can no longer use the `^`(2, nd) derivative test from our textbook. Instead we'll have to look directly at eigenvalues of the Hessian matrix (at critical point). Consider the following function... 

 

> h := (x,y,z) -> x*y*z+x^2*y+23/6*x^3-41/4*x^2-x-2/3*y^3+y^2+z^3/9-7/6*z^2+2*z:
'h(x,y,z)' = h(x,y,z);
 

h(x, y, z) = `+`(`*`(x, `*`(y, `*`(z))), `*`(`^`(x, 2), `*`(y)), `*`(`/`(23, 6), `*`(`^`(x, 3))), `-`(`*`(`/`(41, 4), `*`(`^`(x, 2)))), `-`(x), `-`(`*`(`/`(2, 3), `*`(`^`(y, 3)))), `*`(`^`(y, 2)), `*`... (32)
 

 

First, let's compute the Hessian. 

 

> H := Hessian(h(x,y,z),[x,y,z]);
 

Matrix(%id = 18446744078122124822) (33)
 

 

Notice that the gradient of h at the point (1,2,3) is 0. Thus (1,2,3) is a critical point. 

 

> evalVF(Gradient(h(x,y,z),[x,y,z]),<1,2,3>);
 

Vector[column](%id = 18446744078122125782) (34)
 

 

Next. let's plug in the critical point 1, 2, 3 and see what type of point this is. 

 

> H1 := subs({x=1,y=2,z=3},H); "Eigenvalues" = evalf(LinearAlgebra[Eigenvalues](H1));
 

 

Matrix(%id = 18446744078122126022)
Eigenvalues (35)
 

 

Since we have a mixture of positive and negative eigenvalues, 1, 2, 3is a saddle point. 

 

 

Now let's verify that (0,1,1) is a critical point and see what type it is.
 

> "Gradient" = evalVF(Gradient(h(x,y,z),[x,y,z]),<0,1,1>); H2 := subs({x=0,y=1,z=1},H); "Eigenvalues" = evalf(LinearAlgebra[Eigenvalues](H2));
 

 

 

Gradient
Matrix(%id = 18446744078127651406)
Eigenvalues (36)
 

 

Keeping in mind that the imaginary parts of the numbers are round-off error (`^`(10, -9) and `^`(10, -10) are very very small), all of the eigenvalues are negative at this critical point. This is a relative maximum. 

 

 

Again, let's verify that (2,-1,0) is a critical point and see what type it is. 

 

> "Gradient" = evalVF(Gradient(h(x,y,z),[x,y,z]),<2,-1,0>);

H2 := subs({x=2,y=-1,z=0},H);

"Eigenvalues" = evalf(LinearAlgebra[Eigenvalues](H2));
 

 

 

Gradient
Matrix(%id = 18446744078120398710)
Eigenvalues (37)
 

 

Again, keeping in mind that the imaginary parts of the numbers are round-off error, we have two positive eigenvalues and one negative eigenvalue. This is a saddle point.