7.6 Analyzing the Phase Plane

Mathematica script by Chris Parrish,
cparrish@sewanee.edu

Sources and references for some of these problems include
James Stewart, "Calculus: Concepts and Contexts," Second Edition, Brooks/Cole, 2001
Deborah Hughes-Hallett, Andrew M. Gleason, et. al., "Calculus," Second Edition, John Wiley & Sons, 1998
Robert Fraga, ed., "Calculus Problems for a New Century," The Mathematical Association of America, 1993

Citrus Tree Parasites

Hughes-Hallett, Gleason, et al, Text Example 9.9, pages 553ff

Solve the first equation.

In[567]:=

Clear[p1,p2,null1,null2,c,d]

Solve[50 - 2.5 p1 - 2 p2 == 0,p2]

Out[568]=

Take null1 to be this function.

In[569]:=

null1[p1_] := 0.5 (50 - 2.5 p1)

Now solve the second equation.

In[570]:=

Solve[90 - 150 p1 - 6 p2 == 0,p2]

Out[570]=

Take null2 to be that function.

In[571]:=

null2[p1_] := 5 (3 - 5 p1)

Plot both of these functions to see the regions determined by the null clines.

In[572]:=

c = 0; d = 20;    (* c <= p1 <= d *)

Plot[{null1[p1],null2[p1]},{p1,c,d},
PlotStyle->Red,
PlotRange->{0,25},
PlotLabel->"Null Clines",
AxesLabel->{"p1","p2"}];

What happens in the regions in between the null clines?
Just test for the signs of these derivatives at a few telling points.

In[574]:=

Clear[dp1,dp2,signs,point1,point2,point3]

dp1[p1_,p2_] := 0.001 p1 (50 - 2.5 p1 - 2 p2)
dp2[p1_,p2_] := 0.001 p2 (90 - 150 p1 - 6 p2)

signs[{a_,b_}] := Sign[ {dp1[a,b],dp2[a,b]} ]

point1 = {.1,.1};
point2 = {1,1};
point3 = {20,25};

signs[point1]
signs[point2]
signs[point3]

Out[581]=

Out[582]=

Out[583]=

Let's try to plot a bit of the corresponding phase plane together with a segment of a null cline.

In[584]:=

In[586]:=

Clear[f,p1,p2,a,b,c,d,field,cline,nullcline]

f[p1_,p2_] := p2 (90 - 150 p1 - 6 p2)/(p1 (50 - 2.5 p1 - 2 p2));

a = 0.1; b = 2.1;    (* a <= p1 <= b *)
c = 0.1; d = 2.1;    (* c <= p2 <= d *)

field = PlotVectorField[{1,f[p1,p2]},
{p1,a,b},{p2,c,d},
PlotLabel -> "dy/dx = dp2/dp1",
Axes -> True,
AxesLabel -> {"p1","p2"},
PlotPoints -> 20,
PlotRange->{c,d},
Prolog -> ManganeseBlue,
DisplayFunction->Identity];

cline[p1_] := (90 - 150 p1)/6

nullcline = Plot[cline[p1],{p1,a,b},
PlotStyle->Red,
PlotRange->{c,d},
DisplayFunction->Identity];

Show[{field,nullcline},
DisplayFunction->\$DisplayFunction];

Stretch out the above picture to (W: 504, H: 451) in order to make the arrows come into focus.
Compare the resulting drawing with Figure 9.60, page 555, of Hughes-Hallett, Gleason, et al.

US-Soviet Arms Race, 1945-1960

Hughes-Hallett, Gleason, et al, Exercise 9.9.10, pages 558-9

In[597]:=

Clear[dus,dsov,sov,us,f,a,b,c,d,field,cline,nullcline]

dsov[sov_,us_] := -0.45 sov          + 10.5;
dus[sov_,us_]  :=  8.2  sov - 0.8 us - 142;

f[sov_,us_] := dus[sov,us]/dsov[sov,us];

a = 22; b = 24;    (* a <= sov <= b *)
c = 61; d = 63;    (* c <= us <= d *)

Solve for the nullclines.

In[603]:=

null1[us_] := 10.5/0.45
null2[us_] := (0.8 us + 142)/8.2

Plot these functions to see the regions determined by the null clines.

In[605]:=

Clear[null1,null2,us,sov,a,b,c,d]

null2[sov_] := (8.2 sov - 142)/0.8

a = 22; b = 24;    (* a <= sov <= b *)
c = 61; d = 63;    (* c <= us <= d *)

Plot[null2[sov],{sov,a,b},
PlotStyle->Red,
PlotRange->{c,d},
AspectRatio->1,
Epilog->{Red, Line[{{23.33,c},{23.33,d}}]},
PlotLabel->"Null Clines",
AxesLabel->{"us","sov"}];

What happens in the regions in between the null clines?
Just test for the signs of these derivatives at a few telling points.

In[610]:=

Clear[dp1,dp2,signs,point1,point2,point3,point3]

dsov[sov_,us_] := -0.45 sov          + 10.5;
dus[sov_,us_]  :=  8.2  sov - 0.8 us - 142;

signs[{a_,b_}] := Sign[ {dsov[a,b],dus[a,b]} ]

point1 = {23.1,61.66};  (* "west" *)
point2 = {23.3,61};     (* "south" *)
point3 = {23.5,61.66};  (* "east" *)
point4 = {23.4,63};     (* "north" *)

Print["west:",signs[point1],",  south:", signs[point2],",  east:",signs[point3],",  north:", signs[point4]]

The first number in each one of the four pairs of numbers indicates whether the Soviets will increase or decrease their arms spending. The second number indicates whether the US will increase or decrease its arms spending.

Let's try to plot a bit of the corresponding phase plane together with segments of nullclines.

In[619]:=

Clear[dus,dsov,sov,us,f,a,b,c,d,null1,null2,field,cline,nullcline]

dsov[sov_,us_] := -0.45 sov          + 10.5;
dus[sov_,us_]  :=  8.2  sov - 0.8 us - 142;

f[sov_,us_] := dus[sov,us]/dsov[sov,us];

a = 22; b = 24;    (* a <= sov <= b *)
c = 61; d = 63;    (* c <= us <= d *)

null2[sov_] := (8.2 sov - 142)/0.8

nullcline = Plot[null2[sov],{sov,a,b},
PlotStyle->Red,
PlotRange->{c,d},
AspectRatio->1,
DisplayFunction->Identity];

field = PlotVectorField[{dsov[sov,us],dus[sov,us]},
{sov,a,b},{us,c,d},
PlotLabel -> "dus/dsov = dus[sov,us]/dsov[sov,us]",
Axes -> True,
AxesLabel -> {"sov","us"},
PlotPoints -> 20,
PlotRange->{c,d},
Prolog -> ManganeseBlue,
Epilog->{Red, Line[{{23.33,c},{23.33,d}}]},
DisplayFunction->Identity];

Show[{field,nullcline},
DisplayFunction->\$DisplayFunction];

Created by Mathematica  (April 25, 2004)