UC SCI 111, 2 October 2002
Notebook 7
Applications of maxima/minima

Initializations

Note that the BasicRegression package is not included in the standard Mathematica distribution. It is available from the P-drive on the UC network. To install it locally, it's easiest to copy the entire RMCPackages directory to /AddOns/Applications (inside your Mathematica directory).

Needs["Graphics`ImplicitPlot`"];
Needs["RMCPackages`BasicRegression`"]

Off[General::spell1];
Grad[f_, vars_List] := Map[D[f, #] &, vars]

Mathematica Prerequisites

Arithmetic over lists.

A list in Mathematica is simply an ordered collection of items. The syntax for the list containing the letters a, b, c, and d is:

{a, b, c, d}

Most simple arithmetic commands operate over lists. This is handy.

{a, b, c, d}^2

{a, b, c, d} + {e, f, g, h}

When Mathematica is used to analyze data collected in the real world, the data is generally stored as a list. We can sum the items in a list as follows:

Apply[Plus, {a, b, c, d}]

The Gradient

Recall that the gradient of a function of many variables is a vector comprised of the partial derivatives of the function. The command for generating gradients is provided in this notebook:

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

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 presented not as the graphs of functions, but as the set of solutions for an equation. For example, can we plot an ellipse knowing only that it is determined by the equation 4  (x^2) + (y^2) = 4? After loading the Mathematica package Graphics`ImplicitPlot` (which should have occurred automatically at the beginning of this session), the command ImplicitPlot is available for this purpose.  For example, here is a plot of this ellipse.

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 assignments, 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 Hue[.2] will produce a yellow plot, and Hue[.7]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 simplest way to do this is to give names to the individual plots that you wish to superimpose, and then display them together with the Show command. Here is an example using 3-D graphics:

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

Note that the order in which p1 and p2 appear in the Show command matters. If the two plots have different settings for some common plot option, the option's value will be determined by the first plot listed. Unsatisfactory results are sometimes remedied by switching the arguments inside Show.

It is possible to produce only the final (superimposed) graphic, and to not display the individual plots of which it is comprised. To do this is to add the option setting DisplayFunction->Identity to the individual plots, which will prevent 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 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}, PlotStyle -> Hue[1], 
          DisplayFunction -> Identity];
Show[p1, p2, DisplayFunction -> $DisplayFunction]

Exercise 1.
(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 yellow. Experiment to determine an appropriate viewing rectangle. 
(c) Superimpose the plots.

Optimization example: least squares regression

Below you will find a table showing the world record time (in seconds) for running a mile, and the year in which the record was broken. The times are given in such a way that they are fairly easy to read, yet still expressed in seconds. For instance, in 1911 the record time was 4 minutes and 15.4 seconds. Click anywhere in the table and Enter it into Mathematica.

miledata = ({{1911, 60* 4  +  15.4}, {1913, 60 *4  +  14.6}, {1915, 60 *4  +  12.6}, {1923, 60 *4  +  10.4}, {1931, 60 *4  +  9.2}, {1933, 60 *4  +  7.2}, {1934, 60 *4  +  6.8}, {1937, 60 *4  +  6.4}, {1942, 60 *4  +  6.2}, {1942, 60 *4  +  6.2}, {1942, 60 *4  +  4.6}, {1943, 60 *4  +  2.6}, {1944, 60 *4  +  1.6}, {1945, 60 *4  +  1.4}, {1954, 60 *3  +  59.4}, {1954, 60 *3  +  58.0}, {1957, 60 *3  +  57.2}, {1958, 60 *3  +  54.5}, {1962, 60 *3  +  54.4}, {1964, 60 *3  +  54.1}, {1965, 60 *3  +  53.6}, {1966, 60 *3  +  51.3}, {1967, 60 *3  +  51.1}, {1975, 60 *3  +  51.0}, {1975, 60 *3  +  49.4}, {1979, 60 *3  +  49.0}, {1980, 60 *3  +  48.9}, {1981, 60 *3  +  48.8}, {1981, 60 *3  +  48.7}, {1981, 60 *3  +  47.6}, {1985, 60 *3  +  46.5}});

A good picture is worth a thousand table entries, so let's take a look at these data:

s = ScatterPlot[miledata]

As there is a clear linear trend in the data, it is sensible to find the equation of the line that  best fits them. Since a line y = a + b x is determined by just two numbers, the y-intercept a, and the slope b, our task is to find the two values of these variables that optimize the fit of the line to the given data. But how to we measure the goodness of fit? The most common method is the method of least squares. For any line y = a + b x in the plane, we measure the sum of the squares of the vertical distances from the data points to the line. The line for which this sum is least is optimal.

The illustration below shows an animation: just double click on any frame to start it. It shows a smaller set of just three data points, and several candidates for "best fitting" line (in blue). For each line the vertical distances from point to line are shown (dashed black lines), along with their squares (in red). We seek the line for which the sum of the areas of the squares is least. (Delete the graphics when you are done looking at them by selecting their grouping bracket and choosing Cut in the Edit menu).

Do[SquaredResidualPlot[{1, 2, 3}, {1.6, 1.7, 3.1}, t - 1, 
        1 + 2 (t - 1), PlotRange -> {{0, 4}, {0, 4}}, 
        AxesOrigin -> {0, 0}], {t,  .9, 1.1,  .025}]

So for each data point ((Subscript[x, i]), (Subscript[y, i])) we compute the square of the vertical distance from the point to the line: ((((Subscript[y, i]) - (a + b (Subscript[x, i]))))^2). We then sum these quantities, and regard them as a function of a and  b. The graphic above suggests that this function should have a global minimum. The value ((Subscript[a, 0]), (Subscript[b, 0])) that minimizes the function gives the best line: y = (Subscript[a, 0]) + (Subscript[b, 0]) x. Let's carry out this program, keeping in mind that the data points remain constant throughout.

Here are the data given as two lists. We name them x and y.

x = miledata[[All, 1]]

y = miledata[[All, 2]]

Here is a list of the squared vertical distances from the data points to any line y = a + b x:

Clear[a, b, r, f];
r = (((y - (a + b x)))^2)

f = Apply[Plus, r]

The expression for the function f(a, b) can be simplified if all terms are expanded and collected:

f = Expand[f]

Clear[x, y]

Exercise 2.
(a) Find the two partial derivatives of the function f(a, b) given above.
(b) Find the point ((Subscript[a, 0]), (Subscript[b, 0])) that minimizes f.
(c) Superimpose a plot of the line y = (Subscript[a, 0]) + (Subscript[b, 0]) x with the scatterplot s produced earlier. The purpose is to visually confirm the correctness of your work in part b.
(d) Use the formula for the "least squares" line to predict the world record time for running the mile in the year 2000.

Optimization with constraint: Lagrange's multiplier method

Introduction

Now we want to optimize a function under a certain constraint.  Consider for example the function

f(x, 
        y) = (((x - 0.2))^2) + (((y - 0.5))^2) + (x - 0.2) (y - 
              0.5)

Here is contour plot of f together with the ellipse (x^2) + 2  (y^2) = 1:

g1 = 
      ContourPlot[(((x - 0.2))^2) + (((y - 0.5))^2) +  (x - 
                0.2) (y - 0.5), {x, -1.2, 1.2}, {y, -1, 1}, 
        DisplayFunction -> Identity,  Contours -> 20]; 
    g2 = ImplicitPlot[(x^2) + 2 (y^2) == 1, {x, -1, 1}, 
        DisplayFunction -> Identity, PlotStyle -> Hue[ .2]];
\(\(Show[g1, g2, DisplayFunction -> $DisplayFunction, 
	
        AspectRatio\  -> \ Automatic];\)\)

From the picture we can clearly see where the minimum and maximum of f are located inside and on the boundary of the ellipse. The minimum is of course (0.2, 0.5), which is a critical point of f. The maximum is not a critical point of f. Instead it is precisely the maximum of f when we restrict ourselves to the boundary of the ellipse. More precisely, it is at the point where the yellow line and contour lines are tangent. This is equivalent to saying that the gradients of  f and g = (x^(2 )) + 2  (y^2) are parallel. Let us first define f, g and their partial derivatives:

f = (((x - 0.2))^2) + (((y - 0.5))^2) + (x - 0.2) (y - 
                0.5);
g = (x^2) + 2  (y^2);
fx = D[f, x]; fy = D[f, y];
gx = D[g, x]; gy = D[g, y];

The fact that at the extremal point the gradients are parallel means that there exists a number λ such that (Subscript[f, x]) = λ (Subscript[g, x]) and (Subscript[f, y]) = λ (Subscript[g, y]) at that point. Or, said another way, the gradient of f is λ times the gradient of g. Moreover, g = 1 at our extremal point. This gives us three equations for the three unknowns x,  y, λ which we can try to solve.

NSolve[{g == 1, fx == λ gx, fy == λ gy}, {x, 
          y, λ}] // TableForm

Among the solutions we see two real ones . Clearly one of them is the place where f  is maximal. Can you explain why the other point occurs?

The method we have just employed is called Lagrange's multiplier method, which is widely used in optimization. In general, if we must optimize the function f(x, y) under the constraint g (x, y) = 0 we first form the auxiliary function h(x, y, λ)  =  
        f(x, y)  -  λ g(x, y) then solve the simultaneous system (Subscript[h, x]) = 0,  (Subscript[h, y]) = 0,  (Subscript[h, λ]) = 0.  Here is the way to implement these ideas with minimal typing:

h = f - λ*(g - 1);
NSolve[Grad[h, {x, y, λ}] == {0, 0, 0}, {x, y, λ}] // 
      TableForm

An application: the factory

As an application, imagine a new factory that plans to manufacture widgets. The plant manager must decide how to allocate funds to pay for two things: workers and materials. Now for the sake of argument, suppose that if the plant has x workers and y kilograms of materials, then f(x, y) = 300  (x^0.4) (y^0.6) widgets can be made per day. This is a so-called Cobb-Douglas function, common in economic analysis.

Exercise 3.
Make a contour plot of f over an appropriate domain, and explain why the shape of the function behaves as it should for the context of this problem. For instance, you might consider the extreme cases where either x or y is near 0.

Suppose now that the plant manager has a constraint: she can only spend €1000 per day. Now workers cost €60 per day, and materials cost €2 per kilogram.

Exercise 4.
(a) Express the above constraint as an equation in x and y. 
(b) Plot the "curve" determined by this equation.
(c) Show the constraint curve superimposed on the contour plot for f. Producing a good picture may require you to make some adjustments to the domains of the two plots.

Exercise 5.
If the number of widgets produced each day is to be optimized, roughly how many workers, and how many kilograms of materials should be purchased using the available €1000? Use the plot from the last exercise to make your estimate.

Exercise 6.
Obtain an exact solution to the plant manager's dilemma using Lagrange's multiplier method.

Here's another problem that will yield to the Lagrange multiplier method:

Exercise 7.
Find the point on the hyperbola x y = 1 closest to the point (2, 3).
Hint: distance is minimized precisely when the square of the distance is minimized. Minimizing the square of the distance will simplify some calculations.

An interesting example: Arithmetic and geometric means.

Suppose we are given two positive numbers a, b. The arithmetic mean of these numbers is defined by ((a + b)/2). The geometric mean of a and b is defined by (Sqrt[ab]). It is a well-known fact that the geometric mean is always smaller than the arithmetic mean. The reason is that ((( (Sqrt[a])  -  (Sqrt[b])))^2) ≥ 0. Expanding this we find a-2(Sqrt[ab])+b≥0. Hence a + b ≥ 2 (Sqrt[ a b]), from which our assertion follows.

We can also derive this result using Lagrange multipliers. Let x = (a/(a + b)) and y = (b/(a + b)). The inequality (Sqrt[a b]) ≤ ((a + b)/2), after division on both sides by a + b, comes down to (Sqrt[x y]) ≤ (1/2), or better yet x y ≤ (1/4), under the restriction x + y = 1.

So letting

g = x + y

we maximize

f = x*y

under the restriction g = 1. We introduce the Lagrange multiplier and solve the system by

h=f-λ (g - 1);
Solve[
      Grad[h, {x, y, λ}] == {0, 0, 0}, {x, y, λ}]

So the maximum of x y is attained when x = y = (1/2). We conclude that x y ≤ (1/4) for all positive x,  y with x + y = 1. We have shown the geometric mean of two numbers a and b cannot exceed their arthemetic mean.

Now we play the same game with three numbers a, b, c. We would like to show that ((a b c)^(1/3)) ≤ ((a + b + c)/3). On the left we have the geometric mean of three numbers, on the right you see the arithmetic mean of these numbers. In a similar way as above we can reduce this to the optimization problem for x y z under the restriction x + y + z = 1.

Exercise 8.
(a) Solve this problem using Lagrange's multiplier method.
(b) Determine the maximum of (Subscript[x, 1]) (Subscript[x, 2])... (Subscript[x, n ])under the restriction (Subscript[x, 1]) + (Subscript[x, 2]) +  ...  + 
          (Subscript[x, n])=1  using Lagrange's method and pencil and paper.