Chapter 3

Problems #1-9, pp. 135-139

Problem 1, p. 135

The three-year animal problem. Part a) asks us to draw the state diagram for this situation.


			A--->> .8 -->>	B--->> .9 -->>		C
			^			|			|
			|			|			|
			|_______.8_________|			|
			|_________________3.5________________|

Part a) asks us to construct the leslie matrix:

A = {{0, .8, 3.5}, {.8, 0, 0}, {0, .9, 0}}A  // MatrixForm

The Eigenvalue[] function computes the eigenvalues of A:

eigs = Eigenvalues[A]

The Abs[] function computes the norms of the eigenvalues of A:

Abs[eigs]

To compute the dominant eigenvalue, we need to find the modulus of each eigenvalue.  The dominant eigenvalue is the positive real one, namely 1.51697.  This tells us that the animal, in all its states (including seeds) is growing at a rate of  51.697% (i.e. (1.51697-1)100%).  The Abs[] command finds the absolute value of real numbers and the modulus of complex numbers.

Part d). Track the population for 10 years, given a start of 100 females in the first year, and enough males to keep all the females as pregnant as possible!

Caution:  Matrix multiplication is symbolized with a period.  A space of an asterix is between matrices is an operation, but not matrix multiplication.  The MatrixPower[] command is an efficient way to comput a matrix product.  

Matrix multiplication between a vector and a matrix is very forgiving in Mathematica.  Technically we should be transposing the row  vectors to make them column vector, but Mathematica interprets them the only way that makes sense.  Further, a vector as a list (not a list containing one list) is acceptable.

Clear[X]; X[0] = {100, 0, 0}X[t_] := MatrixPower[A, t].X[0]

After 100 years we have the following distribution. (Presumably the model is not valid for this length of time, but we are verifying the mathematics here).

X[10]

Finally we compute the fraction in each class.  Compare this to the eigenvector associated with the dominate eigenvector (scaled to have sum equal to one).

X[10]
        /((∑_(i = 1)^3) X[10][[i
          ]])

Next we compute the stage distribution vector over a ten year period.  Each row of the table is a year, the columns correspond to the stages.

TableForm[Table[X[t], {t, 0, 10}]]

Part e). The eigenvectors of A are found next.

m = Eigenvectors[A]

The eigenvector associated with the dominant eigenvalue is the first of this list.

m[[1]]

We divide by the sum of the entries of the eigenvector.

m[[1
        ]]/((∑_(i = 1)^3)
           m[[1
            ]][[i
          ]])

This tells us that the eventually 54% of the animal is in first year state, 29% is in second year, and 17% is in the terminal state.  If we check back on our calculations after 100 years, things appear to have reached this state:

X[100]
        /((∑_(i = 1)^3) X[100][[i
          ]])

Problem 2, p. 136

In this problem we're trying to demonstrate that a model is invalid.

A = {{12, 26, 6}, {.3, .92, 0}, {0, .18, 0.67}}A  // MatrixForm

X[0] = {1696, 485, 82}TableForm[Table[X[t], {t, 0, 5}]]

Obviously we'll all soon be smothered by these invading creatures! What's with the eigenvalues? We seem to be growing at a rather astronomical rate:

It's clear that this model is wrong, however, as we started with one 485 mature trees, and after 250 years have 332702 old trees (which could only have come from the mature population). So, unless the model includes a massive transplantation effect, about which we were uninformed, it is completely wrong!

Eigenvalues[A]

One guess as to their problem is that there's a decimal point missing in the first row:

A = {{.12, .26, .6}, {.3, .92, 0}, {0, .18, 0.67}}A  // MatrixForm

Results then look much more sensible: an 8% growth rate, rather than 1126% per 50 years!

Eigenvalues[A]TableForm[Table[X[t], {t, 0, 15}]]

Problem 3, p. 136

Compute a few traces and determinants.

\(A = {{1, 3}, {0, 2}}; 
A\  // MatrixForm\)Det[A]Sum[ A[[i,  i]],  {i,  2} ]

\(A = {{1, 2, 3}, {2, 3, 1}, {1, 2, 1}}; 
A\  // MatrixForm\)Det[A]Sum[ A[[i,  i]],  {i,  3} ]

\(A = {{1, 2, 3}, {2, 3, 1}, {3, 5, 4}}; 
A\  // MatrixForm\)Det[A]Sum[ A[[i,  i]],  {i,  3} ]

The third is singular; the others aren't.

Problem 4, p. 137

Finding eigenvalues and eigenvectors.

<< LinearAlgebra`Orthogonalization`\(A = {{3, 2}, {\(-1\), 0}}; 
A\  // MatrixForm\)Det[A]

Since the  matrix has a non-zero determinant, it has two non-zero eigenvalues:

Det[A - L*IdentityMatrix[2]]solns = Solve[Det[A - L*IdentityMatrix[2]] == 0, L]

\(lam1 = 1; 
lam2 = 2; \)

A - 1*IdentityMatrix[2] // MatrixForm

Hence {-1,1} is an eigenvector  of eigenvalue 1 (corresponding to an eigenvalue of 0 in the shifted matrix):

Eigenvalues[A - 1*IdentityMatrix[2]]Eigenvectors[A - 1*IdentityMatrix[2]]v1 = Eigenvectors[A - 2*IdentityMatrix[2]][[1]]v1 = Normalize[v1]

A - 2*IdentityMatrix[2]

Eigenvalues[A - 2*IdentityMatrix[2]]Eigenvectors[A - 2*IdentityMatrix[2]]v2 = Eigenvectors[A - 2*IdentityMatrix[2]][[2]]v2 = Normalize[v2]

Hence {-2,1} is an eigenvector of eigenvalue 2 (corresponding to an eigenvalue of 0 in the shifted matrix):

\(\(m = Transpose[Eigenvectors[A - 2*IdentityMatrix[2]]]
\)\)

s = LinearSolve[m, {1, 3}]

These are the coefficients {a,b} of the linear combination for the vector {1,3}:
			{1,3} = 7*{-1,1}+-4*{-2,1}

Here are the coordinates of the vector {1,3}

s = LinearSolve[Transpose[{v1, v2}], {1.0, 3.0}]

Now we need to verify that the table equations are equivalent:

X[0] = {1, 3}TableForm[Table[X[t], {t, 0, 2}]]TableForm[Table[s[[1]]*lam1^n*v1 + s[[2]]*lam2^n*v2, {n, 0, 2}]]

We want to prove by induction that

	(A^n)X(0) = a  ((Subscript[λ, 1])^n) (Subscript[v, 1]) + b ((Subscript[λ, 2])^n) (Subscript[v, 2])
	
We've already verified it in the anchoring cases above. Now we consider it true for all integers less than n. Hence

	(A^(n - 1))X(0) = a  ((Subscript[λ, 1])^(n - 1)) (Subscript[v, 1]) + b  ((Subscript[λ, 2])^(n - 1)) (Subscript[v, 2])

We now multiply on the left and right by A, and get

	(A^n)X(0) = a  ((Subscript[λ, 1])^(n - 1)) A(Subscript[v, 1]) + b  ((Subscript[λ, 2])^(n - 1)) A(Subscript[v, 2])
		   = a  ((Subscript[λ, 1])^(n - 1)) ((Subscript[λ, 1]) (Subscript[v, 1])) + b ((Subscript[λ, 2])^(n - 1)) ((Subscript[λ, 2]) (Subscript[v, 2]))
		    = a ((Subscript[λ, 1])^n) (Subscript[v, 1]) + b ((Subscript[λ, 2])^n)(Subscript[v, 2])
		    
Done! QED.

Problem 5, p. 138

Show that if U and V are eigenvectors of the matrix A, associated with eigenvalue λ, then
cU+dV is also an eigenvector of A with eigenvalue λ. 

Simply multiply cU+dV by A, and check the result:

A(cU+dV) = AcU+AdV= cAU+dAV= cλU+dλV= λcU+λdV= λ(cU+dV)

QED.

Problem 6, p. 138

A copy machine in two states, working or broken. Part a is to draw the state diagram, which I will suppress. Part b. is to formulate the transition matrix:

\(A = {{.7, .5}, {.3, .5}}; 
A\  // MatrixForm\)

Part c. We want to look at the probabilities for the next month or so:

X[0] = {1, 0}TableForm[
      Table[X[t]
            /((∑_(i = 1)^2) X[t][[i
              ]]), {t, 0, 30}]]

Probability working tomorrow: 0.7; day 2: .64; day 7: 0.625005, which is essentially the limit (.625). This is the answer to part d:

Eigenvalues[A]eigs = Eigenvectors[A]eigs[[1
        ]]/((∑_(i = 1)^2)
           eigs[[1
            ]][[i
          ]])

Problem 7, p. 138

A copy machine in three states, working, broken but fixable, broken for good. Part a is to draw the state diagram, which I will suppress (broken but fixable is obviously absorbing). Part b. is to formulate the transition matrix:

\(A = Transpose[{{.699, .3, .001}, {.498, .5, .002}, {0, 0, 1}}]; 

    A\  // MatrixForm\)

Part c. For fun, let's look at the probabilities for the next month or so, given that we start with a working machine:

X[0] = {1, 0, 0}TableForm[
      Table[X[t]
            /((∑_(i = 1)^2) X[t][[i
              ]]), {t, 0, 30}]]

Part c. Asks us to compute the fundamental matrix and interpret the results.

TakeColumns[TakeRows[A,{1,2}],{1,2}]
should select the sub-matrix, but it doesn't work.

\(<< \ LinearAlgebra`MatrixManipulation`; 
A\  // MatrixForm\)Fun = {{.699, .498}, {.3, .5}}Fundamental = Inverse[IdentityMatrix[2] - Fun]Fundamental  // MatrixForm

A machine starting in good repair will spend  (on average)  454.545 days in good repair, 272.727 days broken but repairable, and will then be toast. Results are only slightly different if the machine starts broken: then it will spend (on average) 452.727 days in good repair, 273.636 days in the broken but repairable state, then will be toast. Total average lifetime of a machine:

t1 = Sum[ Fundamental[[i, 1]],  {i,  2} ]t2 = Sum[ Fundamental[[i, 2]],  {i,  2} ]

Problem 8, p. 138

Show that if

\(TraditionalForm\`\(												p\_x\)\) = ((Subscript[l, x + 1])/(Subscript[l, x]))	and \(TraditionalForm\`\(		q\_x\)\) = ((Subscript[d, x])/(Subscript[l, x]))

then  (Subscript[q, x]) = 1 - (Subscript[p, x]).

	On page 134, we see that
	
 (Subscript[d, x]) =(Subscript[l, x])-(Subscript[l, x + 1]). (actually the book has it backwards). Thus
 
 \(TraditionalForm\`\(		q\_x\)\) = (((Subscript[l, x])-(Subscript[l, x + 1]))/((Subscript[l, x])))
 so when we sum 
 
 \(TraditionalForm\`\(	q\_x\)\)+(Subscript[p, x]) = ((Subscript[l, x])/(Subscript[l, x]))=1,
 
 QED.

Problem 9, p. 138

We are to calculate some probabilities, using the life tables.

Part a. The probability that a newborn survives to age 5. This is 5p0, which we compute as
((Subscript[l, 0 + 5])/(Subscript[l, 0]))=((Subscript[l, 5])/(Subscript[l, 0]))=(99164/100000)=

.99164

Part b. The probability that a 4-year old survives to age 5. This is  (Subscript[p, 4]), which is stated as .99974, but which Jack actually calculated to show that it should have been rounded to .99975:

99164/99189.0

Part c. The probability that a newborn dies between the ages of 5 and 6. This is  the sum of the number dead in the fifth year divided by the initial cohort:

22/100000*1.01 - 99142/99164.0

Part d. The probability that a newborn dies between the ages of 5 and 10. This is  the sum of the number dead in the fifth through ninth years divided by the initial cohort:

(22 + 18 + 16 + 14 + 13)/100000*1.0.99164*.000837

Part e. The probability that a 5-year old survives to age 10. This is

1.0 - (22 + 18 + 16 + 14 + 13)/99164

or 5p5:

99081/99164.0

Part f. The probability that a 5-year old dies between the ages of 10 and 13. This is

(15 + 15 + 17)/99164.0.999163*(1  -  99034/99081)99081/99164.0  -  99034/99164.0

Project 3.2, p. 140

Part 1 asks us to construct the state diagram, which I'll leave to your imaginations. In Part 2,

A = {{.3, .2, .1}, {.6, .5, .1}, {.1, .3, .8}}A  // MatrixForm

The Eigenvalue[] function computes the eigenvalues of A:

eigs = Eigenvalues[A]

Clear[X]; X[0] = {.1, .1, .8}X[t_] := MatrixPower[A, t].X[0]

After 2 and 10 years we have the following distributions:

X[2]X[10]

Finally we compute the fraction in each class.  Compare this to the eigenvector associated with the dominate eigenvector (scaled to have sum equal to one).

m = Eigenvectors[A]

The eigenvector associated with the dominant eigenvalue is the first of this list.
We divide by the sum of the entries of the eigenvector.

m[[1
        ]]/((∑_(i = 1)^3)
           m[[1
            ]][[i
          ]])

This tells us that eventually, at any given time, 16.3% of the areas will be bare dunes; 30.2% will be grass; and 53.5% will be oak. The limiting matrix will have this vector for each of its three columns.