Goals
After reading this section of notes, you should
know what is a discrete-time model,
know what is a linear difference equation,
know how to use linear algebra techniques to analyze a linear difference equation, and
be aware of some biomathematical applications of linear difference equations.
Background
Let \(N(t)\) denote the number of a single isolated and homogeneous population at time \(t\). A common question is how does the population change over time? Early on, we argued that a general model for this is
\(N(t+\Delta t) = N(t) + G(N(t),\Delta t)\)
which essentially states that the change (increase or decrease) in the population from \(t\) to \(t+\Delta t\) is determined by the current population \(N(t)\) incremented by some “growth” function \(G(N(t),\Delta t)\) which could be positive or negative. A first order approximation for \(G\) is
\(G(N(t),\Delta t) = r \Delta t N(t)\)
where we do not make any assumptions about the sign of \(r\). This leads to
\(N(t+\Delta t) = N(t) + r \Delta t N(t)\)
Previously, we rewrote the last expression as
\(\frac{N(t+\Delta t) - N(t)}{\Delta t} = r N(t)\)
and took the limit as \(\Delta t \rightarrow 0\) which lead to the exponential growth model \(\frac{dN}{dt}=rN\). This involves making an assumption. Specifically, we assume that time \(t\) can be reasonably treated as a continuous variable. There are many situations where this an appropriate assumption. On the other hand, there are situations in which it is more reasonable to view time as a discrete variable. Let’s now make this assumption in the context of our population growth model. To denote this distinction mathematically, we write \(N(t)=N_{t}\), then our simple first-order population growth model becomes
\(N_{t+\Delta t} = (1+r\Delta t)N_{t}\)
If we increment over a fixed-size time step for \(\Delta t\), then we can write the previous equation as
\(N_{t+1} = \lambda N_{t}\)
This is an example of a discrete time model. In general, a discrete time model is a sequence \(x_{0}\), \(x_{1}\), \(\ldots\), \(x_{n}\), \(x_{n+1}\), \(\ldots\), where
\(x_{i}\in \mathbb{R}^{p}\) for all \(i=0,1,\ldots\) and for some \(p\in \mathbb{Z}^{+}\), and
\(x_{i+1}=f(x_{i},x_{i-1},x_{i-2},\ldots,x_{i-k})\) for all \(i=1,2,\ldots\) and some fixed \(k\geq 0\).
Typically, we treat the \(x_{i}\) as unknowns and refer to the expression \(x_{i+1}=f(x_{i},x_{i-1},x_{i-2},\ldots,x_{i-k})\) as a difference equation. A solution of the difference equation is any sequence \(x_{0},x_{1},\ldots\) that satisfies the recurrence defined by the difference equation.
If \(k=0\), then we have
\(x_{i+1}=f(x_{i})\)
and the discrete time model is said to be first-order.
In our example model \(N_{t+1} = \lambda N_{t}\) we have \(f(N_{t})=\lambda N_{t}\) with \(\lambda\) a constant. It is obviously a first-order model. Moreover, it is linear. In general, if a discrete time model can be written as
\(x_{t+1}= A x_{t}\)
where \(x_{t}\in \mathbb{R}^{n}\) for all \(t\) and \(A\) is an \(n\times n\) matrix then we say it is a linear difference equation.
Analysis of Linear Difference Equations
Suppose we have a linear difference equation
\(x_{t+1}= A x_{t}\)
where \(x_{t}\in \mathbb{R}^{n}\) for all \(t\) and \(A\) is an \(n\times n\) matrix, and let \(x_{0}\) be some given initial vector. Then,
\[ \begin{align} x_{1} &= Ax_{0}, \\ x_{2} &= Ax_{1} \\ &= A(Ax_{0}) \\ &= A^2x_{0}, \end{align} \]
and generally,
\[ \begin{align} x_{i+1} &= Ax_{i} \\ &= A(Ax_{i-1}) = A^2x_{i-1} \\ &= A(A(Ax_{i-2})) = A^3x_{i-2} \\ &= \ldots \\ &= A^{k+1}x_{i-k} \\ &= \ldots \\ &= A^{i+1}x_{0} \end{align} \]
Thus, the solution to a linear difference equation \(x_{t+1}=Ax_{t}\) is a sequence \(x_{0}\), \(x_{1}=Ax_{0}\), \(x_{2}=A^2x_{0}\), \(\ldots\), \(x_{n}=A^{n}x_{0}\), \(\ldots\). Notice that to compute each term in the sequence, we must compute the \(n\)-fold matrix power \(A^n\). Typically, this is not convenient to compute. Our goal is to show that there is another approach that can be taken.
Eigenvalue Analysis of Linear Difference Equations
Let \(A\) be an \(n\times n\) matrix and suppose that \(A\) has a complete set of eigenvectors. That is, suppose that \(A\) has \(n\) linearly independent eigenvectors \(v_{1},v_{2},\ldots, v_{n}\) with corresponding eigenvalues \(\lambda_{1},\lambda_{2}\ldots,\lambda_{n}\). Now, let \(x_0\) be a vector in \(\mathbb{R}^{n}\). We may write
\(x_{0}=c_{1}v_{1}+c_{2}v_{2}+\cdots +c_{n}v_{n}\)
where the \(c_{i}\) are uniquely determined scalars. Further, if \(\lambda\), \(v\) is an eigenvalue-eigenvector pair for \(A\), then
\(A^n v = \lambda^n v,\)
an identity you may derive as an exercise. Using this, we see that
\(x_{i} = A^{i}x_{0}=c_{1}\lambda_{1}^i v_{1} + c_{2}\lambda_{2}^i v_{2} + \cdots + c_{n}\lambda_{n}^i v_{n}\)
and this provides another way to express a solution to a linear difference equation. The benefit is that this does not require us to compute matrix powers, only powers of scalars.
Example
Consider the linear difference equation \(x_{t+1}=Ax_{t}\) with
\(A = \left(\begin{array}{cc} 2 & 3 \\ 2 & 1 \end{array}\right),\)
the eigenvalues of \(A\) are
One can also see that if \(\lambda_{1}=4\) and \(\lambda_{2}=-1\), then \(v_{1}=(3,2)^{T}\) and \(v_{2}=(1,-1)^{T}\) are corresponding eigenvectors. Thus, the general solution to this example linear difference equation is
\(x_{i} = c_{1}(4)^{i}\left(\begin{array}{c} 3 \\2 \end{array} \right) + c_{2}(-1)^i\left(\begin{array}{c} 1 \\ -1 \end{array} \right)\)
Suppose for example that \(x_{0}=(1,1)^{T}\). To find the solution with this particular initial condition we must find scalars \(c_{1}\) and \(c_{2}\) so that
\(\left(\begin{array}{c} 1 \\1 \end{array} \right) = c_{1}\left(\begin{array}{c} 3 \\ 2 \end{array} \right) + c_{2}\left(\begin{array}{c} 1 \\ -1 \end{array} \right)\)
this can be done by solving the linear system
\(\left(\begin{array}{cc} 3 & 1 \\ 2 & -1 \end{array}\right)\left(\begin{array}{c} c_{1} \\ c_{2} \end{array} \right) = \left(\begin{array}{c} 1 \\1 \end{array} \right)\)
which is
So, the solution to \(x_{t+1}=Ax_{t}\) with
\(A = \left(\begin{array}{cc} 2 & 3 \\ 2 & 1 \end{array}\right),\)
and \(x_{0}=(1,1)^{T}\) is
\(x_{t}=\frac{2}{5}(4)^t \left(\begin{array}{c} 3 \\ 2 \end{array} \right) - \frac{1}{5}(-1)^t \left(\begin{array}{c} 1 \\ -1 \end{array} \right) = (4)^t \left(\begin{array}{c} \frac{6}{5} \\ \frac{4}{5} \end{array} \right) + (-1)^t \left(\begin{array}{c} -\frac{1}{5} \\ \frac{1}{5} \end{array} \right)\)
An interesting question about the dynamics of discrete time models is, what happens to the system as \(t\rightarrow \infty\). Note that in the expression
\(x_{i} = c_{1}\lambda_{1}^i v_{1} + c_{2}\lambda_{2}^i v_{2} + \cdots + c_{n}\lambda_{n}^i v_{n}\)
if all of the eigenvalues satisfy \(|\lambda_{i}|<1\), then \(x_{i} \rightarrow 0\) as \(i \rightarrow \infty\). Thus, if \(x_{t+1}=Ax_{t}\) is a linear difference equation, we expect \(x_{t} \rightarrow 0\) in the long-term whenever all of the eigenvalues of \(A\) have absolute value less than 1.
We can apply the result we just derived to our simple population growth model \(N_{t+1}=\lambda N_{t}\) by observing that it is a one-dimensional linear difference equation. Therefore, the general solution is \(N_{t}=c\lambda^t N_{0}\) and then the population will die out if \(|\lambda| < 1\). On the other hand, what will happen if \(|\lambda| > 1\)? In the one-dimensional case it is pretty clear, the magnitude of the population will grow without bound. What happens in the general case \(x_{t+1}=Ax_{t}\) whenever \(A\) is an \(n\times n\) matrix?
To answer this question, suppose that our eigenvalues are ordered so that \(|\lambda_1| > |\lambda_2| > \cdots > |\lambda_{n}|\) and then rewrite
\[ \begin{align} x_{t} &= c_{1}(\lambda_{1})^tv_{1} + c_{2}(\lambda_{2})^tv_{2} + \cdots + c_{n}(\lambda_{n})^tv_{n} \\ &= \lambda_{1}^t\left(c_{1} v_{1} + c_{2}\left(\frac{\lambda_{2}}{\lambda_{1}}\right)^t v_{2} + \cdots + c_{n}\left( \frac{\lambda_{n}}{\lambda_{1}}\right)^t v_{n}\right) \end{align} \]
Now, for \(j=2,\ldots n\) we have that \(\left| \frac{\lambda_{j}}{\lambda_{1}}\right| < 1\). So, we expect that as \(t\rightarrow \infty\) each of these terms will go to zero and \(x_{t}\) will asymptotically approach the line determined by the vector \(v_{1}\). For a matrix \(A\), we call the eigenvalue of \(A\) that is largest in absolute value the dominant eigenvalue. Thus, the dominant eigenvalue and corresponding eigenvector determine the long-term dynamics of a discrete dynamical system determined by a linear difference equation.
Example
Consider again the linear difference equation \(x_{t+1}=Ax_{t}\) with
\(A = \left(\begin{array}{cc} 2 & 3 \\ 2 & 1 \end{array}\right)\)
This matrix has dominant eigenvalue \(\lambda_{1}=4\) with corresponding eigenvector \(v_{1}=(3,2)^{T}\). The following figure shows the points \(x_{t}\) in the plane determined by the sequence \(x_{t}=A^{t}x_{0}\) together with the line through the vector \(v_{1}\).
Note that we have plotted on the log scale because there is exponential growth in the magnitude of points in the sequence of vectors \(x_{t}\). This plot clearly shows how the asymptotic dynamics of the discrete time model “follows” the direction of the eigenvector corresponding with the dominant eigenvalue.
Note: Of course it may be the case that a particular \(n\times n\) matrix \(A\) does not have \(n\) distinct eigenvalues. In such cases it is still possible to analyze the linear difference equation \(x_{t+1}=Ax_{t}\) but the linear algebra becomes slightly more complicated. We leave such a discussion to one of the references such as (Allen 2007) or (de Vries et al. 2006).
Application: Age-Structured Models
An important application of linear difference equations as a discrete time model for population growth is provided by the so-called Leslie matrix model. The utility of this model is that it allows us to incorporate additional biologically interesting structure into a population model. By “structure” we mean a feature or features that subsets of a population may or may not have but that plays a distinguishing role in the ecology of the species whose population growth we want to model.
For example, many insects progress through distinct, discrete developmental stages such as egg, larva, pupa, adult. The reproductive behavior and survival rate is different within each of these different developmental stages. Thus, there is a clear age-structure involved in the growth of such species. How does one incorporate such additional structure into a mathematical model for the growth of a population? It is clear that age is a temporal phenomenon but aging may occur on a different time-scale than the overall change in the size of a population. This makes it challenging to incorporate age-structure into a continuous-time model although this can be done using partial differential equations. On the other hand, it is fairly simple to develop discrete time models that include age-structure or other similar population structures. This is done classically using the Leslie matrix model.
Here is the general idea. To start with the simplest case we assume that we have a single population that is closed to migration and we only model females. The population is divided up into discrete groups, say for example into \(N\) groups. Denote by \((x_{i})_{t}\) the number of individuals in group \(i\) at time \(t\). We assume that group 1 corresponds to new offspring. Thus, new individuals are “born into” group 1. For each of the groups \(i=1,2,\ldots, N\), there is a chance of progressing from group \(i\) to group \(i+1\) determined by a value \(s_{i}\). Further, for \(i=1,2,\ldots, N\) we let \(b_{i}\geq 0\) denote the average number of newborn females produced by one female in the \(i\)-th age group that survive through the time interval in which they were born. Then,
\((x_{1})_{t+1} = b_{1}(x_{1})_{t} + b_{2}(x_{2})_{t} + \cdots + b_{N}(x_{N})_{t}\)
\((x_{i+1})_{t+1} = s_{i}(x_{i})_{t}\), for \(i=1,\ldots,N-1\)
We can rewrite this in matrix-vector notation as
\[ X_{t+1} = \left(\begin{array}{c} (x_{1})_{t+1}\\ (x_{2})_{t+1} \\ (x_{3})_{t+1} \\ \vdots \\ (x_{N})_{t+1}\end{array}\right) = \left(\begin{array}{ccccc} b_{1} & b_{2} & \cdots & b_{N-1} & b_{N} \\ s_{1} & 0 & \cdots & 0 & 0 \\ 0 & s_{2} & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & s_{N-1} & 0 \end{array}\right)\left(\begin{array}{c} (x_{1})_{t}\\ (x_{2})_{t} \\ (x_{3})_{t} \\ \vdots \\ (x_{N})_{t} \end{array}\right) = LX_{t} \]
where we call the matrix \(L\) the Leslie matrix. Here is a diagram that illustrates the Leslie matrix model idea in the case where \(N=4\):
Show code
a_graph <-
create_graph() %>%
add_n_nodes(n=4,label = c("1","2","3","4"),
node_aes = node_aes(shape="circle",color = "black")) %>%
add_edges_w_string(edges = "1->2 2->3 3->4 4->1 3->1 2->1 1->1") %>%
select_edges_by_edge_id(edges = 1:7) %>%
set_edge_attrs_ws(
edge_attr = color,
value = "black")
render_graph(a_graph, layout = "circle")
This diagram shows the relationships among the (four in the example) different age groups.
Before we proceed to consider the analysis of the dynamics of a Leslie matrix model, note that for any Leslie matrix \(L\), all of the entries will be nonnegative.
Analyzing the Leslie Matrix Model
Given an \(n\times n\) Leslie matrix \(L\), a typical problem is to determine if its coefficients are such that \(L\) has a strictly dominant eigenvalue \(\lambda_{1}\). That is, a unique eigenvalue \(\lambda_{1}\) such that \(|\lambda_{1}| > |\lambda_{j}|\) for \(j=2,3,\ldots, n\). As we have seen in a previous section, if this is the case then \(\lambda_{1}\) and its corresponding eigenvector \(v_{1}\) determine the long-term dynamics of the Leslie model. In the context of Leslie models, an eigenvector \(v_{1}\) associated with a strictly dominant eigenvalue \(\lambda_{1}\) is called a stable age distribution.
We suspect from our previous analysis of linear difference equations that if \(\lambda_{1}\) is a strictly dominant eigenvalue, then
if \(|\lambda_{1}| < 1\), then the overall population size will decrease over time, and
if \(|\lambda_{1}| > 1\), then the population increases over time in the direction of the stable age distribution \(v_{1}\).
Note that if \(L\) is a Leslie matrix with a strictly dominant eigenvalue \(\lambda_{1}\), then the entries of the stable age distribution \(v_{1}\) will be positive.
Much theory has been developed for a general analysis of Leslie matrix models, see section 1.7 of (Allen 2007) for the details. Here we will simply consider some examples.
Example 1
Suppose that we have a population with two age groups where only one of the two ages is reproducing. Then the Leslie matrix model has
\(L = \left(\begin{array}{cc} 0 & b_{2} \\ s_{1} & 0 \end{array} \right)\)
where \(b_{2}\) and \(s_{1}\) are positive real values. The eigenvalues for this \(L\) are obtained by solving the quadratic polynomial
\(\lambda^{2} - 0 \lambda + (-b_{2}s_{1}) = 0\)
and thus are \(\lambda_{1,2}=\pm\sqrt{b_{2}s_{1}}\)
In this case, there is no strictly dominant eigenvalue because \(|\lambda_{1,2}|=\sqrt{b_{2}s_{1}}\) and there is not a unique eigenvalue with the largest magnitude.
Example 2
Suppose that we have a population with two age groups where both are reproducing. Then the Leslie matrix model has
\(L = \left(\begin{array}{cc} b_{1} & b_{2} \\ s_{1} & 0 \end{array} \right)\)
where \(b_{1}\), \(b_{2}\), and \(s_{1}\) are positive real values. The eigenvalues for this \(L\) are obtained by solving the quadratic polynomial
\(\lambda^{2} - b_{1} \lambda + (-b_{2}s_{1}) = 0\)
and thus are
\(\lambda_{1,2}=\frac{b_{1} \pm \sqrt{b_{1}^2 + 4b_{2}s_{1}}}{2}\)
Thus, there is a strictly dominant eigenvalue \(\lambda_{1} = \frac{b_{1} + \sqrt{b_{1}^2 + 4b_{2}s_{1}}}{2}\) and we have that \(\lambda_{1} > 0\). Furthermore,
\(\lambda_{1} = \frac{b_{1} + \sqrt{b_{1}^2 + 4b_{2}s_{1}}}{2} > \frac{b_{1} + \sqrt{b_{1}^2}}{2} = b_{1}\)
from which we see that if \(b_{1} > 1\), then \(\lambda_{1} > 1\). Note however that this is not the only situation in which we can have \(\lambda_{1} > 1\).
Let’s find the eigenvector \(v_{1}\) corresponding to \(\lambda_{1} = \frac{b_{1} + \sqrt{b_{1}^2 + 4b_{2}s_{1}}}{2}\). Such a vector will satisfy
\[ (L-\lambda_{1} I)v_{1} = \left(\begin{array}{cc} \frac{b_{1}-\sqrt{b_{1}^2+4b_{2}s_{1}}}{2} & b_{2} \\ s_{1} & \frac{-b_{1}-\sqrt{b_{1}^2+4b_{2}s_{1}}}{2} \end{array} \right)\left(\begin{array}{c} x \\ y \end{array} \right) = \left(\begin{array}{c} 0 \\ 0 \end{array}\right) \]
and you can show (for homework) that
\(v_{1} = \left(\begin{array}{c} 1 \\ \frac{s_{1}}{\lambda_{1}} \end{array}\right) = \left(\begin{array}{c} 1 \\ \frac{2s_{1}}{b_{1} + \sqrt{b_{1}^2 + 4b_{2}s_{1}}} \end{array}\right)\)
satisfies the above equation. Thus,
\(L = \left(\begin{array}{cc} b_{1} & b_{2} \\ s_{1} & 0 \end{array} \right)\)
has strictly dominant eigenvalue \(\lambda_{1} = \frac{b_{1} + \sqrt{b_{1}^2 + 4b_{2}s_{1}}}{2}\) with stable age distribution
\(v_{1} = \left(\begin{array}{c} 1 \\ \frac{s_{1}}{\lambda_{1}} \end{array}\right) = \left(\begin{array}{c} 1 \\ \frac{2s_{1}}{b_{1} + \sqrt{b_{1}^2 + 4b_{2}s_{1}}} \end{array}\right)\)
Example 3
Consider the following three-stage Leslie matrix
\[ L = \left(\begin{array}{ccc} 0 & 9 & 12\\ \frac{1}{3} & 0 & 0 \\ 0 & \frac{1}{2} & 0 \end{array}\right) \]
We can compute the eigenvalues and eigenvectors for \(L\) with R:
eigen() decomposition
$values
[1] 2+0i -1+0i -1-0i
$vectors
[,1] [,2] [,3]
[1,] -0.98556187+0i 0.9370426+0i 0.9370426+0i
[2,] -0.16426031+0i -0.3123475-0i -0.3123475+0i
[3,] -0.04106508+0i 0.1561738+0i 0.1561738-0iFrom this, we see that \(L\) has a strictly dominant eigenvalue \(\lambda_{1}=2\) with corresponding stable stage distribution -0.985561874823818+0i, -0.164260312470636+0i, -0.0410650781176591+0i. Let’s simulate the dynamics of \(X_{t+1}=LX_{t}\) for many times-steps and compare the last values with \(v_{1}\):
[,1]
[1,] 2.929681e+31
[2,] 4.882802e+30
[3,] 1.220701e+30We see that the population grows over time. Do the population values grow in the direction of \(v_{1}\)? In order to determine this, we can determine the distance between our last iterate and the line through \(v_{1}\), or equivalently, we can compute the angle between the two vectors.
[1] 3.141593+0iIf the angle between two vectors is 0 or \(\pi\) then they lie on the same line. We see that our two vectros are essentially parallel. We can confirm this by normalizing each of the vectors:
[1] -0.98556187+0i -0.16426031+0i -0.04106508+0i [,1]
[1,] 0.98556187
[2,] 0.16426031
[3,] 0.04106508These vectors are obviously nearly parallel. More examples and problems related to Leslie matrix models will be considered in the homework exercises.
Conclusion
We’ve introduced linear difference equations as an example of discrete time systems and derived the role that the eigenvalues and eigenvectors play in analyzing the dynamics of such systems. Notice that much of this analysis closely parallels the analysis of continuous-time linear autonomous systems. An important example of linear difference equations in biomathematics is the Leslie matrix model approach to representing the population dynamics for a structured population. Next, we will study another class of discrete time systems by looking at first-order nonlinear difference equations. As we will see, even very simple nonlinear difference equations can exhibit very complex dynamics.
Further Reading
For more on discrete time models and linear difference equations in particular, see chapter 1 from either (Allen 2007) or (de Vries et al. 2006). Both of these books also contain references to further reading on the topics of difference equations and their applications within biomathematics.An excellent general reference for structured population models is (Cushing 1998).