II_11_Markov_matrices_Projections_and_Fourier_series.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ This notebook is part of lecture 24 *Markov matrices and Fourier series* in the OCW MIT course 18.06 by Prof Gilbert Strang [1]\n",
    "+ Created by me, Dr Juan H Klopper\n",
    "    + Head of Acute Care Surgery\n",
    "    + Groote Schuur Hospital\n",
    "    + University Cape Town\n",
    "    + <a href=\"mailto:juan.klopper@uct.ac.za\">Email me with your thoughts, comments, suggestions and corrections</a> \n",
    "<a rel=\"license\" href=\"http://creativecommons.org/licenses/by-nc/4.0/\"><img alt=\"Creative Commons Licence\" style=\"border-width:0\" src=\"https://i.creativecommons.org/l/by-nc/4.0/88x31.png\" /></a><br /><span xmlns:dct=\"http://purl.org/dc/terms/\" href=\"http://purl.org/dc/dcmitype/InteractiveResource\" property=\"dct:title\" rel=\"dct:type\">Linear Algebra OCW MIT18.06</span> <span xmlns:cc=\"http://creativecommons.org/ns#\" property=\"cc:attributionName\">IPython notebook [2] study notes by Dr Juan H Klopper</span> is licensed under a <a rel=\"license\" href=\"http://creativecommons.org/licenses/by-nc/4.0/\">Creative Commons Attribution-NonCommercial 4.0 International License</a>.\n",
    "\n",
    "+ [1] <a href=\"http://ocw.mit.edu/courses/mathematics/18-06sc-linear-algebra-fall-2011/index.htm\">OCW MIT 18.06</a>\n",
    "+ [2] Fernando Pérez, Brian E. Granger, IPython: A System for Interactive Scientific Computing, Computing in Science and Engineering, vol. 9, no. 3, pp. 21-29, May/June 2007, doi:10.1109/MCSE.2007.53. URL: http://ipython.org"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Arvo:400,700,400italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=PT+Mono' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Shadows+Into+Light' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Philosopher:400,700,400italic,700italic' rel='stylesheet' type='text/css'>\n",
       "\n",
       "<style>\n",
       "\n",
       "@font-face {\n",
       "    font-family: \"Computer Modern\";\n",
       "    src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "}\n",
       "\n",
       "#notebook_panel { /* main background */\n",
       "    background: #ddd;\n",
       "    color: #000000;\n",
       "}\n",
       "\n",
       "\n",
       "\n",
       "/* Formatting for header cells */\n",
       ".text_cell_render h1 {\n",
       "    font-family: 'Philosopher', sans-serif;\n",
       "    font-weight: 400;\n",
       "    font-size: 2.2em;\n",
       "    line-height: 100%;\n",
       "    color: rgb(0, 80, 120);\n",
       "    margin-bottom: 0.1em;\n",
       "    margin-top: 0.1em;\n",
       "    display: block;\n",
       "}\t\n",
       ".text_cell_render h2 {\n",
       "    font-family: 'Philosopher', serif;\n",
       "    font-weight: 400;\n",
       "    font-size: 1.9em;\n",
       "    line-height: 100%;\n",
       "    color: rgb(200,100,0);\n",
       "    margin-bottom: 0.1em;\n",
       "    margin-top: 0.1em;\n",
       "    display: block;\n",
       "}\t\n",
       "\n",
       ".text_cell_render h3 {\n",
       "    font-family: 'Philosopher', serif;\n",
       "    margin-top:12px;\n",
       "    margin-bottom: 3px;\n",
       "    font-style: italic;\n",
       "    color: rgb(94,127,192);\n",
       "}\n",
       "\n",
       ".text_cell_render h4 {\n",
       "    font-family: 'Philosopher', serif;\n",
       "}\n",
       "\n",
       ".text_cell_render h5 {\n",
       "    font-family: 'Alegreya Sans', sans-serif;\n",
       "    font-weight: 300;\n",
       "    font-size: 16pt;\n",
       "    color: grey;\n",
       "    font-style: italic;\n",
       "    margin-bottom: .1em;\n",
       "    margin-top: 0.1em;\n",
       "    display: block;\n",
       "}\n",
       "\n",
       ".text_cell_render h6 {\n",
       "    font-family: 'PT Mono', sans-serif;\n",
       "    font-weight: 300;\n",
       "    font-size: 10pt;\n",
       "    color: grey;\n",
       "    margin-bottom: 1px;\n",
       "    margin-top: 1px;\n",
       "}\n",
       "\n",
       ".CodeMirror{\n",
       "        font-family: \"PT Mono\";\n",
       "        font-size: 100%;\n",
       "}\n",
       "\n",
       "</style>\n",
       "\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.core.display import HTML, Image\n",
    "css_file = 'style.css'\n",
    "HTML(open(css_file, 'r').read())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Arvo:400,700,400italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=PT+Mono' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Shadows+Into+Light' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Philosopher:400,700,400italic,700italic' rel='stylesheet' type='text/css'>\n",
       "\n",
       "<style>\n",
       "\n",
       "@font-face {\n",
       "    font-family: \"Computer Modern\";\n",
       "    src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "}\n",
       "\n",
       "#notebook_panel { /* main background */\n",
       "    background: #ddd;\n",
       "    color: #000000;\n",
       "}\n",
       "\n",
       "\n",
       "\n",
       "/* Formatting for header cells */\n",
       ".text_cell_render h1 {\n",
       "    font-family: 'Philosopher', sans-serif;\n",
       "    font-weight: 400;\n",
       "    font-size: 2.2em;\n",
       "    line-height: 100%;\n",
       "    color: rgb(0, 80, 120);\n",
       "    margin-bottom: 0.1em;\n",
       "    margin-top: 0.1em;\n",
       "    display: block;\n",
       "}\t\n",
       ".text_cell_render h2 {\n",
       "    font-family: 'Philosopher', serif;\n",
       "    font-weight: 400;\n",
       "    font-size: 1.9em;\n",
       "    line-height: 100%;\n",
       "    color: rgb(200,100,0);\n",
       "    margin-bottom: 0.1em;\n",
       "    margin-top: 0.1em;\n",
       "    display: block;\n",
       "}\t\n",
       "\n",
       ".text_cell_render h3 {\n",
       "    font-family: 'Philosopher', serif;\n",
       "    margin-top:12px;\n",
       "    margin-bottom: 3px;\n",
       "    font-style: italic;\n",
       "    color: rgb(94,127,192);\n",
       "}\n",
       "\n",
       ".text_cell_render h4 {\n",
       "    font-family: 'Philosopher', serif;\n",
       "}\n",
       "\n",
       ".text_cell_render h5 {\n",
       "    font-family: 'Alegreya Sans', sans-serif;\n",
       "    font-weight: 300;\n",
       "    font-size: 16pt;\n",
       "    color: grey;\n",
       "    font-style: italic;\n",
       "    margin-bottom: .1em;\n",
       "    margin-top: 0.1em;\n",
       "    display: block;\n",
       "}\n",
       "\n",
       ".text_cell_render h6 {\n",
       "    font-family: 'PT Mono', sans-serif;\n",
       "    font-weight: 300;\n",
       "    font-size: 10pt;\n",
       "    color: grey;\n",
       "    margin-bottom: 1px;\n",
       "    margin-top: 1px;\n",
       "}\n",
       "\n",
       ".CodeMirror{\n",
       "        font-family: \"PT Mono\";\n",
       "        font-size: 100%;\n",
       "}\n",
       "\n",
       "</style>\n",
       "\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.core.display import HTML, Image\n",
    "css_file = 'style.css'\n",
    "HTML(open(css_file, 'r').read())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sympy import init_printing, Matrix, symbols, eye\n",
    "from warnings import filterwarnings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "init_printing(use_latex = 'mathjax')\n",
    "filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Markov matrices and steady state\n",
    "# Projections and Fourier series"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Markov matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Consider the follow Markov matrix (transition matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left[\\begin{matrix}0.1 & 0.01 & 0.3\\\\0.2 & 0.99 & 0.3\\\\0.7 & 0 & 0.4\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡0.1  0.01  0.3⎤\n",
       "⎢              ⎥\n",
       "⎢0.2  0.99  0.3⎥\n",
       "⎢              ⎥\n",
       "⎣0.7   0    0.4⎦"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = Matrix([[0.1, 0.01, 0.3], [0.2, 0.99, 0.3], [0.7, 0, 0.4]])\n",
    "A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* All the column entries add to 1 (also true for powers of the matrix)\n",
    "* All entries &#8805; 0 (also true for powers of the matrix)\n",
    "* There will be an eigenvalue of 1\n",
    "* All other eigenvalues will have an absolute value of at most 1\n",
    "* For difference equations we will have the following\n",
    "$$ {u}_{k}={A}^{k}{u}_{0}={c}_{1}{\\lambda}_{1}^{k}{x}_{1}+{c}_{1}{\\lambda}_{2}^{k}{x}_{2}+{\\dots} $$\n",
    "* Thus for the eigenvalues other than 1, successive powers will lead to those terms approaching zero and a steady state being reached by the term with the eigenvalue of 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\begin{Bmatrix}1 : 1, & \\frac{49}{200} + \\frac{3 \\sqrt{1009}}{200} : 1, & - \\frac{3 \\sqrt{1009}}{200} + \\frac{49}{200} : 1\\end{Bmatrix}$$"
      ],
      "text/plain": [
       "⎧                ______           ______         ⎫\n",
       "⎪       49   3⋅╲╱ 1009        3⋅╲╱ 1009     49   ⎪\n",
       "⎨1: 1, ─── + ──────────: 1, - ────────── + ───: 1⎬\n",
       "⎪      200      200              200       200   ⎪\n",
       "⎩                                                ⎭"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.eigenvals()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* All the components of the eigenvector of the eigenvalue 1 is positive"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\begin{bmatrix}\\begin{pmatrix}1.0, & 1, & \\begin{bmatrix}\\left[\\begin{matrix}0.857142857142857\\\\47.1428571428571\\\\1.0\\end{matrix}\\right]\\end{bmatrix}\\end{pmatrix}, & \\begin{pmatrix}0.721471405228058, & 1, & \\begin{bmatrix}\\left[\\begin{matrix}0.459244864611511\\\\-1.45924486461151\\\\1.0\\end{matrix}\\right]\\end{bmatrix}\\end{pmatrix}, & \\begin{pmatrix}-0.231471405228058, & 1, & \\begin{bmatrix}\\left[\\begin{matrix}-0.902102007468654\\\\-0.0978979925313461\\\\1.0\\end{matrix}\\right]\\end{bmatrix}\\end{pmatrix}\\end{bmatrix}$$"
      ],
      "text/plain": [
       "⎡⎛1.0, 1, ⎡⎡0.857142857142857⎤⎤⎞, ⎛0.721471405228058, 1, ⎡⎡0.459244864611511⎤⎤\n",
       "⎢⎜        ⎢⎢                 ⎥⎥⎟  ⎜                      ⎢⎢                 ⎥⎥\n",
       "⎢⎜        ⎢⎢47.1428571428571 ⎥⎥⎟  ⎜                      ⎢⎢-1.45924486461151⎥⎥\n",
       "⎢⎜        ⎢⎢                 ⎥⎥⎟  ⎜                      ⎢⎢                 ⎥⎥\n",
       "⎣⎝        ⎣⎣       1.0       ⎦⎦⎠  ⎝                      ⎣⎣       1.0       ⎦⎦\n",
       "\n",
       "⎞, ⎛-0.231471405228058, 1, ⎡⎡-0.902102007468654 ⎤⎤⎞⎤\n",
       "⎟  ⎜                       ⎢⎢                   ⎥⎥⎟⎥\n",
       "⎟  ⎜                       ⎢⎢-0.0978979925313461⎥⎥⎟⎥\n",
       "⎟  ⎜                       ⎢⎢                   ⎥⎥⎟⎥\n",
       "⎠  ⎝                       ⎣⎣        1.0        ⎦⎦⎠⎦"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.eigenvects()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* If the matrix A - 1&#955; is singular, then the eigenvalue is 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left[\\begin{matrix}-0.9 & 0.01 & 0.3\\\\0.2 & -0.01 & 0.3\\\\0.7 & 0 & -0.6\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡-0.9  0.01   0.3 ⎤\n",
       "⎢                 ⎥\n",
       "⎢0.2   -0.01  0.3 ⎥\n",
       "⎢                 ⎥\n",
       "⎣0.7     0    -0.6⎦"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A - eye(3) # A minus the identity matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$-3.03576608295941 \\cdot 10^{-18}$$"
      ],
      "text/plain": [
       "-3.03576608295941e-18"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(A - eye(3)).det() # A computer peculiarity to inidcate 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* The sum of the entries in each column is now zero\n",
    "* We would like a proof involving this (sum of column entries equal zero) as an assumption to give a singular matrix, without having to calculate the determinant\n",
    "* It is easy to see that the rows are linearly dependent and that would give proof of a singular matrix\n",
    "* Look also at the nullspace of A<sup>T</sup>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\begin{bmatrix}\\end{bmatrix}$$"
      ],
      "text/plain": [
       "[]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A_1t = A - eye(3)\n",
    "(A_1t.transpose()).nullspace()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\begin{bmatrix}\\end{bmatrix}$$"
      ],
      "text/plain": [
       "[]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(A.transpose()).nullspace()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\begin{bmatrix}\\end{bmatrix}$$"
      ],
      "text/plain": [
       "[]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.nullspace()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\begin{pmatrix}\\left[\\begin{matrix}1.0 & 0 & 0\\\\0 & 1.0 & 0\\\\0 & 0 & 1.0\\end{matrix}\\right], & \\begin{bmatrix}0, & 1, & 2\\end{bmatrix}\\end{pmatrix}$$"
      ],
      "text/plain": [
       "⎛⎡1.0   0    0 ⎤, [0, 1, 2]⎞\n",
       "⎜⎢             ⎥           ⎟\n",
       "⎜⎢ 0   1.0   0 ⎥           ⎟\n",
       "⎜⎢             ⎥           ⎟\n",
       "⎝⎣ 0    0   1.0⎦           ⎠"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(A - eye(3)).rref()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Note how the eigenvalues of A and A<sup>T</sup> are the same"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.eigenvals() == A.transpose().eigenvals()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* The proof of this lies in the fact that we calculate the eigenvalue(s) by the following equation\n",
    "$$ \\left|{A}-{\\lambda}{I}\\right|={0} $$\n",
    "* The eigenvalue(s) must also solve this equation\n",
    "$$ \\left|{A}^{T}-{\\lambda}{I}\\right|={0} $$\n",
    "* Since &#955;I=&#955;I<sup>T</sup> and both of these equations eqal the same value (0), the eigenvalue(s) must be equal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example of a Markov matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Consider the population of two (isolated) states (*c* and *m*)\n",
    "* Over time there is movement between the two (no loss or entry from outside the system)\n",
    "* We consider the following system fo movement (difference equation)\n",
    "$$ {u}_{k+1}={A}{u}_{k} $$\n",
    "* We can create the following matrix equation\n",
    "$$ \\begin{bmatrix} { u }_{ c } \\\\ { u }_{ m } \\end{bmatrix}_{k+1}=\\begin{bmatrix} 0.9 & 0.2 \\\\ 0.1 & 0.8 \\end{bmatrix}\\begin{bmatrix} { u }_{ c } \\\\ { u }_{ m } \\end{bmatrix}_{k} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left[\\begin{matrix}0.9 & 0.2\\\\0.1 & 0.8\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡0.9  0.2⎤\n",
       "⎢        ⎥\n",
       "⎣0.1  0.8⎦"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = Matrix([[0.9, 0.2], [0.1, 0.8]])\n",
    "A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Let's start at time, *t*=0 (or in our case *k*=0), with the whole population in state *m*\n",
    "* Because the system is closed, the total population stays at 1000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "u0 = Matrix([0, 1000])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* After 1 time period *k*+1 will be given by the following"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left[\\begin{matrix}200.0\\\\800.0\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡200.0⎤\n",
       "⎢     ⎥\n",
       "⎣800.0⎦"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u1 = A * u0\n",
    "u1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left[\\begin{matrix}340.0\\\\660.0\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡340.0⎤\n",
       "⎢     ⎥\n",
       "⎣660.0⎦"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u2 = A * u1\n",
    "u2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Where will we end up?\n",
    "* Let's look at the eigenvalues\n",
    "* One will be one and the other must be the trace minus 1 = 0.7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\begin{Bmatrix}\\frac{7}{10} : 1, & 1 : 1\\end{Bmatrix}$$"
      ],
      "text/plain": [
       "{7/10: 1, 1: 1}"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.eigenvals()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\begin{bmatrix}\\begin{pmatrix}0.7, & 1, & \\begin{bmatrix}\\left[\\begin{matrix}-1.0\\\\1.0\\end{matrix}\\right]\\end{bmatrix}\\end{pmatrix}, & \\begin{pmatrix}1.0, & 1, & \\begin{bmatrix}\\left[\\begin{matrix}2.0\\\\1.0\\end{matrix}\\right]\\end{bmatrix}\\end{pmatrix}\\end{bmatrix}$$"
      ],
      "text/plain": [
       "⎡⎛0.7, 1, ⎡⎡-1.0⎤⎤⎞, ⎛1.0, 1, ⎡⎡2.0⎤⎤⎞⎤\n",
       "⎢⎜        ⎢⎢    ⎥⎥⎟  ⎜        ⎢⎢   ⎥⎥⎟⎥\n",
       "⎣⎝        ⎣⎣1.0 ⎦⎦⎠  ⎝        ⎣⎣1.0⎦⎦⎠⎦"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.eigenvects()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* When can look at what happens after a hundred steps using a small *while* loop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left[\\begin{matrix}666.666666666668\\\\333.333333333334\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡666.666666666668⎤\n",
       "⎢                ⎥\n",
       "⎣333.333333333334⎦"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u = Matrix([0, 1000])\n",
    "i = 0\n",
    "\n",
    "while i < 101:\n",
    "    u = A * u\n",
    "    i += 1\n",
    "u"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* It seems that we are at about &#8532; vs &#8531;\n",
    "* It looks suspiciously like the eigenvector of the eigenvalue 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Remember the equation for the difference equation we used above?\n",
    "$$ {u}_{k}={A}^{k}{u}_{0}={c}_{1}{\\lambda}_{1}^{k}{x}_{1}+{c}_{1}{\\lambda}_{2}^{k}{x}_{2}+{\\dots} $$\n",
    "* We only have two terms, so this will become somewhat shortened"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ {u}_{k}={A}^{k}{u}_{0}={c}_{1}{\\lambda}_{1}^{k}{x}_{1}+{c}_{1}{\\lambda}_{2}^{k}{x}_{2} \\\\ {u}_{k}={c}_{1}{1}^{k}\\begin{bmatrix}2\\\\1\\end{bmatrix}+{c}_{1}{\\frac{7}{10}}^{k}{\\begin{bmatrix}-1\\\\1\\end{bmatrix}} \\\\ {u}_{0}={c}_{1}{1}^{0}\\begin{bmatrix}2\\\\1\\end{bmatrix}+{c}_{1}{\\left({\\frac{7}{10}}\\right)}^{0}{\\begin{bmatrix}-1\\\\1\\end{bmatrix}} \\\\ {u}_{0}={c}_{1}\\begin{bmatrix}2\\\\1\\end{bmatrix}+{c}_{1}{\\begin{bmatrix}-1\\\\1\\end{bmatrix}}=\\begin{bmatrix}0\\\\1000\\end{bmatrix} \\\\  $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\begin{pmatrix}\\left[\\begin{matrix}1 & 0 & \\frac{1000}{3}\\\\0 & 1 & \\frac{2000}{3}\\end{matrix}\\right], & \\begin{bmatrix}0, & 1\\end{bmatrix}\\end{pmatrix}$$"
      ],
      "text/plain": [
       "⎛⎡1  0  1000/3⎤, [0, 1]⎞\n",
       "⎜⎢            ⎥        ⎟\n",
       "⎝⎣0  1  2000/3⎦        ⎠"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result = Matrix([[2, -1, 0], [1, 1, 1000]])\n",
    "result.rref()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* So we have *c*<sub>1</sub>=1000&#247;3 and *c*<sub>2</sub>=2000&#247;3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* This results in the final solution\n",
    "$$ {u}_{k}=\\frac{1000}{3}{1}^{k}\\begin{bmatrix}2\\\\1\\end{bmatrix}+\\frac{2000}{3}{\\left({\\frac{7}{10}}\\right)}^{k}{\\begin{bmatrix}-1\\\\1\\end{bmatrix}} $$\n",
    "* Now we can work out the solution for any time step *k* and even see what happens at time approaches infinity\n",
    "* We note that the second expression disappears as *k* approaches infinity, which represents the steady state"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Projections and Fourier series"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Consider projections with **orthonormal** basis *q*<sub>1</sub>, *q*<sub>2</sub>, ..., *q*<sub>n</sub>\n",
    "* Any vector **v** can be expressed (expanded) as an combination of this basis\n",
    "$$ {v}={x}_{1}{q}_{1}+{x}_{2}{q}_{2}+\\dots+{x}_{n}{q}_{n} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Because the basis *q*<sub>i</sub> were chosen to be orthogonal taking the dot product of *q*<sub>1</sub> on both sides will cancel out all *q*<sub>&#8800;1</sub> factors\n",
    "$$ {q}_{1}^{T}{v}={x}_{1}{q}_{1}^{T}{q}_{1}+{0}+{0}+\\dots \\\\ {q}_{1}^{T}{v}={x}_{1} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* We can see this clearly in matrix form below\n",
    "$$ \\begin{bmatrix} \\vdots  & \\vdots  & \\vdots  \\\\ { q }_{ 1 } & \\cdots  & { q }_{ n } \\\\ \\vdots  & \\vdots  & \\vdots  \\end{bmatrix}\\begin{bmatrix} { x }_{ 1 } \\\\ \\vdots  \\\\ { x }_{ n } \\end{bmatrix}=v\\\\ Qx=v\\\\ x={ Q }^{ -1 }v\\\\ \\because \\quad { Q }^{ -1 }={ Q }^{ T }\\\\ x={ Q }^{ T }v $$\n",
    "* The first component of **x**, *x*<sub>1</sub> is the the first row of Q<sup>T</sup> times **v**\n",
    "* The first row of Q<sup>T</sup> is just q<sub>1</sub>, with a bunch of zeros and just ends as we had above\n",
    "$$ {q}_{1}^{T}{v}={x}_{1} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* The key here was to choose orthonormal basis vectors\n",
    "* Now for how this relates to Fourier series\n",
    "* We might want to write a function as an **infinite** series of expressions\n",
    "$$ f\\left( x \\right)={ a }_{ 0 }+{ a }_{ 1 }\\cos { x } +{ b }_{ 1 }\\sin { x } +{ a }_{ 2 }\\cos { 2x } +{ b }_{ 2 }\\sin { 2x } +\\dots  $$\n",
    "* The idea here is that there is still something orthogonal in each of these expressions (sine and cosine)\n",
    "* Joseph Fourier realized that he could work in function space\n",
    "* The vectors are now functions in an infinitely large space\n",
    "* The basis vectors are *cos*(*x*), *sin*(*x*), *cos*(2*x*)..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* What is the inner (dot) product of functions that make them orthogonal?\n",
    "* For vectors is was the following\n",
    "$$ {v}^{T}{w}={v}_{1}{w}_{1}+{v}_{2}{w}_{2}+\\dots+{v}_{n}{w}_{n} $$\n",
    "* For functions *f* and *g* the analogue is to multiply them and the analogue of all the addition is integration\n",
    "$$ { f }^{ T }g=\\int { f\\left( x \\right)g\\left( x \\right) } dx $$\n",
    "* Now we need to know the lower and upper limit of integration\n",
    "* We note that the sine and cosine functions are periodic and repeat every 2&#960; (these are periodic functions)\n",
    "$$ \\because \\quad f\\left( x \\right)=f\\left( x+2\\pi  \\right)\\\\ { f }^{ T }g=\\int _{ 0 }^{ 2\\pi  }{ f\\left( x \\right) g\\left( x \\right)  } dx $$\n",
    "* Just like the inner product of *pairs* gave zero because of orthogonality we have the same here\n",
    "$$ \\int _{ 0 }^{ 2\\pi  }{ \\sin { x } \\cos { x }  } dx\\\\ u\\left( x \\right) =\\sin { x } \\\\ \\frac { du }{ dx } =\\cos { x } \\\\ \\cos { x } dx=du\\\\ { u }_{ 2 }\\left( 2\\pi  \\right) =\\sin { 2\\pi  } =0\\\\ { u }_{ 1 }\\left( 0 \\right) =\\sin { 0 } =0\\\\ \\int _{ 0 }^{ 0 }{ u } du=0 $$\n",
    "* With some trigonometric identities we can show the same for all the other *pairs*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* A Fourier series is then an expression of a function expanded on this orthonormal basis pairs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* How do we get the coefficients then?\n",
    "* Same as above, i.e. taking the inner product of both sides with *cos*(*x*)\n",
    "$$ \\int _{ 0 }^{ 2\\pi  }{ f\\left( x \\right) \\cos { x }  } dx={ a }_{ 1 }\\int _{ 0 }^{ 2\\pi  }{ { \\left( \\cos { x }  \\right)  }^{ 2 } } dx+0+0+\\dots =\\pi \\\\ { a }_{ 1 }=\\frac { 1 }{ \\pi  } \\int _{ 0 }^{ 2\\pi  }{ f\\left( x \\right) \\cos { x }  } dx $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example problem (Markov matrixes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example problem 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* A particle jumps between positions A and B with the following probabilities\n",
    "    * A to A (stays in A) probability is 0.6\n",
    "    * A to B probability is 0.4 (so all states from A add to 1.0)\n",
    "    * B to B probability is 0.8\n",
    "    * B to A probability is 0.2\n",
    "* If the particle starts in position A, what is the probability that it will be at A and B after the first step, *n*-steps, and &#8734;-steps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\begin{pmatrix}\\left[\\begin{matrix}0.6 & 0.2\\\\0.4 & 0.8\\end{matrix}\\right], & \\left[\\begin{matrix}1\\\\0\\end{matrix}\\right]\\end{pmatrix}$$"
      ],
      "text/plain": [
       "⎛⎡0.6  0.2⎤, ⎡1⎤⎞\n",
       "⎜⎢        ⎥  ⎢ ⎥⎟\n",
       "⎝⎣0.4  0.8⎦  ⎣0⎦⎠"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_vect = Matrix([1, 0])\n",
    "A = Matrix([[0.6, 0.2], [0.4, 0.8]]) # Look carefully at what goes where\n",
    "A, x_vect # Displaying the two matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* After a single step we will have the following"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left[\\begin{matrix}0.6\\\\0.4\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡0.6⎤\n",
       "⎢   ⎥\n",
       "⎣0.4⎦"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A * x_vect"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* A probability of 0.6 of being in position A and a probability of 0.4 of being in position B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Now look at the following trend\n",
    "$$ { p }_{ 1 }=A{ p }_{ 0 }\\\\ { p }_{ 2 }=A{ p }_{ 1 }=A\\left( { Ap }_{ 0 } \\right) ={ A }^{ 2 }{ p }_{ 0 }\\\\ { p }_{ 3 }={ A }^{ 3 }{ p }_{ 0 }\\\\ \\vdots \\\\ { p }_{ n }={ A }^{ n }{ p }_{ 0 } $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* To take the power of a matrix in python is easy, but let's remind ourselves that we need to wotk with eignevalues and eigenvectors\n",
    "* Because all the column entries add to 1, we are dealing with a Markov matrix\n",
    "* One of the eigenvalues will be 1, the trace is 1.4, therefor the other eigenvalue is 0.4\n",
    "* Remember diagonalization?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\begin{pmatrix}\\left[\\begin{matrix}-1.0 & 1.0\\\\1.0 & 2.0\\end{matrix}\\right], & \\left[\\begin{matrix}0.4 & 0\\\\0 & 1.0\\end{matrix}\\right]\\end{pmatrix}$$"
      ],
      "text/plain": [
       "⎛⎡-1.0  1.0⎤, ⎡0.4   0 ⎤⎞\n",
       "⎜⎢         ⎥  ⎢        ⎥⎟\n",
       "⎝⎣1.0   2.0⎦  ⎣ 0   1.0⎦⎠"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S, D = A.diagonalize()\n",
    "S, D"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Note the matrix of eigenvalues (diagonal) *D*\n",
    "* Note the matrix of eigenvectors\n",
    "* So for the eigenvalue of 1 we have the following eigenvector (steady state)\n",
    "$$ \\begin{bmatrix}1\\\\2\\end{bmatrix} $$\n",
    "* For the other eigenvalue we have the following eigenvector (decay)\n",
    "$$ \\begin{bmatrix}-1\\\\0\\end{bmatrix} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ {u}_{k}={A}^{k}{u}_{0}={c}_{1}{\\lambda}_{1}^{k}{x}_{1}+{c}_{1}{\\lambda}_{2}^{k}{x}_{2} \\\\ {u}_{k}={c}_{1}{1}^{k}\\begin{bmatrix}1\\\\2\\end{bmatrix}+{c}_{1}{\\frac{4}{10}}^{k}{\\begin{bmatrix}-1\\\\1\\end{bmatrix}} \\\\ {u}_{0}={c}_{1}{1}^{0}\\begin{bmatrix}1\\\\2\\end{bmatrix}+{c}_{1}{\\left({\\frac{4}{10}}\\right)}^{0}{\\begin{bmatrix}-1\\\\1\\end{bmatrix}} \\\\ {u}_{0}={c}_{1}\\begin{bmatrix}1\\\\2\\end{bmatrix}+{c}_{1}{\\begin{bmatrix}-1\\\\1\\end{bmatrix}}=\\begin{bmatrix}1\\\\0\\end{bmatrix} \\\\  $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Solving for the constants we can create an augmented matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\begin{pmatrix}\\left[\\begin{matrix}1 & 0 & \\frac{1}{3}\\\\0 & 1 & - \\frac{2}{3}\\end{matrix}\\right], & \\begin{bmatrix}0, & 1\\end{bmatrix}\\end{pmatrix}$$"
      ],
      "text/plain": [
       "⎛⎡1  0  1/3 ⎤, [0, 1]⎞\n",
       "⎜⎢          ⎥        ⎟\n",
       "⎝⎣0  1  -2/3⎦        ⎠"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u = Matrix([[1, -1, 1], [2, 1, 0]])\n",
    "u.rref()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* This gives us the following solution\n",
    "$$ {u}_{n}=\\frac{1}{3}\\left({1}^{n}\\right)\\begin{bmatrix}1\\\\2\\end{bmatrix}+\\frac{-2}{3}{\\left({\\frac{4}{10}}\\right)}^{n}\\begin{bmatrix}-1\\\\1\\end{bmatrix} \\\\  { p\\left( A \\right)  }_{ n }=\\frac { 1 }{ 3 } +\\frac { -2 }{ 3 } { \\left( \\frac { 4 }{ 10 }  \\right)  }^{ n }\\left( -1 \\right) \\\\ { p\\left( B \\right)  }_{ n }=\\frac { 1 }{ 3 } \\left( 2 \\right) +\\frac { -2 }{ 3 } { \\left( \\frac { 4 }{ 10 }  \\right)  }^{ n }\\left( 1 \\right) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* So at infinity it will be at position A with a probability of &#8531; and at position B with a probability of &#8532;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left[\\begin{matrix}0.333333333362792\\\\0.666666666725585\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡0.333333333362792⎤\n",
       "⎢                 ⎥\n",
       "⎣0.666666666725585⎦"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(A ** 1000000) * x_vect # Just to show how easy it is just to take the power of A in python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}