Math 236  •  1 May 2001

Spotted owls and flying squirrels: 
a discrete linear dynamical system
and its analysis using 
eigenvalues and eigenvectors

Initialization cells: utility print and plot functions

The functions in the cells below are utilities for printing and plotting sequences.  These cells are initialization cells and should have been evaluated automatically when you opened this notebook—provided you agreed then to let them be evaluated.  If you did not, be sure to evaluate them now, before you proceed.

Off[General::spell1];
PrintVectorSequence::usage  =  "Print2VectorSequence[seq, a, b] \
prints a table, with labeled columns, of {a, seq[a][[1\
]], seq[b][[2]]}, ... \
{b, seq[b][[1]], seq[b]\
[[2]]}."; 
PlotVectorSequence::usage  =  "PlotVectorSequence[vector, first, \
number] up to dim 5";
On[General::spell1];

\(PrintVectorSequence[seq_, \ a_, \ b_]\  := \ 
      Module[{k}, \ 
        TableForm[Table[Flatten[{k, \ x[k]}], \ {k, \ a, \ b}], \ 
          TableSpacing\  -> \ {1, \ 6}, \ 
          TableHeadings\  -> \ {None, \ {"k", \ \*""\!\(\(x\_k\)\&\(\(\(\\\)\(\\ \)\) -> \(\(\\\)\(\\ \)\)\)\)[[1]]"", \ \*""\!\(\(x\_k\)\&\(\(\(\\\)\(\\ \)\) -> \(\(\\\)\(\\ \)\)\)\)[[2]]""}}]]\)

Needs["Graphics`Colors`"]

The following function was written by Frank Wattenberg.

PlotVectorSequence[vect_,  a_,  n_]  :=  
      Module[{dim,  k,  j,  largest,  curve,  dots,  CurveColor},  
        dim  =  Length[vect[1]];  
largest  =  
          Max[Table[vect[k][[j]],  {k,  a,  n},  {j,  1,  
                dim}]];  
dots[j_]  :=  
          Map[Point,  
            Table[{k,  vect[k][[j]]},  {k,  a,  
                n}]];  
curve[j_]  :=  
          Line[Table[{k,  vect[k][[j]]},  {k,  a,  
                n}]];  
CurveColor[1]  =  Red;  
        CurveColor[2]  =  Green;  CurveColor[3]  =  Blue;  
        CurveColor[4]  =  Black;  
        CurveColor[5]  =  Purple;  
Show[
          Graphics[{PointSize[0.02],  
              Map[Join,  Table[dots[j],  {j,  1,  dim}]]}],  
          Map[Graphics,  
            Table[{CurveColor[j],  curve[j]},  {j,  1,  dim}]],  
          AspectRatio  ->  1/GoldenRatio,  
          Axes  ->  Automatic,  
          PlotRange  ->  {0, largest}]]

The model

Time will be discrete (t = 0, 1, 2, 3, ...), not continuous.

t = time (years) from initial time
o (t) = quantity of spotted owls at time t
s (t) = quantity of flying squirrels at time t

Initial data: 
			o(0)  =  (Subscript[o, 0]) (given number)
			s(t)  = (Subscript[s, 0])  (given number)
Dynamics:
			   \(o(
              t\  + 1)\  = \ \ \ \ \ \ \ \ \ \ \ \ \ 0.4\ \(o(t)\)\  + \ 
              0.3\ \(s(t)\), 

          s(t\  + \ 1)\  = \ \ \(-0.325\)\ \(o(t)\)\  + \ 1.2\ \(s(t)\)\)

State vector at time t:	   (Overscript[x(t), ->]) = (({{o(t)}, {s(t)}}))

Matrix form of system:	(Overscript[(x(t  +  1)),  ->])  = (({{0.4, 0.3}, {-0.325, 1.2}})) (Overscript[(x(t)),  ->])  = A (Overscript[(x(t)),  ->])

A = (({{0.4, 0.3}, {-0.325, 1.2}}));

a  =  Rationalize[A];  a  //  MatrixForm

Examples of three different initial states

Example 1

Clear[x0,  x]

\(x0\  = \ {60. , \ 130. }; \ \ \ \ \ \ \ \  (*\ initial\ state\ *) 

    x[0]\  = \ x0;\)
x[t_]  :=  a  .  x[t  -  1]

Another way to generate successive state vectors:  x[t_] := MatrixPower[A, t] . x0

x[1]

{x[1], x[2],  x[3]}

Table[x[t],  {t,  0,  10}]  //  TableForm

The function PrintVectorSequence, defined in an initialization cell near the start of this notebook, prints a table that includes the times and labels the columns:

PrintVectorSequence[x,  0,  10]

What's happening?

x[1]  -  x[0]
((x[1]  -  x[0])/x)[0]

x[1]  ==  1.05 x[0]

x[2]  ==  1.05 x[1]

Both populations grow by 5% each year.

Visualize as follows.

PlotVectorSequence is a function written by Frank Wattenberg's that appears in the section "Initialization cells" at the start of this notebook.

First plot both populations (vertical axis) against time (horizontal axis):

PlotVectorSequence[x,  0,  20];

Visualize the joint behavior via a parametric plot in phase space—squirrels (vertical axis) vs. owls (horizontal axis):

ListPlot[
        Table[x[t],  {t,  0,  10}],      PlotStyle  ->  {PointSize[0.025]},  
        AxesLabel  ->  {"owls",  "squirrels"}];

Slope there gives ratios of components of state vectors.  Look at ratios of the two populations:

o[
          t_]  :=  x[
            t][[1]];
s[
          t_]  :=  x[
            t][[2]];
Table[o[t]/s[t],  {t,  0,  10}]

Compare:

x0[[1]]/x0[[2]]

A  .  x0
1.05 x0

A  .  x0  ==  1.05 x0

x0 = {60., 130.} is an eigenvector of A.
1.05 is the eigenvalue of A associated with this eigenvector.

Example 2

Clear[x0,  x]

\(x0\  = \ {200. , \ 100. }; \ \ \  (*\ new\ initial\ condition\ *) 

    x[0]\  = \ x0;\)
x[t_]  :=  a  .  x[t  -  1]

PrintVectorSequence[x,  0,  10]

What's happening?

x[1]  -  x[0]
((x[1]  -  x[0])/x)[0]

Both populations decline by 45 % per year.

x[1]  ==  0.55  x[0]

x[2]  ==  0.55 x[1]

PlotVectorSequence[x,  0,  10];

ListPlot[
        Table[x[t],  {t,  0,  10}],     PlotStyle  ->  {PointSize[0.025]},  
        AxesLabel  ->  {"owls",  "squirrels"}];

Look at ratios of the two populations:

o[
          t_]  :=  x[
            t][[1]];
s[
          t_]  :=  x[
            t][[2]];
Table[o[t]/s[t],  {t,  0,  10}]

Compare:

x0[[1]]/x0[[2]]

A . x0
0.55 x0

A  .  x0  ==  0.55 x0

x0 = {200., 100.} is an eigenvector of A.
0.55 is the eigenvalue of A associated with this eigenvector.

Example 3

Clear[x0,  x]

x0  =  {100. ,  175. };
x[0]  =  x0;
x[t_]  :=  a  .  x[t  -  1]

PrintVectorSequence[x,  0,  10]

What's happening?

x[1]  -  x[0]
((x[1]  -  x[0])/x0)

x[2]  -  x[1]
((x[2]  -  x[1])/x)[1]

No such simple relationship as in examples 1 and 2!

PlotVectorSequence[x,  0,  10];

ListPlot[
        Table[x[t],  {t,  0,  10}],     PlotStyle  ->  {PointSize[0.02]},  
        AxesLabel  ->  {"owls",  "squirrels"}];

o[
          t_]  :=  x[
            t][[1]];
s[
          t_]  :=  x[
            t][[2]];
Table[o[t]/s[t],  {t,  0,  20}]

Why is that the asymptotic (long-run) ratio?

Eigenvalue/eigenvector analysis of Example 3

From Example 1 and Example 2, the coefficient matrix A has the following two eigenvectors and associated eigenvalues:

v1  =   {60. ,  130. };  λ1  =  1.05;
v2  =   {200. ,  100. };  λ2  =  0.55;

A (Overscript[(Subscript[v, 1]), ->]) =  (Subscript[λ, 1]) (Overscript[(Subscript[v, 1]), ->]),
A (Overscript[(Subscript[v, 2]) , ->])=  (Subscript[λ, 2]) (Overscript[(Subscript[v, 2]), ->])

The two eigenvectors are linearly independent, hence form a basis of (ℝ^2).  Write the initial state vector in Example 3 as a linear combination (Overscript[((Subscript[x, 0])),  ->])   =  
        (Subscript[c, 1]) (Overscript[((Subscript[v, 1])),  ->])   +  
          (Subscript[c, 2]) (Overscript[((Subscript[v, 2])),  ->]) of this basis:

soln  =  
      First[ Solve[
          x0  ==  c1 v1  +  c2 v2,  {c1,  c2}] ]

Thus eigenvalue/eigenvector analysis yields the explicit formula for (Overscript[(x(t)),  ->])  here:
			(Overscript[(x(
              0)),  ->])   =  (Overscript[((Subscript[x, 0])),  ->]) ,
			\(TraditionalForm\`\(\(x(
              t)\)\& -> \)\  = \ \(c\_1\ A\^t\ \(\(v\_1\)\& -> \)\  + \ c\_2\ A\^t\ \(\(v\_2\)\& -> \)

        \ \ \ \ \ \ \ \ \ \(\(=\)\(\ \)\(c\_1\ λ\_1\^t\ \(\(v\_1\)\& -> \)\  + \ c\_2\ λ\_2\ \(\(v\_2\)\& -> \)

          \(\(=\)\(\ \)\(\((1.05)\)\^t\ c\_1\ \(\(v\_1\)\& -> \)\  + \ \((0.55)\)\^t\ c\_2\ \(\(v\_2\)\& -> \)\)\)\)\)\)\)

Look at the second (vector) term:

Table[(((0.55))^t)  c2 v2  /.  soln,  {t,  0,  20}]

Since 0  <  0.55  <  1, then (lim_( t  ->  ∞)) 
          (0.55^t)  =  0 and so 
	(lim_( t  ->  ∞)) 
          (0.55^t) (Subscript[c, 2]) (Overscript[((Subscript[v, 2])),  ->])   =  (Overscript[0,  ->]) .
Hence for large t,
	\(TraditionalForm\`\(\(x(
              t)\)\& -> \)\  ≈ \ \ \((1.05)\)\^t\ c\_1\ \(\(v\_1\)\& -> \)\).

In other words, as t increases, the coefficient (((0.55))^t) becomes very small and hence so does the second term of \(TraditionalForm\`\(\(x(
              t)\)\& -> \)\  = \ \((1.05)\)\^t\ c\_1\ \(\(v\_1\)\& -> \)\  + \ \((0.55)\)\^t\ c\_2\ \(\(v\_2\)\& -> \)\).  Thus the behavior of the two populations is really determined by the first term—the one involving the "dominant" eigenvalue 1.1.

Now we can account for the long-run behavior of the ratio between the two populations sizes.  Recall the ratios:

Table[o[t]/s[t],  {t,  0,  20}]

Now we know:
	\(TraditionalForm\`\(o(t)\)\/\(s(t)\)\  = \ \(\((\(\(x(t)\)\& -> \))\)\_1\/\((\(\(x(t)\)\& -> \))\)\_2\  ≈ \ \(\(1.05\^t\) \(c\_1\) \(\(\(\ \)\(\((v\_1)\)\_1\)\)\& -> \)\)\/\(\(1.05\^t\) \(c\_1\) \(\(\(\ \)\(\((v\_1)\)\_2\)\)\& -> \)\)\  = \ \(\(\(\ \)\(\((v\_1)\)\_1\)\)\& -> \)\/\(\(\ \)\(\((v\_1)\)\_2\)\)\)\)

The ratio \(TraditionalForm\`\(\(\(\ \)\(\((v\_1)\)\_1\)\)\& -> \)\/\(\(\ \)\(\((v\_1)\)\_2\)\)\) is:

v1[[1]]/v1[[2]]

And that ratio is the asymptotic, or long-run, ratio—the limit of the ratio o(t)/s(t)—found experimentally above.

Finding eigenvalues and eigenvectors with Mathematica

Eigensystem[a]

(a  -  λ IdentityMatrix[2])  .  {60, 
          130}  /.  λ  ->  21/20

(a  -  λ IdentityMatrix[2])  .  {200, 
          100}  /.  λ  ->  11/20

If you want just the eigenvalues:

Eigenvalues[a]
Eigenvalues[A]

Reference

The dynamical system studied here appears as Exercise 5, Section 5.6, of David C. Lay, Linear Algebra and Its Applications, 2nd ed.Addison-Wesley, Reading, Mass., 2000, page 346.