Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I want to form the function $h=f-\lambda_{1}g_{1}-\lambda_{2}g_{2}$ where $f$ is the function to optimize subject to the constraints $g_{1}=0$ and $g_{2}=0$ so that I can solve the first partial derivatives with respect to $\lambda_{1}$ and $\lambda_{2}$. Can someone get me started using $f(x,y,z)=xy+yz$ subject to the constraints $x^2+y^2-2=0$ and $x^2+z^2-2=0$?

share|improve this question
Can someone explain how to solve the system of equations that follow from all of the partial derivatives set equal to zero? This is what I have so far: – Logan Nov 11 '13 at 23:16
Do you want just to optimize the target function subject to the stated constraints? If so, then why not just Minimize[{x y + y z, x^2 + y^2 - 2 == 0, x^2 + z^2 - 2 == 0}, {x, y}] and similarly with Maximize? Or do you insist on explicitly implementing the Lagrange method? – murray Jan 29 '14 at 14:53
up vote 26 down vote accepted

We define the function f and multiple constraint functions g1, g2:

f[x_, y_, z_] := x y + y z
g1[x_, y_] := x^2 + y^2 - 2
g2[x_, z_] := x^2 + z^2 - 2

then, in order to find necessary conditions for constrained extrema we introduce the Lagrange function h with Lagrange multipliers λ1 and λ2:

h[x_, y_, z_, λ1_, λ2_] := f[x, y, z] - λ1 g1[x, y] - λ2 g2[x, z]

Now we solve an appropriate system of equations satisfying necesary conditions (i.e. vanishing of all first derivatives of h):

  Column[ pts = {x, y, z} /. 
          FullSimplify @ Solve[ D[h[x, y, z, λ1, λ2], #] == 0 & /@ {x, y, z, λ1, λ2}, 
                                      {x, y, z, λ1, λ2}], Frame -> All]]

enter image description here

A bit nicer way of finding all the solutions uses Grad - a new function in Mathematica 9 for vector analysis:

{x, y, z} /. Solve[ Grad[ h @@ #, #] == 0, #]& @ {x, y, z, λ1, λ2} // FullSimplify

The above table contains all critical points of the Lagrange function h. For sufficient conditions one can use Maximize and Minimize, e.g.:

FullSimplify @  ToRadicals @ 
Maximize[{f[x, y, z], g1[x, y] == 0, g2[x, z] == 0}, {x, y, z}]
{1 + Sqrt[2], {x -> -(1/Sqrt[2 + Sqrt[2]]), 
               y -> -Sqrt[1 + 1/Sqrt[2]],
               z -> -Sqrt[1 + 1/Sqrt[2]]}}

We add a graphics with contours of constrained minima and maxima, the contraint functions ass well as all critical points of h:

  ContourPlot3D[{ f[x, y, z] ==  1 + Sqrt[2], 
                  f[x, y, z] == -1 - Sqrt[2], 
                  g1[x, y] == 0, g2[x, z] == 0}, 
                  {x, -2.3, 2.3}, {y, -2.3, 2.3}, {z, -2.3, 2.3}, 
                  ContourStyle -> {Directive[Cyan, Opacity[0.5]], 
                                   Directive[Green, Opacity[0.5]], 
                                   Directive[Orange, Opacity[0.15]], 
                                   Directive[Orange, Opacity[0.15]]}, Mesh -> None], 
  Graphics3D[{Magenta, PointSize[0.015], Point[pts]}]]

enter image description here

On the cyan surfaces we have maxima, on the green ones - minima and the solutions of the necessary conditions are denoted with the magenta points lying on the tube constraints.

share|improve this answer
Wow, awesome. That looks great. Thanks. – Logan Nov 12 '13 at 6:56

Another possible way (using a hammer to kill a fly perhaps...) with the VariationalMethods package

<< VariationalMethods`

f[x_, y_, z_] := x y + y z
g1[x_, y_] := x^2 + y^2 - 2
g2[x_, z_] := x^2 + z^2 - 2

eqs = 
   f[x[t], y[t], z[t]] - (λ1[t] g1[x[t], y[t]] + λ2[t] g2[x[t], z[t]]), 
   {x[t], y[t], z[t], λ1[t], λ2[t]}, t] /. x_[t] -> x;

See the resulting equations:


(* y-2 x (λ1+λ2)==0
   x+z-2 y λ1==0
   y-2 z λ2==0
   2-x^2-z^2==0 *)

And solve as the in the other answers!

share|improve this answer
gradient[g_, vars_] :=  Table[D[g@@vars, vars[[j]]], {j, 1, Length[vars]}]

system1[lstConst_, vars_] := Join[ Join@@ 
Table[gradient[lstConst[[j]], vars], {j, 1, Length[lstConst]}],      

system2[f_, lstConst_, vars_, lambda_] := Join[ gradient[f, vars] - 
Sum[ lambda[[j]]*gradient[lstConst[[j]], vars], {j, 1, 
Length[lstConst]}],Table[lstConst[[j]]@@vars, {j, 1, Length[lstConst]}]] ;     

criticalPointsSystem1[lstConst_, vars_] :=   Solve[system1[lstConst, vars] == 
 Table[ 0, {j, 1, (Length[vars] + 1)*Length[lstConst]}], 
vars] /. {(x_ -> y_) -> y} ;

criticalPointsSystem2[f_, lstConst_, vars_, lambda_] :=   
Map[ Function [x, Take[x, Length[vars]]], 
Solve[system2[f, lstConst, vars, lambda] == 
Table[0, {j, 1, Length[vars] + Length[lambda]}],
      Join[vars, lambda]]] /. {(x_ -> y_) -> y};

criticalPointsLagrangeM[f_, lstConst_, vars_, lambda_] := 
Join[criticalPointsSystem1[lstConst, vars], 
criticalPointsSystem2[f, lstConst, vars, lambda]];

optimizeByLagrangeM[f_, lstConst_, vars_, lambda_, type_] := 
Which[ToUpperCase[type] == "MINIMIZE",Min[Map[Function[x, f @@ x], 
criticalPointsLagrangeM[f, lstConst, vars, lambda]]],
ToUpperCase[type] == "MAXIMIZE", Max[Map[Function[x, f@@x],
criticalPointsLagrangeM[f, lstConst, vars, lambda]]],True, 
Print["The given type of optimization problem is not supported"]];

f[x_, y_] := x; g[x_, y_] := y^2 + x^4 - x^3; (* test  *)

optimizeByLagrangeM[f, {g}, {x, y}, {\[Lambda]}, "MiNimize"]

optimizeByLagrangeM[f, {g}, {x, y}, {\[Lambda] }, "Maximize"]
share|improve this answer
Hi ! Please, refer to the help centre and learn how to properly format your code and edit your post. – Sektor Mar 24 '15 at 19:25
Done. Thanks Sektor. – TAWA Mar 24 '15 at 20:10
Thank you ! Hope you stick around to contribute to the community :) ! + 1 BTW :D – Sektor Mar 24 '15 at 20:10

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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