Graphical-Numerical
Optimization Methods
and Lagrange Multipliers

copyright © 2000 by Paul Green and Jonathan Rosenberg

 

 

Constrained Extremal Problems in Two Variables

In this notebook, we will examine the problem of finding the extreme values of a function on a bounded region.  We will start by finding the extreme values of the function  on the region .  Extreme values can occur either at critical points of f interior to the region, or along the boundary. What we must do, therefore, is evaluate f at those critical points that satisfy the inequality defining the region, and compare those values to the maximum and minimum along the boundary. We will find the latter by using the method of Lagrange multipliers.  We begin by defining the functions f and g in MATLAB.

 

syms x y z

f=4-x^2-y^2+x^4+y^4

g=exp(x^2)+2*exp(y^2)  

 

f =

4-x^2-y^2+x^4+y^4

g =

exp(x^2)+2*exp(y^2)  

 

Next, we find the critical points of f. For reasons that will become apparent later, we include the derivative of f with respect to z, even though it is identically 0 since f does not depend on z.

 

gradf=jacobian(f,[x,y,z])

[xfcrit,yfcrit]=solve(gradf(1),gradf(2));

[xfcrit,yfcrit]  

 

gradf =

[ -2*x+4*x^3, -2*y+4*y^3,          0]

ans =

[            0,            0]

[  1/2*2^(1/2),            0]

[ -1/2*2^(1/2),            0]

[            0,  1/2*2^(1/2)]

[            0, -1/2*2^(1/2)]

[  1/2*2^(1/2),  1/2*2^(1/2)]

[ -1/2*2^(1/2), -1/2*2^(1/2)]  

 

 

Next, we evaluate g at the critical points to determine which of them are in the region.

 

gfun=inline(vectorize(g))

double(gfun(xfcrit,yfcrit))  

 

gfun =

     Inline function:

     gfun(x,y) = exp(x.^2)+2.*exp(y.^2)

ans =

    3.0000

    3.6487

    3.6487

    4.2974

    4.2974

    4.9462

    4.9462  

 

We see that only the first three critical points are in the region. Consequently we are only interested in the values of f at those points.

 

ffun=inline(vectorize(f))

ffun(xfcrit([1,2,3]),yfcrit([1,2,3]))  

 

ffun =

     Inline function:

     ffun(x,y) = 4-x.^2-y.^2+x.^4+y.^4

ans =

[    4]

[ 15/4]

[ 15/4]  

 

We must now address the problem of finding the extreme values of f on the boundary curve. The theory of Lagrange multipliers tells us that the extreme values will be found at points where the gradients of f and g are parallel. This will be the case precisely when the cross-product of the gradients of f and g is zero. It is for this reason that we have defined the gradient of f, and will define the gradient of g, to be formally three-dimensional. We observe that since the gradients of f and g lie in the plane (spanned by i and j), only the third component of the cross-product,

(df/dx)(dg/dy) - (df/dy)(dg/dx),

will ever be different from zero. This gives us a system of two equations, the solutions of which will give all possible locations for the extreme values of f on the boundary.

 

gradg=jacobian(g,[x,y,z])

gradcross=cross(gradf,gradg);

lagr=gradcross(3)  

 

gradg =

[ 2*x*exp(x^2), 4*y*exp(y^2),            0]

lagr =

4*(-2*x+4*x^3)*y*exp(y^2)-2*(-2*y+4*y^3)*x*exp(x^2)  

 

[xboundcrit,yboundcrit]=solve(g-4,lagr)  

 

xboundcrit =

0

yboundcrit =

.63676142165505314316845091584897  

 

MATLAB appears to have found a constrained critical point. However it cannot possibly be the only one, since f must take both a maximum and a minimum along the curve g=4. We shall have to resort to graphical-numerical methods. We begin by using genplot to  plot the loci of the two equations we are trying to solve.

 

genplot(g, -2:.1:2,-2:.1:2,'contour',[4,4],'r'); hold on;

genplot(lagr, -2:.05:2,-2:.05:2,'contour',[0,0],'b'); hold off;  

 

 

We see now that the critical point found by solve is one of eight. Four of them are on the axes, and can be easily be found analytically. The other four will require the use of the root finder newton2d. Let us deal with the first four.  We simply have to set x or y to zero and solve the equation g=4 for the other. This gives

 

xaxiscrits=solve(subs(g-4,y,0));

[xaxiscrits,[0;0]]

yaxiscrits=solve(subs(g-4,x,0));

[[0;0],yaxiscrits]  

 

ans =

[  log(2)^(1/2),             0]

[ -log(2)^(1/2),             0]

ans =

[               0,  log(3/2)^(1/2)]

[               0, -log(3/2)^(1/2)]  

The numerical values of these critical points are:

double([xaxiscrits,[0;0]])

double([[0;0],yaxiscrits])  

 

ans =

    0.8326         0

   -0.8326         0

ans =

         0    0.6368

         0   -0.6368  

and the corresponding f-values are:

 

double(ffun(xaxiscrits,[0;0]))

double(ffun([0;0],yaxiscrits))   

 

ans =

    3.7873

    3.7873

ans =

    3.7589

    3.7589  

 

The remaining four constrained critical points seem to be near .

 

[xb1,yb1]=newton2d(g-4,lagr,.5,.5)  

 

xb1 =

    0.6000

yb1 =

    0.4994  

[xb2,yb2]=newton2d(g-4,lagr,-.5,.5) 

 

xb2 =

   -0.6000

yb2 =

    0.4994  

 

[xb3,yb3]=newton2d(g-4,lagr,-.5,-.5) 

 

xb3 =

   -0.6000

yb3 =

   -0.4994  

 

[xb4,yb4]=newton2d(g-4,lagr,.5,-.5) 

 

xb4 =

    0.6000

yb4 =

   -0.4994  

ffun([xb1,xb2,xb3,xb4],[yb1,yb2,yb3,yb4])  

 

ans =

    3.5824    3.5824    3.5824    3.5824  

 

We now see that the maximum value of f in the region of interest is 4, at the interior critical point (0,0), and the minimum value is 3.5824, taken at the four points we just found.  We can check this graphically by superimposing a contour plot of f on a picture of the bounding curve.  Note that we select only those contours of f with values close the ones we want, in order to get the best possible picture, and have added the colorbar command to add a legend indicating what colors correspond to what f-values.

 

genplot(g, -2:.1:2,-2:.1:2,'contour',[4,4],'k'); hold on;

genplot(f, -2:.05:2,-2:.05:2,'contour',[3:.1:4.2]); hold off;

colorbar

 

Problem:  Find the minimum and maximum values of the function    on the region defined by .

 

A Physical Application

Here's an example from quantum mechanics that illustrates how the Lagrange multiplier method can be used.  Consider a particle of mass m free to move along the x-axis, sitting in a "potential well" given by a potential energy function V(x) that tends to +infinity as |x| tends to infinity. Then the "ground state" of the particle is given by a "wave function" phi(x), where phi(x)2 gives the probability of finding the function near x.  Since the total probability of finding the function somewhere is 1, we have the constraint that the integral of phi(x)2 must be 1. A basic principle of physics states that phi(x) will be the function minimizing the energy E(phi), which is given by the integral of

(h2/2m)phi'(x)2 + V(x)phi(x)2,

where h is Planck's constant. Usually it's not possible to find the function exactly, but a convenient way of approximating it is to guess a form for phi(x) involving some parameters, and then adjust the parameters to minimize E(phi) subject to the constraint that the integral of phi(x)2 is = 1. For simplicity let's take h=m=1.  Say for example that V(x) = x4/2.  Let's choose a reasonable form for a function phi(x) that dies rapidly at infinity. By symmetry, phi(x) has to be an even function of x, so let's try:

syms a b c  

phi=a*exp(-x^2)*(1+b*x^2+c*x^4)  

 

phi =

a*exp(-x^2)*(1+b*x^2+c*x^4)  

 

g=simple(int(phi^2,-inf,inf))  

 

g =

1/512*a^2*2^(1/2)*pi^(1/2)*(256+128*b+96*c+48*b^2+105*c^2+120*b*c)  

So our constraint is g=1.  Now let's compute the energy function. (For convenience we're leaving out the factor of ½ in both terms inside the integral.)

 

phiprime=simple(diff(phi))

energy=simple(int(phiprime^2+x^4*phi^2,-inf,inf)) 

 

phiprime =

-2*a*x*exp(-x^2)*(1+b*x^2+c*x^4-b-2*c*x^2)

energy =

1/8192*a^2*2^(1/2)*pi^(1/2)*(13995*c^2+10248*c*b+4864-128*b+3472*b^2-1248*c)  

Now we solve the Lagrange multipler equations.  This time we have a function E(a,b,c) that we need to minimize subject to a constraint g(a,b,c) = 1, so the equations are

grad E - lambda grad g = 0.

syms lam 

eg=jacobian(energy,[a,b,c]);

gg=jacobian(g,[a,b,c]);

lagr=eg-lam*gg;  

 

[asol,bsol,csol,lsol]=solve(g-1,lagr(1),lagr(2),lagr(3),a,b,c,lam);

double([asol,bsol,csol,lsol])  

 

ans =

   0.5971 - 0.0000i  -7.5249             4.8791 + 0.0000i  16.7155         

  -0.5971 + 0.0000i  -7.5249             4.8791 + 0.0000i  16.7155         

   0.7902 + 0.0000i   0.5222 + 0.0000i  -0.0577             1.0611         

  -0.7902 - 0.0000i   0.5222 + 0.0000i  -0.0577             1.0611         

   0.7177 - 0.0000i  -3.1700 + 0.0000i  -0.2087             7.5359         

  -0.7177 + 0.0000i  -3.1700 + 0.0000i  -0.2087             7.5359            

Clearly there's some small imaginary round-off error here.  We can get rid of it by taking the real parts:

real(double([asol,bsol,csol,lsol]))  

 

ans =

    0.5971   -7.5249    4.8791   16.7155

   -0.5971   -7.5249    4.8791   16.7155

    0.7902    0.5222   -0.0577    1.0611

   -0.7902    0.5222   -0.0577    1.0611

    0.7177   -3.1700   -0.2087    7.5359

   -0.7177   -3.1700   -0.2087    7.5359  

Now we need to look for the minimum value of the energy:

efun=inline(vectorize(energy));  

real(double(efun(asol,bsol,csol)))  

 

ans =

   16.7155

   16.7155

    1.0611

    1.0611

    7.5359

    7.5359  

The minimum energy for this choice of the form of phi(x) is thus 1.0611.  We can plot the probability density:

bestphi=subs(phi,[a,b,c],[asol(3),bsol(3),csol(3)]);

ezplot(bestphi^2); title('probability density')  

 

 

Additional Problems

1. Find the minimum value of the function f(x, y) = 4 - x2 - y2 + x4 + y4 subject to the constraint g(x, y) = x2 - 2x + y2 - 3 <= 0. You will find many interior critical points and many solutions to the Lagrange multiplier equations.

2. Find and plot the function phi(x) (of the same form as in the example above) minimizing the energy for the potential V(x) = (x2 + x4)/2, with the same constraint that the integral of phi(x)2 =1, and compare the result with the example above.  Again take h=m=1.