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.
f=4-x^2-y^2+x^4+y^4
g=exp(x^2)+2*exp(y^2)
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.
[xfcrit,yfcrit]=solve(gradf(1),gradf(2));
[xfcrit,yfcrit]
[ -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.
double(gfun(xfcrit,yfcrit))
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(xfcrit([1,2,3]),yfcrit([1,2,3]))
Inline function:
ffun(x,y) =
4-x.^2-y.^2+x.^4+y.^4
ans =
[ 4]
[ 15/4]
[ 15/4]
(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.
gradcross=cross(gradf,gradg);
lagr=gradcross(3)
[ 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)
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]
[ 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([[0;0],yaxiscrits])
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))
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)
0.6000
yb1 =
0.4994
[xb2,yb2]=newton2d(g-4,lagr,-.5,.5)
-0.6000
yb2 =
0.4994
[xb3,yb3]=newton2d(g-4,lagr,-.5,-.5)
-0.6000
yb3 =
-0.4994
[xb4,yb4]=newton2d(g-4,lagr,.5,-.5)
0.6000
yb4 =
-0.4994
ffun([xb1,xb2,xb3,xb4],[yb1,yb2,yb3,yb4])
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
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
phi=a*exp(-x^2)*(1+b*x^2+c*x^4)
a*exp(-x^2)*(1+b*x^2+c*x^4)
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.)
energy=simple(int(phiprime^2+x^4*phi^2,-inf,inf))
-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.
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])
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]))
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)))
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')
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.