7.4 Models of Population Growth

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

Population of Switzerland

Hughes-Hallett, Gleason, et al, Exercise 9.7.1, page 536

In[1083]:=

<<Graphics`Colors` <<Graphics`PlotField`

In[1085]:=

Clear[p, p0, t] <br /> RowBox[{RowBox[{p0, =, 6.6}], ;}] RowBox[{RowBox[{RowBox[{r, =, 0.002}] ... d", AxesLabel {"t (yrs since 1988)", "population (millions)"}] ;

[Graphics:HTMLFiles/7.4_populationGrowth_3.gif]

Spread of Information by Mass Media

Hughes-Hallett, Gleason, et al, Exercise 9.7.5a, page 537

Let's choose some notation:
  population size  = m
  number of informed people = i
  number of people who know initially = i_0
Then, di/dt = k (m - i).

In[1090]:=

?DSolve

DSolve[eqn, y, x] solves a differential equation for the function y, with independent variable ... tial equations. DSolve[eqn, y, {x1, x2, ... }] solves a partial differential equation. More…

In[1091]:=

Clear[i,t,k,m]

DSolve[ i'[t] == k (m - i[t]), i[t],t]

Out[1092]=

{{i[t] m + ^(-k t) C[1]}}

We know that i[0] = i_0, so that
   i_0 = i[0] = m + c;
thus,
  c = i_0 - m.
Choose units so that m = 1 and k = 1 and i_0 = 0.

In[1093]:=

m = 1 ; k = 1 ; i0 = 0 ; <br /> i[t_] := m - (m - i0) ^(-k t)  Plot[i[t], {t,  ... y Mass Media", AxesLabel {"t", "proportion informed"}] ;

[Graphics:HTMLFiles/7.4_populationGrowth_7.gif]

Spread of Information by Word of Mouth

Hughes-Hallett, Gleason, et al, Exercise 9.7.5b, page 537

Let's choose some notation:
  population size  = m
  number of informed people = i
  number of people who know initially = i_0
Then, di/dt = k i (m - i).

In[1098]:=

Clear[i,t,k,m]

DSolve[ i'[t] == k i[t] (m - i[t]), i[t],t]

Out[1099]=

{{i[t]  (^(k m t) m)/(^(k m t) - ^(m C[1]))}}

Here is a large-scale view of the slope field.
Take k = 1.

In[1100]:=

Clear[f,i,t,k,m,a,b,c,d];

f[i_,t_] := i (m - i);

m = 1;
a = 0; b = 1;    (* a <= t <= b *)
c = 0; d = 1;    (* c <= i <= d *)

init1 = 0.0;
init2 = 0.05;
init3 = 0.75;

pts = {{0,init1},{0,init2},{0,init3}};

field = PlotVectorField[{1,f[i,t]},
                        {t,a,b},{i,c,d},
                        PlotLabel -> "Spread of Information by Word of Mouth",
                        Axes -> True,
                        AxesLabel -> {"t","i"},
                        PlotPoints -> 20,
                        Prolog -> ManganeseBlue,
                        Epilog -> {Red,PointSize[0.02],
                                   Map[Point,pts]}];

[Graphics:HTMLFiles/7.4_populationGrowth_9.gif]

Let's plot the solutions which pass through the red dots.
Since m = 1, we have
  i[t] = 1/(1 + c Exp[-kt])
Hence, c = (1 / i[0]  - 1)

In[1110]:=

Clear[i1,i2,i3,t]
m = 1;
k = 1;

c1 = 1/init1 - 1;
c2 = 1/init2 - 1;
c3 = 1/init3 - 1;

i1[t_] := 1/(1 + c1 Exp[-k t])
i2[t_] := 1/(1 + c2 Exp[-k t])
i3[t_] := 1/(1 + c3 Exp[-k t])

plots = Plot[{i1[t],i2[t],i3[t]},{t,a,b},
             PlotStyle->Indigo,
             PlotRange->{0,1},
             DisplayFunction->Identity];
             
Show[{field,plots},
     DisplayFunction->$DisplayFunction];

                                      1 Power :: infy : Infinite expression  --- encountered. More…                                      0.`

[Graphics:HTMLFiles/7.4_populationGrowth_11.gif]

In[1122]:=

Clear[i, t, m, i0, k] ; <br /> i[t_] := m/(1 + ((m - i0)/i0) ^(-k m t))

When is the information spreading the most rapidly?
di/dt is the rate of spread of information

In[1124]:=

i'[t]

Out[1124]=

(^(-k m t) k m^2 (-i0 + m))/(i0 (1 + (^(-k m t) (-i0 + m))/i0)^2)

When is this function a maximum? -- when its derivative is zero.

In[1125]:=

i''[t]

Out[1125]=

m ((2 ^(-2 k m t) k^2 m^2 (-i0 + m)^2)/(i0^2 (1 + (^(-k m t) (-i0 + m))/i0)^3) - (^(-k m t) k^2 m^2 (-i0 + m))/(i0 (1 + (^(-k m t) (-i0 + m))/i0)^2))

Let's simplify this expression.

In[1126]:=

Simplify[i''[t]]

Out[1126]=

-(^(k m t) i0 k^2 (i0 + ^(k m t) i0 - m) m^3 (-i0 + m))/((-1 + ^(k m t)) i0 + m)^3

When is this zero?

(1) when i0 = 0
(2) when k = 0
(3) when m = i0
(4) when  i0  +  i0  Exp[kmt] - m = 0
(5) when m = 0

Case (4) is the most interesting

In[1127]:=

Solve[i0 + i0 Exp[k m t] - m == 0, t]

Solve :: ifun : Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information.  More…

Out[1127]=

{{tLog[-(i0 - m)/i0]/(k m)}}

Gause's Yeast

Hughes-Hallett, Gleason, et al, Exercise 9.7.6, page 537

In[1128]:=

Clear[data];

data = {{0,0.37},{10,8.87},{18,10.66},{23,12.5},{34,13.27},{42,12.87},{47,12.7}};

dots = ListPlot[data,
                PlotStyle->{MarsOrange,PointSize[0.02]}];

[Graphics:HTMLFiles/7.4_populationGrowth_18.gif]

Let's take L = 13.

In[1131]:=

L = 13;

Exercise 9.7.6b
Estimate k.

In[1132]:=


Clear[p,k];

p[0] = 0.37;
p[10] = 8.87;

estimatedDerivativeAtZero = (p[10] - p[0])/10;

k = estimatedDerivativeAtZero / (p[0] (1 - p[0]/L))

Out[1136]=

2.3646

Exercise 9.7.6c

As a check, look at p[18] and compare your answer with the value of 10.66 from the table.

In[1137]:=

k = .28;

a = (L - p[0])/p[0];

p[t_] := L/(1 + a Exp[- k t]);

c = 0;
d = 40;

plot = Plot[p[t],{t,c,d},
     PlotStyle->Red,
     PlotLabel->"p[t]",
     PlotRange->{0,13},
     AxesLabel->{"t (hrs)","p (num cells)"}];
     
p[18]

[Graphics:HTMLFiles/7.4_populationGrowth_20.gif]

Out[1144]=

10.6472

How well does our solution fit the original data?

In[1145]:=

Show[dots,plot];

[Graphics:HTMLFiles/7.4_populationGrowth_22.gif]

Threshold Population

Hughes-Hallett, Gleason, et al, Exercise 9.7.17a, page 537

Graph the derivative,dp/dt = deriv[p], against p:
    dp/dt = deriv[p] = p (p - 6)

In[1146]:=

Clear[deriv,p,c,d];

deriv[p_] := p (p - 6);

c = -1; d = 8;    (* c <= p <= d *)

Plot[deriv[p],{p,c,d},
     PlotStyle->Red,
     PlotLabel->"dp/dt = p (p - 6)",
     AxesLabel->{"p (population)","dp/dt"}];

[Graphics:HTMLFiles/7.4_populationGrowth_23.gif]

Hughes-Hallett, Gleason, et al, Exercise 9.7.17b, page 537

Find the general solution.

In[1150]:=

DSolve[p'[t] == p[t] (p[t] - 6), p[t], t]

Out[1150]=

{{p[t]  -6/(-1 + ^(6 t + 6 C[1]))}}

Define it, and then check that this definition satisfies the differential equation.

In[1151]:=

Clear[p,t,c,d,const,p0];

p[t_] := 6/(1 + const Exp[6t]);

Simplify[p'[t] == p[t] (p[t] - 6)]

Out[1153]=

True

Find find the solution satisfying p[0] = 5:
         p0 =  p[0] = 6/(1 + const)  => const = (6 / p0) - 1
If p0 = 5, take const = -6/5.

In[1154]:=

p[t_] := 6/(1 + const Exp[6t]);

const = 6/p0 - 1;
p0 = 5;

c = 0; d = 1;    (* c <= t <= d *)

diesOut = Plot[p[t],{t,c,d},
               PlotStyle->Red,
               PlotRange->{0,8},
               PlotLabel->"soln[t] = 6/(1 + 1/5 Exp[6t])",
               AxesLabel->{"t (time)",None}];

[Graphics:HTMLFiles/7.4_populationGrowth_26.gif]

Hughes-Hallett, Gleason, et al, Exercise 9.7.17c, page 537

Find soln for generic p[0].
Try it for p[0] = 8; so that const = -1/4.
We have to shorten the domain of the variable t to get a decent picture,
because this function blows up at t = 0.2.

In[1159]:=

Clear[p,t,c,d,const,p0];

p[t_] := 6/(1 - 1/4 Exp[6t]);

c = 0; d = 0.2;    (* c <= t <= d *)

grows = Plot[p[t],{t,c,d},
             PlotStyle->Red,
             PlotRange->{0,30},
             PlotLabel->"soln[t] = 6/(1 - 1/4 Exp[6t])",
             AxesLabel->{"t (time)",None}];

[Graphics:HTMLFiles/7.4_populationGrowth_27.gif]

Hughes-Hallett, Gleason, et al, Exercise 9.7.17d, page 537

In[1163]:=

Needs["Graphics`PlotField`"]
Clear[f,p,t,a,b,c,d];

(* dp/dt = f[p,t] = p (p - 6) *)

f[p_,t_] := p (p - 6);

a = 0; b = 8;    (* a <= t <= b *)
c = 0; d = 8;    (* c <= p <= d *)

pts = {};

field = PlotVectorField[{1,f[p,t]},
                        {t,a,b},{p,c,d},
                        PlotLabel -> "dp/dt = p (p - 6)",
                        Axes -> True,
                        AxesLabel -> {"t","p"},
                        PlotPoints -> 20,
                        Prolog -> ManganeseBlue,
                        Epilog -> {Red,PointSize[0.02],
                                   Map[Point,pts]}];

[Graphics:HTMLFiles/7.4_populationGrowth_28.gif]

Let's try to plot everything in the same picture.

In[1171]:=

Clear[p1,t,c1,d1,const,field,grows,diesOut];

p1[t_] := 6/(1 + 1/5 Exp[6t]);

c1 = 0; d1 = 1;    (* c <= t <= d *)

diesOut = Plot[p1[t],{t,c1,d1},
               PlotStyle->Red,
               PlotRange->{0,8},
               DisplayFunction->Identity];
               
               
Clear[p2,t,c2,d2,const];
             
p2[t_] := 6/(1 - 1/4 Exp[6t]);

c2 = 0; d2 = 0.1;    (* c <= t <= d *)

grows = Plot[p2[t],{t,c2,d2},
             PlotStyle->Red,
             PlotRange->{0,10},
             DisplayFunction->Identity];

Clear[f,p,t,a,b,c,d];

f[p_,t_] := p (p - 6);

a = 0; b = 1;    (* a <= t <= b *)
c = 0; d = 10;    (* c <= p <= d *)

field = PlotVectorField[{1,f[p,t]},
                        {t,a,b},{p,c,d},
                        PlotLabel -> "dp/dt = p (p - 6)",
                        Axes -> True,
                        AxesLabel -> {"t","p"},
                        PlotPoints -> 10,
                        Prolog -> ManganeseBlue,
                        AspectRatio->0.5,
                        DisplayFunction->Identity];
                                          
Show[{field,grows,diesOut},
      DisplayFunction->$DisplayFunction];

[Graphics:HTMLFiles/7.4_populationGrowth_29.gif]

Could do a bit better with those funky arrowheads, but you get the idea.

Sustainable Yield

Hughes-Hallett, Gleason, et al, Exercises 9.7.19-21, pages 540-541

Let's look at the function which describes the rate of growth
of the fish population when there is no fishing.

In[1189]:=

Clear[r,p,c,d,k];

r[p_] := p (2 - k p);

k = 0.01;          (* growth rate parameter *)
c = 0; d = 300;    (* c <= p <= d *)

rate = Plot[r[p],{p,c,d},
            PlotStyle->Red,
            PlotLabel->"dp/dt = p (2 - k p)",
            AxesLabel->{"p (population size)","dp/dt"}];

[Graphics:HTMLFiles/7.4_populationGrowth_30.gif]

In the absence of fishing, the fish population should grow
when there are fewer than 200 fish, and it should decrease
when there are more than 200 fish. Let's see if we can get
a slope field to tell that story.

In[1194]:=

Needs["Graphics`PlotField`"]
Clear[rate,p,t,a,b,c,d];

rate[p_,t_] := p (2 - k p);

k = 0.01;              (* growth rate parameter *)
removalRate = 75;      (* fish per year *)
a = 0; b = 1;          (* a <= t <= b *)
c = 199; d = 200.5;    (* c <= p <= d *)

pts = {};

field = PlotVectorField[{1,rate[p,t]},
                        {t,a,b},{p,c,d},
                        PlotLabel -> "Growth of Fish Population\n   dp/dt = p (2 - k p)",
                        Axes -> True,
                        AxesLabel -> {"t","p"},
                        PlotPoints -> 20,
                        AspectRatio->1.0,
                        Prolog -> ManganeseBlue,
                        Epilog -> {Red,PointSize[0.02],
                                   Map[Point,pts]}];

[Graphics:HTMLFiles/7.4_populationGrowth_31.gif]

Now we allow fishing. Fishermen (fisherfolk) remove 75 fish per year,
and this contributes to the rate of change in the growth of the fish population.

In[1203]:=

Clear[dp,p,c,d,k];

dp[p_] := p (2 - k p) - removalRate;

removalRate = 75;
k = 0.01;          (* growth rate parameter *)
c = 0; d = 300;    (* c <= t <= d *)

rateWithFishing = Plot[dp[p],{p,c,d},
                       PlotStyle->Red,
                       PlotLabel->"dp/dt = p (2 - k p) - removalRate",
                       AxesLabel->{"p (population size)","dp/dt"}];

[Graphics:HTMLFiles/7.4_populationGrowth_32.gif]

Oops! -- big effect! Populations that start with fewer than 50 fish will be driven to extinction, while populations with more fish when the fishing season opens will tend towards the stable size of 150 fish -- lower than the 200 we saw earlier.

Here is a view of the associated slope field.

In[1209]:=

Clear[r2,p,t,a,b,c,d];

r2[p_,t_] := p (2 - k p) - removalRate;

k = 0.01;               (* growth rate parameter *)
removalRate = 75;       (* fish per year *)
a = 0;   b = 6;         (* a <= t <= b *)
c = 20; d = 160;        (* c <= p <= d *)

pts = {{0,40},{0,60},{0,140},{0,160}};

field = PlotVectorField[{1,r2[p,t]},
                        {t,a,b},{p,c,d},
                        PlotLabel -> "Growth under Fishing Pressure\ndp/dt = p (2 - k p) - removalRate",
                        Axes -> True,
                        AxesLabel -> {"t","p"},
                        PlotPoints -> 12,
                        AspectRatio->1.0,
                        Prolog -> ManganeseBlue,
                        Epilog -> {Red,PointSize[0.02],
                                   Map[Point,pts]}];

[Graphics:HTMLFiles/7.4_populationGrowth_33.gif]

Now let's plot some particular solutions.

Here, p[0] = 40, 60, and 160.

In[1217]:=

Clear[p,t,p0,k,p40,p60,p140,p160,c,d]

p0 = 60;                   (* p0 = initial population size *)
k = 1 - 100/(p0 - 50);

k60  = 1 - 100/(60 - 50);
k140 = 1 - 100/(140 - 50);
k160 = 1 - 100/(160 - 50);

p60[t_]  := 50 + 100/(1 - k60 Exp[-t])
p140[t_] := 50 + 100/(1 - k140 Exp[-t])
p160[t_] := 50 + 100/(1 - k160 Exp[-t])

c = 0; d = 6;    (* c <= t <= d *)

plots = Plot[{p60[t],p140[t],p160[t]},{t,c,d},
             PlotStyle->Red,
             PlotLabel->"p[0] = 60, 140, and 160",
             PlotRange->{40,160},
             AxesLabel->{"t (time)","p (population size)"}];

[Graphics:HTMLFiles/7.4_populationGrowth_34.gif]

Here is an example in which the population dies out; set p[0] = 40.

In[1228]:=

c = 0; d = 1.2;    (* c <= t <= d *)

k40 = 1 - 100/(40 - 50);

p40[t_] := 50 + 100/(1 - k40 Exp[-t])

pop40 = Plot[p40[t],{t,c,d},
             PlotStyle->Red,
             PlotRange->{0,50},
             PlotLabel->"p[0] = 40",
             AxesLabel->{"t (time)","p (population size)"}];

[Graphics:HTMLFiles/7.4_populationGrowth_35.gif]

Can we get all of this into one picture?

In[1232]:=

Show[{field,plots,pop40},
     DisplayFunction->$DisplayFunction];

[Graphics:HTMLFiles/7.4_populationGrowth_36.gif]

Now consider the effect of a harvesting rate of H = 100 fish per year.

In[1233]:=

Clear[dp,p,c,d,k];

dp[p_] := p (2 - k p) - removalRate;

removalRate = 100;
k = 0.01;          (* growth rate parameter *)
c = 0; d = 200;    (* c <= t <= d *)

rateWithFishing = Plot[dp[p],{p,c,d},
                       PlotStyle->Red,
                       PlotLabel->"dp/dt = p (2 - k p) - removalRate",
                       AxesLabel->{"p (population size)","dp/dt"}];

[Graphics:HTMLFiles/7.4_populationGrowth_37.gif]

There is a single root at p = 100, and for any other value of p the derivative dp/dt is negative. Thus, if the fishing rate is 100 fish per year, a population with more than 100 fish will approach an equilibrium population size of 100 fish. Any smaller population will die out.

In[1239]:=

Needs["Graphics`PlotField`"]
Clear[rate,p,t,a,b,c,d];

rate[p_,t_] := p (2 - k p) - removalRate;

k = 0.01;              (* growth rate parameter *)
removalRate = 100;     (* fish per year *)
a = 0; b = 10;          (* a <= t <= b *)
c = 90; d = 110;        (* c <= p <= d *)

pts = {};

field = PlotVectorField[{1,rate[p,t]},
                        {t,a,b},{p,c,d},
                        PlotLabel -> "Growth of Fish Population\n   dp/dt = p (2 - k p) - 100",
                        Axes -> True,
                        AxesLabel -> {"t","p"},
                        PlotPoints -> 20,
                        AspectRatio->1.0,
                        Prolog -> ManganeseBlue,
                        Epilog -> {Red,PointSize[0.02],
                                   Map[Point,pts]}];

[Graphics:HTMLFiles/7.4_populationGrowth_38.gif]

Finally, we consider the effect of a harvesting rate of H = 200 fish per year.

In[1248]:=

Clear[dp,p,c,d,k];

dp[p_] := p (2 - k p) - removalRate;

removalRate = 200;
k = 0.01;          (* growth rate parameter *)
c = 0; d = 200;    (* c <= t <= d *)

rateWithFishing = Plot[dp[p],{p,c,d},
                       PlotStyle->Red,
                       PlotLabel->"dp/dt = p (2 - k p) - removalRate",
                       AxesLabel->{"p (population size)","dp/dt"}];

[Graphics:HTMLFiles/7.4_populationGrowth_39.gif]

The derivative dp/dt is always negative! Thus, if the fishing rate is 200 fish per year, any population of fish will die out, no matter how large the initial population.

In[1254]:=

Needs["Graphics`PlotField`"]
Clear[rate,p,t,a,b,c,d];

rate[p_,t_] := p (2 - k p) - removalRate;

k = 0.01;              (* growth rate parameter *)
removalRate = 200;     (* fish per year *)
a = 0; b = 10;          (* a <= t <= b *)
c = 90; d = 110;        (* c <= p <= d *)

pts = {};

field = PlotVectorField[{1,rate[p,t]},
                        {t,a,b},{p,c,d},
                        PlotLabel -> "Growth of Fish Population\ndp/dt = p (2 - k p) - 200",
                        Axes -> True,
                        AxesLabel -> {"t","p"},
                        PlotPoints -> 10,
                        AspectRatio->1.0,
                        Prolog -> ManganeseBlue,
                        Epilog -> {Red,PointSize[0.02],
                                   Map[Point,pts]}];

[Graphics:HTMLFiles/7.4_populationGrowth_40.gif]


Created by Mathematica  (April 25, 2004)