7
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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}] $\endgroup$ – Dr. belisarius Oct 16 '14 at 23:30
6
$\begingroup$

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]

Blockquote

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]

Blockquote

| improve this answer | |
$\endgroup$
  • $\begingroup$ Replacing NSolve with Solve results in an amusingly large output. $\endgroup$ – DumpsterDoofus Oct 17 '14 at 0:32
  • $\begingroup$ @DumpsterDoofus It is same with Solve[eq1, Reals] // N $\endgroup$ – Junho Lee Oct 17 '14 at 1:19
  • $\begingroup$ 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. $\endgroup$ – DumpsterDoofus Oct 17 '14 at 1:27
  • $\begingroup$ @DumpsterDoofus This is not proper for hand-write solving, so he was explained with "my calc IV students to solve with Mathematica" $\endgroup$ – Junho Lee Oct 17 '14 at 1:33
  • $\begingroup$ 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. $\endgroup$ – ChazS Oct 17 '14 at 2:25
5
$\begingroup$

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}

| improve this answer | |
$\endgroup$
4
$\begingroup$
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}]

enter image description here

| improve this answer | |
$\endgroup$
  • $\begingroup$ 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} ] $\endgroup$ – ChazS Oct 17 '14 at 2:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.