Not able to get a solution to a Lagrange Multiplier Equation

Question

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[{
  gradf == L gradg + M gradh,
  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.

Answer 1

These are your definded gradient of functions

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[{
   gradf == L gradg + M gradh,
   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]

Answer 2

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[{
  scaledgradf == L gradg + M gradh,
  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}

Answer 3

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][[1]] > 0,
    y];

ptsMax = {y /. #[[2]], #[[1]]} & /@
  (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 /. #[[2]], #[[1]]} & /@
  (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}]