# 16 In-Class Assignment: Linear Dynamical Systems¶

## Agenda for today’s class (80 minutes)¶

%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-1-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'


## 1. Epidemic Dynamics - Discrete Case¶

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

A = np.matrix([[0.95, 0.04, 0, 0],[0.05, 0.83, 0, 0],[0, 0.1, 1, 0],[0,0.03,0,1]])
sym.Matrix(A)


Do this: If we start with $$x(0) = (1, 0, 0, 0)$$ for day 0. Use the for loop to find the distribution of the four groups after 50 days.

x0 = np.matrix([[1],[0],[0],[0]])
x  = x0
for i in range(50):
x = A*x
print(x)


Do this: Write a program to apply the above transformation matrix for 200 iterations and plot the results.

#Put your answer to the above question here


## 1. Epidemic Dynamics - Continuous Model¶

Instead of using the discrete markov model, we can also use a continuous model with ordinary differential equations.

For example, we have that

$\dot{x}_1 = {dx_1(t)\over dt} = -0.05x_1(t)+ 0.04 x_2(t)$

It means that the changes in the susceptible group depends on susceptible and infected individuals. It increase because of the recovered people from infected ones and it decreases because of the infection.

Similarly, we have the equations for all three groups. $$$\dot{x}_2 = {dx_2(t)\over dt} = 0.05x_1(t)-0.17 x_2(t) \\ \dot{x}_3 = {dx_3(t)\over dt}= 0.1 x_2(t) \\ \dot{x}_4 = {dx_4(t)\over dt} = 0.03 x_2(t)$$$

Do this: We can write it as system of ODEs as $$$\dot{x}(t) = Bx(t)$$$$Write down the matrix$$B$ in numpy.matrix

# Put your answer to the above question here.


Do this: Plot all the distribution for 200 days. Then compare it with the discrete version.

x0 = np.matrix([[1],[0],[0],[0]])
n = 200
x_all = np.matrix(np.zeros((4,n)))

### Your code starts here ###

### Your code ends here ###
for i in range(4):
plt.plot(x_all[i].T)

D2, C2 = np.linalg.eig(B)
np.allclose(C2*np.diag(D2)*C2**(-1), B)
C = C2


## 3. Population Dynamics¶

In this section, we consider the distribution of a population at different ages. Let $$x(t)$$ be a 100-vector with $$x_i(t)$$ denoting the number of people with agent $$i-1$$.

The birth rate is given in the vector $$b$$, where $$b_i$$ is the average number of births per person with age $$i-1$$. $$b_i=0$$ for $$i<13$$ and $$i>50$$.

The death rate is given by the vector $$d$$, where $$d_i$$ is the portion of those aged $$i-1$$ who dies this year.

Then we can model the population with the dynamic systems. We consider the 0-year old first. It includes all newborn from all ages, so we have $$$x_1(t+1) = b^\top x(t)$$$$Then we consider all other ages for$$i>1$$.$$$$x_i(t+1) = (1-d_i)x_{i-1}(t)$$$

population = np.array([
3.94415,       3.97807,       4.09693,       4.11904,       4.06317,       4.05686,       4.06638,       4.03058,
4.04649,       4.14835,       4.17254,       4.11442,       4.10624,       4.11801,       4.16598,       4.24282,
4.31614,       4.39529,       4.50085,       4.58523,       4.51913,       4.35429,       4.26464,       4.19857,
4.24936,       4.26235,       4.15231,       4.24887,       4.21525,       4.22308,       4.28567,       3.97022,
3.98685,       3.88015,       3.83922,       3.95643,       3.80209,       3.93445,       4.12188,       4.36480,
4.38327,       4.11498,       4.07610,       4.10511,       4.21150,       4.50887,       4.51976,       4.53526,
4.53880,       4.60590,       4.66029,       4.46463,       4.50085,       4.38035,       4.29200,       4.25471,
4.03751,       3.93639,       3.79493,       3.64127,       3.62113,       3.49260,       3.56318,       3.48388,
2.65713,       2.68076,       2.63914,       2.64936,       2.32367,       2.14232,       2.04312,       1.94932,
1.86427,       1.73696,       1.68449,       1.62008,       1.47107,       1.45533,       1.40012,       1.37119,
1.30851,       1.21287,       1.16142,       1.07481,       0.98572,       0.91472,       0.81421,       0.71291,
0.64062,       0.53800,       0.43556,       0.34499,       0.28139,       0.21698,       0.16944,       0.12972,
0.09522,       0.06814,       0.04590,       0.03227])
d = np.array([
0.00623,       0.00044,       0.00027,       0.00020,       0.00016,       0.00012,       0.00011,       0.00011,
0.00012,       0.00011,       0.00010,       0.00013,       0.00013,       0.00015,       0.00020,       0.00025,
0.00037,       0.00047,       0.00064,       0.00071,       0.00076,       0.00087,       0.00087,       0.00088,
0.00094,       0.00092,       0.00095,       0.00093,       0.00099,       0.00101,       0.00103,       0.00109,
0.00110,       0.00114,       0.00115,       0.00120,       0.00131,       0.00137,       0.00146,       0.00156,
0.00162,       0.00185,       0.00201,       0.00216,       0.00243,       0.00258,       0.00298,       0.00325,
0.00351,       0.00387,       0.00413,       0.00454,       0.00494,       0.00533,       0.00571,       0.00602,
0.00670,       0.00710,       0.00769,       0.00828,       0.00860,       0.00932,       0.00998,       0.01101,
0.01250,       0.01282,       0.01404,       0.01515,       0.01687,       0.01830,       0.01967,       0.02133,
0.02347,       0.02562,       0.02800,       0.03083,       0.03441,       0.03711,       0.04126,       0.04448,
0.04964,       0.05539,       0.06149,       0.06803,       0.07673,       0.08561,       0.09540,       0.10636,
0.11802,       0.13385,       0.15250,       0.16491,       0.18738,       0.20757,       0.22688,       0.25196,
0.27422,       0.29239,       0.32560,       0.34157])
b = np.array([
0.00000,       0.00000,       0.00000,       0.00000,       0.00000,       0.00000,       0.00000,       0.00000,
0.00000,       0.00000,       0.00020,       0.00020,       0.00020,       0.00020,       0.00020,       0.01710,
0.01710,       0.01710,       0.01710,       0.01710,       0.04500,       0.04500,       0.04500,       0.04500,
0.04500,       0.05415,       0.05415,       0.05415,       0.05415,       0.05415,       0.04825,       0.04825,
0.04825,       0.04825,       0.04825,       0.02250,       0.02250,       0.02250,       0.02250,       0.02250,
0.00510,       0.00510,       0.00510,       0.00510,       0.00510,       0.00035,       0.00035,       0.00035,
0.00035,       0.00035,       0.00000,       0.00000,       0.00000,       0.00000,       0.00000,       0.00000,
0.00000,       0.00000,       0.00000,       0.00000,       0.00000,       0.00000,       0.00000,       0.00000,
0.00000,       0.00000,       0.00000,       0.00000,       0.00000,       0.00000,       0.00000,       0.00000,
0.00000,       0.00000,       0.00000,       0.00000,       0.00000,       0.00000,       0.00000,       0.00000,
0.00000,       0.00000,       0.00000,       0.00000,       0.00000,       0.00000,       0.00000,       0.00000,
0.00000,       0.00000,       0.00000,       0.00000,       0.00000,       0.00000,       0.00000,       0.00000,
0.00000,       0.00000,       0.00000,       0.00000])

plt.subplot(1,3,1)
plt.plot(population)
plt.title('population in 2010')
plt.subplot(1,3,2)
plt.plot(b)
plt.title('birth rate in 2010')
plt.subplot(1,3,3)
plt.plot(d)
plt.title('death rate in 2010')

plt.show()


Do this: Find the $$100\times 100$$ matrix A2 such that $$$x(t+1)=Ax(t)$$$

A2 =
print(A2)


Do this: Plot the population distribution at year 2020.

Written by Dr. Dirk Colbry, Michigan State University