The MATLAB Notebook v1.5.2

Population Dynamics

We are going to look at two models for population growth of a species. The first is a standard exponential growth/decay model that describes quite well the population of a species becoming extinct, or the short-term behavior of a population growing in an unchecked fashion. The second, more realistic model, describes the growth of a species subject to constraints of space, food supply, and competitors/predators.

Exponential Growth/Decay

We assume that the species starts with an initial population P₀. The population after n time units is denoted P_n. Suppose that in each time interval, the population increases or decreases by a fixed proportion of its value at the beginning of the interval. Thus P_n = P_n-₁ + r P_n-₁, n ³ 1. The constant r represents the difference between the birth rate and the death rate. The population increases if r is positive, decreases if r is negative, and remains fixed if r = 0.

Here is a simple M-file that will compute the population at stage n, given the population at the previous stage and the rate r.

function X = itseq(f,Xinit,n,r)

% computing an iterative sequence of values

X = zeros(n + 1, 1);

X(1) = Xinit;

for i = 1:n

X(i + 1) = f(X(i), r);

end

In fact, this is a simple program for computing iteratively the values of a sequence

a_n = f(a_n-₁), n ³ 1, provided you have previously entered the formula for the function f and the initial value of the sequence a₀. Note the extra parameter r is built into the algorithm.

Now let's use the program to compute two populations at five year intervals for different values of r:

r = 0.1; Xinit = 100; f = inline('x*(1 + r)', 'x', 'r');

X = itseq(f, Xinit, 100, r);

format long; X(1:5:101)

ans =

1.0e+006 *

0.00010000000000

0.00016105100000

0.00025937424601

0.00041772481694

0.00067274999493

0.00108347059434

0.00174494022689

0.00281024368481

0.00452592555682

0.00728904836851

0.01173908528797

0.01890591424713

0.03044816395414

0.04903707252979

0.07897469567994

0.12718953713951

0.20484002145855

0.32989690295921

0.53130226118483

0.85566760466078

1.37806123398224

r = -0.1; X = itseq(f, Xinit, 100, r);

X(1:5:101)

ans =

1.0e+002 *

1.00000000000000

0.59049000000000

0.34867844010000

0.20589113209465

0.12157665459057

0.07178979876919

0.04239115827522

0.02503155504993

0.01478088294143

0.00872796356809

0.00515377520732

0.00304325272217

0.00179701029991

0.00106111661200

0.00062657874822

0.00036998848504

0.00021847450053

0.00012900700782

0.00007617734805

0.00004498196225

0.00002656139889

In the first case, the population is growing rapidly; in the second, decaying rapidly. In fact, it is clear from the model that, for any n, the quotient P_n/P_n₊₁= (1 + r), and therefore it follows that P_n = P₀(1 + r)ⁿ, n ³ 0. This accounts for the expression exponential growth/decay. The model predicts a population growth without bound (for growing populations), and is therefore not realistic. Our next model allows for a check on the population caused by limited space, limited food supply, competitors, and predators.

Logistic Growth

The previous model assumes that the relative change in population is constant, that is

(P_n₊₁ – P_n)/P_n = r.

Now let's build in a term that holds down the growth, namely

(P_n₊₁ – P_n)/P_n = r – uP_n.

We shall simplify matters by assuming that u = 1 + r , so that our recursion relation becomes

P_n₊₁ = u P_n(1 – P_n),

where u is a positive constant. In this model, the population P is constrained to lie between 0 and 1, and should be interpreted as a percentage of a maximum possible population in the environment in question. So let us set up the function we will use in the iterative procedure:

clear f; f = inline('u*x*(1 - x)', 'x', 'u');

Now let's compute a few examples; and use plot to display the results.

u = 0.5; Xinit = 0.5; X = itseq(f, Xinit, 20, u); plot(X)

u = 1; X = itseq(f, Xinit, 20, u); plot(X)

u = 1.5; X = itseq(f, Xinit, 20, u); plot(X)

u = 3.4; X = itseq(f, Xinit, 20, u); plot(X)

In the first computation, we have used our iterative program to compute the population density for 20 time intervals, assuming a logistic growth constant u = 0.5, and an initial population density of 50%. The population seems to be dying out. In the remaining examples, we kept the initial population density at 50%; the only thing we varied was the logistic growth constant. In the second example, with a growth constant u = 1, once again the population is dying out—although more slowly. In the third example, with a growth constant of 1.5 the population seems to be stabilizing at 33.3...%. Finally, in the last example, with a constant of 3.4 the population seems to oscillate between densities of approximately 45% and 84%.

These examples illustrate the remarkable features of the logistic population dynamics model. This model has been studied for more than 150 years, its origins lying in an analysis by the Belgian mathematician Verhulst. Here are some of the facts associated with this model. We will corroborate some of them with MATLAB. In particular, we shall use bar as well as plot to display some of the data.

(1) The logistic constant cannot be larger than 4.

In order for the model to work, the output at any point must be between 0 and 1. But the parabola ux(1 - x), for 0 £ x £ 1, has its maximum height when x = 1/2, where its value is u/4. To keep that number between 0 and 1, we must restrict u £ 4. Here is what happens if u is bigger than 4:

u = 4.5; Xinit = 0.9; X = itseq(f, Xinit, 10, u)

X =

1.0e+072 *

0.00000000000000

-0.00000000000000

-3.49103403458070

(2) If 0 £ u £ 1, the population density tends to zero for any initial configuration.

X = itseq(f, 0.99, 100, 0.8); X(101)

ans =

1.939524024691387e-012

X = itseq(f, 0.75, 20, 1);

bar(X)

(3) If 1 < u £ 3, the population will stabilize at density 1 -1/u for any initial density other than zero.

The third of the original four examples corroborates the assertion (with u = 1.5 and 1 -1/u = 1/3). In the following examples, we set u = 2, 2.5 and 3 respectively, so that 1 -1/u equals 0.5, 0.6 and 0.666..., respectively. The convergence in the last computation is rather slow (as one might expect from a boundary case—or bifurcation point).

X = itseq(f, 0.25, 100, 2); X(101)

ans =

0.50000000000000

X = itseq(f, 0.75, 100, 2); X(101)

ans =

0.50000000000000

X = itseq(f, 0.5, 20, 2.5);

plot(X)

X = itseq(f, 0.75, 100, 3);

bar(X); axis([0 100 0 0.8])

(4) If 3 < u < 3.56994..., then there is a periodic cycle.

The theory is quite subtle. For a fuller explanation, the reader may consult Encounters with Chaos, by Denny Gulick (McGraw-Hill), section 1.5. In fact there is a sequence

u₀ = 3 < u₁ = 1 + Ö6 < u₂ < u₃ < ... < 4,

such that: between u₀ and u₁ there is a cycle of period 2; between u₁ and u₂ there is cycle of period 4; and in general, between u_k and u_k₊₁ there is a cycle of period 2^k⁺¹. In fact one knows that, at least for small k, one has the approximation u_k₊₁ » 1 + Ö(3+u_k). So

u1 = 1 + sqrt(6)

u1 =

3.44948974278318

u2approx = 1 + sqrt(3 + u1)

u2approx =

3.53958456106175

This explains the oscillatory behavior we saw in the last of the original four examples (with u₀ < u = 3.4 < u₁). Here is the behavior for u₁ < u = 3.5 < u₂. The command bar is particularly effective here for spotting the cycle of order 4.

X = itseq(f, 0.75, 100, 3.5);

bar(X); axis([0 100 0 0.9])

(5) There is a value u < 4 beyond which—chaos!

It is possible to prove that the sequence u_k tends to a limit u_¥. The value of u_¥, sometimes called the Feigenbaum parameter, is aproximately 3.56994... Let's see what happens if we use a value of u between the Feigenbaum parameter and 4.

X = itseq(f, 0.75, 100, 3.7); plot(X)

This is an example of what mathematicians call a chaotic phenomenon! It is not random—the sequence was generated by a precise, fixed mathematical procedure, but the results manifest no discernible pattern. Chaotic phenomena are unpredictable, but with modern methods (including computer analysis), mathematicians have been able to identify certain patterns of behavior in chaotic phenomena. For example, the last figure suggests the possibility of unstable periodic cycles and other recurring phenomena. Indeed a great deal of information is known. The aforementioned book by Gulick is a fine reference, as well as the source of an excellent bibliography on the subject.

Rerunning the Model with SIMULINK

The logistic growth model that we have been exploring lends itself particularly well to simulation using SIMULINK. Here is a simple SIMULINK model that corresponds to the above calculations:

Let's briefly explain how this works. If you ignore the Discrete Pulse Generator block and the Sum block in the lower left for a moment, this model implements the equation

x at next time = u x (1 – x) at old time,

which is the equation for the logistic model. The Scope block displays a plot of x as a function of (discrete) time. However, we need somehow to build in the initial condition for x. The simplest way to do this is as illustrated here: we add to the right-hand side a discrete pulse that is the initial value of x at time t = 0 and is 0 thereafter. Since the model is discrete, you can achieve this by setting the period of the Discrete Pulse Generator block to something longer than the length of the simulation, and setting the width of the pulse to 1 and the amplitude of the pulse to the initial value of x. The outputs from the model in the two interesting cases of u = 3.4 and u = 3.7 are shown here:

Output with u = 3.4

Output with u = 3.7

In the first case of u = 3.4, the periodic behavior is clearly visible. On the other hand, when u = 3.7, we get chaotic behavior.