UC SCI 111, vM+. 
Notebook 13

Needs["Graphics`Colors`"];
Off[General::spell1];
<<  Graphics`ImplicitPlot`

Mathematica: plotting curves revisited.

ImplicitPlot.

Sometimes we want to make a picture of a curve that's not the graph of a function. We have seen before that one way of doing this uses Mathematica's command ParametricPlot. But this only helps when you have a parametrization of the curve. Often you know an equation defining the curve but not a parametrization. In this situation the command ImplicitPlot can help. Let's demonstrate it by the all-famous example.

ImplicitPlot[x^2 + y^2 == 1, {x, -1, 1}];

Exercise 4.
A parametrization of this circle C is given by 

	x = \(TraditionalForm\`\(2\ t\)\/\(\(\ \)\(\(t\^\(\(2\)\(\ \)\)\)\(+\)\(\ \)\(1\)\(\ \)\)\)\)		y =\(TraditionalForm\`\(\(1\)\(-\)\(\ \)\(t\^2\)\(\ \)\)\/\(\(\ \)\(\(t\^\(\(2\)\(\ \)\)\)\(+\)\(\ \)\(1\)\(\ \)\)\)\)	

a) Plot this parametrization using ParametricPlot.

Fix the point  P=(0,1)  on  C. Consider (Subscript[l, t]), the line that passes through P with slope -t. In general, (Subscript[l, t]) will pass C  in one other point, (Subscript[Q, t]).

b) Prove that the parametrization above is indeed a parametrization of C by computing the coordinates of (Subscript[Q, t]).

This method is called linear projection from a given point.

Optimization with and without constraints.

Two variable functions.

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

f = (x + 1) (y - 1) Exp[-x^2 - y^2]

Let us now plot its three dimensional graph. Here we have given the interesting range of x- and y-coordinates. Of course, when studying another function, you have to decide which range of values will give the interesting features of the graph.

Plot3D[f, {x, -3, 3}, {y, -3, 3}, PlotPoints -> 30, 
      AspectRatio -> 1];

So now we seem to have an idea of where the (local) maxima and minima are located. We can determine their positions by determining the common zeros of the partial derivatives of our function. First we determine the derivatives,

fx := D[f, x]

fy := D[f, y]

Now we solve the local maxima/minima problem.

extrema = NSolve[{fx == 0, fy == 0}, {x, y}]

Sometimes Mathematica tells us that it may not have found all solutions when using NSolve. So to make sure that we have indeed found all local maxima/minima it is better to be sure and consider the graphs carefully. Looking at a so-called contourplot may also help,

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

Can you decide where the local maxima/minima should all be? To determine the function values we give the following command

f  /. extrema

There is one local minimum which is not so clear in the contour plot. To have a better look we make a contourplot near this point,

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

Can you explain why this plot gives a much better picture of the local minimum?

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

g := x y Sin[x^2 + y^2]/(1 + x^2 y^2)

And plot its three dimensional graph

Plot3D[g, {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, -3, 3}, {y, -3, 3}, Contours -> 25, 
      ContourLines -> False, PlotPoints -> 40];

You can probably understand that NSolve is not able to determine all local extrema for us and we have to resort to some specialized methods to determine individual extrema up to some precision. We shall not do this here however.

To grit your teeth you can now determine (local) extrema for the following functions

	f(x,y)  =  x y (1-x-y)
	
	f(x,y) = (((x^4)+(y^4)-10 (x^2)+20 y)/(10 + 
            (x^2) + (y^2)))
	
Use the plot commands given above for f to get a picture of its graph. Also solve the extrema problem with NSolve.  It is not always easy  to get the right picture. So you have to use a lot of common sense as well.

Exercise 5.
For each of the two functions above prepare a sketch (on paper) of the level lines of the function and the location of the (local) extrema and saddle points. Also give the coordinates of these points.

Least Squares Method.

In many laboratory experiments the set of measurements provides us with a finite set of data which can be plotted in a two-dimensional graph. Here is an example.

data = {{1, 3.2}, {2.2, 4}, {2.9, 4.9}, {3.6, 5.5}, {4.3, 6}, {5.1, 
            6.4}, {6, 8}, {6.5, 8.8}, {7.2, 9}, {8, 9.5}};

dataplot = 
        ListPlot[data, PlotStyle -> AbsolutePointSize[5]];

We are given that these data correspond to a linear function y=ax+b. The reason that they do not lie on a straight line is the error which is inevitable in measurements. What we like to do is find a function y=ax+b which fits as closely to the data as possible. Denote the x-coordinates of the data by (Subscript[x, 1]), (Subscript[x, 2]), (Subscript[x, 3]), .... and the y-coordinates by (Subscript[y, 1]), (Subscript[y, 2]), (Subscript[y, 3]), ....We choose a, b in such a way that the sum of all   \(TraditionalForm\`\(\((y\_i\  - a\ x\_\(\(i\)\(\ \)\) - b)\)\(\ \)\)\^2\) is minimal. This is called the method of least squares. Let us work this out in our example.

s = 
      Expand[Sum[(((data[[i, 2]] - a data[[i, 1]] - b))^2), {i, 
            10}]]

Now we determine the stationary point

solutions = NSolve[{D[s, a] == 0, D[s, b] == 0}, {a, b}]

So there is one stationary point. Since we know that there exists a global minimum the point found must be this minimum. Thus we have found our optimal linear function.

f = a x + b  /. solutions

Here is a combined plot of the measurements and the linear fit.

Show[dataplot, 
        Plot[f, {x, 0, 8}, DisplayFunction  ->  Identity ]];

It is interesting to note that we found exactly one stationary point. However, this becomes evident if we looks at the partial derivatives,

D[s, a]

D[s, b]

They are linear functions in a, b. To find the stationary points we must solve 2 linear equations in 2 variables. Usually this has a unique solution.

Exercise 6.
a) In the following command a list of 100 data is generated. Using the same steps as above, determine a linear fit for these data using the method of least squares.

data = Table[{Random[Real, {0, 10}], 0}, {k, 100}]; Clear[k];
For[k = 1, k <= 100, k++, 
      data[[k, 2]] = 
        2*data[[k, 1]] + 3 + 
          Tan[Random[Real, {-1.4, 1.4}]]]

b) It would be interesting to have a formula for a linear fit of a set of data ((Subscript[x, 1]), (Subscript[y, 1])), ((Subscript[x, 2]), (Subscript[y, 2])),((Subscript[x, 3]), (Subscript[y, 3])),.........  In the following command we define a set of four arbitrary data. Determine a linear fit by using the method of least squares. See if you can guess a formula for the linear fit of n  arbitrary data.

data = {{x1, y1}, {x2, y2}, {x3, y3}, {x4, y4}};

In Mathematica the very general function Fit performs linear fits by using the method of least squares in many different situations.

c) Read the help function for Fit and use it to perform a linear least squares approximation on the first set of data.

The distance between two line segments.

Let p and q be two points in (ℝ^3). The straight line through p and q is, by definition, the set

	{ (1-t)p + t q: t ϵ ℝ }	=	{ t(q-p) + p : t ϵ ℝ }
	
By restricting the parameter t tho the interval [0,1] we obtain the line segment pq between p and q. Note that the direction of the line is determined by the vector q-p. N.B: recall from Exercise 8.2.2 that directions correspond to unit vectors, so the direction is given by

	((q - p)/(|| q - p ||))
	
As a warming up we ask the following question: which point on this line is closest to the origin? This geometric question can be reformulated analytically as: for what value of t is || t(q-p) + p || minimal? As will be evident soon, it is slightly easier to avoid the square root and look for a minimum of

	(1/2) || t ( q-p ) +  p \(TraditionalForm\`\( || \^\(\(\ \)\(2\)\)\)\)	=	(1/2)( t ( q-p ) + p )·( t ( q-p ) + p )
	
Indeed, the derivative of this function of t is given by

	(q-p) · (t (q-p) + p )
	
Verify this! Hint: recall the general properties of inner products from Kreyszig, section 8.2.
By putting this expression equal to zero we find the answer

	t = - \(TraditionalForm\`\(\((\ q - \ p\ )\)\ ·\ p\)\/\(\(||\)\(\ \)\(q - \ p\)\(\ \)\( || \^\(\(\ \)\(2\)\)\)\)\)
	
But in addition to this explicit expression, one can derive a geometrical property from the condition

	(q-p) · ( t (q-p) + p ) = 0

by interpreting q-p as, modulo normalization, the direction of the line an t (q-p ) + p as, again modulo normalization, the direction of the line that connects the origin to the point on the original line: the fact that the inner product of these two directions is zero means that the two lines are orthogonal, that is perpendicular to each other, see Kreyszig, section 8.2. In words: the point closest to the origin can be found by gliding a point along the line and looking when the connection between this point and the origin is perpendicular to the line.

Note that we take so far for granted that the extremum is a minimum, as this seems evident from the geometry; in fact you can easily verify it analytically an you are invited to do so.

Now let v and w be two other points in (ℝ^3). The problem posed in the following exercise is to determine the distance between the line segments pq and vw, which is by definition the minimal distance between two points, one on each segment.

Exercise 7.
Use  tϵ  [0,1] to parameterize pq and u ϵ  [0,1] to parameterize vw. Let

	f(t,u) = (1/2)(distance)^2
	
where the distance is between the point with parameter t on pq and the point with parameter u on vw.
a) Write down f(t,u) explicitly.

To consider concrete examples we must choose points p, q, v and w. We exploit the random number generator to do this.

r  :=  Random[Integer, {0, 100}]
{p, q, v, w}  =  
        Table[{r, r, r}, {4}];
\(\(Show[
Graphics3D[{
PointSize[ .02], 
Point /@ {p, q, v, 
                w}, 
Red, 
Line[{p, 
                q}], 
Blue, 
Line[{v, 
                w}]
}], 
Boxed\  -> \ 
          False
];\)
	\)

b) Define the distance function f(t,u) in Mathematica and make a 3-dimensional plot of its graph for (t,u) ϵ [0,1]x[0,1]. Repeat this for several choices of p,q,v and w. Also make contour plots and mark where extrema occur. Investigate  whether an extremum is a maximum or a minimum by considering the graph. Do you ever find a maximum in the interior of [0,1] × [0,1]?
c) There are three essentially different configurations for the two points that yield the minimal distance. Explain what these three types are. To do so, you need to make explicit what is meant by essentially.

Finally, let us remove the restriction that the points lie between p and q, respectively, between v and w. In other words, let us remove the restriction

	0 ⩽ t ⩽ 1,	0 ⩽ u ⩽ 1,
	
when considering f(t,u).

d)  Show that f(t,u) has a unique critical point whenever the two lines are not parallel. Also show that the critical point is a minimum. 
Hint: the Cauchy-Schwarz inequality

	| a·b | ⩽ ||a||·||b||
	
is strict unless a is a multiple of b.
e) Show that the line through the two points that yield the shortest distance is perpendicular to both the line through p and q and the line through v and w. In which of the three configurations of part c) is this also true?

Lagrange's multiplier method.

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)
	
in the region (x^2)+2(y^2)≤1.  Here is a contourplot of f together with the graph of 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 -> RGBColor[1, 1, 0]];
\(\(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 stationary point of f. The maximum is not a stationary 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 contourlines 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 fx=λ*gx, fy=λ*gy at that point. 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, λ}]

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 Langrange'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,  g = 0.

Exercise 8.
Determine the maximum and minimum of the function f(x, y) = (x^2) - 2 (y^2)when x,y are restricted to the circle (((x - 1))^2) + (((y - 1))^2) = 2.

Several constraints.

Using Lagrange's method we can also solve problems with more than two variables and several constraints. As an example consider the ellipse 3 (((x - 3))^2) + (((y - 5))^2) = 1 and the ellipse 2 (((x - 1))^2) + (((y - 1))^2) =1.    Here is a picture ,

\(\(ImplicitPlot[{3 \((x - 3)\)^2 + \((y - 5)\)^2 == 
            1, 
		2 \((x - 1)\)^2 + \((y - 1)\)^2 == 1}, {x, 0, 4}, {y, 0, 
          6}, 
	AspectRatio\  -> \ 1, \ PlotPoints\  -> \ 50];\)\)

We like to determine a point ((Subscript[x, 1]), (Subscript[y, 1])) on the first ellipse and ((Subscript[x, 2]), (Subscript[y, 2])) on the second ellipse such that their distance is minimal. We now have an optimization problem with four variables (Subscript[x, 1]), (Subscript[y, 1]), (Subscript[x, 2]), (Subscript[y, 2 ])and two constraints 
	
	3 ((((Subscript[x, 1]) - 3))^2) + ((((Subscript[y, 1]) - 5))^2) = 1 	    2 ((((Subscript[x, 2]) - 1))^2) + ((((Subscript[y, 2]) - 1))^2) = 1.  
	
The function to be minimized is

d = (((x1 - x2))^2) + (((y1 - y2))^2)

We introduce the constraint functions

f = 3 (((x1 - 3))^2) + (((y1 - 5))^2); 
    g = 2 (((x2 - 1))^2) + (((y2 - 1))^2);

To each constraint we associate a Lagrange multiplier, for example λ and μ. Define

h = d - λ f - μ g;

We have to deal with four partial derivatives since there are four variables,

hx1 = D[h, x1]; hx2 = D[h, x2]; hy1 = D[h, y1]; 
    hy2 = D[h, y2];

Since the problem is too complicated for NSolve we use FindRoot with the approximate solution (2.5,4) (1.5,1.5) which can be inferred from the picture. We simply guessed the value for the multipliers

FindRoot[{hx1 == 0, hx2 == 0, hy1 == 0, hy2 == 0, f == 1, g == 1}, {x1, 
        2.5}, {y1, 4}, {x2, 2}, {y2, 2}, {λ, 1}, {μ, 
        1}]

Exercise 9.
Consider the hyperbolas xy=1 and the ellipse (((x - 4))^2) + 2 (((y - 4))^2) = 1. Show that they have no point of intersection. Determine the point P on the hyperbola and Q on the ellipse such that the distance between P and Q is minimal.

Arithmetic and geometric mean.

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 geometric 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[ab]) from which our assertion follows.

We can also show this result using Lagrange multipliers. Let x=a/(a+b) and y=b/(a+b). The inequality (Sqrt[ab])  ≤(a + b)/2 after division on both sides by a+b ,comes down to (Sqrt[xy]) ≤1/2 or better yet, xy ≤ (((1/2))^2), 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 by

h=f-λ (g - 1);hx = D[h, x];hy = D[h, y];

Now solve the Lagrange multiplier system,

Solve[{hx == 0, hy == 0, g == 1}, {x, y, λ}]

So the maximum of xy is attained when x=y=1/2. We conclude that xy ≤ (((1/2))^(2 )) for all positive x,y with x+y=1.

Now we play the same game with three numbers a,b,c. We like to show that (abc^(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 three numbers. In a similar way as above we can reduce this to the optimization problem for xyz under the restriction x+y+z=1.

Exercise 10.
a) Solve this optimization problem by using Lagrange multipliers.
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.