{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**This is a modified version of Chapter 9 notebook**"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"] add https://github.com/VMLS-book/VMLS.jl"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [],
"source": [
"using VMLS\n",
"using LinearAlgebra"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Chapter 9\n",
"# Linear dynamical systems\n",
"### 9.1 Linear dynamical systems\n",
"Let’s simulate a time-invariant linear dynamic system \n",
"$$\n",
"x_{t+1} = Ax_t, t = 1, . . . , T,\n",
"$$\n",
"with dynamics matrix\n",
"$$\n",
"A =\n",
"\\left[\\matrix{\n",
"0.97 & 0.10 & −0.05 \\cr\n",
"−0.3 & 0.99 & 0.05 \\cr\n",
"0.01 & −0.04 & 0.96\n",
"}\\right]\n",
"$$\n",
"\n",
"and initial state $x_1 = (1, 0,−1)$. We store the state trajectory in the $n × T$ matrix\n",
"`state_traj`, with the $i$th column $x_t$. We plot the result in Figure 9.1."
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x_1 = [1,0,-1]; # initial state\n",
"n = length(x_1); T = 50;\n",
"A = [ 0.97 0.10 -0.05 ; -0.3 0.99 0.05 ; 0.01 -0.04 0.96 ]\n",
"state_traj = [x_1 zeros(n,T-1) ];\n",
"for t=1:T-1 # Dynamics recursion\n",
"state_traj[:,t+1] = A*state_traj[:,t];\n",
"end\n",
"using Plots\n",
"plot(1:T, state_traj', xlabel = \"t\", \n",
" label = [\"(x_t)_1\", \"(x_t)_2\", \"(x_t)_3\"],\n",
" size=[300,200],legend=false)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Figure 9.1** Linear dynamical system simulation.\n",
"\n",
"### 9.2 Population dynamics\n",
"We can create a population dynamics matrix with just one simple line of Julia. The following code predicts the 2020 population distribution in the US using the data of Section [9.2](https://web.stanford.edu/~boyd/vmls/vmls.pdf#section.9.2) of VMLS, which are available through the VMLS function `population_data`. The result is shown in Figure 9.2."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"D = population_data();\n",
"b = D[\"birth_rate\"];\n",
"d = D[\"death_rate\"];"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"100"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"length(b)"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(b,size=[300,200])"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"100"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"length(d)"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(d,size=[300,200])"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = [b'; diagonal(1 .- d[1:end-1]) zeros(length(d)-1)];\n",
"x = D[\"population\"];\n",
"plot(x,size=[300,200])"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Import 3 100-vectors: population, birth_rate, death_rate\n",
"\n",
"\n",
"for k = 1:10\n",
"global x\n",
"x = A*x;\n",
"end;\n",
"using Plots\n",
"plot!(x, legend=false, xlabel = \"Age\", \n",
" ylabel = \"Population (millions)\",size=[300,200])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Figure 9.2** Predicted age distribution in the US in 2020.\n",
"\n",
"\n",
"Note the keyword `global` in the for-loop. Without this statement, the scope of the\n",
"variable x created by the assignment `x = A*x` would be local to the for-loop, i.e.,\n",
"this variable does not exist outside the loop and is different from the `x` outside the\n",
"loop.\n",
"\n",
"**9.3 Epidemic dynamics**\n",
"Let’s implement the simulation of the epidemic dynamics from VMLS §[9.3](https://web.stanford.edu/\\\\%7Eboyd/vmls/vmls.pdf#section.9.3). The plot is in figure 9.3."
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"T = 210;\n",
"A = [ 0.95 0.04 0 0 ; 0.05 0.85 0 0 ; 0 0.10 1 0 ; 0 0.01 0 1 ];\n",
"x_1 = [1,0,0,0];\n",
"state_traj = [x_1 zeros(4,T-1) ]; # State trajectory\n",
"for t=1:T-1 # Dynamics recursion\n",
"state_traj[:,t+1] = A*state_traj[:,t];\n",
"end\n",
"using Plots\n",
"plot(1:T, state_traj', xlabel = \"Time t\", \n",
" label = [\"Susceptible\", \"Infected\", \"Recovered\", \"Deceased\"],\n",
" size=[400,200],legend=false)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Figure 9.3** Simulation of epidemic dynamics."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**9.4 Motion of a mass**\n",
"Let’s simulate the discretized model of the motion of a mass in §[9.4](https://web.stanford.edu/~boyd/vmls/vmls.pdf#section.9.4) of VMLS. See figure 9.4."
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"h = 0.01; m = 1; eta = 1;\n",
"A = [ 1 h ; 0 1-h*eta/m ];\n",
"B = [ 0 ; h/m ];\n",
"x1 = [0,0];\n",
"K = 600; # simulate for K*h = 6 seconds\n",
"f = zeros(K); f[50:99] .= 1.0; f[100:139] .= -1.3;\n",
"X = [x1 zeros(2,K-1)];\n",
"for k=1:K-1\n",
"X[:,k+1] = A* X[:,k] + B*f[k]\n",
"end\n",
"using Plots\n",
"plot(X[1,:], xlabel=\"k\", ylabel=\"Position\",\n",
" size=[400,200],legend=false )"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(X[2,:], xlabel=\"k\", ylabel=\"Velocity\", size=[400,200],legend=false )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Figure 9.4** Simulation of a mass moving along a line: position (top) and\n",
"velocity (bottom)."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.5.2",
"language": "julia",
"name": "julia-1.5"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.5.2"
}
},
"nbformat": 4,
"nbformat_minor": 4
}