Project Help Notebook by Keith Stroyan

Help NoteBook Copyright 1998 by Academic Press, Inc. - All rights reserved.
Modify this Notebook to prepare your project report.

Project 33: The Local Second Derivative Test for Max-min & Locally Stable Equilibria

by A. Student and B. Pupil

with technical help from Da LabMonitor and Help Notebook

Initialization

Off[General::spell];
Off[General::spell1];
Off[InterpolatingFunction::dmval];
Off[ParametricPlot::pptr];
Off[Graphics::gprim];
Off[Graphics::gptn];
Off[NDSolve::ndsz];

Needs["Graphics`Colors`"]
	
colorQ[a_]:=(Head[a]===RGBColor);
grayQ[a_]:=(Head[a]===GrayLevel);
thickQ[a_]:=(Head[a]===Thickness);
absthickQ[a_]:=(Head[a]===AbsoluteThickness);
notruleQ[u_] := !(Head[u]===Rule);
ruleQ[u_] := (Head[u]===Rule);

Norm::usage = "Norm[expr] returns the Euclidean length of the
		vector, 'expr'.";
Norm[a_?VectorQ] := Sqrt[a.a];

tol = .5*10^-6;

Vector2D[tail_,tip_,opts___]:=
	Module[{ntail,ntip,sopts,Punit,up,arrow},
	ntail = N[tail];
	ntip = N[tip];
	
	sopts = Join[Select[{opts},colorQ],Select[{opts},
				thickQ],Select[{opts},absthickQ],Select[{opts},grayQ]];

If[Chop[ntail.ntail+ntip.ntip,tol] == 0,
			(*too small to plot or not number*)
			Return[Graphics[Point[{0,0}]]];
	,						
If[Chop[(ntip-ntail).(ntip-ntail),tol] == 0,
		Return[Graphics[{Join[sopts,{Line[{ntail,ntip}]}]}]];
	,
	Punit = (.075 Norm[ntip-ntail]/Norm[ntip-ntail]) (ntip-ntail) ;
	up = 0.15 {-Punit[[2]],Punit[[1]]};
	arrow = Line[{ntail,ntip,ntip-Punit+up,ntip,ntip-Punit-up}];
	]];
Return[Graphics[{Join[{AbsoluteThickness[2]},Join[sopts,{arrow}]]}]];
];

		
uArrow[vec_,opts___]:= 
	Module[{a,b},
		If[Norm[N[vec]]<1,b = vec,b=vec/Norm[vec]];
				a=ToExpression["Tail"]/.Select[{opts},ruleQ]/.
					{ToExpression["Tail"]->{0,0}};
				Return[Vector2D[a,a+b,opts]]
				];

flow2D[{f_,g_},{x_,y_,t_},initVals_,tFinal_,dt_,opns___] :=
 Module[{soln,arr,pt,crv,tt},
soln = Flatten[Map[
	NDSolve[{x[0]==#[[1]],y[0]==#[[2]],
	x'[t] == f[x[t],y[t]],y'[t] == g[x[t],y[t]]},
	{x[t],y[t]},{t,0,tFinal}]&,
										initVals],1];
											
pt[tt_] := Map[Graphics[{Blue,AbsolutePointSize[4],Point[{#[[1]],#[[2]]}]}]&,
{x[t],y[t]}/.soln/.t->tt];
										
crv[tt_] := ParametricPlot[Evaluate[{x[t],y[t]}/.soln],{t,0,tt},opns,
			AspectRatio->1,DisplayFunction->Identity];
			
Show[{crv[dt/100],pt[0]},opns,AspectRatio->1,
DisplayFunction->$DisplayFunction];
Do[Show[{crv[tt],pt[tt]},DisplayFunction->$DisplayFunction];,
{tt,dt,tFinal,dt}];
	Return[soln];	
];

vectorFlow2D[{f_,g_},{x_,y_,t_},initVals_,tFinal_,dt_,opns___] :=
 Module[{soln,arr,pt,crv,tt},
soln = Flatten[Map[
	NDSolve[{x[0]==#[[1]],y[0]==#[[2]],
	x'[t] == f[x[t],y[t]],y'[t] == g[x[t],y[t]]},
	{x[t],y[t]},{t,0,tFinal}]&,
										initVals],1];
								
	arr[tt_] := Map[uArrow[{f[#[[1]],#[[2]]],g[#[[1]],#[[2]]]},Tail->#,Red]&,
	{x[t],y[t]}/.soln/.t->tt];
	
pt[tt_] := Map[Graphics[{Blue,AbsolutePointSize[4],Point[{#[[1]],#[[2]]}]}]&,
{x[t],y[t]}/.soln/.t->tt];
										
crv[tt_] := ParametricPlot[Evaluate[{x[t],y[t]}/.soln],{t,0,tt},opns,
			AspectRatio->1,DisplayFunction->Identity];
			
Show[{crv[dt/100],pt[0],arr[dt/100]},opns,AspectRatio->1,
DisplayFunction->$DisplayFunction];
Do[Show[{crv[tt],arr[tt],pt[tt]},DisplayFunction->$DisplayFunction];,
{tt,dt,tFinal,dt}];
	Return[soln];	
];

SetOptions[Graphics,Axes->False,AxesOrigin->{0,0},AspectRatio->Automatic];
SetOptions[Graphics3D,ViewPoint->{5,2,1},Boxed->False,AspectRatio->Automatic];
SetOptions[ParametricPlot3D,ViewPoint->{5,2,1},AspectRatio->Automatic];
SetOptions[Plot3D,ViewPoint->{5,2,1},AspectRatio->Automatic];
SetOptions[ParametricPlot,AspectRatio->Automatic];
SetOptions[Plot,AspectRatio->Automatic];

colorQ[a_]:=(Head[a]===RGBColor);
grayQ[a_]:=(Head[a]===GrayLevel);
thickQ[a_]:=(Head[a]===Thickness);
absthickQ[a_]:=(Head[a]===AbsoluteThickness);
notruleQ[u_] := !(Head[u]===Rule);
ruleQ[u_] := (Head[u]===Rule);

Norm::usage = "Norm[expr] returns the Euclidean length of the
		vector, 'expr'.";

Norm[a_?VectorQ] := Sqrt[a.a];

Angle[a_?VectorQ,b_?VectorQ] := ArcCos[(a.b)/(Norm[a] Norm[b])]

Proj[a_?VectorQ,b_?VectorQ] := ((a.b)/(b.b))b

tol = .5*10^-6;

OrthUp[a_?VectorQ] :=
	Module[{b,c,d,f},
		d = N[a/Norm[a]];
		If[Chop[d[[3]]^2 - 1,tol] == 0, b={1,0,0} ,
		b = {0,0,1} - d[[3]] d; b = b/Norm[b];];
	Return[b]];
	
OrthSet[a_?VectorQ] :=
	Module[{b,c,d},
		d = N[a/Norm[a]];
		If[Chop[d[[1]],tol]==0,b={1,0,0},
			 If[Chop[d[[2]],tol]==0,b={0,1,0},
			  b={-d[[2]],d[[1]],0}]];
			b=b/Norm[b];c = Cross[d,b];
	Return[{b,c/Norm[c],d}]];

				
Rotate[a_?VectorQ,theta_] :=
	Module[{u,v,r},
		u = OrthSet[a];
		v = Transpose[u];
		r = N[{{Cos[theta],-Sin[theta],0},{Sin[theta],Cos[theta],0},{0,0,1}}];
	Return[v.r.u]];

Vector2D[tail_,tip_,opts___]:=
	Module[{ntail,ntip,sopts,Punit,up,arrow},
	ntail = N[tail];
	ntip = N[tip];
	
	sopts = Join[Select[{opts},colorQ],Select[{opts},
				thickQ],Select[{opts},absthickQ],Select[{opts},grayQ]];

If[Chop[ntail.ntail+ntip.ntip,tol] == 0,
			(*too small to plot or not number*)
			Return[Graphics[Point[{0,0}]]];
	,						
If[Chop[(ntip-ntail).(ntip-ntail),tol] == 0,
		Return[Graphics[{Join[sopts,{Line[{ntail,ntip}]}]}]];
	,
	Punit = 0.075 (ntip-ntail) ;
	up = 0.15 {-Punit[[2]],Punit[[1]]};
	arrow = Line[{ntail,ntip,ntip-Punit+up,ntip,ntip-Punit-up}];
	]];
Return[Graphics[{Join[{AbsoluteThickness[2]},Join[sopts,{arrow}]]}]];
];

Vector3D[tail_,tip_,opts___]:=
	Module[{ntail,ntip,sopts,Punit,up,u2,u3,u4,rotmat,arrow},
	ntail = N[tail];
	ntip = N[tip];
	
	sopts = Join[Select[{opts},colorQ],Select[{opts},
				thickQ],Select[{opts},absthickQ],Select[{opts},grayQ]];

If[Chop[ntail.ntail+ntip.ntip,tol] == 0,
			(*too small to plot or not number*)
			Return[Graphics3D[Point[{0,0,0}]]];					
	,		
If[Chop[(ntip-ntail).(ntip-ntail),tol] == 0,
		Return[Graphics3D[{Join[sopts,{Line[{ntail,ntip}]}]}]];
	,
	Punit = (ntip-ntail);
	rotmat = Rotate[Punit,N[Pi/2]];
	up = OrthUp[Punit];
	Punit = .075 Punit;
	up = .15 Norm[Punit] up;
	u2 = rotmat.up;
	u3 = rotmat.u2;
	u4 = rotmat.u3;
	arrow = Line[{ntail,ntip,
				ntip - Punit + up,ntip,
				ntip - Punit + u2,ntip,
				ntip - Punit + u3,ntip,
				ntip - Punit + u4}];
	]];
Return[Graphics3D[{Join[{AbsoluteThickness[2]},Join[sopts,{arrow}]]}]];
];

aRCtAN[y_,x_] := Module[{v},
	If[ Chop[N[x],tol] == 0 ,If[ N[y] < 0, v = -N[Pi/2], v = N[Pi/2]],
		v = N[ArcTan[y/x]];
		If[ x < 0 && y < 0, v = N[Pi -  v],v]];
	Return[v]];

Arc2D[X_,Y_,opts___] := Graphics[{opts,
		Disk[{0,0},1,{Min[
							aRCtAN[
						N[X[[2]]],N[X[[1]]]
								]
							,
							aRCtAN[
							N[Y[[2]]],N[Y[[1]]]
								]
						]
						,
					Max[
							aRCtAN[
						N[X[[2]]],N[X[[1]]]
								]
							,
							aRCtAN[
							N[Y[[2]]],N[Y[[1]]]
								]
						]}
				]}
				];
Arc3D[X_,Y_] := 
	Module[{t,U,V},
		U[V_]:= N[(1/Sqrt[V.V]) V];
		Return[Graphics3D[
		Polygon[Append[Table[U[t X + (1-t) Y],{t,0,1,0.1}],{0,0,0}]]
				]
		]
	];

(*Zero Arrow causes trouble. Put in Vector?*)

Vector[tail_,tip_,opts___]:=
	If[Length[N[tail]]==2,Return[Vector2D[tail,tip,opts]],
		Return[Vector3D[tail,tip,opts]]];
		
Arrow[vec_,opts___]:= 
	Module[{a,b},
				If[Length[N[vec]]==2,
				a=ToExpression["Tail"]/.Select[{opts},ruleQ]/.
					{ToExpression["Tail"]->{0,0}};
				Return[Vector2D[a,a+vec,opts]],
				a=ToExpression["Tail"]/.Select[{opts},ruleQ]/.
					{ToExpression["Tail"]->{0,0,0}};
				Return[Vector3D[a,a+vec,opts]]]
				];

Frame2D[delta_] := 
      Graphics[{
            {Green,Line[{{-delta/4,0},{delta,0}}]},
                    Text["x",{1.01 delta,0}],
            {Green,Line[{{0,-delta/4},{0,delta}}]},
                    Text["y",{0, 1.01 delta}]}];

Frame3D[delta_] := 
      Graphics3D[{
            {Green,Line[{{-delta/4,0,0},{delta,0,0}}]},
                    Text["x",{1.01 delta,0,0}],
            {Green,Line[{{0,-delta/4,0},{0,delta,0}}]},
                    Text["y",{0, 1.01 delta,0}],
            {Green,Line[{{0,0,-delta/4},{0,0,delta}}]},
                    Text["z",{0,0,1.01 delta}]}];

View2D[stuff_] := Show[stuff,Frame2D[1]];
View3D[stuff_] := Show[stuff,Frame3D[1]];

My Project Introduction

Write this last.  Begin with the basic question and state your main finding clearly, but with few technical details.  You might even modify the title of this section to make it more interesting.

Beginning Explorations or Computer Experiments

We begin with a computer explanation of steepest ascent described in the print project.  We have included some fussy code to help you see how the metohd works and some exercises to check both your computer and mathematical understanding.

Ascent and Descent by the Gradient of Height

Suppose the height of a graph over a region of the (x,y)-plane is given by a function,

z = h[x,y]

for example, enter the next function:

Clear[x, y, z, h, f, g, gradh];
h[x_, y_]  =  
      x y (E^(-(1/2) ((x^2) + (y^2))))

Theorem 18.12 of the text tells us that the gradient vector

∇h = (({{∂(h/∂)x}, {∂(h/∂)y}}))

when evaluated above a particular  (x,y), points in the (x,y)-direction that will make  h[x,y]  increase fastest if we change  x  and   y  in that direction.

For example, the function h[x,y] entered above has the following general symbolic gradient (enter it now):

f[x_, y_]  = Factor[ Evaluate[D[h[x, y], x]]];
g[x_, y_] = Factor[ Evaluate[D[h[x, y], y]]];
gradh[x_, y_] = {f[x, y], g[x, y]};
MatrixForm[gradh[x, y]]

At the point (x,y) = (2,3), h[x,y]  will change at the rate |∇h[2,3]|  if we change  (x,y)  in the direction  ∇h[2,3],

Print["Rate of fastest increase in height = ", 
        (Sqrt[gradh[2, 3] . gradh[2, 3]]), " ≈ ", 
        (Sqrt[gradh[2. , 3. ] . gradh[2. , 3. ]]), 
        " units up per unit over in direction = ", 
        MatrixForm[gradh[2, 3]], " ≈ ", 
        MatrixForm[gradh[2. , 3. ]]];

This information is more revealing when dispalyed on a topographical map of our mountain:

topo = 
        ContourPlot[h[x, y], {x, -4, 4}, {y, -4, 4}, 
          FrameLabel -> {"x", "y"}, RotateLabel -> False, 
          PlotPoints -> 30, DisplayFunction -> Identity];
topo = 
        Show[topo, Arrow[ 10 gradh[1/2, 3], Tail -> {1/2, 3}, Blue], 
          Graphics[{AbsolutePointSize[4], Blue, Point[{1/2, 3}]}], 
          Arrow[10 gradh[-2, 0], Tail -> {-2, 1}, Orange], 
          Graphics[{AbsolutePointSize[4], Orange, Point[{-2, 1}]}], 
          DisplayFunction -> $DisplayFunction];
Print["Gradients at (1/2,3) and (-2,1)."];

We can describe the path of steepest ascent by saying a "nearsighted" mountain climber moves at the (x,y)- velocity given by the gradient,

(dX/dt) =  ∇h[X]

If we give a starting value of the climb, X[(Subscript[t, 0])]  =  (Subscript[X, 0]) = {(Subscript[x, 0]),  (Subscript[y, 0])}, then the initial value problem that says, "Start at (Subscript[X, 0]), and climb in the direction of fastest ascent at forward (x,y)-speed |∇h[X]|." is the mathematical problem:

Find X[t] = {x[t],y[t]} so that

X[(Subscript[t, 0])]  =  (Subscript[X, 0])  and  (dX/dt) =  ∇h[X]

Clear[soln];
tFinal = 40; 
    soln =  NDSolve[
        {x'[t]  ==  f[x[t], y[t]],  y'[t]  ==  g[x[t], y[t]],  
          x[0] ==   .5, y[0]  ==  3}, {x[t], y[t]}, {t, 0, tFinal}];
{{x1[t_], y1[t_]}} = {x[t], y[t]} /. soln;
curve1  =  
        ParametricPlot[{x1[t], y1[t]}, {t, 0, tFinal}, 
          PlotStyle -> {Thickness[0.01], Blue}, PlotRange -> All, 
          DisplayFunction -> Identity];
Show[topo, curve1, DisplayFunction -> $DisplayFunction];

Computing Exercise:

Add the fastest ascent solution path to the topo graphic beginning at the orange point (-2,1).

??

Motion in 3-D

The derivative  (dX/dt)  gives the velocity of the (x,y)-motion, but not the climber's velocity on the  3D  mountain.  The height varies by  z = h[x,y], so if we have a parametric curve  {x[t],y[t]}  in (x,y)-space, the  3D  curve is given by

XYZ[t] = (({{x[t]}, {y[t]}, {z[t]}}))

where  z[t] = h[x[t],y[t]]

\(\(D3Path\  = \ 
        ParametricPlot3D[{x1[t], y1[t], h[x1[t], y1[t]]}, {t, 0, tFinal}, 
	
          Ticks -> None, Boxed -> False, 
          Prolog -> {AbsoluteThickness[4]}, DisplayFunction -> Identity];
      \)\)
hill  =  
        Plot3D[h[x, y] -  .01, {x, -4, 4}, {y, -4, 4}, 
          Ticks -> False, AxesLabel -> {"x", "y", "z"}, 
          Mesh -> False, BoxRatios -> {1, 1, 1}, PlotPoints -> 40, 
          DisplayFunction -> Identity];
Show[hill, D3Path, 
        Graphics3D[
          {AbsolutePointSize[4], DarkGreen, Point[{1/2, 3, h[1/2, 3]}]}], 
        ViewPoint -> {5, 3, 1}, PlotRange -> All, 
        DisplayFunction -> $DisplayFunction];

Computer Exercise

Give a 3D view of the fastest ascent solution path beginning at the (x,y)-point (-2,1).

Exercise: 3D Velocity and Speed

Show that the velocity of the 3-D motion of the fastest ascent climber is

V[t] = (({{(dx/dt)}, {(dy/dt)}, {∇h[x[t], y[t]].(({{(dx/dt)}, {(dy/dt)}}))}}))

and the speed of the 3D motion is

\(\@\(\(|\)\(∇h[x[t], y[t]]\)\( | \^2\)\(+\(\(|\)\(∇h[x[t], y[t]]\)\( | \^4\)\)\)\)\)

In particular, show that the rate of vertical ascent is at the rate \(\(\(|\)\(∇h[x[t], y[t]]\)\( | \^2\)\)\) up per unit time for the fastest ascent motion.

Hint: If we move in the (x,y)-direction of  ∇h[x[t],y[t]] what rate up per unit over do we climb?  Notice that we move over at speed   ∇h[x[t],y[t]]  for the fastest ascent motion.

Note

The remainder of the project concentrates on analyzing the 2D curves  (x[t],y[t]) for gradient ascent and descent.  This analysis will show when paths near a critical point are attracted to it or repelled from it.

Exercise: Descent

The fastest"nearsighted" way down a mountain of height  h[x,y] moves at the (x,y)- velocity given by,

(dX/dt) =  -∇h[X]

a) Explain why.

b) Explain why this is "nearsighted."

If we give a starting value of the climb, X[(Subscript[t, 0])]  =  (Subscript[X, 0]) = {(Subscript[x, 0]),  (Subscript[y, 0])}, then the initial value problem that says, "Start at (Subscript[X, 0]), and move in the direction of fastest descent at forward (x,y)-speed |∇h[X]|." is the mathematical problem:

Find X[t] = {x[t],y[t]} so that

X[(Subscript[t, 0])]  =  (Subscript[X, 0])  and  (dX/dt) =  -∇h[X]

c) Plot the path of steepest descent starting at (1/2,3)  on the topographical map of the function  h[x,y]  above.

d)  Explain why the descent paths are the same (x,y)-paths as steepest ascent paths, but moving the opposite direction.  In particular, if all ascent paths in a small neighborhood of a point ((Subscript[x, c]), (Subscript[y, c])) move toward the point, what is the behavior of descent paths in this neighborhood?

Dynamic Analysis

We begin with some hints needed in your project report.  See the print project and section above for help completing them.

Beginning Hint: Moving Up Hill

Suppose the height of a graph over a region of the (x,y)-plane is given by a function,

z = h[x,y]

Also suppose parametric (x,y)-paths are given by the gradient dynamical system (or gradient differential equation)

(dX/dt) =  ∇h[X]

and any choice of initial values, where  X=(({{x}, {y}}))  and  ∇h[X] = (({{∂(h/∂)x[x, y]}, {∂(h/∂)y[x, y]}})).

Explain why all such paths flow away from the location of minima of  h[x,y]  and toward the location of maxima of h[x,y].

Hint: Getting Stuck

Critical points for maxima and minima in the interior of a domain of definition of  h[x,y]  are those points (Subscript[X, c]) that satisfy

∇h[(Subscript[X, c])]  =  0  or  ({{∂(h/∂)x[(Subscript[x, c]),(Subscript[y, c])]=0}, {∂(h/∂)y[(Subscript[x, c]),(Subscript[y, c])]=0}})

a) Explain why these points are equilibria of the gradient dynamical system above.

b) What is the solution path  X[t]  of the gradient dynamical system if the initial value of the path is   X[0] = (Subscript[X, c]), a critical value?

Hint: Attraction, Repulsion, and Indecision

Suppose (Subscript[X, c]) satisfies  ∇h[(Subscript[X, c])]  =  0.

a)  Suppose there is a small neighborhood of (Subscript[X, c]) so that every solution path of the gradient dynamical system that starts at a point in that neighborhood flows toward the point.  This happens for the following system near (1,1).

Clear[h, f, g, x, y, t];
h[x_, y_]  =  
        x y (E^(-(1/2) ((x^2) + (y^2))));
f[x_, y_]  =  Evaluate[D[h[x, y], x]];
g[x_, y_] =  Evaluate[D[h[x, y], y]];
initconds = 
        Flatten[Table[{x, y}, {x, - .2, 1.8,  .4}, 
            {y, - .2, 1.8,  .4}], 1];
tfinal = 10;
dt = 2;
Print["Flow solution of the differential equations:"];
Print["dx = ", f[x, y], " dt"];
Print["dy = ", g[x, y], " dt"];
flow2D[{f, g}, {x, y, t}, initconds, tfinal, dt, 
        AxesLabel -> {"x", "y"}, 
        PlotRange -> {{-1, 2}, {-1, 2}}];

Explain why this means (Subscript[X, c]) is the (x,y)-location of a maximum for  h[x,y].

b)  Suppose there is a small neighborhood of (Subscript[X, c]) so that every solution path of the gradient dynamical system that starts at a point in that neighborhood flows away from the point.  This happens for the following system near (1,-1).

Clear[h, f, g, x, y, t];
h[x_, y_]  =  
        x y (E^(-(1/2) ((x^2) + (y^2))));
f[x_, y_]  =  Evaluate[D[h[x, y], x]];
g[x_, y_] =  Evaluate[D[h[x, y], y]];
initconds = 
        Flatten[Table[{x, y}, {x,  .6, 1.4,  .2}, 
            {y, -1.4, - .6,  .2}], 1];
tfinal = 10;
dt = 2;
Print["Flow solution of the differential equations:"];
Print["dx = ", f[x, y], " dt"];
Print["dy = ", g[x, y], " dt"];
flow2D[{f, g}, {x, y, t}, initconds, tfinal, dt, 
        AxesLabel -> {"x", "y"}, 
        PlotRange -> {{-1, 2}, {-2, 1}}];

Explain why this means (Subscript[X, c]) is the (x,y)-location of a minimum for  h[x,y].

c)  Suppose that for every small neighborhood of (Subscript[X, c]) so that some solution paths of the gradient dynamical system that start at a point in that neighborhood flows away from the point, but other solution paths flow toward (Subscript[X, c]).  This happens for the following system near (0,0).

Clear[h, f, g, x, y, t];
h[x_, y_]  =  
        x y (E^(-(1/2) ((x^2) + (y^2))));
f[x_, y_]  =  Evaluate[D[h[x, y], x]];
g[x_, y_] =  Evaluate[D[h[x, y], y]];
initconds = 
        Flatten[Table[{x, y}, {x, - .4,  .4,  .2}, 
            {y, - .4,  .4,  .2}], 1];
tfinal = 10;
dt = 2;
Print["Flow solution of the differential equations:"];
Print["dx = ", f[x, y], " dt"];
Print["dy = ", g[x, y], " dt"];
flow2D[{f, g}, {x, y, t}, initconds, tfinal, dt, 
        AxesLabel -> {"x", "y"}, 
        PlotRange -> {{-1, 2}, {-2, 1}}];

Explain why this means (Subscript[X, c]) is NOT the (x,y)-location of a maximum nor a minimum for  h[x,y].

Hint 33.2 The Second Derivative Test

First, read Theorems 24.1, 2, 3, 4, Procedure 24.5 and Theorem  24.6 in the text.    Then review the basic types of linear systems.  The text Figure 24.9 shows the basic local flows.  When the characertistic values are both positive, the equilibrium is a repellor.  When the characteristic values are both negative, the equilibrium is an attractor.  If one characteristic root is positive and the other is negative, the equilibrium neither attracts nor repells.  Gradient dynamical systems do not have complex characteristic roots.

Work Hint 33.2 of the print project and type your solution here in your report.  Conclude this part of your report with a summary of cases like:

Condition for a Minimum: Two Positive Characteristic Roots

The critical point ((Subscript[x, c]), (Subscript[y, c])) satisfying  ∇h[(Subscript[x, c]), (Subscript[y, c])] = 0  is a local minimum of the height function and a repelling equilibrium of the gradient dynamical system provided:

(∂^2) (h/∂)(x^2)[(Subscript[x, c]), (Subscript[y, c])](∂^2) (h/∂)(y^2)[(Subscript[x, c]), (Subscript[y, c])]-((((∂^2) (h/(∂x ∂y))[(Subscript[x, c]), (Subscript[y, c])]))^2)>0

and

(∂^2) (h/∂)(x^2)[(Subscript[x, c]), (Subscript[y, c])]+(∂^2) (h/∂)(y^2)[(Subscript[x, c]), (Subscript[y, c])]>0

Condition for a Maximum: Two Negative Characteristic Roots

Condition for a Saddle: Opposite Sign Characteristic Roots

Condition for Failure of the Test

Programming the test

You could even program the results of the test by using If[.] commands.

?? If

For example, test various cases of the following program when  u  and  v  are numbers.

u = 1;
v = 0;
If[u == 0, "First zero", 
      If[u < 0, "First negative", 
        If[v == 0, "Second zero", 
          If[v < 0, "Second negative", "Second positive"]]]]

Finding Critical Points or Gradient Equilibria

Following is a Mathematica solution of the classification of critical points.

Define the Differential Equations

Clear[f, g, h, x, y];
h[x_, y_]  =  
        x y (E^(-(1/2) ((x^2) + (y^2))));
f  = Factor[D[h[x, y], x]]
g = Factor[ D[h[x, y], y]]

Find the Equilibria

Solve[{f == 0, g == 0}, {x, y}]

If the algebraic Solve[.] command does not work, you will have to use FindRoot[.] and experiment with initial guesses until you find all solutions.

FindRoot[{f == 0, g == 0}, {x,  .7}, {y,  .7}]

The Symbolic Matrix of Partials

aX  =  D[f,  x]; 
aY  =  D[f,  y]; 
bX  =  D[g,  x]; 
bY  =  D[g,  y]; 
Print[MatrixForm[{{aX,  aY},  {bX,  bY}}]]

The Specific Matrix at (0,0):

x = 0;
y = 0;
Print[MatrixForm[{{aX,  aY},  {bX,  bY}}]];
charPoly = (aX - r) (bY - r) - aY bX

The Second Partial Derivative Test Expressions at (0,0):

(∂^2) (h/∂)(x^2)[(Subscript[x, c]), (Subscript[y, c])](∂^2) (h/∂)(y^2)[(Subscript[x, c]), (Subscript[y, c])]-((((∂^2) (h/(∂x ∂y))[(Subscript[x, c]), (Subscript[y, c])]))^2)

x = 0;
y = 0;
Sign[aX bY  -  bX aY]

(∂^2) (h/∂)(x^2)[(Subscript[x, c]), (Subscript[y, c])]+(∂^2) (h/∂)(y^2)[(Subscript[x, c]), (Subscript[y, c])]

Sign[aX + bY]

The results of the Local Second Derivative Test are ___

Find the Characteristic Roots for the (0,0)-equilibrium

N[Solve[charPoly == 0, r]]

Roots of opposite sign mean ___???___

Plot3D[h[x + dx, y + dy], {dx, -1, 1}, {dy, -1, 1}, 
        AxesLabel -> {"x", "y", "h"}, Ticks -> None];

The Second Partial Derivative Test Expressions at (1,1):

(∂^2) (h/∂)(x^2)[(Subscript[x, c]), (Subscript[y, c])](∂^2) (h/∂)(y^2)[(Subscript[x, c]), (Subscript[y, c])]-((((∂^2) (h/(∂x ∂y))[(Subscript[x, c]), (Subscript[y, c])]))^2)

x = 1;
y = 1;
Sign[aX bY  -  bX aY]

(∂^2) (h/∂)(x^2)[(Subscript[x, c]), (Subscript[y, c])]+(∂^2) (h/∂)(y^2)[(Subscript[x, c]), (Subscript[y, c])]

Sign[aX + bY]

The results of the Local Second Derivative Test are ___

Characteristic roots at (1,1):

The characteristic equation and roots at another equilibrium are:

x = 1;
y = 1;
Print[MatrixForm[{{aX,  aY},  {bX,  bY}}]];
charPoly = (aX - r) (bY - r) - aY bX;
Print[charPoly, " = 0"];
N[Solve[charPoly == 0, r]]

Both roots negative means __???__

Plot3D[h[x + dx, y + dy], {dx, -1, 1}, {dy, -1, 1}, 
        AxesLabel -> {"x", "y", "h"}, Ticks -> None];

The Second Partial Derivative Test Expressions at (1,-1):

(∂^2) (h/∂)(x^2)[(Subscript[x, c]), (Subscript[y, c])](∂^2) (h/∂)(y^2)[(Subscript[x, c]), (Subscript[y, c])]-((((∂^2) (h/(∂x ∂y))[(Subscript[x, c]), (Subscript[y, c])]))^2)

x = 1;
y = -1;
Sign[aX bY  -  bX aY]

(∂^2) (h/∂)(x^2)[(Subscript[x, c]), (Subscript[y, c])]+(∂^2) (h/∂)(y^2)[(Subscript[x, c]), (Subscript[y, c])]

Sign[aX + bY]

The results of the Local Second Derivative Test are ___

Characteristic roots at (1,-1):

The characteristic equation and roots at another equilibrium are:

x = 1;
y = -1;
Print[MatrixForm[{{aX,  aY},  {bX,  bY}}]];
charPoly = (aX - r) (bY - r) - aY bX;
Print[charPoly, " = 0"];
N[Solve[charPoly == 0, r]]

Both roots positive means __???___

Plot3D[h[x + dx, y + dy], {dx, -1, 1}, {dy, -1, 1}, 
        AxesLabel -> {"x", "y", "h"}, Ticks -> None];

The Second Partial Derivative Test Expressions at (-1,1):

(∂^2) (h/∂)(x^2)[(Subscript[x, c]), (Subscript[y, c])](∂^2) (h/∂)(y^2)[(Subscript[x, c]), (Subscript[y, c])]-((((∂^2) (h/(∂x ∂y))[(Subscript[x, c]), (Subscript[y, c])]))^2)

x = -1;
y = 1;
Sign[aX bY  -  bX aY]

(∂^2) (h/∂)(x^2)[(Subscript[x, c]), (Subscript[y, c])]+(∂^2) (h/∂)(y^2)[(Subscript[x, c]), (Subscript[y, c])]

Sign[aX + bY]

The results of the Local Second Derivative Test are ___

Characteristic roots at (-1,1):

The characteristic equation and roots at another equilibrium are:

x = -1;
y = 1;
Print[MatrixForm[{{aX,  aY},  {bX,  bY}}]];
charPoly = (aX - r) (bY - r) - aY bX;
Print[charPoly, " = 0"];
N[Solve[charPoly == 0, r]]

Both roots positive means __???___

Plot3D[h[x + dx, y + dy], {dx, -1, 1}, {dy, -1, 1}, 
        AxesLabel -> {"x", "y", "h"}, Ticks -> None];

The Second Partial Derivative Test Expressions at (-1,-1):

(∂^2) (h/∂)(x^2)[(Subscript[x, c]), (Subscript[y, c])](∂^2) (h/∂)(y^2)[(Subscript[x, c]), (Subscript[y, c])]-((((∂^2) (h/(∂x ∂y))[(Subscript[x, c]), (Subscript[y, c])]))^2)

x = -1;
y = -1;
Sign[aX bY  -  bX aY]

(∂^2) (h/∂)(x^2)[(Subscript[x, c]), (Subscript[y, c])]+(∂^2) (h/∂)(y^2)[(Subscript[x, c]), (Subscript[y, c])]

Sign[aX + bY]

The results of the Local Second Derivative Test are ___

Characteristic roots at (-1,-1):

The characteristic equation and roots at another equilibrium are:

x = -1;
y = -1;
Print[MatrixForm[{{aX,  aY},  {bX,  bY}}]];
charPoly = (aX - r) (bY - r) - aY bX;
Print[charPoly, " = 0"];
N[Solve[charPoly == 0, r]]

Plot3D[h[x + dx, y + dy], {dx, -1, 1}, {dy, -1, 1}, 
        AxesLabel -> {"x", "y", "h"}, Ticks -> None];

Hint 33.4 Local Max-min

Solve hint 33.4 of the print project.  Specifically, find and classify all the critical points of the function  h[x,y] = 3 ((x^2) + 3 (y^2)) (e^(-((x^2) + (y^2))))

Clear[x, y, h];
h[x_, y_] = 
        4 ((x^2) + 3 (y^2)) (E^(-((x^2) + (y^2))));
      
Plot3D[h[x, y], {x, -2, 2}, {y, -2, 2}, 
        AxesLabel -> {"x", "y", "h"}, PlotPoints -> 40, 
        Mesh -> False, ViewPoint -> {5, 1, 2}];