Math 121 Lab 5:
Finding extrema of a function

© Copyright by T.I. Lakoba
The loop used in Exercise 2 was improved following a suggestion by Elijah Jackson (Fall 2022). 
Last modified  2022-11-16

In this Lab, you will use and compare two different methods of finding  coordinates of critical points of a function of two variables. In Exercise 1 you will use a method covered in the textbook and in the lectures. In Exercise 2 you will consider an alternative method, which is a prototype of a large group of practical methods of finding extrema of functions of several variables.

Clear["`*"]

Exercise 1

Here you will first estimate the location of critical points of a surface by inspection of its contour plot. Then you will attempt to find those points by using the general method described in Section 14.7 of the textbook.

(a)   Define a function given below. 

Technical notes:

1.  Make sure to use the correct syntax of defining a function of two arguments. As a reminder, see and follow the Important technical notes 2 and 3 in Exercise 2(a) of Lab 4. 
 
2.  When defining the function, copy and paste the expression on the right-hand side. (For some reason, Mathematica then appends apostrophes after 0.3 and 0.15; simply delete them.) It is very important that you use the correct expression for the function here. This is why it is okay to copy and paste this expression as opposed to entering it yourself - in this particular case. If you miss a parenthesis or enter an incorrect coefficient somewhere, your function may look completely different than the one intended here, and this will derail the rest of your work in this Exercise. 

3. Use the correct “e” symbol - the bold/fancy one from the pallette; this was explained in Exercise 1(a) of Lab 4. Alternatively, the capital “E” or “Exp[...]” will work just fine. 


   f(x,y)=  (e^(-2 (x^2)+(y^2)-0.3` (y^4)))  *( 1+0.15` *(Cos[6 y]-Cos[9 (x+y)]+Sin[5 x]+Sin[7 (x-y)])* Sin[x*y] )

Make its contour plot using the range  [-1.5, 1.5]  for  x  and  [-2.5, 2.5] for  y. 

As I also explained in Exercise 2(a) of Lab 4, after you have defined a function, avoid copy-pasting its expression. That is, use the function name,  f[x,y],  not its expression, in the ContourPlot command. 
Here and below, if you copy-paste an expression where you should have used a function name, your score will be HARSHLY reduced.


Use 20 contours per plot and don’t forget to set PlotRange -> All (which you should do always).  Specify the option ContourShading->False; this will make your contours appear on a white background. Give a name to your plot.

Looking at this plot and using Notes from Sec. 14.7 about contour plots of critical points, determine the following.
   (i)  What critical points (local extrema and saddles) does this surface have? Where (i.e. for what values of  x  and  y) approximately are these points located?  If any of the types of critical points is absent for this surface, say so.  (Note, however, that you should find at least one saddle point.) Round your estimates to one decimal place. 
   (ii)  Near what point(s) does this surface appear to have the steepest slope?

Note that you may make a 3D plot of the surface to check your answers.

(b)   In this part, you will first use Theorem 2 from Sec. 14.7 to attempt to find all critical points of the surface. Then you will use Theorem 3 (a.k.a. Second Derivative Test) to determine their types (i.e., local extrema, saddles, or neither of these). 

First, define the partial derivatives of  f(x,y) denoting them as  fx(x,y) and  fy(x,y).  Look up the syntax of the derivative either in Exercise 2(a) of Lab 4 or under “Help” in the main menu. 
Again, follow the Important technical notes 2 and 3 in Exercise 2(a) of Lab 4.

Next, to find the critical points, use the command  FindRoot[{fx[x,y]==0,fy[x,y]==0},{x,x_initial},{y,y_initial}]  to numerically solve the simultaneous equations  fx(x,y)=0  and  fy(x,y)=0. Here  x_initial  and  y_initial  are the coordinates of the initial guess that you will provide (see below). An initial guess is required for an algorithm such as FindRoot to try to solve those equations. 

The initial guesses that you need to use are listed in the large Answer Box below. They fall into two categories:

   Category (i):  These three pairs of numbers are close to the estimates of the critical point locations that you found by inspection in part (a);
   
   Category (ii):  Points  Q1=(1.1, 0.05),  Q2=(1.1, -0.05), Q3=(0.2, 1.3),  and  Q4=(0.1, -1.2).
   Note:   Points Q1 and Q2 are intentionally chosen to be far from the estimated locations of the critical points. These points will be used later in Exercise 2.  
   Points Q3 and Q4 are close to the estimates of the extrema, yet are slightly different from the initial guesses in Category (i). 
   
   Technical note:   FindRoot is a somewhat “touchy” command in the following sense. If you accidentally forget to type two “=”s (i.e. “==”) and type only one “=” (i.e., type FindRoot[{fx[x,y]=0,...  instead of  FindRoot[{fx[x,y]==0,... ), Mathematica will generate some error message and will stop working correctly. The only way to overcome this is to go back to the very top of this Lab, execute the Clear[“`*”] command, and then execute all commands that you have defined above. 

Record the results which Mathematica gives for the coordinates of the critical points in each of these cases. Round the answers to 2 decimal places.  State whether you think each answer found by Mathematica is a real (i.e., actual) critical point or a spurious numerical artifact.  Here is how you can tell if an answer is spurious: 
For the initial guesses supplied in Category (i): 
   - if the answer is far from your initial guess, the answer must be spurious.  
For the guesses supplied in Category (ii): 
   - if changing the last decimal place of either coordinate of the initial guess changes the answer (also in either coordinate) considerably (e.g., by 1.0 or so, or more), then that answer is spurious. E.g., if initial guess (0.1, -1.2) yields a root (6.3, -5.4) while initial guess (0.2, -1.2) yields a root (6.4, -4.1), then this root is spurious. 
   
   Technical note:  For some of your answers, Mathematica will generate an error message “(MessageTemplate[FindRoot, cvmit, "Failed to converge to the requested accuracy or precision within 100 iterations.", , , , , Local])
This is normal; just move on, paying attention to the above instructions about real and spurious roots. Any other error message, however, is not normal, and you should check your work and perhaps resort to using the Clear[“`*”] command mentioned in the previous Technical note.

You may try to use other points as initial guesses if you wish and see what happens, but please record only the answers to (i) and (ii).


Did you find the answers found by Mathematica to be sensitive to the initial guess?  
	An answer is sensitive if a small change in the initial guess leads to a substantially larger change in the answer; it is insensitive otherwise. So, to answer the above question, look at the answer that  FindRoot  found for the initial guess  Q3  and compare it with the answer found for one of the initial guesses in  Category (i); then do a similar comparison for Q4.

(c)   Now, for each of the real (i.e., not spurious) critical points which Mathematica has found, apply the  second derivatives test to determine whether this point is a maximum, minimum, or a saddle.   Show your work.  Keep in mind the note about not copy-pasting, stated in part (a) of this Exercise. 

Suggestion:  Define the quantity denoted as D in Theorem 3 of Sec. 14.7 to be a function of x and y (again, follow the Important technical notes 2 and 3 in Exercise 2(a) of Lab 4).  Then evaluate it at the critical points (use their values rounded to 2 decimal places).   However, do not name that quantity D, because in Mathematica that symbol is reserved for a built-in function.  Also, do not attempt to simplify the expressions for any of the partial derivatives and for the quantity D (which, as per above, is not named “D”), because Mathematica will take forever to simplify them.

Clear["`*"]
(* Just in case, clear out the workspace. 

Note that the symbols  
 "(*" and "*)"  mark the beginning and end of a comment in Mathematica.
I will use them to place clarifying comments after those commands that require clarification or explanation.

*)

In the next Exercise, you will show that one can always find a local extremum if one follows the direction of the gradient. Given the relation between the direction of the gradient and the direction of the steepest ascent along a surface, this method can be called the “method of steepest ascent”. This name will be further explained as we proceed. 


Let us begin with the following
Example

Define the same function as in Exercise 1 (copy-pasting the expression below is okay in this particular case):

   f(x,y)=  (e^(-2 (x^2)+(y^2)-0.3` (y^4)))  *( 1+0.15` *(Cos[6 y]-Cos[9 (x+y)]+Sin[5 x]+Sin[7 (x-y)])* Sin[x y] )

f[x_,y_]=Exp[-(2*x^2-y^2+0.3*y^4)]*(1+0.15*(Sin[5*x]+Cos[6*y]+Sin[7*(x-y)]-Cos[9*(x+y)])*Sin[(x*y)])

This Example has two parts. In the first part, I will explain the idea behind the method of steepest ascent. In the second part I will present an efficient implementation of that idea.

Part 1 of Example  -  Idea

Select a "starting point" (this name will become clear as you proceed) to be  P=(1, 2).  The starting point is the analog of the initial guess used in part (b) of Exercise 1. 
Then plot this point along with the contour plot of  f(x,y) :

OldPointX=1; 
OldPointY=2; 
OldPoint=Point[{OldPointX,OldPointY}]
(* 
The command  Point  tells Mathematica that the object you have defined is a point in the plane.
*)

w1=ContourPlot[f[x,y],{x,-1.2,1.2},{y,-2.3,2.3},PlotRange->All,ContourShading->False,Contours->20];
(*
  This is the contour plot.
*)
pts=List[OldPoint];  
(*
This initiates a list of points to be used below.
*)
w2=Graphics[{PointSize[0.05],pts}];
(*
This plots the above list of points. For now, we have only one point in this list, but there will be more later.
*)
Show[w1,w2]
(*
This plots the contour plot and the list of points together.
*)

Technical note:  If the contour plot that you see here looks any different than the one you obtained in Exercise 1, then go back to that Exercise and redefine your f[x,y] there, making sure that you obtain the same contour plot as here. Then redo those of your steps in Exercise 1 that require correction.

Now, from Exercise 1, you know that the function has two local maxima.  You may suspect that if you keep moving in the direction of positive slope on the surface, you will eventually reach one of those maxima.  It is also reasonable to suspect that if you move in the direction of the steepest slope (i.e., the steepest ascent), you will need the least amount of time to reach that maximum. 

As you remember, the direction of the steepest ascent (i.e., the fastest increase of the function) occurs in the direction of the gradient.  Therefore, let us move into the direction of the gradient by the amount

                                                 |(Overscript[▽, →])f | dt
                                                 
per step, where  |(Overscript[▽, →])f |  is the length of the gradient and  dt is some small number denoting the size of the step. Let us take:

dt=0.2
(* This rather large value is chosen for illustration purposes only. Later we will set  dt  to a smaller value. *)

Define the components of the gradient, and then execute the subsequent commands. The purpose of each command is explained after it in the form of a comment.

fx[x_,y_]=D[f[x,y],x];
fy[x_,y_]=D[f[x,y],y];
(* 
These are the partial derivatives of f[x,y], the same ones that you were supposed to define in Exercise 1. They are redefined here because we have cleared the workspace after Exercise 1.  
*)

NewPointX=OldPointX+fx[OldPointX,OldPointY]*dt;
NewPointY=OldPointY+fy[OldPointX,OldPointY]*dt;
NewPoint=Point[{NewPointX,NewPointY}]
(* 
Add the vector   (Overscript[▽, →])f * dt  to the "old point" to obtain the "new point". Thus, you have started at OldPoint and "stepped" in the direction of the steepest ascent.
*)

pts=Append[pts,NewPoint]
(*
Append the "new point" to the existing list of points.
*)

w3=Graphics[{PointSize[0.05],pts}];
(*
Prepare to plot the augmented list of points.
*)

Show[w1,w3]
(*
Show the starting point and the new point together with the contour plot.
*)

The old and new points are quite close, with the new point being a bit closer to the peak.

Now repeat the entire process:

OldPointX=NewPointX;
OldPointY=NewPointY;
(* 
The "new point" at the previous step becomes the "old point" at the current step.
*)

NewPointX=OldPointX+fx[OldPointX,OldPointY]*dt;
NewPointY=OldPointY+fy[OldPointX,OldPointY]*dt;
NewPoint=Point[{NewPointX,NewPointY}]
(* 
This finds the "new point" at the current step as we continue to move in the direction of the gradient.
*)

pts=Append[pts,NewPoint];
(*
Append the "new point" to the existing list of points.
*)

wnew=Graphics[{PointSize[0.05],pts}];

Show[w1,wnew]

Your next "new point" should be closer to the peak than at the previous iteration. Return the cursor into the command box above and keep pressing "Shift+Enter" to see your iterations move towards the peak. This happens for early iterations. Soon, however, when the point reaches some vicinity of the pinnacle (i.e., the point with  x ≈ -0.02  and  y ≈  1.30, which you found in Exercise 1(b)), the iterations begin to overshoot the pinnacle and never appear to settle on it. (You should also look at the coordinates of your point - they appear right below the command box and right above the plot - and notice how they change.) This overshooting behavior occurs because your step  dt  has been chosen too large. In the next part of this Example you will see how to efficiently perform many iterations with a smaller  dt.

You may now treat yourself to the 3D view of the result of your work:

pSurf3D=Plot3D[f[x,y],{x,-1.2,1.2},{y,-2.3,2.3},AxesLabel->{"X","Y"," "},PlotRange->{0,2.7},Mesh->False];
(*
Plot the surface in 3D.
*)

pLast3D=Graphics3D[{PointSize[0.05],Point[{NewPointX,NewPointY,f[NewPointX,NewPointY]}]}];
(* 
Plot the point obtained at the last iteration, in 3D.
*)

Show[pSurf3D,pLast3D]

To complete this example, let us see if we can find the maximum by numerically solving the simultaneous equations  (Subscript[f, x])(x,y)=0,  (Subscript[f, y])(x,y)=0  with the initial guess  (1, 2)  (which is the starting point we used in this Example):

FindRoot[ {fx[x,y] == 0, fy[x,y] == 0}, {x,1},{y,2} ]

As you see, the answer Mathematica gives has nothing to do with the peak. Worse yet, even if you start rather close to the peak (at one of the points of the period-3 orbit found above):

FindRoot[ {fx[x,y] == 0, fy[x,y] == 0}, {x,0.27},{y,1.32} ]

the numerical answer is still not relevant to the peak!

To conclude from this example, you have shown that the method of "steepest ascent" is capable of leading you to a relative maximum. On the other hand, solving the equation   
				(Overscript[▽, →])f  = (Overscript[0, ->])
numerically will not, in general, give you the peak's coordinates.

FYI:
This observation is behind the idea of a large family of widely used numerical methods called "methods of steepest descent" (which are analogous to those of "steepest ascent" except that they find relative minima of functions by moving in the opposite direction of the gradient).  They are studied in courses on Numerical Methods or Numerical Analysis.

Part 2 of Example  -  Efficient Implementation

In order to see your iterations move towards the peak in Part 1 of the Example, you had to perform a somewhat time-consuming process of placing your cursor into a command box, pressing “Shift+Enter”, observing a point being added to the plot, and then repeating this process a few times. While in this manner you can clearly see what happens after each iteration, this process requires repetitive effort. 
     Also, if you “Delete All Output” and then restore it via “Evaluate Notebook”, you (and I!) will have to go through the entire process again in order to restore your final plot. Fortunately, there is an automated method to perform several similar iterations at once: it is called “looping”. 
     Below I show you an example involving one of Mathematica’s loops -- the “For” loop. Namely, the single command box below performs 10 iterations and generates a single plot at the end showing their result.

Clear["wnew"]
OldPointX=1; 
OldPointY=2; 
OldPoint=Point[{OldPointX,OldPointY}];
(* 
Use the same starting point as above.
*)

dt=0.05;
(* 
We are now choosing a smaller dt in the hopes that iterations with smaller steps will converge on the pinnacle. 
*)

pts=List[OldPoint];
(* 
Define the step size and initiate the list of points. 
*)

NewPointX=OldPointX+fx[OldPointX,OldPointY]*dt;
NewPointY=OldPointY+fy[OldPointX,OldPointY]*dt;
NewPoint=Point[{NewPointX,NewPointY}];
(* 
Outside the "loop", obtain the first "new point".
*)

For[i=0,i<10,i++,
OldPointX=NewPointX;
OldPointY=NewPointY;
NewPointX=OldPointX+fx[OldPointX,OldPointY]*dt;
NewPointY=OldPointY+fy[OldPointX,OldPointY]*dt;
NewPoint=Point[{NewPointX,NewPointY}];
pts=Append[pts,NewPoint];
wnew=Graphics[{PointSize[0.05],pts}];
]
(*
Execute commands in a loop 10 times.
The parameter "i" does not play any role other than just counting the number of iterations.
*)

NewPoint
(* Display the coordinates of the last point of the iterations.
*)

Show[w1,wnew]
(*
Show the results of all the 10 iterations combined.
*)

End of Example 
___________________________________

Two obvious questions still remain:

   (i)  As you know from Exercise 1, the surface in question has two peaks. What determines which of these two peaks you will end up at?

   (ii)  In addition to the peaks, you were also supposed to find (at least one) saddle point of this surface.  Then, under what conditions (if any) will the method of steepest ascent lead you to the saddle point and keep you there?

You will examine these questions in Exercise 2.

Exercise 2

(a)   Repeat the calculations of the above Example starting at the point  Q1=(1.1, 0.05) (this is one of the initial guesses used in part (b) of Exercise 1).  Throughout this Exercise, use  dt = 0.05. Show your work and produce a plot showing all your iterations (i.e., you need to plot the entire list of your iterations, not just the last point on that list).  

Your goal is to use as many iterations as needed for the trajectory to settle at a vicinity of one of the extrema (see the next paragraph for the precise definition of “vicinity”). You will have to experiment to find the minimum required number of iterations. (You can, of course, put in “10000” for that number, which should be enough for the trajectory to settle, but then your loop will take a long time to run, and I will subtract points for that.)  

Use the “For” loop described in Part 2 of the Example before this Exercise. As mentioned above, you will need to decide by experimentation how many iterations you need to do in order to achieve the following:
- When each coordinate of your solution is rounded to 2 decimal places, the 2nd decimal place must not change with subsequent iterations.  
Example:  Suppose that iterations 789, 790, 791 yield solutions  0.317, 0.314, 0.308. When rounded to 2 decimals, these are 0.32, 0.31, 0.31. Thus, you would stop at iteration 791. 

Technical notes:

1.  You should introduce only minor changes to the “For” loop, as described in the first two paragraphs of this Exercise. In particular, you should not attempt to change its syntax. 

2.  As I have stressed several times already, you should follow the Important technical note 3 of Exercise 2(a) of Lab 4. If you neglect to explicitly indicate arguments of the functions that you use in the “For” loop, Mathematica may freeze, causing you to waste time reviving it, and you may even risk losing all your work.

3.  Give a name to your last plot.

Side note:  In “real life” problems, one does not experiment with the number of iterations. Instead, one programs into one’s code a “termination condition”. Such a condition compares the results of two or more consecutive iterations and stops the loop if those results are found to be sufficiently close of one another (say, within 10^(-12)). If the iterations have indeed produced several consecutive close-by results, then one says that they have found a critical point. Alternatively, iterations may spend some short time near a critical point and then move on elsewhere; in this case they fail to find that particular critical point.

(b)    Repeat part (a) of this Exercise for the starting point being  Q2=(1.1,-0.05).  
Then plot the results of your iterations together with the combined plot obtained in part (a).

(c)    Record the x- and  y-coordinates of your last iterations in parts (a) and (b). Round them to two decimal places. Also record the coordinates of the true local extrema found by Mathematica in Exercise 1(b). (These are the answers which you were supposed to call “real” in the large answer box at the end of 1(b); they are also the ones which agree with the local extrema found in 1(a).)

Have you been able to find the saddle point using the method of steepest ascent? (Please see the Side note at the end of part (a) regarding the meaning of the word “find” here.)
Whether you answer `yes’ or `no’, make sure to describe the behavior of the trajectory near the saddle point.

Examine the last plot in part (b). 
	What can you observe about the direction of the trajectories relative to the level curves? 
	State a fact from Sec. 14.6 that explains this observation.