COMPUTER PROJECT 1

YOUR NAME HERE

Submit to BLACKBOARD with the title:
LastNameFirstInitialProject1.nb,
for example,
GoodmanFProject1.nb

HOMEMADE CONTOUR PLOTS

In this part of the project, we find how to draw a countour of a function of two variables by solving a differential equation.  Let's take the function
f[x, y]  =  (x^3) - 2 x y + (y^3)  
and work on plotting a contour such as
  (x^3) - 2 x y + (y^3)  =  22 
 which runs through the point (1, 3).
 
 The idea presented in class is that the gradient of the function is always perpendicular to the contour lines.  Let's enter the function and compute the gradient:

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

Clear[gradf,  m,  n]
m[x_,  y_]  = (∂_x) f[x, y]
n[x_,  y_]  =  (∂_y) f[x, y]
gradf[x_, y_]  :=  {m[x, y],  n[x, y]}

Now the tangent direction to a contour line is a direction perpendicular to the gradient.  Following the discussion in the lecture, define a vector
 velocity[jx, y]
 with the following properties:
 (1)  velocity[x,y] is perpendicular to gradf[x, y].
 (2)  velocity[x, y]  points in the direction obtained by rotating gradf[x,y] by 90 degrees clockwise in the plane.
 (3)  velocity[x,y] is a unit vector.
 
In the following cell, write in a formula for velocity[x, y].
Note:  The length of a vector v  =  {a,  b}  is computed in Mathematica by
Norm[v] or Norm[{a,  b}]

velocity[x_,  y_] :=  ENTER FORMULA HERE

We extract the two components of the velocity vector field for convenience:

v1[x_,  
        y_]  :=  velocity[x,  y][[1]]
v2[x_,  
        y_]  :=  velocity[x, y][[2]]

The differential equation for a level curve through a point (a,  b)  is

x' (t)  =  v1[ x,  y]
y'[t]  =  v2[x,  y]
x[0]  =  a
y[0]  =  b

This instructs a "blind particle" to start out at (a, b) and, if it finds itself at some position (x, y)  at some time t, to continue with vector velocity velocity[x, y].  In other words, the particle is instructed to move with speed 1 in the direction given by the velocity vector velocity[x,y].

The next cell defines a Mathematica function which for the input [a, b], solves the differential equation for the particle starting at position (a,  b).

Clear[curve]
curve[a_, b_]  :=  
      NDSolve[
{
x'[t]  == 
            v1[x[t], y[t]],  
y'[t]  ==  
            v2[x[t], y[t]],  
x[0]  ==  
            a,  
y[0]  == 
            b
}, 
{x,  
          y},    
{t,  -100,  100}]

Let's compute curve[1,  3].

curve[1,  3]

The result is an "InterpolatingFunction". That is just a function given by a table of numerical data.   To see the result we have to plot it as a parametric curve.  When you evaluate the next input cell, you should get a curve which looks like this.  If you do not get such a curve, then you must have got the wrong formula for the velocity vector, and you had better go back to fix your
 work.

ParametricPlot[
      Evaluate[{x[t],  y[t]} /. curve[1, 3]],  {t,  -6, 9},  
      AspectRatio ->  1]

Let's compare this with the output of Mathematica's  ContourPlot[] function.

ContourPlot[
      f[x, y], {x,  -5,  5},  {y,  -5, 
        5},   Contours -> {22},  
      ContourShading ->  False]

Another interesting contour is that through (1/3,  1/3). You should get the following picture:

ParametricPlot[
      Evaluate[{x[t],  y[t]} /. curve[1/3, 1/3]],  {t,  -6, 9},  
      AspectRatio ->  1]

The following cell defines a function countour[ ] which calculates the contour corresponding to the starting point (a, b).

Clear[  contour]
contour[{a_, b_}]  := 
      ParametricPlot[
        Evaluate[{x[t],  y[t]} /. curve[a, b]],  {t,  -20, 20},  
        AspectRatio ->  1,  
        PlotRange ->  {{-8, 8},  {-8, 8}}, 
        DisplayFunction ->  Identity]

The following cell lists a bunch of starting points and then calculates and displays the contour for each of these starting points.  The result is a home made contour plot for the function f[x, y].

Clear[ startPoints,  
      plots]

    startPoints  =  {
{1, 1},  
{2, 
          2}, 
{3,  3}, 
{4, 
          4},  
{1/3,  
          1/3}, 
{-1/2, -1/
            2},  
{-1,  -1}, 
{-2,  -2}, 
{-3, -3}, 
{-4,  -4}
}

    plots  =  Map[ contour,  startPoints];

    Show[plots,  
      DisplayFunction ->  $DisplayFunction]

Two questions to investigate.  These questions will require you to do further computations,  by adapting the computations above;

(1)   Try to find out for which values of  a the contour through (a,  a)  (like for a  =  1/3)  is a closed loop and for which values of a the contour is not a closed loop (like for a  =  2).    To do this, you might want to copy the definition of countour[ ], change the plot range to something more appropriate (smaller), and draw one or more new contour plots with different starting points.   Be careful to avoid the starting points
(0,0) and (2/3, 2/3), since the gradient is zero at those points!

(2)   All of the non-closed loop contours seem to be asymptotic to some line.  What is that line?    To do this, you might want to copy the definition of countour[ ], change the plot range to something more appropriate (larger), and draw one or more new contour plots with different starting points.

DO YOU WORK FOR THESE TWO QUSTIONS BELOW HERE:

CRITICAL POINTS

The function 
f[x, y]  =  (x^3) - 2 x y + (y^3)
has two critical points as shown by the folloiwng computation:

Clear[x, y, f]
f[x_, y_] := (x^3) - 2 x y + (y^3);
Clear[gradf,  m,  n]
m[x_,  y_]  = (∂_x) f[x, y]
n[x_,  y_]  =  (∂_y) f[x, y]
gradf[x_, y_]  :=  {m[x, y],  n[x, y]};

We solve for critical points.  Only two of the four solutions are real (as opposed to complex valued).

Solve[
      gradf[x, y] == 0,  {x,  
        y}]

:Let's look at contour plots and ordinary plots near the two critical points.

ContourPlot[
      f[x, y],  {x,   .6,   .73},  {y,   .6,   .73}]
Plot3D[f[x, y],  {x,   .6,   .73},  {y,   .6,   .73}]

ContourPlot[
      f[x,  y],  {x,  - .1,   .1},  {y,  - .1,   .1}]
Plot3D[
      f[x, y], {x,  - .1,   .1},  {y,  - .1,   .1} ]

Two questions:

(1)  Describe the behavior of the function f[x, y] near each of the two critical points.  Are the critical points maxima, minima, or saddle points?   This question requires a verbal answer but no computation.

(2)  Compare the behavior of contours which start near the critical point (0,0), but at points
(a, a)
(a, -a)
(-a, a)
(-a, -a)
for a small.
This question requires additional computations AND a verbal answer.

(3)  How do contours behave which start near the critical point (2/3, 2/3)?  Do all contours starting near (2/3, 2/3) behave the same?
This question requires additional computations AND a verbal answer.

DO YOUR WORK FOR THESE QUESTIONS BELOW HERE: