16 Pre-Class Assignment: Linear Dynamical Systems#

Readings for this topic (Recommended in bold)#

Boyd Chapter 9 pg 163-173

Goals for today’s pre-class assignment#

Linear Dynamical Systems
Markov Models
Ordinary Differential Equations

1. Linear Dynamical Systems#

A linear dynamical system is a simple model of how a system changes with time. These systems can be represented by the following “dynamics” or “update equation”:

\[x_{(t+1)} = A_tx_t\]

Where $t$ is an integer representing the progress of time and $A_t$ is an $n \times n$ matrices called the dynamics matrices. Often the above matrix does not change with $t$. In this case the system is called “time-invariant”.

We have seen a few “time-invariant” examples in class.

✅ DO THIS: Review Chapter 9 in the Boyd and Vandenberghe text and become familiar with the contents and the basic terminology.

2. Markov Models#

This is not the first time we have used Dynamical Linear Systems.

✅ DO THIS: Review Markov models in 10–Eigenproblems_pre-class-assignment.ipynb. See how this is a special type of linear dynamical systems that work with state probabilities.

Example#

The dynamics of infection and the spread of an epidemic can be modeled as a linear dynamical system.

We count the fraction of the population in the following four groups:

Susceptible: the individuals can be infected next day
Infected: the infected individuals
Recovered (and immune): recovered individuals from the disease and will not be infected again
Deceased: the individuals died from the disease

We denote the fractions of these four groups in $x(t)$. For example $x(t)=(0.8,0.1,0.05,0.05)$ means that at day $t$, 80% of the population are susceptible, 10% are infected, 5% are recovered and immune, and 5% died.

We choose a simple model here. After each day,

5% of the susceptible individuals will get infected
95% of the susceptible individuals will not get infected, but remain susceptible
3% of infected inviduals will die
10% of infected inviduals will recover and be immune to the disease
4% of infected inviduals will recover but will not be immune to the disease
83% of the infected inviduals will remain infected
All recovered and immune individuals will remain recovered and immune
All deceased individuals will remain deceased

✅ Do this: Write the dynamics matrix $P$ for the above markov model. Hint: The columns of the matrix should add to 1.

# Put your matrix here
import numpy as np

✅ Do this: Recall that a steady state vector of a Markov model is a vector $\vec{x}$ such that:

$P\vec{x} = \vec{x}$, i.e., $\vec{x}$ is an eigenvector corresponding to the eigenvalue $1$
$\sum_{i = 1}^{n}x_i = 1$ and $x_i \ge 0$ for all $i$, i.e., the entries are nonnegative and add up to $1$.

Find all of the steady state vectors for the dynamics matrix you found in the previous question. Interpret the meaning of these steady state vectors. Come to class ready to discuss your findings. Hint: There may be more than two steady state vectors for the dynamics matrix $P$.

# Put any code for the above problem here

Put your answers for the above problem here

3. Ordinary Differential Equations#

Ordinary Differential Equations (ODEs) are yet another form of linear dynamical systems and are a scientific model used in a wide range of problems of the basic form:

$$\dot{x} = A x$$#

These are equations such that the instantaneous rate of change in $x$ (i.e. $\dot{x}$ is the derivative of $x$) is dependent on $x$. Many systems can be modeled with these types of equations.

Here is a quick video that introduces the concepts of Differential Equations. The following is a good review of general ODEs.

from IPython.display import YouTubeVideo
YouTubeVideo("8QeCQn7uxnE",width=640,height=360, cc_load_policy=True)
# https://www.youtube.com/watch?v=8QeCQn7uxnE

Now consider an ODE as a system of linear equations:

\[\dot{x_t} = A x_t\]

Based on the current $x$ vector at time $t$ and the matrix $A$, we can calculate the derivative $\dot{x}$ at time $t$. Once we know the derivative, we can increment the time by some small amount $dt$ and calculate a new value of $x$ as follows:

\[x_{t+1} = x_t + \dot{x_t}dt\]

Then we can do the exact sequence of calculations again for $t+2$. The following function has the transition matrix ($A$), the starting state vector ($x_0$) and a number of time steps ($N$) and uses the above equations to calculate each state and return $x$ at each time step:

The following code generates a trajectory of points starting from x_0, applying the matrix $A$ to get $x_1$ and then applying $A$ again to see how the system progresses from the start state.

%matplotlib inline
import matplotlib.pylab as plt
import numpy as np
import sympy as sym
sym.init_printing()

def traj(A, x, n):
    dt = 0.01
    x_all = np.matrix(np.zeros((len(x),n)))   # Store all points on the trajectory
    for i in range(n):  
        x_dot = A*x         # First we transform x into the derrivative
        x = x + x_dot*dt    # Then we estimate x based on the previous value and a small increment of time.
        x_all[:,i] = x[:,0] 
    return x_all

For example the following code uses the matrix $A= \begin{bmatrix}1 & 1 \\ 1 & -2\end{bmatrix}$ and the starting point $(1,1)$ over 50 timesteps to get a graph:

A = np.matrix([[1,1],[1,-2]])
x0 = np.matrix([[1],[1]])

x_all = traj(A, x0, 50)
plt.scatter(np.asarray(x_all[0,:]),np.asarray(x_all[1,:]))

plt.scatter(list(x0[0,:]),list(x0[1,:])) #Plot the start point as a reference

✅ Do this: Consider the dynamical system $\dot{x_t} = A x_t$ with dynamics matrix

\[\begin{split}A= \begin{bmatrix}2 & 3 \\ 4 & -2\end{bmatrix}.\end{split}\]

Using the traj() function, plot the trajectory of the system over $50$ timesteps from each of the starting locations $(1.5,1)$, $(-1.5,-1)$, $(-1,2)$, $(1,-2)$ and $(2,-2)$.

A = np.matrix([[2,3],[4,-2]])
x0 = np.matrix([[1.5, -1.5, -1, 1, 2],[1, -1, 2, -2, -2]])

# Put your code here

✅ Do this: Now, consider the dynamical system $\dot{x_t} = A x_t$ with dynamics matrix

\[\begin{split}A= \begin{bmatrix}6 & -1 \\ 1 & 4\end{bmatrix}.\end{split}\]

Again, using the traj() function, plot the trajectory of the system over $50$ timesteps from each of the starting locations $(1.5,1)$, $(-1.5,-1)$, $(-1,2)$, $(1,-2)$ and $(2,-2)$.

# Put your code here

✅ Do this: Finally, consider the dynamical system $\dot{x_t} = A x_t$ with dynamics matrix

\[\begin{split}A= \begin{bmatrix}5 & 2 \\ -4 & 1\end{bmatrix}.\end{split}\]

Again, using the traj() function, plot the trajectory of the system over $50$ timesteps from each of the starting locations $(1.5,1)$, $(-1.5,-1)$, $(-1,2)$, $(1,-2)$ and $(2,-2)$.

# Put your code here

Written by Dr. Dirk Colbry, Michigan State University
This work is licensed under a Creative Commons Attribution-NonCommercial 4.0 International License.

###STARTFOOTER###

CMSE/MTH 314 Course Website

16 Pre-Class Assignment: Linear Dynamical Systems

Contents