Lagrange Multipliers with Maple

Copyright 2001, 2002 by James F. Hurley, University of Connecticut Department of Mathematics, Unit 3009, Storrs CT 06269-3009. All rights reserved.

The interactive version of this worksheet is in MSB 203 in the Math Lab folder inside the Server1 volume.


The method of Lagrange multipliers applies to constrained optimization problems, in which the object is find the maximum and minimum of a function
f : R^3 R on a level surface g ( x , y , z ) = k in R^3 . Maple's built-in routine for solving systems of equations is often helpful for such problems, because Lagrange's method involves solving a system of equations

(*) f ( x , y , z ) = lambda g ( x , y , z ),
g ( x , y , z ) = k


in the variables
x , y , z , and lambda . Since both gradient vectors have three coordinates, this is a system of four equations in four unknowns. Hand solution of such a system usually requires some insight in eliminating variables, and the complexity of the associated algebra can be daunting. Maple's solve command can't solve every such system, but it can handle virtually all exercises in texts, whose authors designed them to be solvable by hand.

Example 1 . As of June, 2001, FedEx Corporation's size limitations on rectangular boxes are maximum dimension (called the length of the box) at most 119 inches, and girth (length plus the perimeter of the ends) at most 165 inches. What is the box of maximum volume that FedEx will accept for shipment?

Solution. Let z be the maximum dimension of the box in inches, and let x and y denote the other two dimensions in inches. Then the problem is to


maximize
V = f ( x , y , z ) = xyz

subject to


0
z 119, 0 x z , 0 y z , and 2 x + 2 y + z 165.

To get the most out of the constrained variables x , y and z , replace the last inequality by the corresponding equation . That defines the level surface g ( x , y , z ) = 165 for the function g , where

g(x,y,z) = 2*x+2*y+z .


This problem is feasible for hand solution, but Maple provides a handy way to check the solution
x = y = 27.5, z = 55 from hand implementation of the Lagrange-multiplier method. The routine below uses Maple's partial differentiation commands diff(f(x, y), x) , diff(f(x, y), y) and diff(f(x, y), z) to produce the system of equations (*) above. (This capability is valuable in case of a discrepancy between a hand-computed answer and an answer book: you can check all the partial differentiation in the hand solution.) Note that since this routine doesn't put the partial derivatives together into the gradient vector, invoking the linalg package is not necessary. Execute the code to see Maple's answers.

> f := (x, y, z) -> x*y*z;
g := (x, y, z) -> 2*x + 2*y + z;
f_x := (x, y, z) -> diff(f(x, y, z), x):
f_y := (x, y, z) -> diff(f(x, y, z), y):
f_z := (x, y, z) -> diff(f(x, y, z), z):
evalf(f_x(x, y, z));
evalf(f_y(x, y, z));
evalf(f_z(x, y, z));
g_x := (x, y, z) -> diff(g(x, y, z), x):
g_y := (x, y, z) -> diff(g(x, y, z), y):
g_z := (x, y, z) -> diff(g(x, y, z), z):
evalf(g_x(x, y, z));
evalf(g_y(x, y, z));
evalf(g_z(x, y, z));

f := proc (x, y, z) options operator, arrow; x*y*z ...

g := proc (x, y, z) options operator, arrow; 2*x+2*...

y*z

x*z

x*y

2.

2.

1.


The next routine asks Maple to solve the system of equations (*) for the
f and g in this problem.

> solve( { evalf(f_x(x, y, z))= lambda*evalf(g_x(x, y, z)),
evalf(f_y(x, y, z))= lambda*evalf(g_y(x, y, z)),
evalf(f_z(x, y, z))= lambda*evalf(g_z(x, y, z)),
g(x, y, z) = 165});

{x = 0., y = 0., lambda = 0., z = 165.}, {x = 0., l...
{x = 0., y = 0., lambda = 0., z = 165.}, {x = 0., l...


Notice that Maple finds all solutions to the system, even those for which a variable is 0, which of course are not relevant in designing an actual box. Only the last solution has all positive dimensions, and it does agree with the solution obtained by hand calculation. The following example checks the hand solution to an example in the forthcoming second edition of Multivariable Calculus by the author of this worksheet.

Example 2 . Find the maximum and minimum of f ( x , y ) = 2*x-y subject to the constraint x^2+y^2 <= 16 .

Solution. First, find the extreme values on the bounding circle x^2+y^2 = 16 of the feasible region, which is the circular disk of radius 4 centered at the origin. To finish the problem, investigate the behavior of f inside the feasible region. Proceeding as before, re-code the formulas for f and g , and then execute the edited routine, which of course involves just the two variables x and y .

> f := (x, y) -> 2*x - y ;
g := (x, y) -> x^2 + y^2 - 16;
f_x := (x, y, z) -> diff(f(x, y), x):
f_y := (x, y, z) -> diff(f(x, y), y):
evalf(f_x(x, y));
evalf(f_y(x, y));
g_x := (x, y) -> diff(g(x, y), x):
g_y := (x, y) -> diff(g(x, y), y):
evalf(g_x(x, y));
evalf(g_y(x, y));

f := proc (x, y) options operator, arrow; 2*x-y end...

g := proc (x, y) options operator, arrow; x^2+y^2-1...

2.

-1.

2.*x

2.*y

Then, have Maple solve the system of equations

f ( x , y ) = lambda g ( x , y ),
g(x,y) = 0 .


Execute the code below to see Maple's solutions to this system.

> (ans1, ans2) := solve( {
evalf(f_x(x, y))= lambda*evalf(g_x(x, y)),
evalf(f_y(x, y))= lambda*evalf(g_y(x, y)),
g(x, y) = 0
} );

ans1, ans2 := {y = -1.788854382, x = 3.577708764, l...
ans1, ans2 := {y = -1.788854382, x = 3.577708764, l...

> subs(ans1, f(x, y));

8.944271910

> subs(ans2, f(x, y));

-8.944271910


From this, Maple's answers for extreme-value points and the corresponding maximum and minimum of
f subject to the two constraints are clear. How do they compare to the values found by hand:

minimum
-4*sqrt(5) at x = -8/sqrt(5) , y = 4/sqrt(5) ,
maximum
4*sqrt(5) at the point with coordinates that are the negatives of those values?

Maple converts those numbers to floating-point numbers that agree with the output from its
solve command except in the last decimal place. (That discrepancy is the result of floating-point roundoff error.)

> evalf([-8/sqrt(5), 4/sqrt(5), -4*sqrt(5)]);

[-3.577708765, 1.788854382, -8.944271912]

To identify the maximum and minimum of f on the entire closed disk , note that the gradient is never 0, so there can be no local maximum or minimum inside the circle. Thus, the values found above are the actual maximum and minimum of f on the entire feasible region.

The concluding example below (Example 10.7, Section 4.10 of the current edition of
Multivariable Calculus ) illustrates the Lagrange-multiplier method for optimizing a function subject to more than one constraint. .

Example 3 . Find the maximum and minimum of f(x,y,z) = x*y+x*z on the curve of intersection of the cylinders x^2+y^2 = 1 and xz = 1 .

Solution . Proceeding as before, re-code the formulas in Exampl 1 for f(x,y,z) = x*y+x*z and g(x,y,z) = x^2+y^2-1 , add the second constraint function h ( x , y , z ) = x*z-1 to the first routine, and then execute the edited routine:

> f:= (x, y, z) -> x*y + x*z;
g := (x, y, z) -> x^2 + y^2 - 1;
h := (x, y, z) -> x*z - 1;
f_x := (x, y, z) -> diff(f(x, y, z), x):
f_y := (x, y, z) -> diff(f(x, y, z), y):
f_z := (x, y, z) -> diff(f(x, y, z), z):
evalf(f_x(x, y, z));
evalf(f_y(x, y, z));
evalf(f_z(x, y, z));
g_x := (x, y, z) -> diff(g(x, y, z), x):
g_y := (x, y, z) -> diff(g(x, y, z), y):
g_z := (x, y, z) -> diff(g(x, y, z), z):
evalf(g_x(x, y, z));
evalf(g_y(x, y, z));
evalf(g_z(x, y, z));
h_x := (x, y, z) -> diff(h(x, y, z), x):
h_y := (x, y, z) -> diff(h(x, y, z), y):
h_z := (x, y, z) -> diff(h(x, y, z), z):
evalf(h_x(x, y, z));
evalf(h_y(x, y, z));
evalf(h_z(x, y, z));

f := proc (x, y, z) options operator, arrow; x*y+x*...

g := proc (x, y, z) options operator, arrow; x^2+y^...

h := proc (x, y, z) options operator, arrow; x*z-1 ...

y+z

x

x

2.*x

2.*y

0.

z

0.

x


Then, use the following routine to have Maple solve the system of four equations in the four variables
x , y , z and lambda . (If Maple gives an error message about multiple assignments, try adding or deleting answers: here we initially guess that there may be as many as four, which turns out to be accurate.)

f ( x , y , z ) = lambda g ( x , y , z ) + mu g ( x , y , z )
g(x,y,z) = 0 , h(x,y,z) = 1.

> (ans1, ans2, ans3, ans4) := solve( {
evalf(f_x(x, y, z))= lambda*evalf(g_x(x, y, z)) + mu*evalf(h_x(x, y, z)),
evalf(f_y(x, y, z))= lambda*evalf(g_y(x, y, z)) + mu*evalf(h_y(x, y, z)),
evalf(f_z(x, y, z))= lambda*evalf(g_z(x, y, z)) + mu*evalf(h_z(x, y, z)), g(x, y, z) = 0, h(x, y, z) = 0
} );

ans1, ans2, ans3, ans4 := {mu = 1., y = -.707106781...
ans1, ans2, ans3, ans4 := {mu = 1., y = -.707106781...
ans1, ans2, ans3, ans4 := {mu = 1., y = -.707106781...
ans1, ans2, ans3, ans4 := {mu = 1., y = -.707106781...
ans1, ans2, ans3, ans4 := {mu = 1., y = -.707106781...

> subs(ans1, f(x, y, z));

1.500000000

> subs(ans2, f(x, y, z));

1.500000000

> subs(ans3, f(x, y, z));

.4999999998

> subs(ans4, f(x, y, z));

.4999999998


From this, Maple's answers for extreme-value points and the corresponding maximum
f subject to the two constraints are clear. The minimum value appears so close to 0.5 that the discrepancy appears due to floating-point roundoff error.

How do these answers compare to the values found by hand? As you can check (or read in the text), hand computation gives minimum
1/2 at x = sqrt(1/2) , y = -sqrt(1/2) , z = sqrt(2) and also at (x, y, z) = (-1/sqrt(2), 1/sqrt(2), -sqrt(2)) ; the same hand work yields maximum value 3/2 at the points (x, y, z) = (1/sqrt(2), 1/sqrt(2), sqrt(2)) and (x, y, z) = (-1/sqrt(2), -1/sqrt(2), -sqrt(2)) . As the following shows, Maple converts those numbers to floating-point numbers that substantially agree with the output from the solve command.

> evalf([-1/sqrt(2), 1/sqrt(2), -sqrt(2), 1/2]);

[-.7071067810, .7071067810, -1.414213562, .50000000...