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 41: Cooperation & Competition Between Species

by A. Student and B. Pupil

with technical help from Da LabMonitor and Help Notebook

Special Mathematica Functions: flow2D[.] & vectorFlow2D[.]

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];	
];

BeginPackage["Graphics`flow2D`"]
flow2D::usage =
"flow2D[{f,g},{x,y},initialValues,tFinal,dt,(options)]
	shows the flow solution of dx = f[x,y] dt, dy = g[x,y] dt 
	starting with the initial (x,y)-values = initialValues = {{a,b},{c,d},...}
	flowing up to time tFinal with animation frames every dt time units."

Begin["`Private`"]
stuff goes here 
End[]
EndPackage[]

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

The NoteBook Flow2D.nb contains a step-by-step explanation of the Mathematica solution of thecompetition equations

(dx/dt)  =  x (1 - (x/7) - (y/7))

(dy/dt)  =  y (1 - (x/10) - (y/5))

Then it gives the following computation where all the steps are combined in one program.

The Flow of Text Example 22.4: Fierce Competition

This computation takes a long time, because there are a lot of initial conditions and a lot of time steps.

Clear[f, g, x, y, t];
f[x_, y_] := x (1 - (x/10) - (y/5));
g[x_, y_] := y (1 - (x/7) - (y/7));
initconds = 
        Flatten[Table[{x, y}, {x, - .1, 9.9, 2}, {y, - .1, 9.9, 2}], 
          1];
tfinal = 10;
dt = 1;
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, 10}, {-1, 10}}];

To show the flow as a real time animation, highlight the cell containing the graphs and press the <Command> + <y>  keys at the same time.

Notice that the species  x  and  y  do not coexist.  One or the other ends up extinct in the long term.  Explain why the flow above shows this.

Practice Exercise

Modify the computation above to solve the weak competition equations:

(dx/dt)  =  x (1 - (x/7) - (y/7))

(dy/dt)  =  y (1 - (x/10) - (y/5))

Notice the different limiting behavior?  Why does this flow show that these species coexist?

Carrying Capacity and The Meaning of the Equations

See HINT 41.3 and explain the biological meaning of the equations

(dx/dt)  =  x (1 - (x/h) - b  (y/k))

(dy/dt)  =  y (1 - c  (x/h) - (y/k))

Mathematica Experiments (HINT 41.4)

Clear[f, g, x, y, t];
f[x_, y_] := x (1 - (x/h) - b  (y/k));
g[x_, y_] := y (1 - c  (x/h) - (y/k));
h = 7;
k = 10;
b = 2/3;
c = 3/4;
initconds = 
        Flatten[Table[{x, y}, {x, - .1, 9.9, 2}, {y, - .1, 9.9, 2}], 
          1];
tfinal = 10;
dt = 1;
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, 10}, {-1, 10}}];

Mathematica Experiments (HINT 41.4 - alternative formulation)

It simplifies the mathematical analysis if we change variables to

u = (x/h)

v = (y/k)

so that the carrying capacities are taken out of the equations.  One unit of  u  is worth  one unit of  v, even if  x is elephants and  y  is mice.

Clear[f, g, u, v, t];
f[u_, v_] := u (1 - u - b v);
g[u_, v_] := v (1 - c u - v);
b = 2/3;
c = 3/4;
initconds = 
        Flatten[Table[{u, v}, {u, - .1, 1.9,  .4}, 
            {v, - .1, 1.9,  .4}], 1];
tfinal = 10;
dt = 1;
Print["Flow solution of the differential equations:"];
Print["du = ", f[u, v], " dt"];
Print["dv = ", g[u, v], " dt"];
flow2D[{f, g}, {u, v, t}, initconds, tfinal, dt, 
        AxesLabel -> {"u", "v"}, 
        PlotRange -> {{- .1, 2}, {- .1, 2}}];

Verifying the Conjecture on the Role of  b  and  c

Once you have experimented a little and can see what the roles of  b  and  c  are in competition, prove your conjecture as in Hint 41.6

Hint 41.6

Hint 41.7

Read the print project hint and compare the text approach from the following excerpt from the text NoteBook LocStability.nb

A Symbolic Computation: Competing Species

The computations throughout this section are done symbolically.  In the Ecology project on competing species you can show geometrically that stability of a positive equilibrium for both of the two "species"  x  and  y  of the equations below depends on the competition terms  b  and  c  being less than one.  (The parameters  h  and  k  are the "carrying capacities" of the environment for the species when they are alone in the environment.)
 	
	dx/dt = m x(1 - x/h - b y/k);
	dy/dt = n y (1 - c x/h - y/k);
			
	The parameter  b  gives a competition multiplier for  y  relative to  y's  carrying capacity in competition with  x  in the  x-equation.  Likewise,  c  gives the relative competitiveness of  x  in the  y-equation.
	
		The stability for this problem can be done as follows, but it is quite a lot easier to solve geometrically, as you can see in the competing species project.

The Positive Equilibrium Point

This computation will find the equilibrium points to be used throughout the computations.

Clear[x, y, h, k, b, c];
equns = {1 == (x/h) + ((b y)/k), 1 == (y/k) + ((c x)/h)}

equil = Solve[equns, {x, y}]

equil = Simplify[equil]

Question

When are both of these points positive? (The parameters are all positive.)

The Matrix of Partial Derivatives

This step will give the matrix of the partial derivatives symbolically.

f = m x(1 - x/h - b y/k);
g = n y(1 - c x/h - y/k);
aX = D[f,x];
aY = D[f,y];
bX = D[g,x];
bY = D[g,y];
matrixOfPartials = {{aX,aY},{bX,bY}};
Print[TableForm[matrixOfPartials]]

The Matrix at the Positive Equilibrium

This step will give the matrix of the partial derivatives at the interior equilibrium point:

x = (h - b h)/(1 - b c);
y = (k - c k)/(1 - b c);
Print[TableForm[Simplify[{{aX,aY},{bX,bY}}]]]

Notice the common factor in the denominator of all these terms

Multiply this matrix by the common factor

TableForm[Simplify[(???){{aX,aY},{bX,bY}}]]

The Characteristic Polynomial at Equilibrium

The next step gives the characteristic equation for the equilibrium.  We can simplify it somewhat by multiplying by the common term in all the denominators.

charPoly = ((aX - r)(bY - r) - aY bX)
simplePoly = Collect[Collect[Together[(b c - 1)charPoly],r]/(b c -1),r]

The General Characteristic Roots

This step will give the characteristic roots for the given equilibrium point which can be used to decide what type of equilibrium point it is.

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

When are the characteristic roots opposite signs?

Roots of a quadratic equation

 r^2 + B r + C == 0
 
have the same sign if  4 C  is positive and the discriminant is positive, then Sqrt[B^2 - 4 C] < B.  Our equation in terms of the partial derivative coefficients is

(aX - r)(bY - r) - bX aY = r^2 - (aX + bY) + (aX bY - bX aY) == 0

so in terms of the collected parameters above,

	B = -(aX + bY)
	C = (aX bY - bX aY)

C > 0?

Solve[aX bY == aY bX,{b,c}]

Simplify[(aX bY - aY bX)(1 - b c)^2/(m n)]

When is the product (1 - b)(1 - c)(1 - b c) > 0?
Why is C > 0  or  (aX bY - aY bX) > 0  under the same circumstances?

Solve[(aX + bY)^2 == 4(aX bY - bX aY),{b,c}]
Solve[(m+n)^2 == 4 m n(1-b c),{b,c}]

Notice that the simplified solution gives
	
		b c = -(m - n)^2/(4 m n)
		
and recall that  b  and  c  are positive fractional competition rates.

Cooperation between Species

Mathematica Experiments (HINT 41.9)

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

Analytical Computation of The Equilibrium

Analytical Computation of Stability

The Biological Meaning of  f > 1  or  g > 1