Exponential and logistic growth 

If you click on !!! (the menu bar at the top) maple will execute all commands.  

You can also enter each command as you go along  reading this document. You may enter your own commands at any time and place in this document. 

 

> restart; -1
 

> with(plots); -1
 

> `assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
`assign`(logi, proc (r, K) option operator; description
 

Error, `:=` unexpected
 

>
 

Exponential growth 

The projection function for logistic growth is  

> `assign`(F, proc (P) options operator, arrow; `+`(P, `*`(r, `*`(P))) end proc); 1
 

proc (P) options operator, arrow; `+`(P, `*`(r, `*`(P))) end proc (1)
 

where r>-1 is a constant. When r>0, F(P)>P and the population is increasing; when -1<r<0,  F(P)<P and the populationen is decreasing. 

Here are the graphs of the projection functions F(P)=P+rP  for r=0.05 (increasing population) and -0.05 (decreasing population)  together with the diagonal F(P)=P (constant population) 

> plot([`*`(`+`(1, 0.5e-1), `*`(P)), `*`(`+`(1, -0.5e-1), `*`(P)), P], P = 0 .. 100, labels = [F, F(P)], color = [green, green, red]); 1
 

Plot_2d
 

The projection function for the increasing population lies above the projection function for the constant population, and the projectione function for the decreasing population lies below it. 

We may illustrate the population dymanics by cobwebs 

> display([cobweb(proc (P) options operator, arrow; `*`(`+`(1, 0.5e-1), `*`(P)) end proc, 20, 32, 0, 100), cobweb(proc (P) options operator, arrow; `*`(.95, `*`(P)) end proc, 90, 32, 0, 100)], labels = ...
display([cobweb(proc (P) options operator, arrow; `*`(`+`(1, 0.5e-1), `*`(P)) end proc, 20, 32, 0, 100), cobweb(proc (P) options operator, arrow; `*`(.95, `*`(P)) end proc, 90, 32, 0, 100)], labels = ...
 

Plot_2d
 

The growth function ΔF(P)=F(P)-P=rP is increasing/decreasing with increasing population density 

> plot([`+`(`-`(`*`(0.5e-1, `*`(P)))), `+`(`*`(0.5e-1, `*`(P)))], P = 0 .. 100, labels = [P, `ΔF`(P)]); 1
 

Plot_2d
 

The population orbit with P0=2 when r=0.05 

> dynplotfilm(proc (P) options operator, arrow; `*`(`+`(1, 0.5e-1), `*`(P)) end proc, 2); 1
 

Plot_2d
 

Population orbits with different initial population sizes P0 when r=0.05 

> display(seq(dynplot(proc (P) options operator, arrow; `*`(`+`(1, 0.5e-1), `*`(P)) end proc, `+`(2, `*`(1000, `*`(u))), 50), u = 0 .. 10), labels = [t, P[t]], insequence = false); 1
display(seq(dynplot(proc (P) options operator, arrow; `*`(`+`(1, 0.5e-1), `*`(P)) end proc, `+`(2, `*`(1000, `*`(u))), 50), u = 0 .. 10), labels = [t, P[t]], insequence = false); 1
 

Plot_2d
 

The population orbit with P0=2 when r=-0.05 

> dynplotfilm(proc (P) options operator, arrow; `*`(.95, `*`(P)) end proc, 100); 1
 

Plot_2d
 

Population orbits with different initial population sizes P0 when r=-0.05 

> display(seq(dynplot(proc (P) options operator, arrow; `*`(.95, `*`(P)) end proc, `+`(100, `*`(100, `*`(u))), 50), u = 0 .. 10), labels = [t, P[t]], insequence = false); 1
display(seq(dynplot(proc (P) options operator, arrow; `*`(.95, `*`(P)) end proc, `+`(100, `*`(100, `*`(u))), 50), u = 0 .. 10), labels = [t, P[t]], insequence = false); 1
 

Plot_2d
 

Logistic growth 

The projection function is 

> `assign`(F, proc (P) options operator, arrow; `+`(P, `*`(r, `*`(`+`(1, `-`(`/`(`*`(P), `*`(K)))), `*`(P)))) end proc); 1
 

proc (P) options operator, arrow; `+`(P, `*`(r, `*`(`+`(1, `-`(`/`(`*`(P), `*`(K)))), `*`(P)))) end proc (2)
 

Here is the gaph of the projection function   

> (logi(0.5e-1, 100))(P); 1
 

`*`(`+`(1.05, `-`(`*`(0.5000000000e-3, `*`(P)))), `*`(P)) (3)
 

with r=0.05 and carrying capacity K=100 

> display(plot([(logi(0.5e-1, 100))(P), P], P = 0 .. 2200, y = 0 .. 600, labels = [P,
display(plot([(logi(0.5e-1, 100))(P), P], P = 0 .. 2200, y = 0 .. 600, labels = [P,
 

Plot_2d
 

The population is increasing, F(P)>P, when 0<P<K, and decreasing, F(P)<P, when P>K.  

This may be a little difficult to see on the above graph. It becomes a little more clear if we zoom in on the left part of the graph 

> display(plot([(logi(0.5e-1, 100))(P), P], P = 0 .. 110, y = 0 .. 110, labels = [P,
display(plot([(logi(0.5e-1, 100))(P), P], P = 0 .. 110, y = 0 .. 110, labels = [P,
 

Plot_2d
 

Another way of saying this is that the growth function ΔF(P)=F(P)-P is positive for 0<P<K and negative for P>K. 

Here is the graph of the growth function 

The growth function ΔF(P)=F(P)-P  is maximal  when the population densiyt P=K/2 is half of the carrying capacity and  negative beyond the carrying capacity K=100. 

> display([plot(`+`((logi(0.5e-1, 100))(P), `-`(P)), P = 0 .. 110, labels = [P, `ΔF`(P)]), implicitplot([P = 50, y = 1.25], P = 0 .. 110, y = 0 .. 1.3, color = [blue, green])]); 1
display([plot(`+`((logi(0.5e-1, 100))(P), `-`(P)), P = 0 .. 110, labels = [P, `ΔF`(P)]), implicitplot([P = 50, y = 1.25], P = 0 .. 110, y = 0 .. 1.3, color = [blue, green])]); 1
 

Plot_2d
 

The growth function ΔF(P)  is maximal  when the population density P=K/2 is half of the carrying capacity. The maximal value of the growth function is called Maximal Sustainable Yield or MSY.  

The maximal value of the growth function ΔF(P) is achieved when P=K/2 and it is  MSY = max ΔF(P) =1/4*K*r. In this particular case, MSY = 1.25. 

>
 

Here are the population orbits at two initial population sizes P0=2 and P0=130 

> display([dynplotfilm(logi(0.5e-1, 100), 2), dynplotfilm(logi(0.5e-1, 100), 130)]); 1
 

Plot_2d
 

Population orbits with several different initial population sizes P0 

> display(seq(dynplot(logi(0.5e-1, 100), `+`(2, `*`(5, `*`(u))), 120), u = 0 .. 10), labels = [P,
display(seq(dynplot(logi(0.5e-1, 100), `+`(2, `*`(5, `*`(u))), 120), u = 0 .. 10), labels = [P,
 

Plot_2d
 

When r<2 there are two equilibria. The equilibrium 0 is unstable and the equilibrium K is stable. 

For instance when r=1.8 we can see this by solving the equation F(P)=P for P and computing the derivative F'(P) at the two solutions 

> stab1(logi(1.8, 100)); 1
 

stab1(proc (P) options operator, arrow; `*`(`+`(1, `*`(1.8, `*`(`+`(1, `-`(`/`(`*`(P), `*`(100))))))), `*`(P)) end proc) (4)
 

or by drawing cobweb diagrams 

> display(cobweb(logi(1.8, 100), 2, 15, 0, 120), cobweb(logi(1.8, 100), 120, 10, 0, 120)); 1
 

Plot_2d
 

The stable equilibrium determines the long term behaviou becasuse all population orbits (from not too big initial populations) converge to the carrying capacity 

> seq(floor(iterate(logi(1.8, 100), P0, 200)), P0 = [1, 5, 6, 150, 250]); 1
 

100, 100, 100, 100, `+`(`-`(infinity)) (5)
 

Let us now consider what happens when r>2. Here are four examples with increasing parameter r. 

> dynplotfilm(logi(2.2, 100), 2); 1
 

Plot_2d
 

> dynplotfilm(logi(2.5, 100), 2); 1
 

Plot_2d
 

> dynplotfilm(logi(2.6, 100), 2); 1
 

Plot_2d
 

> dynplotfilm(logi(2.7, 100), 2); 1
 

Plot_2d
 

Logistic growth with varying parameter r - stable cycles and chaos emerge 

> display(seq(dynplot(logi(`+`(1, `*`(0.5e-1, `*`(s))), 100), 50, 100, 0, 150), s = 0 .. 40), insequence = true); 1
 

Plot_2d
 

When r is slightly bigger than 2 stable cycles occur.  

We look for 2-cycles in logi(2.2,100). 

> `assign`(F, logi(2.2, 100)); 1
 

proc (P) options operator, arrow; `*`(`+`(1, `*`(2.2, `*`(`+`(1, `-`(`/`(`*`(P), `*`(100))))))), `*`(P)) end proc (6)
 

The population dynamics with P0=5 indicates that there is a 2-cycle 

> dynplot(F, 5, 100); 1
 

Plot_2d
 

We can find the 2-cycle by solving the equation F(F(P))=P for P. We also compute the derivative (FF)'(P) at each solution to determine stability. 

> stab2(F); 1
 

stab2(F) (7)
 

Alternatively, we can use the cobweb diagram for FF 

> display(cobweb(proc (P) options operator, arrow; F(F(P)) end proc, 40, 20, 0, 130), cobweb(proc (P) options operator, arrow; F(F(P)) end proc, 50, 20, 0, 130)); 1
 

Plot_2d
 

It looks as if there are stable 2-cycles through 70 and 115 and unstable 2-cycles at 0 and 100. 

The stable 2-cycles determine the long term behavior because population orbits will converge to them.  

We compute the 100th generation F^100(P0)  with four different initial populations 

> seq(iterate(proc (P) options operator, arrow; F(F(P)) end proc, P0, 100), P0 = [40, 50, 80, 120]); 1
 

116.2844349, 74.62465608, 74.62465608, 116.2844349 (8)
 

and check that we indeed have a 2-cycle 

> orbit(F, %[1], 6); 1
 

116.2844349, 74.62465608, 116.2844349, 74.62465608, 116.2844349, 74.62465608 (9)
 

and that it is stable 

> evalf(subs(P = 74.62465608, diff(F(F(P)), P))); 1
 

.1600000127 (10)
 

> evalf(subs(P = 116.2844349, diff(F(F(P)), P))); 1
 

.160000011 (11)
 

We next look for 4-cycles in logi(2.5,100); 

> `assign`(F, logi(2.5, 100)); 1
 

proc (P) options operator, arrow; `*`(`+`(1, `*`(2.5, `*`(`+`(1, `-`(`/`(`*`(P), `*`(100))))))), `*`(P)) end proc (12)
 

> dynplot(F, 2, 50); 1
 

Plot_2d
 

There seems to be a stable 4 cycle near 53. 

Any orbit near 53 will converge to it. 

> iterate(F, 55, 200); 1
 

53.59475563 (13)
 

We check that this is indeed a 4-cycle 

> orbit(F, %, 4); 1
 

53.59475563, 115.7716989, 70.12378950, 122.4996169 (14)
 

and that it is stable 

> subs(P = 53.59475563, diff(F(F(F(F(P)))), P)); 1
 

-0.30500006e-1 (15)
 

We could also try to compute it directly or use cobweb 

> stab4(logi(2.5, 100)); 1
 

stab4(proc (P) options operator, arrow; `*`(`+`(1, `*`(2.5, `*`(`+`(1, `-`(`/`(`*`(P), `*`(100))))))), `*`(P)) end proc) (16)
 

> display(cobweb(proc (P) options operator, arrow; F(F(F(F(P)))) end proc, 20, 20, 0, 130), cobweb(proc (P) options operator, arrow; F(F(F(F(P)))) end proc, 45, 20, 0, 130), cobweb(proc (P) options oper...
display(cobweb(proc (P) options operator, arrow; F(F(F(F(P)))) end proc, 20, 20, 0, 130), cobweb(proc (P) options operator, arrow; F(F(F(F(P)))) end proc, 45, 20, 0, 130), cobweb(proc (P) options oper...
 

Plot_2d
 

We approximate the stable cycles by orbits 

> seq(iterate(proc (P) options operator, arrow; F(F(F(F(P)))) end proc, P0, 100), P0 = [20, 45, 82, 110]); 1
 

120.0000000, 70.12378950, 53.59475563, 122.4996169 (17)
 

> orbit(F, 70.12378950, 8); 1
 

70.12378950, 122.4996169, 53.59475563, 115.7716989, 70.12378950, 122.4996169, 53.59475563, 115.7716989
70.12378950, 122.4996169, 53.59475563, 115.7716989, 70.12378950, 122.4996169, 53.59475563, 115.7716989
(18)
 

This is the 4-cycle we found before. We have now also found also another cycle, which is a 2-cycle 

> orbit(logi(2.5, 100), 60.00000000, 8); 1
 

60.00000000, 120.0000000, 60.00000000, 120.0000000, 60.00000000, 120.0000000, 60.00000000, 120.0000000
60.00000000, 120.0000000, 60.00000000, 120.0000000, 60.00000000, 120.0000000, 60.00000000, 120.0000000
(19)
 

but it is unstable 

> evalf(subs(P = 60, diff(F(F(F(F(P)))), P))); 1
 

1.562500000 (20)
 

> evalf(subs(P = 60, diff(F(F(P)), P))); 1
 

-1.250000000 (21)
 

The unstability can be illustrated by orbits through points close to P0=20 

> seq(iterate(F, P0, 100), P0 = [19.995, 20, 20.005]); 1
 

122.4996169, 120.0000000, 115.7716989 (22)
 

Here is the cobweb diagram with r=2.5: 

> cobweb(logi(2.5, 100), 50, 1000, 0, 130); 1
 

Plot_2d
 

and here is the cobweb diagram with r=2.6: 

> cobweb(logi(2.6, 100), 50, 1000, 0, 130); 1
 

Plot_2d
 

The stable orbits are indicated in the bifurcation diagram 

> pointplot([seq(seq([`+`(1.8, `*`(0.5e-2, `*`(i))), orbit(logi(`+`(1.8, `*`(0.5e-2, `*`(i))), 100), 80, 250)[j]], j = 230 .. 250), i = 1 .. 240)], style = point, color = blue); 1
pointplot([seq(seq([`+`(1.8, `*`(0.5e-2, `*`(i))), orbit(logi(`+`(1.8, `*`(0.5e-2, `*`(i))), 100), 80, 250)[j]], j = 230 .. 250), i = 1 .. 240)], style = point, color = blue); 1
 

Plot_2d
 

>
 

Logistic growth with harvest - try to change the harvest 

> dynplotfilm(`+`(logi(.5, 100), `-`(10.5)), 30); 1
 

Plot_2d
 

> dynplot(`+`(logi(.5, 100), `-`(4)), 20, 100); 1
 

Plot_2d
 

Maximal Sustainable Yield MSY 

> `assign`(F, logi(.8, 200)); 1
 

proc (P) options operator, arrow; `*`(`+`(1, `*`(.8, `*`(`+`(1, `-`(`/`(`*`(P), `*`(200))))))), `*`(P)) end proc (23)
 

> `assign`(MSY, `+`(`*`(.8, `*`(`*`(`/`(1, 4), 200))))); 1
 

40.0 (24)
 

An equilibrium  for logistic growth with harvest H is a  population size P so that growth ΔF(P) equals harvest, H.  

The equilibria can be determined from the intersection points between the graph of the growth function and the horizontal line at level H 

> plot([`+`(F(P), `-`(P)), MSY, `+`(`*`(.1, `*`(MSY))), `+`(`*`(.5, `*`(MSY))), `+`(`*`(.8, `*`(MSY)))], P = 0 .. 200); 1
 

Plot_2d
 

Let us consider the case where the harvest is 

> `assign`(H, `+`(`*`(.5, `*`(MSY)))); 1
 

20.00 (25)
 

The equilibria are 

> stab1(proc (P) options operator, arrow; `+`(F(P), `-`(H)) end proc); 1
 

stab1(proc (P) options operator, arrow; `+`(F(P), `-`(H)) end proc) (26)
 

The smallest equilibrium P1=29 is unstable, the biggest P2=170 is stable. This means that if the population size is less than 29, the population will go extinct; if the population is more than 29, it will in the long term stabilize at 

population size 170. 

Here are the population orbits from P0=28, P0=30 and P0=50 

> display([dynplot(proc (P) options operator, arrow; `+`(F(P), `-`(H)) end proc, 30, 20), dynplot(proc (P) options operator, arrow; `+`(F(P), `-`(H)) end proc, 50, 20), dynplot(proc (P) options operator...
display([dynplot(proc (P) options operator, arrow; `+`(F(P), `-`(H)) end proc, 30, 20), dynplot(proc (P) options operator, arrow; `+`(F(P), `-`(H)) end proc, 50, 20), dynplot(proc (P) options operator...
 

Plot_2d
 

We can also see this from the cobweb diagrams from P0=28 and P0=30 

> display([cobweb(proc (P) options operator, arrow; `+`(F(P), `-`(H)) end proc, 28, 7, 0, 180), cobweb(proc (P) options operator, arrow; `+`(F(P), `-`(H)) end proc, 30, 20, 0, 180)]); 1
 

Plot_2d
 

>