Information for Maple assignment 4


Each student will get individual data for the assignment.

Optimization using Maple


Introduction

Here is the problem:

Theoretical results imply that the function f(x,y,z) = x2+x+3yz has a maximum and a minimum on the (solid) ball defined by g(x,y,z) = x2+y2+z2 ≤ 1. Find these maximum and minimum values by

The critical point of the objective function

We first discuss the problem of finding critical points of the objective function and determining whether or not they are relevant to our problem. A critical point of f(x,y,z) is determined by the condition ∇f(x,y,z) = 0. Since our objective function is quadratic the resulting equations will be linear and will have a unique solution (or, in general, possibly an infinite number of solutions, but the latter case will not occur for data in this lab). In our example ∇f(x,y,z) = <2x+1,3z,3y>, so that the equations to be solved are 2x+1 = 0, 3z = 0, and 3y = 0.

Here is one way to get Maple to solve these (rather trivial, in this example) equations and to test the resulting critical point.

> f:=x^2+x+3*y*z;

                      f := x2 + x + 3 y z

> g:=x^2+y^2+z^2;

                      g := x2 + y2 + z2

> CPeqx:=diff(f,x)=0;

                      CPeqx := 2 x + 1 = 0

> CPeqy:=diff(f,y)=0;

                      CPeqy := 3 z = 0

> CPeqz:=diff(f,z)=0;

                      CPeqz := 3 y = 0

> CPsol:=solve({CPeqx,CPeqy,CPeqz});

                      CPsol := { x = -1/2, y = 0, z = 0 }

> CPgvalue:=subs(CPsol,g);

                       CPgvalue := 1/4

> CPfvalue:=subs(CPsol,f);

                       CPfvalue := -1/4

Comments on the Maple commands

The command CPeqx:=diff(f,x)=0 labels the equation ∂ f / ∂ x = 0 as CPeqx so that it can be referred to later. The = is part of an equation, and the symbols := are used by Maple to assign a value. In this case, the "label" CPeqx is assigned as its value the equation diff(f,x)=0. The solve step finds all solutions of the system of equations (in this case there is only one such solution). Finally, the commands CPgvalue:=... and CPfvalue:=... substitute the critical point into the constraint and objective functions.

In this case, because the value of g at the critical point is less than 1, the critical point lies inside the domain of interest, and CPfvalue is a candidate for an extreme value of the objective function.

Extrema on the boundary and Lagrange multipliers

The boundary of our region of interest is the surface g(x,y,z)= x2+y2+z2=1, a sphere which easy to understand and to visualize. The objective function f(x,y,z)=x2+x+3yz may take its maximum value, its minimum value, or both on this surface, and points at which this can happen may be found by the method of Lagrange multipliers.

The constraint surface together with level surfaces of the objective function

To appreciate the method of Lagrange multipliers it is helpful to visualize the constraint surface together with surfaces on which the objective function is constant: f(x,y,z) = C. The first picture on the right shows the sphere together with the surface on which f = C with C = 1; the latter surface is a hyperboloid of one sheet. The choice C = 1 is not significant here; it shows a typical situation. The hyperboloid divides the sphere into two pieces. On one piece of the sphere, the function f(x,y,z) has values larger than 1, and on the other piece the values of the function are less than 1. Clearly 1 is not an extreme value.

The second picture shows the sphere together with the surface f(x,y,z) = 2 (again a hyperboloid of one sheet). Now the hyperboloid does not cut the sphere into two pieces; it just touches it tangentially at a single point. The value of f(x,y,z) at any points other than that point of tangency must be less than 2, so that 2 must be the maximum value of f(x,y,z) on the sphere.

The first picture above was drawn with the following Maple commands (similar commands were used for the second picture and other pictures below):

pg:=implicitplot3d(x^2+y^2+z^2=1,x=-2..2,y=-2..2,z=-2..2,axes=none,numpoints=1800,color=blue):
p1:=implicitplot3d(x2+x+3*y*z=.5,x=-2..2,y=-2..2,z=-2..2,axes=none,numpoints=1800,color=green);
pt1:=textplot3d([-1.0,2.5,-2,typeset(f(x,y,z)=1)],font=[Helvetica,16],axes=none):
display3d(pg,p1,pt1,orientation=[50,20,0]);

Here the implicitplot3d commands create the pictures of the two surfaces, the textplot3d command creates the label giving the value of C, and the display3d command combines these into the final picture. The numpoints=1800 option increases the number of points sampled from the default of 1000. With only 1000 points, the graphs drawn for the surfaces look very polyhedral—that is, they resemble more closely linear approximations to the true smooth graphs. The orientation option specifies the point of view; it corresponds to angles ϑ = 50, φ = 20, and ψ = 0.

Lagrange multipliers and extreme values

Consider again the surface f(x,y,z)=x2+x+3y = C and its intersection with the constraint surface, the sphere g(x,y,z)=x2+y2+z2=1. As discussed above, if the surfaces are not tangent at an intersection point (such surfaces are called transverse) then there will be values of f greater and less than C on the constraint surface, but this will not happen (or more precisely, may not happen) if the surfaces are tangent. The Lagrange multiplier equations follow from this observation: we look for extreme values among the solutions x, y, and z of the equations
   ∇f(x,y,z)=λ∇g(x,y,z)
Here λ is a real number, and the equality indicates that the tangent planes are parallel because the normal vectors are pointing in the same direction. This one vector equation is three scalar equations.
   g(x,y,z)=1
This guarantees that the solutions we seek will be on the constraint surface.

In this specific case, we want to solve these equations
   2x+1= λ(2x)
   3z=λ(2y)
   3y=λ(2z)
   1=x2+y2+z2

Finding solutions using Maple

Here is one way to get solutions, assuming that Maple has already been given the functions f and g. The first command instructs Maple to consider only real numbers in solving the Lagrange multiplier equations; the need for this command is discussed further under "A cautionary example" below.

> with(RealDomain):

> LMeqx:=diff(f,x)=L*diff(g,x);

                      LMeqx := 2 x + 1 = 2 L x

> LMeqy:=diff(f,y)=L*diff(g,y);

                      LMeqy := 3 z = 2 L y

> LMeqz:=diff(f,z)=L*diff(g,z);

                      LMeqz := 3 y = 2 L z

> LMsols:=solve({LMeqx,LMeqy,LMeqz,g=1});

      LMsols := { L = 3/2, x = 1, y = 0, z = 0 }, 

                { L = 1/2, x = -1, y = 0, z = 0 },

                { L = -3/2, x = -1/5, y = -2 √3 / 5, z =  2 √3 / 5 },

                { L = -3/2, x = -1/5, y = -2 √3 / 5, z =  2 √3 / 5 },

> LMvalues:=seq(subs(LMsols[j],f),j=1..nops({LMsols}));

       LMvalues := 2, 0, -8/5,  -8/5

Comments on the Maple commands

The commands here are similar to those used above to find the critical point of f, and we comment only on what is new. L is used instead of λ. The Lagrange multiplier equations in this example have four solutions; the solve step finds all the solutions. Finally, the command LMvalues:=... substitutes the values found into the objective function. Inside, the nops({LMsols}) instructs Maple to count the number of elements in the set of solutions -- that is, the number of solutions to the Lagrange multiplier equations. Then Maple will substitute each of these solutions, in order, into the objective function.

The solutions and their geometric interpretations

From the above we have learned that the maximum value of the objective function on the sphere is 2, and this is attained at the point (1,0,0). The minimum value is –8/5, attained at the two points (−1/5,−;2√3/5, 2√3/5) and (−;1/5,2√3/5,−;2√3/5). The value f(x,y,z) = 0 is obtained from one solution of the Lagrange multiplier equations but does not correspond to either a maximum or a minimum of the objective function.

Combining the critical point and Lagrange multiplier results

Our analysis of the critical point gave us one candidate for an extreme value of f(x,y,z) in the domain g(x,y,z) ≤ 1; this was the value
−1/4
at the critical point itself. The Lagrange multiplier equations gave us three candidate values:
2, 0, and −8/5,
By inspection of these four values we conclude that the maximum value is 2 and the minimum is −8/5. If the critical point had not fallen inside the domain g ≤ 1 then we would have had to consider only the values from the Lagrange multiplier analysis.

Further comments

A common error

Students working by hand frequently discard two of the solutions found by Maple. These are the solutions with y=0 and z=0 -- perhaps these are thought to be "trivial" solutions which can't give useful answers. Such solutions should not be thrown away. For the objective function x+3xy these solutions are not extreme values but are more complicated.

A cautionary example

Suppose we keep the unit sphere x2+y2+z2=1 as constraint, and consider the objective function F(x,y,z)=x+0.5yz. What follows is true for F(x,y,z)=x+kyz if 0≤k≤1, but keeping track of the results for one value of k will be complicated enough!

The maximum value of this F on the unit sphere is 1, and this value occurs only at the point (1,0,0). The minimum value is −1, which occurs only at the point (−1,0,0). A picture of the constraint surface and the level surface of the objective function which corresponds to the maximum value, the graph of x+0.5yz=1, is shown to the right.

Algebraically, this is a solution to the Lagrange multiplier equations which has both y=0 and z=0. This sort of solution is exactly one which may be thrown away when people compute "by hand".

Here is a sequence of Maple commands similar to those above done. We only display the results of the last command. This should give the values of the objective function which result. The extreme values should be among these "candidates".

> F:=x+0.5*y*z:

> g:=x^2+y^2+z^2:

> LMeqx:=diff(F,x)=L*diff(g,x):

> LMeqy:=diff(F,y)=L*diff(g,y):

> LMeqz:=diff(F,z)=L*diff(g,z):

> LMsols:=solve({LMeqx,LMeqy,LMeqz,g=1}):

> LMvalues:=seq(subs(LMsols[j],x+.5*y*z),j=1..nops({LMsols}));

      LMvalues := 1., -1., 1.250000000, 1.250000000, -1.250000000, -1.250000000}
The values which Maple returns as candidates for extreme values are ±1 and ±1.25.

To the right is a graph of the unit sphere and the surface x+0.5yz=1.25, at an angle similar to what was just shown before. The surface for the candidate "extreme value" does not seem to touch the unit sphere! What happened?

Let's show a bit more information about by looking in detail at the result of one of the intermediate steps, the solve command.

> LMsols:=solve({LMeqx,LMeqy,LMeqz,g=1});

        LMsols := {z = 0., x = 1., y = 0., L = 0.5000000000},
 
           {z = 0., y = 0., x = -1., L = -0.5000000000},

    {x = 2., L = 0.2500000000, z = 1.224744871 I, y = 1.224744871 I},

   {x = 2., L = 0.2500000000, z = -1.224744871 I, y = -1.224744871 I},

   {x = -2., z = 1.224744871 I, y = -1.224744871 I, L = -0.2500000000},

   {x = -2., y = 1.224744871 I, z = -1.224744871 I, L = -0.2500000000}
The third element of LMsols is
     {x = 2., L = 0.2500000000, z = 1.224744871 I, y = 1.224744871 I}
where I is Maple's notation for i, or the square root of −1. If we compute 2.+.5* 1.224744871 I* 1.224744871 I (that's x+.5yz at this "point") then the result, to 10 digit accuracy, will be 1.250000000.

Here is a quote from a Maple help page. It is an important part of the design of Maple which should be known to every user:

By default, Maple performs computations under the assumption that the underlying number system is the complex field.

So, for example, we have:

> solve(x^2+1);

                         I, -I
The Maple command with(RealDomain), which we introduced above, changes this default setting. In the case of computations with Lagrange multipliers such a change is useful. Consider the following commands and responses.
> with(RealDomain):

> solve(x^2+1);

> LMsols:=solve({eqx,eqy,eqz,g=1});

         LMsols := {z = 0., x = 1., y = 0., L = 0.5000000000},

           {z = 0., y = 0., x = -1., L = -0.5000000000}

> values:=seq(subs(LMsols[j],x+.5*y*z),j=1..nops({LMsols}));

                               values := 1., -1.
When the RealDomain package is loaded, Maple declares that x2+1=0 has no roots. It also discards any elements of LMsols which have non-zero imaginary parts, and so it declares as a result of the substitution step, only those values which correspond to real solutions of the Lagrange multiplier equations. We don't get the outputs ±1.25 because the intermediate steps involved imaginary numbers which RealDomain discards. Those solutions won't be geometrically valid.


Maintained by the Math 251 guru, currently (Spring 2015) E. Speer, speer@math.rutgers.edu. Last modified 3/8/2015.