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 39: Voodoo Eats the Mice
(Logistic Growth with Hunting)

by A. Student and B. Pupil

with technical help from Da LabMonitor and Help Notebook

Begin by selecting "Evaluate Notebook" from Kernel/Evaluation on the menu bar at the Mathematica top.  (Say "yes" to evaluate initialization.)

Special Mathematica Functions: flow2D[.] & vectorFlow2D[.] (copy this into your notebook if different, don't look inside -THIS IS PROGRAMMING!)

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 = (.15 Norm[ntip-ntail]/Norm[ntip-ntail]) (ntip-ntail) ;
	up =  {-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]]
				];
				
vectorFlow1D[f_,{y_,t_},initVals_,tFinal_,dt_,opns___] :=
 Module[{soln,arr,pt,pt0,ptt,seg,crv,tt,x},
soln =Flatten[Map[
	NDSolve[{x[0]==0,y[0]==#,
	x'[t] == 1,y'[t] == f[y[t]]},
	{x[t],y[t]},{t,0,tFinal}]&,
										initVals],1];
											
pt[tt_] := Map[Graphics[{Blue,AbsolutePointSize[5],Point[{tFinal,#}]}]&,
y[t]/.soln/.t->tt];

pt0[tt_] := Map[Graphics[{Blue,AbsolutePointSize[5],Point[{0,#}]}]&,
y[t]/.soln/.t->tt];

ptt[tt_] := Map[Graphics[{Blue,AbsolutePointSize[4],Point[{#[[1]],#[[2]]}]}]&,
{x[t],y[t]}/.soln/.t->tt];

seg[tt_] := Map[Graphics[{Green,Line[{{#[[1]],#[[2]]},{tFinal,#[[2]]}}]}]&,
{x[t],y[t]}/.soln/.t->tt];

	arr[tt_] := Map[uArrow[{0,f[#]},Tail->{tFinal,#},Red]&,
	y[t]/.soln/.t->tt];
										
crv[tt_] := Plot[Evaluate[y[t]/.soln],{t,0,tt},opns,
			AspectRatio->1,DisplayFunction->Identity];
			
Show[{crv[dt/100],pt[0],pt0[0],seg[0],arr[0]},opns,
PlotRange->{{0,10},{0,10}},AspectRatio->1,
DisplayFunction->$DisplayFunction];
Do[Show[{crv[tt],pt[tt],ptt[tt],seg[tt],arr[tt]},DisplayFunction->$DisplayFunction];,
{tt,dt,tFinal,dt}];
	Return[soln];	
];

flow1D[f_,{y_,t_},initVals_,tFinal_,dt_,opns___] :=
 Module[{soln,ptt,pt,pt0,seg,crv,tt,x},
soln =Flatten[Map[
	NDSolve[{x[0]==0,y[0]==#,
	x'[t] == 1,y'[t] == f[y[t]]},
	{x[t],y[t]},{t,0,tFinal}]&,
										initVals],1];
											
pt[tt_] := Map[Graphics[{Blue,AbsolutePointSize[5],Point[{tFinal,#}]}]&,
y[t]/.soln/.t->tt];

pt0[tt_] := Map[Graphics[{Blue,AbsolutePointSize[5],Point[{0,#}]}]&,
y[t]/.soln/.t->tt];

ptt[tt_] := Map[Graphics[{Blue,AbsolutePointSize[4],Point[{#[[1]],#[[2]]}]}]&,
{x[t],y[t]}/.soln/.t->tt];

seg[tt_] := Map[Graphics[{Green,Line[{{#[[1]],#[[2]]},{tFinal,#[[2]]}}]}]&,
{x[t],y[t]}/.soln/.t->tt];
;
										
crv[tt_] := Plot[Evaluate[y[t]/.soln],{t,0,tt},opns,
			AspectRatio->1,DisplayFunction->Identity];
			
Show[{crv[dt/100],pt[0],pt0[0],seg[0]},opns,
PlotRange->{{0,10},{0,10}},AspectRatio->1,
DisplayFunction->$DisplayFunction];
Do[Show[{crv[tt],pt[tt],ptt[tt],seg[tt]},DisplayFunction->$DisplayFunction];,
{tt,dt,tFinal,dt}];
	Return[soln];	
];

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.

Exponential Growth

Constant per capita fertility ...

(dx/dt) =  r x

With initial population  (Subscript[x, 0]),

x[0] = (Subscript[x, 0])
(dx/dt) =  r x	⟺	x[t] = (Subscript[x, 0]) (e^(r t))

Why solve this with Mathematica?  We know it Knee Jerk!

Finding the rate from two values

Clear[x, x0, r, t];
eMice[t_] = 
        x0 (E^(r t));
Print["Mice as a function of t in general = ", 
        eMice[t]];
x0 = 1000;
Print["Mice as a function of t in with 1000 initial = ", 
        eMice[t]];
Solve[eMice[1/3] == 1100, r]
NSolve[eMice[1/3] == 1100, r]

A week's worth:

Clear[x, x0, r, t];
eMice[t_] = 
        x0 (E^(r t));
Print["Mice as a function of t in general = ", 
        eMice[t]];
x0 = 1000;
r = 0.285931;
Print["Mice as a function of t in with these parameters = ", 
        eMice[t]];
Plot[eMice[t], {t, 0, 7}, 
        PlotStyle -> {AbsoluteThickness[3], Red}, 
        PlotLabel -> "Mice for a week."];

Mice for a Year

eMice[365]

How big is this? Say compared with the mass of the earth or a volume one square mile on the base and how high???

Logistic Growth

The constant fertility (exponential) model predicts unrealistic growth ... so we introduce a cut-off term  (1-x/c) in the growth model.

(dx/dt) = r x (1 - (x/c))

bla, bla, bla.  We graphed this in class and found that the crucial fact was to find out where  (dx/dt) = 0, (dx/dt) < 0,(dx/dt) > 0.  We can be silly and solve this equation with Mathematica,

Clear[x, r, c, h];
growthRate[x_] = r x (1 - (x/c));
Solve[growthRate[x] == 0, x]

Solution of the Differential Equation with Mathematica

Clear[x, r, c, h];
growthRate[x_] = r x (1 - (x/c));
x0 = 1000;
r = 0.285931;
c = 2500;
tFinal = 90;
solution = 
        NDSolve[{x'[t] == growthRate[x[t]], x[0] == x0}, 
          x[t], {t, 0, tFinal}];
x[t_] = x[t] /. solution[[1]];
plot1 = 
        Plot[x[t], {t, 0, tFinal}, PlotRange -> {0, 5000}, 
          PlotStyle -> {AbsoluteThickness[3], Red}, 
          PlotLabel -> "Logistic mice for a summer."];

Clear[x, t];
x0 = 5000; 
    solution = 
      NDSolve[{x'[t] == growthRate[x[t]], x[0] == x0}, 
        x[t], {t, 0, tFinal}];
x[t_] = x[t] /. solution[[1]];
plot2 = 
        Plot[x[t], {t, 0, tFinal}, PlotRange -> {0, 5000}, 
          PlotStyle -> {AbsoluteThickness[3], Blue}, 
          PlotLabel -> "Logistic mice for a summer."];

Show[plot1, plot2];

Lots of initial conditions for  (dx/dt) = f[x]

Thanks to the flow1D[.] package hidden in the initialization, we can animate many initial populations at once.

Clear[t, r, c, h, x, f];
growthRate[x_] := r x (1 - x/c);
r = 0.285;
c = 2500;
initXs = Table[x, {x, 0, 5000, 1000}];
tFinal = 5;
dt = 0.25;
flow1D[growthRate, {x, t}, initXs, tFinal, dt, 
        PlotRange -> {{0, tFinal + dt}, {0, 5000}}];

Linear Hunting

The constant fertility (exponential) model predicts unrealistic growth ... so we introduce a cut-off term  (1-x/c) in the growth model.

(dx/dt) = r x (1 - (x/c)) - h x

bla, bla, bla.  We graphed this in class and found that the crucial fact was to find out where  (dx/dt) = 0, (dx/dt) < 0,(dx/dt) > 0.  We can be silly and solve this equation with Mathematica,

Clear[x, r, c, h];
growthRate[x_] = 
        r x (1 - (x/c)) - h x;
Solve[growthRate[x] == 0, x]

Solution of the Differential Equation with Mathematica

Clear[x, r, c, h];
growthRate[x_] = 
        r x (1 - (x/c)) - h x;
x0 = 1000;
r = 0.285931;
c = 2500;
h = r/5;
xE = c (1 - h/r);
Print["When h = ", 
        h, " the positive equilibrium population, xE = ", 
        xE];
tFinal = 90;
solution = 
        NDSolve[{x'[t] == growthRate[x[t]], x[0] == x0}, 
          x[t], {t, 0, tFinal}];
x[t_] = x[t] /. solution[[1]];
plot1 = 
        Plot[x[t], {t, 0, tFinal}, PlotRange -> {0, 5000}, 
          PlotStyle -> {AbsoluteThickness[3], Red}, 
          PlotLabel -> "Hunted mice for a summer."];

Clear[x, t];
x0 = 5000;
solution = 
        NDSolve[{x'[t] == growthRate[x[t]], x[0] == x0}, 
          x[t], {t, 0, tFinal}];
x[t_] = x[t] /. solution[[1]];
plot2 = 
        Plot[x[t], {t, 0, tFinal}, PlotRange -> {0, 5000}, 
          PlotStyle -> {AbsoluteThickness[3], Blue}, 
          PlotLabel -> "Hunted mice for a summer."];

Show[plot1, plot2];

Lots of initial conditions for  (dx/dt) = f[x]

Thanks to the flow1D[.] package hidden in the initialization, we can animate many initial populations at once.

Clear[t, r, c, h, x, growthRate];
growthRate[x_] = r x (1 - (x/c)) - h x;
x0 = 1000;
r = 0.285931;
c = 2500;
h = r/2;
xE = c (1 - h/r);
initXs = Table[x, {x, 0, 5000, 1000}];
tFinal = 5;
dt = 0.25;
Print["When h = ", 
        h, " the positive equilibrium population, xE = ", 
        xE];
flow1D[growthRate, {x, t}, initXs, tFinal, dt, 
        PlotRange -> {{0, tFinal + dt}, {0, 5000}}];

Another expression for growth with hunting

(dx/dt) = r x (1  -  (x/c))  -  
          h x = (r - h) x  - (r/c) (x^2)

Can be written

(dx/dt)=(r - h) x(
              1  - (x/((1 - h/r) c))) ...   or something like that what does this tell you about h = r or h > r???

Non-Linear Hunting

The constant fertility (exponential) model predicts unrealistic growth and the linear hunting model lets Voodoo eat ALL the mice, so we try nonlinear hunting where Vodoo is MORE successful at LARGE populations, but less successful at small ones.

(dx/dt) = r x (1 - (x/c)) - h (x^2)

bla, bla, bla.  We graphed this in class and found that the crucial fact was to find out where  (dx/dt) = 0, (dx/dt) < 0,(dx/dt) > 0.  We can be silly and solve this equation with Mathematica,

Clear[x, r, c, h];
growthRate[x_] = 
        r x (1 - (x/c)) - h (x^2);
Solve[growthRate[x] == 0, x]

It might help you understand if you re-write the growth law in terms of the new equilibrium point.

Simplify[(c/(1 + ((h c)/r)))] == 
      Simplify[((c r)/(c h + r))]

Why and how can we re-write the nonlinear hunting model as a new logistic equation???  What does this tell us about extinction?  ...

Solution of the Differential Equation with Mathematica

Clear[x, r, c, h, xE];
xE = (c/(1 + ((h c)/r)));
growthRate[x_] = 
        r x (1 - (1/xE) x);
x0 = 1000;
r = 0.285931;
c = 2500;
h = r/3000;
Print["When h = ", 
        h, " the positive equilibrium population, xE = ", 
        xE];
tFinal = 90;
solution = 
        NDSolve[{x'[t] == growthRate[x[t]], x[0] == x0}, 
          x[t], {t, 0, tFinal}];
x[t_] = x[t] /. solution[[1]];
plot1 = 
        Plot[x[t], {t, 0, tFinal}, PlotRange -> {0, 5000}, 
          PlotStyle -> {AbsoluteThickness[3], Red}, 
          PlotLabel -> "Hunted mice for a summer."];

Clear[x, t];
x0 = 5000;
solution = 
        NDSolve[{x'[t] == growthRate[x[t]], x[0] == x0}, 
          x[t], {t, 0, tFinal}];
x[t_] = x[t] /. solution[[1]];
plot2 = 
        Plot[x[t], {t, 0, tFinal}, PlotRange -> {0, 5000}, 
          PlotStyle -> {AbsoluteThickness[3], Blue}, 
          PlotLabel -> "Hunted mice for a summer."];

Show[plot1, plot2];

Lots of initial conditions for  (dx/dt) = f[x]

Thanks to the flow1D[.] package hidden in the initialization, we can animate many initial populations at once.

Clear[x, r, c, h, xE, growthRate];
xE = (c/(1 + ((h c)/r)));
growthRate[x_] = r x (1 - (1/xE) x);
x0 = 1000;
r = 0.285931;
c = 2500;
h = r/2000;
initXs = Table[x, {x, 0, 5000, 1000}];
tFinal = 5;
dt = 0.25;
Print["When h = ", 
        h, " the positive equilibrium population, xE = ", 
        xE];
flow1D[growthRate, {x, t}, initXs, tFinal, dt, 
        PlotRange -> {{0, tFinal + dt}, {0, 5000}}];