Applications of Eigenvalues and Eigenvectors¶
With the concepts related to eigenvalues and eigenvectors in place, we return to examine Discrete Dynamical Systems.
Infectious disease model¶
For the \(SIRS\) model of infectious disease, we had the following discrete dynamical system.
Given an initial condition \(X_0\), we know that \(X_n = A^nX_0\). We are interested in determining the behavior of the system for large \(n\). We might now recognize that this calculation is exactly the same as the Power Method from the previous section, and therefore expect that the sequence of vectors produced should tend toward the eigenvector corresponding to the largest eigenvalue.
In this case the components of the vector have individual meaning, so let’s calculate the first \(20\) iterations and plot \(s_t\), \(i_t\), and \(r_t\) to get a sense of how they are changing. For this calculation we store each vector \(X_t\) as a column in an array named \(\texttt{results}\).
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
A = np.array([[0.95, 0, 0.15],[0.05,0.8,0],[0,0.2,0.85]])
## T is final time
T = 20
## X at time 0
X = np.array([[0.95],[0.05],[0]])
## The first column of results contains the initial values
results = np.copy(X)
for i in range(T):
X = A@X
results = np .hstack((results,X))
## t contains the time indices 0, 1, 2, ..., T
t = np.linspace(0,T,T+1)
## s, i, r values are the rows of the results array
s = results[0,:]
i = results[1,:]
r = results[2,:]
fig,ax = plt.subplots()
## The optional label keyword argument provides text that is used to create a legend
ax.plot(t,s,'b+',label="Susceptible");
ax.plot(t,i,'rx',label="Infectious");
ax.plot(t,r,'g+',label="Recovered");
ax.set_ylim(0,1.1)
ax.grid(True)
ax.legend();
Based on the calculation it appears that the state of the population has reached an equilibrium after 20 weeks. In the equilibrium state, each category of the population, \(S\), \(I\), and \(R\), have as many individuals entering the category as leaving it. In terms of the matrix equation, if \(X\) is the vector that contains the equilibrium values of \(s_t\), \(i_t\), and \(r_t\), then \(X\) must be a solution to the equation \(AX=X\), since \(X_{t-1}=X_t\) when the components are not changing. The equation \(AX=X\) implies that \(X\) is an eigenvector of \(A\) corresponding to an eigenvalue of one.
The equation \(AX=X\) represents a linear system, even though it is not written in the typical form. We can rearrange it to look more familiar by using the tools of matrix algebra. The key is to recognize that \(X\) can be written as \(IX\), where \(I\) is the identity matrix. This gives \(AX=IX\), which can be rewritten as \((A-I)X = 0\). This matrix equation is the type with which we have been working throughout the previous chapters.
Before attempting to solve the system \((A-I)X = 0\), we first recognize it as a Homogeneous System and apply some reasoning based on our knowledge of linear systems. We know that the solution set of a homogeneous system certainly contains the trivial solution \(X=0\). Therefore if there is also a nonzero solution to \((A-I)X=0\), the system must have an infinite number of solutions since linear systems always have zero, exactly one, or an infinite number of solutions. If the system \((A-I)X=0\) does not have a unique solution there must be a free variable, which means that one of the columns of the matrix \(A-I\) will not have a pivot. In this case the direct methods of solution based on \(LU\) or \(QR\) factorizations will fail. We must go back to the basic row operations and calculate the RREF of the augmented matrix \([(A-I)|0]\).
import laguide as lag
I = np.eye(3)
ZERO = np.zeros((3,1))
augmented_matrix = np.hstack((A-I,ZERO))
reduced_matrix = lag.FullRowReduction(augmented_matrix)
print(reduced_matrix)
[[ 1. 0. -3. 0. ]
[ 0. 1. -0.75 0. ]
[ 0. 0. 0. 0. ]]
In the reduced system for the equilibrium values \(s\), \(i\), and \(r\), we can take \(r\) as the free variable and write \(s=3r\) and \(i=0.75r\). For any value of \(r\), a vector of the following form is an eigenvector for \(A-I\), corresponding to the eigenvalue one.
Although this eigenvector is not unique, there is one more constraint on our equilibrium values that we have not yet addressed. The variables \(s\), \(i\), and \(r\) represent percentages of the population in each category, therefore we must also require that \(s+i+r = 1\) in order for the results to be meaningful in the application. This gives the equation \(3r + 0.75r + r = 1\), which gives the unique equilibrium values of \(r = 4/19\), \(s=12/19\), and \(i=3/19\). If we carry out a large number of iterations, we see that the computed values are very close to the theoretical equilibrium values.
## T is final time
T = 100
## X at time 0
X = np.array([[0.95],[0.05],[0]])
for i in range(T):
X = A@X
print("Expected values of s, i, r:")
print("r =",12/19,"i =", 3/19,"i =", 4/19,'\n')
print("Computed values of s, i, r at time ",T,":")
print(X)
Expected values of s, i, r:
r = 0.631578947368421 i = 0.15789473684210525 i = 0.21052631578947367
Computed values of s, i, r at time 100 :
[[0.63157895]
[0.15789474]
[0.21052632]]
The results agree nicely with our analysis, but there are two remaining questions:
Why is it that \(s_t+i_t+r_t = 1\) for every \(t\) in the calculation, even though this constraint does not appear to be enforced?
What are the results if we change the initial population composition given by \(X_0\)?
The answer to the first question is simple, but not immediately obvious. If we add the three equations that define \(s_t\), \(i_t\), and \(r_t\), we find that \(s_t+i_t+r_t = s_{t-1} + i_{t-1} + r_{t-1}\). This means that as long as the three equations are true, the sum \(s_t+i_t+r_t\) does not change as we carry the iteration forward in time. Since we start with the components of \(X_0\) summing to one, the components of \(X_t\) will sum to one for all other \(t\).
In order to answer the second question, we assume that we start with an arbitrary \(X_0\), although we will require that the components of \(X_0\) to be nonnegative and add to one so that the results are meaningful in the application. As we did when considering the Power Method, we write \(X_0\) as a linear combination of \(\{V_1, V_2, V_3\}\), the eigenvectors of \(A\). This makes it easy to express \(X_t = A^tX_0\) in terms of the eigenvectors as well, just as did before.
Based on the results of our calculations, we expect that \(\lambda_1 = 1\) is the largest eigenvalue.
If \(|\lambda_2|< 1\) and \(|\lambda_3|< 1\) then the second and the third terms will approach zero as \(t\) gets large and \(X_t\) will approach the direction of \(V_1\). We use SciPy to verify that \(\lambda_2\) and \(\lambda_3\) are in fact a pair of complex numbers with magnitude less than one.
import scipy.linalg as sla
EVAL, EVEC = sla.eig(A)
print(EVAL,'\n')
## The absolute function finds the absolute value of real numbers, or magnitude of complex numbers.
print(np.absolute(EVAL[1]))
print(np.absolute(EVAL[2]))
[1. +0.j 0.8+0.08660254j 0.8-0.08660254j]
0.8046738469715539
0.8046738469715539
Therefore unless \(c_1 = 0\), the vectors \(X_t\) approach the direction of \(V_1\). The case where \(c_1 = 0\) corresponds to initial conditions with negative entries, which are not relevant in the application. The conclusion is that \(X_t\) will approach the equilibrium state for all meaningful initial configurations of the population.
Exercises¶
Exercise 1: Experiment with a range of initial conditions in the infectious disease model to provide evidence that an equilibrium state is reached for all meaningful initial states.
## Code solution here.
Exercise 2: Perform an analysis similar to that in the example for the following infectious disease model. In this model the rate at which individuals move from the Recovered category to the Susceptible category is less than that in the example. Make a plot similar to that in the example and also calculate the theoretical equilibrium values for \(s\), \(i\), and \(r\).
## Code solution here.