# How can I implement the method of Lagrange multipliers to find constrained extrema?

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$?

-
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 at 14:53

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):

TraditionalForm[
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]]


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:

Show[
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]}]]


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.

-
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 =
EulerEquations[
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:

eqs//TableForm

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


And solve as the in the other answers!

{x,y,z}/.FullSimplify[Solve[eqs,{x,y,z,λ1,λ2}]]//TableForm
`
-