UC SCI 111, 2 Oct 2002
Notebook 6
Optimisation in 2 variables

Off[General::spell1];
Grad[f_, vars_List] := Map[D[f, #] &, vars]
Needs["Graphics`ImplicitPlot`"]

Mathematica Prerequisites

Plotting functions of two variables with Plot3D.

Suppose we have a function f of two variables (enter the cells below):

Clear[f, x, y]

f[x_, y_] := x y (E^(-((x^2) + (y^2))))

We plot f with the function Plot3D, as follows:

Plot3D[f[x, y], {x, -2, 2}, {y, -2, 2}]

By default the front edges of the box running left to right are parallel to the x-axis, the side edges running front to back are parallel to the y-axis, and the vertical edges are parallel to the z-axis. By default Mathematica takes the domain and divides it into a 15×15 grid, and produces the plot by evaluating the function at each of the (15^2) = 225 points on the grid. You will often want greater resolution; this can be achieved with the option PlotPoints. The setting PlotPoints->40, for instance, will cause (40^2) = 1600 points to be sampled. It takes  longer to produce such a graph, but the results are often worth it.

Plot3D[f[x, y], {x, -2, 2}, {y, -2, 2}, 
      PlotPoints -> 40]

Memory management.

3-D graphics consume an enormous amount of your computer's memory, and can cause things to slow down (or worse) if you produce too many in one session. You can practice good memory management by periodically deleting graphics that are no longer needed. Simply go to the 'Kernel' menu, and choose 'Delete All Output' if you want to get rid of everything. You can be more selective by clicking on the cell bracket of a graphic and deleting it when you are finished looking at it.

Plotting functions of two variables with ContourPlot.

Use the command ContourPlot (same syntax as Plot3D) to view the level curves of a function of two variables. A level curve at height z is the collection of points (x, y) in the plane for which z = f(x, y). The ContourPlot command will by default show 10 level curves, and will shade the regions between them according to their heights. Higher regions are lighter, while lower regions are darker.

ContourPlot[f[x, y], {x, -2, 2}, {y, -2, 2}]

The option Contours controls the number of level curves shown. Show more or fewer curves by setting Contours to any number you wish.

ContourPlot[f[x, y], {x, -2, 2}, {y, -2, 2}, Contours -> 20, 
      PlotPoints -> 35]

The option Contours may also be set to a list of z-values. In this case, level curves are shown at the heights you specify. This is used extensively in the next section.

ContourPlot[f[x, y], {x, -2, 2}, {y, -2, 2}, 
      Contours -> {- .05, 0,  .05}, PlotPoints -> 35]

Exercise 1.
Describe in words the level "curve" at height z = 0 for the function f(x, y) = 
        x y (e^(-((x^2) + (y^2)))). Reconcile your answer with the plot above and with the algebraic definition of f.

Plotting curves in the plane with ImplicitPlot.

We have used the Plot command to plot the graphs of functions of one variable. But how can one plot curves that are not the graphs of functions, for example a closed curve such as an ellipse? After loading the Mathematica package Graphics`ImplicitPlot` (which should have occured automatically at the beginning of this session), the command ImplicitPlot is available for this purpose.  For example, here is a plot of the ellipse with equation 4  (x^2) + (y^2) = 4.

ImplicitPlot[4  (x^2) + (y^2) == 4, {x, -2, 2}]

This command takes as its first argument an equation in two variables. Note that equations use the double equal sign (a special syntax is needed since the single equal sign is used for making asignments, a fundamentally different operation).

A handy option for both Plot and ImplicitPlot is PlotStyle, which controls the appearance of the curve. The setting PlotStyle->Hue[1] will produce a red plot, while PlotStyle->Hue[.7] will produce a blue plot.

Superimposing two plots with Show.

It is useful to be able to take two graphics (for instance two plots produced by different plotting commands) and superimpose one upon the other. This is accomplished with the Show command.  The slickest way to do this is to add the option setting DisplayFunction->Identity to the individual plots, which prevents them from being rendered individually. Then add the option setting DisplayFunction->$DisplayFunction to the Show command to render the results. Here is an example using 3-D graphics:

p1 = 
        Plot3D[4  (x^2) + (y^2), {x, -2, 2}, {y, -2, 2}, 
          PlotPoints -> 30, 
          DisplayFunction -> Identity];
p2 = 
        Plot3D[4, {x, -2, 2}, {y, -2, 2}, 
          DisplayFunction -> Identity];
Show[p1, p2, DisplayFunction -> $DisplayFunction]

And here's an example using 2-D graphics:

p1 = 
        ContourPlot[4  (x^2) + (y^2), {x, -2, 2}, {y, -6, 6}, 
          PlotPoints -> 30, 
          DisplayFunction -> Identity];
p2 = 
        ImplicitPlot[x y == 1, {x, -2, 2}, {y, -6, 6}, 
          PlotStyle -> Hue[0.3], 
          DisplayFunction -> Identity];
Show[p1, p2, DisplayFunction -> $DisplayFunction]

Exercise 2.
(a) Make a contour plot of the function f(x, y) = sin(x y)
(b) Make a plot of the curve determined by the equation (x^3) + 3  x y + (y^3) = -2. Render the plot in blue. Experiment to determine an appropriate viewing rectangle. 
(c) Superimpose the plots.

Taking derivatives of functions of many variables with D.

A function of two variables can be differentiated with respect to either variable. It is said to have two partial derivatives. One simply differentiates (using the usual rules) with respect to one variable, treating the second variable as a constant. In Mathematica, the usual differentiation commands still work. But now it is clear why the syntax demands that we specify the variable with respect to which the differentiation occurs.

D[Sin[x y],  x]

D[Sin[x y],  y]

We may also use the BasicInput Palette for such calculations:

(∂_x) Sin[x y]

This is a common notation for partial derivatives. Throughout this notebook, the partial derivative of a function f(x, y) with respect to x will be denoted (∂_x) f(x, y).

The partial derivatives  measure the slope of the function in the x and y directions, respectively. At the top of a peak or the bottom of a valley, both slopes are zero (the plane tangent to the graph of the surface is horizontal). This leads to a guiding principle for optimization: Maximum and minimum values of a differentiable function of two variables can only occur at a point (x, y) where both partial derivatives evaluate to zero.

A convenient notation exists for expressing both partial derivatives simultaneously. We define the gradient of a function f(x, y) to be the ordered two-tuple whose components are the two partial derivatives. This notebook provides the command Grad, which will give the gradient.

Grad[Sin[x y], {x, y}]

Optimization without constraint

Example 1.

In this section we consider functions of two variables and determine their local maxima and minima. We start with a sample function,

Clear[f, x, y];
f[x_, y_] := x y (E^(-((x^2) + (y^2))))

We wish to find the (x, y) coordinates of those points where f is optimized. A plot is a good place to start. For our purposes, a contour plot is most useful.

ContourPlot[f[x, y], {x, -2, 2}, {y, -2, 2}, 
      PlotPoints -> 30]

The plot shows four extrema consisting of two maxima and two minima. A graphical approach is simple: zoom in on the area of the plot near an extremum. Here we choose the maximum in the upper right of the above graphic.

ContourPlot[f[x, y], {x,  .6, 1}, {y,  .6, 1}, 
      PlotPoints -> 30]

We can repeat this procedure as often as necessary to attain the accuracy we desire.

ContourPlot[f[x, y], {x,  .7,  .71}, {y,  .7,  .71}, 
      PlotPoints -> 30]

We conclude that f has a local maximum near the point ( .707,   .707).

There is also an analytical approach. At an extremum, both partial derivatives evaluate to zero. That is, the 3-D graph has horizontal tangents in both the x direction (west to east) and the y direction (south to north). We should be able to find all potential extrema at once by solving the system of equations

(∂_x) f(x, y) = 0
(∂_y) f(x, y) = 0

for the variables x and y. Since we know how to ask Mathematica to solve systems of equations, we embark upon this task.

fx = D[ f[x, y], x]

fy = D[ f[x, y], y]

We first try the Solve command, as this gives exact answers. Is it up to the task?

extrema = Solve[{fx == 0, fy == 0}];
extrema // TableForm

The messages are simply warnings, or disclaimers. Mathematica is unwilling to guarantee the results; we should verify them some other way. Note for example that the solution x = 0, y = 0 does not occur in the list.
Since only the first four rows contain real valued solutions, we restrict our attention to these. The fourth row contains the only solution where both x and y are positive and real. Does this coincide with our graphical solution?

N[(1/(Sqrt[2]))]

Excellent!

We now see that the first four rows coincide precisely with the four extrema we identified in the contour plot. Using the analytical approach we had the two-pronged advantage of finding all extrema in one go, and in attaining perfect accuracy. We can also find the values assumed by f at these points:

f[x, y] /. extrema[[1]]

f[x, y] /. extrema[[2]]

f[x, y] /. extrema[[3]]

f[x, y] /. extrema[[4]]

One last comment is in order. The very last solution to the system was simply y = 0. Since no x value was mentioned, x can assume any value whatever. However, a quick look at the system reveals that y = 0 is not a general solution (check (∂_y) f(x,y)). Mathematica is not infallable, and did warn us that there may be problems. One should attempt to verify all potential solutions graphically or some other way before trusting the results.

Example 2.

To give an example showing that things can become really complicated take the following function:

Clear[g, x, y];
g[x_, 
        y_] := ((x y Sin[(x^2) + (y^2)])/(1 + (x^2) + (y^2)))

Plot3D[g[x, y], {x, -3, 3}, {y, -3, 3}, 
        PlotPoints -> 50];

If you cannot see the local extrema straight away, you can also make a contourplot. In the following plot we turned off drawing the actual contourlines, which in this case gives a nice visual effect.

ContourPlot[g[x, y], {x, -3, 3}, {y, -3, 3}, Contours -> 25, 
        ContourLines -> False, PlotPoints -> 40];

We can visually attempt to find the coordinates of the maximum in the first quadrant that is closest to the origin by zooming.

ContourPlot[g[x, y], {x,  .8, 1.2}, {y,  .8, 1.2}, 
        PlotPoints -> 40];

We note that there is a local maximum near the point ( .94,  .94).

We now attempt an analytical approach.

gx = D[ g[x, y], x]

gy = D[ g[x, y], y]

Aren't you glad you didn't have to do that by hand? But the complexity of these expressions leaves us doubtful that the system of equations will easily admit solution.

Solve[{gx == 0, gy == 0}]

Not good. We must abandon hope of finding precise solutions. We now use the next weapon in our equation solving arsenal.

NSolve[{gx == 0, gy == 0}]

Interesting, but only the first of these solutions tells us a value for x and for y: the point (0, 0). A careful look at the plot of g shows that this point is not an extremum, but a so-called "saddle point" or "mountain pass". It shows that not all solutions of the system of equations need be extrema.

We can tell the command to solve for only x and y by adding a second argument:

NSolve[{gx == 0, gy == 0}, {x, y}]

Drat! Our worst fears are confirmed -these are truely nasty equations. We are now forced to use our final weapon, FindRoot, and seek solutions one by one.

max = FindRoot[{gx == 0, gy == 0}, {x,  .94}, {y,  .94}]

Here is the value of the relative maximum we have just found:

g[x, y] /. max

If we want other extrema, we have to go after them one by one by choosing appropriate initial guesses to put into the FindRoot command above.

Further examples.

Here are two more functions to investigate:

f(x, y) = x y (1 - x - y)

g(x, y)=(((x^4)+(y^4)-10 (x^2)+20 y)/(10 + (x^2) + (y^2)))

Exercise 3.
For each of the two functions above prepare a plot the level curves of the function and the location of the (local) extrema and saddle points. You will need to spend some time finding an appropriate domain. Also give the coordinates of these points.

First order Taylor approximation

First order Taylor approximation: The tangent plane.

Recall that for a function f(x) of one variable, the first order Taylor polynomial at the point x = (Subscript[x, 0]) is  the line

y = f((Subscript[x, 0])) + f' ((Subscript[x, 0])) (x - (Subscript[x, 0]))

This is simply the equation of the tangent line to the graph of f at the point x = (Subscript[x, 0]).

Similarly, for a function f(x, y) of two variables, there is a first order Taylor polynomial at the point ((Subscript[x, 0]), (Subscript[y, 0])) with equation

z = 
      f((Subscript[x, 0]), (Subscript[y, 0])) + (∂_x) f((Subscript[x, 0]), (Subscript[y, 0])) (x - 
              (Subscript[x, 0])) + (∂_y) f((Subscript[x, 0]), (Subscript[y, 0])) (y - 
              (Subscript[y, 0]))

This is simply the equation of the tangent plane to the graph of f at the point ((Subscript[x, 0]), (Subscript[y, 0])).

Consider the following example.

Clear[f, x, y];
f[x_, y_] := (E^(-((x^2) + (y^2))))

p1 = 
      Plot3D[f[x, y], {x, -2, 2}, {y, -2, 2}, 
        PlotPoints -> 30]

We will compute the expression for the tangent plane to f at the point ((1/10), -(2/10)).

fx = D[ f[x, y], x] /. {x -> 1/10, y -> -2/10}

fy = D[ f[x, y], y] /. {x -> 1/10, y -> -2/10}

t[x_, y_] = 
      f[1/10, -2/10] + fx (x - 1/10) + fy (y + 2/10)

We can show the tangent plane together with the function as follows:

p2 = 
        Plot3D[t[x, y], {x, - .3,  .5}, {y, - .4, 0}, 
          DisplayFunction -> Identity];
Show[{p1, p2}]

Exercise 4
Calculate the first order Taylor polynomial for the function f(x, y) = 1 - (((x - 1))^2) - (((y - 1))^2) at the point (0, 0).

drop: Vectors and the dot product.

A brief review of vectors is in order. A vector in the plane can be thought of as a magnitude and a direction. It is typically represented geometrically as an arrow in the plane; the length of the arrow is its magnitude, and the direction in which it points is its direction. It may be placed anywhere in the plane, but it is common to place the "tail" of all vectors at the origin. A common notation for the vector whose "head" is at the point ((Subscript[u, 1]), (Subscript[u, 2])) when its tail is at the origin, is simply {(Subscript[u, 1]), (Subscript[u, 2])} (using curly brackets). The special vector {0, 0} has no direction and a magnitude of 0. It is called the zero vector.

Needs["Graphics`Arrow`"];
u = 
        Graphics[{Hue[1], 
            Arrow[{0, 0}, {2, 1}]}];
o = 
        Graphics[{PointSize[ .02], 
            Point[{0, 0}]}];
\(\(t1 = 
        Graphics[
          Text[\*""(\!\(u\_1\),\!\(u\_2\))"", {2.3, 1}, 
            TextStyle -> {FontSize -> 
                  12}]];\)\)
t2 = 
        Graphics[
          Text["(0,0)", {0,  .2}, 
            TextStyle -> {FontSize -> 
                  12}]];
Show[u, o, t1, t2, AspectRatio -> Automatic, 
        PlotRange -> All];

Exercise 5
A vector has its tail at the point ((Subscript[x, 0]), (Subscript[y, 0])) and its head at the point (x, y). In the standard "curly bracket" notation, find an expression for this vector. Hint: Where will the head of the vector be when its tail is at the origin?

The dot product of the vector {(Subscript[u, 1]), (Subscript[u, 2])} with the vector {(Subscript[v, 1]), (Subscript[v, 2])} is defined to be the number

(Subscript[u, 1]) (Subscript[v, 1]) + (Subscript[u, 2]) (Subscript[v, 2])

Mathematica can calculate dot products; just use the period as the operator:

{1, 2} . {3, 4}

Exercise 6
Using only paper and pencil, show that two nonzero vectors are perpendicular if and only if their dot product is zero. Hint: Consider the picture below with two vectors rooted at the origin, and compute the lengths of the three sides of the resulting triangle. Now think like Pythagoras...

Needs["Graphics`Arrow`"];
u = 
        Graphics[{Hue[1], 
            Arrow[{0, 0}, {1, 2}]}];
v = 
        Graphics[{Hue[1], 
            Arrow[{0, 0}, {4, -2}]}];
o = 
        Graphics[{PointSize[ .02], 
            Point[{0, 0}]}];
\(\(t1 = 
        Graphics[
          Text[\*""(\!\(u\_1\),\!\(u\_2\))"", {1.5, 2}, 
            TextStyle -> {FontSize -> 
                  12}]];\)\)
\(\(t2 = 
        Graphics[
          Text[\*""(\!\(v\_1\),\!\(v\_2\))"", {4, \(-1.6\)}, 
            TextStyle -> {FontSize -> 
                  12}]];\)\)
t3 = 
        Graphics[
          Text["(0,0)", { .5,  .1}, 
            TextStyle -> {FontSize -> 
                  12}]];
Show[u, v, o, t1, t2, t3, AspectRatio -> Automatic, 
        PlotRange -> All];

drop: Geometrical interpretation of the gradient

Consider a differentiable function f(x, y) of two variables, and a point ((Subscript[x, 0]), (Subscript[y, 0])) in its domain. Imagine the contour plot of this function, with a level curve shown passing through the point ((Subscript[x, 0]), (Subscript[y, 0])). One can find the equation of the line in the plane tangent to this level curve at the point ((Subscript[x, 0]), (Subscript[y, 0])) quite easily: The level curve is simply the intersection of the 3-D plot of f with the horizontal plane z = f((Subscript[x, 0]), (Subscript[y, 0])). The intersection of the tangent plane to f at ((Subscript[x, 0]), (Subscript[y, 0])) with this horizontal plane gives the desired tangent line.

Exercise 7
Using pencil and paper, show that the tangent line to the level curve at ((Subscript[x, 0]), (Subscript[y, 0])) is given by (∂_x) f((Subscript[x, 0]), (Subscript[y, 0])) (x - 
                (Subscript[x, 0])) + (∂_y) f((Subscript[x, 0]), (Subscript[y, 0])) (y - (Subscript[y, 0])) = 
        0.

Exercise 8
Use the results of exercises 5, 6, and 7 to show that the gradient vector is perpendicular to the tangent line to the level curve of f at ((Subscript[x, 0]), (Subscript[y, 0]))