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 . We wish to find and classify all of the critical points of .
> | 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); |
First, we should find the critical points. This is done by solving the system of equations: and
> | solve({diff(f(x,y),x)=0,diff(f(x,y),y)=0}); |
(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); |
(2) |
(0,2) is a relative maximum [negative eigenvalues OR positive determinant and negative diagonal entries]
Example: Let . Let's find and classify all of the critical points of .
> | 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); |
First, we should find the critical points. This is done by solving the system of equations: and
> | solve({diff(f(x,y),x)=0,diff(f(x,y),y)=0}); |
(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)); |
(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)); |
(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)); |
(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)); |
(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)); |
(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)); |
(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); |
The following command finds where and (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)})); |
(10) |
Note: 1.224744872 is actually
The following command computes the Hessian matrix of 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)); |
(11) |
The critical point 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)); |
(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); |
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); |
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); |
First, we will plot the plane tangent to 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); |
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 . Maple tracks down the (nearest) solution . 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 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])))); |
(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 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])))); |
(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 . 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])))); |
(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])))); |
(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 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); |
Example: Let's find the points on the level surface 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); |
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); |
(17) |
So 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 . 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 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); |
(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 and 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 -- 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]); |
(19) |
Note: "RootOf()" stands for a number which is a solution of this means or . So we actually have FOUR solutions!
> | plot(Z^10-2,Z=-1.1..1.1); |
The above graph shows that has TWO -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 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; |
(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)); |
(21) |
Therefore, at the critical points, we have D=80 > 0 and . 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)); |
(22) |
Note...
> | z=2/(2^(1/10))^4; |
(23) |
So the four solutions are and .
Next, let's solve the problem again. But this time we'll use Lagrange multipliers. First: our objective function is the distance from the origin (squared) and our constraint is the level surface defined by .
> | 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; |
(24) |
Next, we need to compute the gradient of and the gradient of 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; |
(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>; |
(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))); |
(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 is . Then the roots of are So .
Example: The temperature at a point on a certain metal plate is degrees Fahrenheit. Let begin by defining the function and graphing the surface .
> | 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); |
Now let's find all of the function 's critical points. Define "critPts" to be the list of all REAL critical points (notice that the roots of 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]]; |
(28) |
Next, let's define the discriminant "" and test all of the REAL critical points.
Note: We will assign the discriminant the name "" since "" 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)); |
(29) |
The discriminant is negative for our first critial point , so it is a saddle point. We can see that the discriminant is positive and the second partial of is positive for our last two critical points and . Therefore, these points are local minima.
Suppose that our metal plate is circular with radius 2 and is centered at . 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 ). 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))); |
(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 . The lowest temperature is -2 degrees and this occurs at "critPts[2]" and "critPts[3]" which are the points and .
Example: Let's solve the following optimization problem using Lagrange multipliers: We want to find the min/max values of subject to the constraint . Moreover, we want to find where the min/max values occur and create a plot showing the relevant level curves of and as well as a few gradient vectors.
The objective function is and the constraint equation is .
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}); |
(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 , as a guess for lambda we should use λ=. 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 and . 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); |
The maximum value of constrained to is . This occurs when ) and ).
The minimum value of constrained to is . This occurs when and .
Example: Here we classify critical points of a function of 3 variables. We can no longer use the 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); |
(32) |
First, let's compute the Hessian.
> | H := Hessian(h(x,y,z),[x,y,z]); |
(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>); |
(34) |
Next. let's plug in the critical point and see what type of point this is.
> | H1 := subs({x=1,y=2,z=3},H); "Eigenvalues" = evalf(LinearAlgebra[Eigenvalues](H1)); |
(35) |
Since we have a mixture of positive and negative eigenvalues, is 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)); |
(36) |
Keeping in mind that the imaginary parts of the numbers are round-off error ( and 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)); |
(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.