The Leslie Model for Population Growth

This idea is best illustrated by an example.  Suppose we are studying a population of birds, which has been introduced into a new area where there are few predators.   In order to understand how the population grows,  we divide it into three age groups:  birds of age less than 2 years, birds of age between 2 and 4 years, and birds of age between 4 and 6 years.   Now as a matter of fact, some of these birds live longer than 6 years, but they are so few in number that we will neglect them and assume that no bird lives longer than 6 years.  Let's label our groups 1, 2 and 3.  From extensive empirical data we have that the surivival rate for group 1 is 50%, which the survival rate for group 2 is 25%.  On the other hand, these birds produce no offspring in the first 2 years, an average of 4 offspring in the next two years and 3 offspring in the last two years. Let's convert the percentages into decimal numbers and write the reproductive numbers for each group as an array  a,  and the survival rates for each group as an array b:

a = {0.0, 4, 3};
b = {0.5, 0.25};

Notice that the vector  b  only has two elements.  We didn't bother to write the survival rate for group 3; by our assumptions, it is 0.  Also notice the decimals.  This had the effect of asking Mathematica to treat every number in the arrays as finite precision floating point numbers rather than exact fractions.  This changes how Mathematica tries to do some calculations.   Finally, observe the semicolins.  In this context, they have the effect of suppressing output.  Sometimes this isn't desirable.  Remove them and reclick if you want to see output.

Now for the problem of modelling this population.  Let the population in a given (2-year) time period be the vector x = {x1, x2, x3}, where the first coordinate is the group 1 population, etc.  How would we get the population vector y = {y1, y2, y3} for the next time period?    A little thought leads to the following system:

	y1 = a[[1]]*x1 + a[[2]]*x2 + a[[3]]*x3
	y2 = b[[1]]*x1
	y3 =                    b[[2]]*x2

See the pattern?  This is really a discrete linear dynamical system.    In fact we can write it as
	
			y  =  L . x

Let's build the appropriate matrix   L  from the data   a   and  b.  Such a matrix is called  a  "Leslie  matrix".   It will take several steps:

L = Table[Switch[i-j,1, b[[j]], _,0],{i,3},{j,3}];
L[[1]] = a;

Let's have a look at L:

MatrixForm[L]

Here is a possible  initial population vector:

x0  = Table[100, {i,3}]

Let's see what happens in one time period and then a subsequent one:

L.x0

L.L.x0

L.L.L.x0

Recall what the situation is from our earlier eigenvalues notebook.  Eigenvalues govern the behavior of powers of  L  (provided  L  is diagonalizable, which we'll assume here).  In fact, we have that if  P  diagonalizes  L , i.e.,  L^{-1}.P.L = D, a diagonal matrix,  then for all integers  k

		L^k  =  P . D^k . P^{-1}.
		
Now if there is a  dominant  eigenvalue,  i.e., one that is larger in absolute value than all the others,  then  this eigenvalue  governs growth.  Experiment with the previous two cells (run them out a little farther) and make an estimate as to the growth rate.  Compare this against the following:

Eigensystem[L]

So it appears the population is going to grow.  Suppose that we can control the population in some predictable way which has the effect of removing a precentage of each age group from the population.  To fix an order of events for our model, we take the point of view that each population group is allowed to grow over a time period, and then at the end of the period a certain fraction  h[[i]]  is harvested from the group before  they "graduate" to the next level.  In the case that the percentage is the same for all groups, we say that we have a  "uniform harvest policy."  In any case, let  H be the "harvest matrix",  that is,  H is a diagonal matrix with h[[i]] in the ith diagonal spot.  The population dynamics now change.  Here is the new model:

		y  =  L . x  -  H . L . x  =  (I  -  H) . L . x
		
Thus the new transition matrix from one state to the next is   A  =  (I  -  H) . L  Interestingly enough, as long as  H  is sensible  (all diagonal entries between 0 and 1),  the new matrix  A  is also a Leslie matrix.

Here is a sample uniform harvest vector :

h = Table[.05,{i,3}]

This leads to a harvest matrix:

H = DiagonalMatrix[h];

Next, you can construct the iteration matrix for this population model:

A = (IdentityMatrix[3]-H).L

Here is how you interpret this matrix:  A population is allowed to grow for one year according to the matrix  L.  At the end of that year, harvest takes place, resulting in the new population for the next year.  You now have the necessary data for experiments.  For example, here is how the population fares after one year with the above data:

x1 = A.x0

Here is a way to find the largest eigenvalue of  the resulting iteration   matrix A  for population growth (which you want to be 1 for a stable population).

Max[Abs[Eigenvalues[A]]]

You now have the necessary data for experiments.  For example, you can change

You now have the necessary data for experiments.  For example, to get the long term population distribution for this stragegy, you could find the eigenvector corresponding to this eigenvalue, then normalize it by dividing by the sum of the entries.  Dot it with {1,1,1} to get the sum.  Or just let the population run on for a while and normalize the resulting vector to get a population.  A few the the basic theoretical facts about Leslie matrices are the following:

DEFINITION:  The square matrix  A  is a Leslie matrix if A has the form
		{{a_1,  a_2, a_3, ...   ,  a_{n-1}, a_n},
		  {b_1,  0   ,   0  ,  ...  ,  0           ,  0},
		  {0    , b_2,   0  ,  ...  ,  0           ,  0},
		    ...
		  {0   ,  0   ,    0  , ...   , b_{n-1},   0}},
		  
where the numbers  b_i  are positive and the numbers a_i  are nonnegative.

THEOREM 1.   If  A  is a Leslie matrix with at least one positive entry in the first row, then  A      will have exactly one positive eigenvalue which is a simple (not repeated) eigenvalue of  A .

THEOREM 2.  If  A  is a Leslie matrix with at least two adjacent positive entries in the first row, then the eigenvalue of Theorem 1 will be a dominant eigenvalue, i.e., it is strictly larger than any other eigenvalue in absolute value.

Observe that these theorems fail if the hypotheses are weakened.  Click on the following examples:

a = {{0,0,0},{1,0,0},{0,1,0}};
MatrixForm[a]

Eigenvalues[a]

a = {{0,0,1},{1,0,0},{0,1,0}};
MatrixForm[a]

Eigenvalues[a]