# Not able to get a solution to a Lagrange Multiplier Equation

I assigned several constrained optimization problems to my calc IV students to solve with Mathematica, but this one is being stubborn.

The function to optimize is: f(x,y,z) = y Exp[x-z] with the constraints: 9x^2 + 4y^2 + 36z^2 = 36 and xy + yz =1.

This obviously requires the gradient of the functions and two Lagrange multipliers. After storing the gradients into variables, here's the system of equations needed to solve:

gradf = Grad[y E^(x - z), {x, y, z}];

gradg = Grad[9 x^2 + 4 y^2 + 36 z^2 - 36, {x, y, z}];

gradh = Grad[x y + y z - 1, {x, y, z}];

NSolve[{
9 x^2 + 4 y^2 + 36 z^2 == 36,
x y + y z == 1},
{x, y, z}, Reals
]


(note, the L and M are lambda and mu)

There are four critical values Mathematica should find, but no matter the solve technique I try or options used for the various solve methods, Mathematica never stops running. Built-in commands FindMaximum and FindMaxValues, etc., eventually find the solutions, but the original systems still remains illusive.

Can anyone shed some light on what options to use. Unfortunately Maple wins this round with the "allvalues" command.

Thanks for the help.

• Plot[Evaluate[y Exp[(x - z)] /. # & /@ Solve[9 x^2 + 4 y^2 + 36 z^2 == 36 && x y + y z == 1, {x, y, z}]], {y, -3, 3}] Oct 16 '14 at 23:30

gradf = Grad[y E^(x - z), {x, y, z}]
gradg = Grad[9 x^2 + 4 y^2 + 36 z^2 - 36, {x, y, z}]
gradh = Grad[x y + y z - 1, {x, y, z}]


I think you can use following this, firstly I eliminated L, M.

eq1 = Eliminate[{
9 x^2 + 4 y^2 + 36 z^2 == 36,
x y + y z == 1}, {L, M}];
sol = NSolve[eq1, Reals];

y E^(x - z) /. sol


{-0.0687504, 0.408449, -5.35065, 9.7938}

Visualization 1

I have shown the points of critical value on the valued graphics of y E^(x - z)

fig1 = ContourPlot3D[{9 x^2 + 4 y^2 + 36 z^2 == 36,
x y + y z == 1}, {x, -3, 3}, {y, -3, 3}, {z, -3, 3}, Mesh -> None,
ContourStyle -> Directive[Opacity[0.5]]];
fig2 = Graphics3D[{PointSize[0.03], Red, Point[{x, y, z}] /. sol}];
Show[fig1, fig2] Visualization 2

I made ColorFunction for the function y E^(x - z)

colorf = Module[{rs, f, x, y, z},
rs[v_] := Rescale[v, {0, 1}, {-3, 3}];
f = rs@y  E^(rs@x - rs@z);
Function @@ {{x, y, z}, ColorData["Temperature"][
If[f < 0,
Rescale[Log@(-f), {-4, 4}, {0, 1/2}],
Rescale[Log@(f), {-4, 4}, {1/2, 1}]]]
}];


and added functional Mesh lines as the function values modified with log since the values are exponentially very big .

fig3 = ContourPlot3D[
{9 x^2 + 4 y^2 + 36 z^2 == 36, x y + y z == 1},
{x, -3, 3}, {y, -3, 3}, {z, -3, 3},
Mesh -> {Join[-2^Range[-4, 9, 2], 2^Range[-4, 9, 2]]},
MeshFunctions -> {Function[{x, y, z}, y E^(x - z)]},
ColorFunction -> colorf,
ContourStyle -> Opacity[0.5]];
fig4 = Graphics3D[{PointSize[0.03], Red, Point[{x, y, z}] /. sol}];
Show[fig3, fig4, Boxed -> False, Axes -> False] • Replacing NSolve with Solve results in an amusingly large output. Oct 17 '14 at 0:32
• @DumpsterDoofus It is same with Solve[eq1, Reals] // N Oct 17 '14 at 1:19
• Yup, I was just thinking about how the OP assigned this as a homework assignment for his students, and that his students would probably be terrified by the long Root[...] expressions if they happened to try to symbolically solve it and were expecting something simple. Oct 17 '14 at 1:27
• @DumpsterDoofus This is not proper for hand-write solving, so he was explained with "my calc IV students to solve with Mathematica" Oct 17 '14 at 1:33
• Wow, Junho. Thanks for the insight on the Eliminate command. I have not used that before; time to play. Hard to believe that this immediately and effectively reduced the system to a solvable problem. Wish I had this kind of computing power when I was in school 30 years ago. Oct 17 '14 at 2:25

NSolve can deal primarily with polynomial systems. In this problem one can easily eliminate the exponential term from the gradient of f. Since stationary point conditions depend only on directions of gradients and not on their magnitudes, one can divide gradients by nonzero functions.

scaledgradf = Grad[y E^(x - z), {x, y, z}]/E^(x - z);
gradg = Grad[9 x^2 + 4 y^2 + 36 z^2 - 36, {x, y, z}];
gradh = Grad[x y + y z - 1, {x, y, z}];

(sol=NSolve[{
9 x^2 + 4 y^2 + 36 z^2 == 36,
x y + y z == 1},
{x, y, z, L, M}, Reals
]);//Timing


{0.107982, Null}

y E^(x - z) /. sol


{-5.35065, -0.0687504, 0.408449, 9.7938}

f[y_] = y Exp[(x - z)] /. Solve[
9 x^2 + 4 y^2 + 36 z^2 == 36 && x y + y z == 1,
{x, y, z}] // Simplify // Quiet;

cons = List @@ Reduce[
Cases[f[y], Sqrt[t_] -> t, Infinity][] > 0,
y];

ptsMax = {y /. #[], #[]} & /@
(Outer[NMaximize[{#1, #2}, y] &,
f[y],
cons, 1] // Flatten[#, 1] &)


{{-0.545867, -0.0687504}, {2.96569, 3.62876}, {-0.452387, -0.12009}, {1.76806, 9.7938}}

ptsMin = {y /. #[], #[]} & /@
(Outer[NMinimize[{#1, #2}, y] &,
f[y],
cons, 1] // Flatten[#, 1] &)


{{-2.96569, -2.41294}, {0.904622, 0.408449}, {-2.15701, -5.35065}, {0.452387, 1.70417}}

Plot[Evaluate[f[y]], {y, -3, 3},
Epilog -> {Red, PointSize[Medium],
Tooltip[Point[#], #] & /@ ptsMax,
Darker[Blue],
Tooltip[Point[#], #] & /@ ptsMin}] • Thanks Bob. I'm working through your code and learning new things. Interesting to see these solutions points while comparing it to my rudimentary plot: ContourPlot3D[{ 9 x^2 + 4 y^2 + 36 z^2 == 36, x y + y z == 1 }, {x, -4, 4}, {y, -4, 4}, {z, -2, 2}, BoxRatios -> {2, 3, 1} ] Oct 17 '14 at 2:28