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
Assignment wrap up

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 th progress of time and $A_t$ are an $n \times n$ matrix 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-invarient” 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
Decreased: 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 immuned, and 5% died.

We choose a simple model here. After each day,

5% of the susceptible individuals will get infected
3% of infected inviduals will die
10% of infected inviduals will recover and immuned to the disease
4% of infected inviduals will recover but not immuned to the disease
83% of the infected inviduals will remain

✅ Do this: Write the dynamics matrix for the above markov linear dynamical system. Come to class ready to discuss the matrix. (hint the columns of the matrix should add to 1).

# Put your matrix here

✅ Do this: Review how we solved for the long term steady state of the markov system. See if you can find these probabilities for your dyamics matrix.

# Put your matrix here

3. Ordinary Differential Equations¶

Ordinary Differential Equations (ODEs) are yet another for 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 is 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)

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 at $\dot{x}$ at time $t$. Once we know the derivative, we can increment the time to 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 all of the $x$ statues:

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()

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
<ipython-input-4-97c04e533e83> in <module>
----> 1 get_ipython().run_line_magic('matplotlib', 'inline')
      2 import matplotlib.pylab as plt
      3 import numpy as np
      4 import sympy as sym
      5 sym.init_printing()

~/REPOS/MTH314_Textbook/MakeTextbook/envs/lib/python3.9/site-packages/IPython/core/interactiveshell.py in run_line_magic(self, magic_name, line, _stack_depth)
   2342                 kwargs['local_ns'] = self.get_local_scope(stack_depth)
   2343             with self.builtin_trap:
-> 2344                 result = fn(*args, **kwargs)
   2345             return result
   2346 

~/REPOS/MTH314_Textbook/MakeTextbook/envs/lib/python3.9/site-packages/decorator.py in fun(*args, **kw)
    230             if not kwsyntax:
    231                 args, kw = fix(args, kw, sig)
--> 232             return caller(func, *(extras + args), **kw)
    233     fun.__name__ = func.__name__
    234     fun.__doc__ = func.__doc__

~/REPOS/MTH314_Textbook/MakeTextbook/envs/lib/python3.9/site-packages/IPython/core/magic.py in <lambda>(f, *a, **k)
    185     # but it's overkill for just that one bit of state.
    186     def magic_deco(arg):
--> 187         call = lambda f, *a, **k: f(*a, **k)
    188 
    189         if callable(arg):

~/REPOS/MTH314_Textbook/MakeTextbook/envs/lib/python3.9/site-packages/IPython/core/magics/pylab.py in matplotlib(self, line)
     97             print("Available matplotlib backends: %s" % backends_list)
     98         else:
---> 99             gui, backend = self.shell.enable_matplotlib(args.gui.lower() if isinstance(args.gui, str) else args.gui)
    100             self._show_matplotlib_backend(args.gui, backend)
    101 

~/REPOS/MTH314_Textbook/MakeTextbook/envs/lib/python3.9/site-packages/IPython/core/interactiveshell.py in enable_matplotlib(self, gui)
   3511         """
   3512         from IPython.core import pylabtools as pt
-> 3513         gui, backend = pt.find_gui_and_backend(gui, self.pylab_gui_select)
   3514 
   3515         if gui != 'inline':

~/REPOS/MTH314_Textbook/MakeTextbook/envs/lib/python3.9/site-packages/IPython/core/pylabtools.py in find_gui_and_backend(gui, gui_select)
    278     """
    279 
--> 280     import matplotlib
    281 
    282     if gui and gui != 'auto':

ModuleNotFoundError: No module named 'matplotlib'

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 (0,0) 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 refernce

✅ Do this: Let $$A= \begin{bmatrix}2 & 3 \\ 4 & -2\end{bmatrix}$$

Write a loop over the points $(1.5,1)$, $(-1.5,-1)$, $(-1,2)$ and plot the results of the traj function:

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: Let $$A= \begin{bmatrix}6 & -1 \\ 1 & 4\end{bmatrix}$$

Write a loop over the points $(1.5,1)$, $(-1.5,-1)$, $(-1,2)$, $(1,-2)$ and $(2,-2)$ and plot the results of the traj function:

# Put your code here

✅ Do this: Let $$A= \begin{bmatrix}5 & 2 \\ -4 & 1\end{bmatrix}$$

Write a loop over the points $(1.5,1)$, $(-1.5,-1)$, $(-1,2)$, $(1,-2)$ and $(2,-2)$ and plot the results of the traj function:

# Put your code here

4. Assignment wrap up¶

✅ Assignment-Specific QUESTION: Where you able to get the ODE code working in the above example. If not, where did you get stuck?