LAB 10: AGE-SPECIFIC POPULATION GROWTH AND MATRIX EIGENVALUES

Mathematics:

Demographic predictions of population growth are based on a number of mathematical models. This project exploits the Leslie matrix model, where the female population is divided into age classes of equal duration. As time progresses, the number of females, within each class, changes because of birth, death, and aging.

Suppose we study n age classes of females over maximum age in T years and denote the number of females in the j-th class at the k-th time instance by xj(k). Each age class is T/n years in duration. It makes sense therefore to study changes of the number of females xj(k) between two observation times in T/n years. All females of the j-th class at the k-th time pass into the (j+1)-th class at the (k+1)-th time, except for those who died during the time interval. Denote bj be the fraction of females that survives in the j-th class, i.e.0  bj  1 and

xj+1(k+1) = bj xj(k),   j = 1,2,…,n-1

The first class at the (k+1)-th time has the number of females x1(k+1) that was born during the time interval by females of all other classes. Denote aj be the average number of daughters born by each female in the j-th class during the time inteval, i.e. aj 0 and

x1(k+1) = a1 x1(k) + a2 x2(k) + … + an xn(k)

The following birth and death parameters were registered from the year 1965 for Canadian females:

 Age [0,5) [5,10) [10,15) [15,20) [20,25) [25,30) [30,35) [35,40) [40,45) [45,50) aj 0.00 0.00024 0.05861 0.28608 0.44791 0.36399 0.22259 0.10457 0.02826 0.00240 bj 0.99651 0.99820 0.99802 0.99729 0.99694 0.99621 0.99460 0.99184 0.98700 ---

In this table, the age classes over T = 50 years are neglected since they rarely have children, i.e. aj = 0 ,  j > n, where n = 10 age classes with duration in T/n = 5 years.

Introducing matrix notations, the Leslie model can be written as the discrete dynamical system x(k+1) = L x(k), where x(k) is age distribution vector at the k-th time and L is the Leslie matrix:

L =

The Leslie matrix L is supposed to have at least two non-zero successive entries aj and aj+1. The matrix L has n eigenvalues 1,…,n with n eigenvectors z1,…,zn, such that Lzj =j zj. It can be shown that there is unique positive eigenvalue + which is dominant, i.e. all other negative and complex eigenvalues have smaller absolute value: |j | < +. Let z+ be the eigenvector for the eigenvalue +. Then, the long-term behaviour of the age distribution of the female population is:

x(k) = c +k z+

The population increases if + > 1, decreases if + < 1 and stabilizes if + = 1.

Objectives:

·         Understand the long-term behaviour of the female population within the Leslie model

·         Exploit dynamical modeling of the population growth within the Leslie model

·         Exploit the power method for finding the dominant eigenvalue and its eigenvector of the matrix L

MATLAB script for dynamical modeling of the population growth within the Leslie model

If Canadian women continued to reproduce and die as they did in 1965, the dynamical growth of female population in the ten age classes every five years would follow this graph:

Every five years, the increase in the number of females at each age class, defined as yj(k+1) = (xj(k+1) – xj(k))/xj(k), will approach to a uniform ratio for large k, which is approximately 0.07622 or 7.622%. In the year of 2040, the histogram of population is distributed according to the graph:

If the vector x(k) is normalized such that x1(k) = 1, the vector x(k) approaches approximately for large k:

 Age class [0,5) [5,10) [10,15) [15,20) [20,25) [25,30) [30,35) [35,40) [40,45) [45,50) Proportion 1.000 0.9259 0.8588 0.7964 0.7380 0.6836 0.6328 0.5848 0.5390 0.4943

Steps in writing the MATLAB script:

1. Define T = 50, n = 10, and dT = L/n. Define vectors a and b for the birth and death parameters.
2. Initialize the Leslie matrix L from given vectors a and b.
3. Define the initial age distribution vector x(0) as vector of ones.
4. Organize a loop and compute the age distribution vector x(k) for k = 1,2,…,20. Save the results into a matrix X.
5. Plot dynamics of each age class in time. Scale the time axis t starting with 1965 for k = 0.
6. Compute the yield vector yj(k) and save the results into a matrix Y.
7. Plot the increase in the number of females of each age class in time t.
8. For the year of 2040, compute the histogram of population at each age class.
9. Normalize the vector x(k) for the year 2040 such that x1(k) = 1. Compare x(k) with the table above.

MATLAB script for finding the dominant eigenvalue and its eigenvector of the Leslie matrix

Suppose that the Leslie matrix L is diagonalizable, i.e. the eigenvectors zj  can be combined into matrix P:

P = [ z1, z2, …, zn ]

and eigenvalues j   are combined into matrix D = diag[1, 2 , …, n] such that L = P D P-1. Then, the age distribution vector x(k) at the k-th time is:

x(k) = Lk x(0)  = P Dk P-1 x(0)

If the unique positive eigenvalue+  is dominant, i.e. all other eigenvalues satisfy: |j  | < +, then the age distribution vector x(k)  has the limiting behavior in the limit k :

x(k) c z+

where z+ is the eigenvector for + and c is a constant.

The dominant eigenvalue + and its eigenvector z+ can be found in the iterative power method based on the asymptotic limit above. Construct a sequence:

x(k)  x(k+1) = L x(k)

where ck  = maxj | L x(k) |j, i.e. the largest entry in the vector x(k+1)  is normalized by 1. Then, the sequence converges to the eigenvector z+, while the normalization constant ck converges to +. The dominant eigenvalue for the problem is + = 1.07622 and the eigenvector z+ is given in the table:

 Age class [0,5) [5,10) [10,15) [15,20) [20,25) [25,30) [30,35) [35,40) [40,45) [45,50) Eigenvector 1.000 0.9259 0.8588 0.7964 0.7380 0.6836 0.6328 0.5848 0.5390 0.4943

Steps in writing the MATLAB script:

1. Initialize steps (1)-(3) as in the first script.
2. Use MATLAB function "[P,D] = eig(L)" to find the matrix of eigenvectors P and the diagonal matrix of eigenvalues D. Identify the unique positive eigenvalue + and confirm that it is dominant. Compare the corresponding eigenvector z+ with the one given above.
3. Write a loop for finding x(k+1), normalizing the maximal element of the vector x(k+1) by one.
4. Run the loop for 20 values and see how close the normalization constant ck  approaches the dominant eigenvalue +.
5. Check how close the vector x(20) approaches the eigenvector z+.

Exploiting the MATLAB script:

1. Design termination of the loop when a required criterion of accuracy of the dominant eigenvalue is reached. For example if the difference |ck+1 – ck | < tol for any given tolerance tol, an infinite loop can be terminated.

QUIZ:

1)      The value of the dominant eigenvalue+ is determined by the net reproduction rate of the population:

R = a1 + a2 b1 + a3 b1 b2 + … + an b1 b2 … bn-1

If R > 1, the population increases. If R < 1, the population decreases. If R = 1, the population stabilizes. Write a function "[R,status] = NetRate(a,b)" that computes the scalar R from the vectors a and b and states in the string status whether the population increases, decreases or stabilizes.

2)      Consider the birth-death parameters for the animal population:

Age

[0,5)

[5,10)

[10,15)

aj

0

4

3

bj

0.5

0.25

---

Find the net population growth for the animal population and state whether it increases, decreases, or stabilizes. Compute the dynamics of the animal population in time and confirm the prediction.