Chapter_2_Elements_of_matrix_theory.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "toc": true
   },
   "source": [
    "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n",
    "<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#Chapter-2---Elements-of-Matrix-Theory\" data-toc-modified-id=\"Chapter-2---Elements-of-Matrix-Theory-1\">Chapter 2 - Elements of Matrix Theory</a></span><ul class=\"toc-item\"><li><span><a href=\"#2.1.2-The-Jordan-Normal-Form\" data-toc-modified-id=\"2.1.2-The-Jordan-Normal-Form-1.1\">2.1.2 The Jordan Normal Form</a></span><ul class=\"toc-item\"><li><span><a href=\"#Example-2.5-Revisiting-the-wireless-sensor-network-example\" data-toc-modified-id=\"Example-2.5-Revisiting-the-wireless-sensor-network-example-1.1.1\">Example 2.5 Revisiting the wireless sensor network example</a></span></li><li><span><a href=\"#NumPy/-SciPy-approach\" data-toc-modified-id=\"NumPy/-SciPy-approach-1.1.2\">NumPy/ SciPy approach</a></span></li><li><span><a href=\"#SymPy-approach\" data-toc-modified-id=\"SymPy-approach-1.1.3\">SymPy approach</a></span></li></ul></li><li><span><a href=\"#2.1.3-Semi-convergence-and-convergence-for-discrete-time-linear-systems\" data-toc-modified-id=\"2.1.3-Semi-convergence-and-convergence-for-discrete-time-linear-systems-1.2\">2.1.3 Semi-convergence and convergence for discrete-time linear systems</a></span><ul class=\"toc-item\"><li><span><a href=\"#Definition-2.6-(Spectrum-and-spectral-radius-of-a-matrix)\" data-toc-modified-id=\"Definition-2.6-(Spectrum-and-spectral-radius-of-a-matrix)-1.2.1\">Definition 2.6 (Spectrum and spectral radius of a matrix)</a></span></li></ul></li><li><span><a href=\"#2.2.1-The-spectral-radius-for-row-stochastic-matrices\" data-toc-modified-id=\"2.2.1-The-spectral-radius-for-row-stochastic-matrices-1.3\">2.2.1 The spectral radius for row-stochastic matrices</a></span><ul class=\"toc-item\"><li><span><a href=\"#Theorem-2.8-(Geršgorin-Disks-Theorem)\" data-toc-modified-id=\"Theorem-2.8-(Geršgorin-Disks-Theorem)-1.3.1\">Theorem 2.8 (Geršgorin Disks Theorem)</a></span></li></ul></li><li><span><a href=\"#2.3.3-Applications-to-matrix-powers-and-averaging-systems\" data-toc-modified-id=\"2.3.3-Applications-to-matrix-powers-and-averaging-systems-1.4\">2.3.3 Applications to matrix powers and averaging systems</a></span><ul class=\"toc-item\"><li><span><a href=\"#Theorem-2.13-(Powers-of-non-negative-matrices-with-a-simple-and-strictly-dominant-eigenvalue)\" data-toc-modified-id=\"Theorem-2.13-(Powers-of-non-negative-matrices-with-a-simple-and-strictly-dominant-eigenvalue)-1.4.1\">Theorem 2.13 (Powers of non-negative matrices with a simple and strictly dominant eigenvalue)</a></span><ul class=\"toc-item\"><li><span><a href=\"#Example-2.14-Wireless-sensor-network\" data-toc-modified-id=\"Example-2.14-Wireless-sensor-network-1.4.1.1\">Example 2.14 Wireless sensor network</a></span></li></ul></li></ul></li><li><span><a href=\"#Exercises-2.18\" data-toc-modified-id=\"Exercises-2.18-1.5\">Exercises 2.18</a></span></li><li><span><a href=\"#Exercises-2.19\" data-toc-modified-id=\"Exercises-2.19-1.6\">Exercises 2.19</a></span></li></ul></li></ul></div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib widget\n",
    "\n",
    "# Import packages\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "import networkx as nx\n",
    "import scipy.linalg as spla\n",
    "from sympy import Matrix\n",
    "\n",
    "# For interactive graphs\n",
    "import ipywidgets as widgets\n",
    "\n",
    "# Import self defined functions\n",
    "import lib  # General library\n",
    "\n",
    "# Settings\n",
    "custom_figsize= (6, 4) # Might need to change this value to fit the figures to your screen\n",
    "custom_figsize_square = (5, 5)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter 2 - Elements of Matrix Theory\n",
    "These Jupyter Notebook scripts contain some examples, visualization and supplements accompanying the book \"Lectures on Network Systems\" by Francesco Bullo http://motion.me.ucsb.edu/book-lns/. These scripts are published with the MIT license. **Make sure to run the first cell above to import all necessary packages and functions and adapt settings in case.** In this script it is necessary to execute cell by cell chronologically due to reocurring examples (excepts for e.g. the Exercises in the end). (Tip: Use the shortcut Shift+Enter to execute each cell). Most of the functions are kept in separate files to keep this script neat."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.1.2 The Jordan Normal Form\n",
    "### Example 2.5 Revisiting the wireless sensor network example\n",
    "The following cells are showing the computation of the Jordan Normal Form $J$, the invertible transformation matrix $T$ and some of its dependencies. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Defining the A matrix again\n",
    "A = np.array([[1/2, 1/2, 0., 0.],\n",
    "              [1/4, 1/4, 1/4, 1/4],\n",
    "              [0., 1/3, 1/3, 1/3],\n",
    "              [0., 1/3, 1/3, 1/3]\n",
    "])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is the possibility to calculate the Jordan Normal Form directly with the package SymPy https://docs.sympy.org/latest/index.html. However, we are determining the Jordan Normal Form via determining the generalized eigenvectors (read more for literature recommendations about generalized eigenvectors in the book) with the SciPy package first to discuss some possibilities and problems with non symbolic toolboxes.\n",
    "\n",
    "### NumPy/ SciPy approach\n",
    "\n",
    "From the documentation of scipy.linalg.eig: *'Solve an ordinary or generalized eigenvalue problem of a square matrix.'*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Right eigenvectors\n",
    "lambdas, eigv = spla.eig(A)\n",
    "\n",
    "# Left eigenvectors\n",
    "lambdas2, eigw = spla.eig(A.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Due to numerical instabilities, the zero values are not reflected and it can be seen, how the expected eigenvalue of 1 is not precise. The zeros can be fixed with:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def correct_close_to_zero(M, tol=1e-12):\n",
    "    M.real[abs(M.real) < tol] = 0.0\n",
    "    if M.imag.any():\n",
    "        M.imag[abs(M.imag) < tol] = 0.0\n",
    "    return M\n",
    "\n",
    "eigv_cor = correct_close_to_zero(eigv)\n",
    "eigw_cor = correct_close_to_zero(eigw)\n",
    "lambdas_cor = correct_close_to_zero(lambdas)\n",
    "lambdas2_cor = correct_close_to_zero(lambdas2)\n",
    "\n",
    "print(\"Right eigenvectors:\")\n",
    "lib.matprint(eigv_cor)\n",
    "print(\"\\n\")\n",
    "print(\"Left eigenvectors:\")\n",
    "lib.matprint(eigw_cor)\n",
    "print(\"\\n\")\n",
    "print(\"Eigenvalues (right):\")\n",
    "lib.matprint(lambdas_cor)\n",
    "print(\"\\n\")\n",
    "print(\"Eigenvalues (left) for matching later:\")\n",
    "lib.matprint(lambdas2_cor)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are two options now for $T^{-1}$: Taking the inverse of the right eigenvectors (which contains again numerical instabilities) or building it from the left eigenvectors, what would include some sorting to match the eigenvalue order from the  right eigenvector (often it is the case, that they are already aligned since calling scipy.linalg.eig twice on a matrix with the same eigenvalues)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = eigv_cor.copy()*-1  # Rescale the eigenvectors to match eigenvalues later\n",
    "# Sorting if necessary, remember to use transpose, since in T^-1 the rows represent the left eigenvectors.\n",
    "Tinv = eigw_cor.T.copy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can simply compute J, when compared, is fairly close to the solution in the book, however, due to numerical intabilities not precise. Further on, the order of the eigenvalues might be different than the on from the book."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "J = correct_close_to_zero(Tinv@A@T)\n",
    "print(\"Jordan Normal Form via SciPy/Numpy:\")\n",
    "lib.matprint(J)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SymPy approach\n",
    "\n",
    "Now we use a symbolic toolbox package SymPy from python as a blackbox. Note, that also here the order of the eigenvalues might be different!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Asym = Matrix(A) # Sympy Matrix toolbox object\n",
    "Tsym, Jsym = Asym.jordan_form()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we can compare them with our previous results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Jordan Normal Form SymPy:\")\n",
    "lib.matprint(np.array(Jsym).astype(np.float64))\n",
    "print(\"Jordan Normal Form SciPy:\")\n",
    "lib.matprint(J)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.1.3 Semi-convergence and convergence for discrete-time linear systems\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Definition 2.6 (Spectrum and spectral radius of a matrix)\n",
    "We display the spectrum of the previous A matrix with the spectrum radius for visualization purpose. Additionally, we also show how the spectrum of a randomly generated matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax213 = plt.subplots(figsize=custom_figsize_square)\n",
    "lib.plot_spectrum(A, ax213);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_M1=8\n",
    "\n",
    "# A unifornmly distributed, positive, row stochastic matrix vs not row stochastic\n",
    "M1 = np.random.uniform(0, 1,(n_M1,n_M1))\n",
    "M1 = M1 / M1.sum(axis=1, keepdims=1) # Row-stochastic\n",
    "M2 = M1 - 0.05 # Not row-stochastic\n",
    "\n",
    "fig, (ax2131, ax2132) = plt.subplots(1,2, figsize=(custom_figsize_square[0]*2, custom_figsize_square[1]))\n",
    "lib.plot_spectrum(M1, ax2131);\n",
    "lib.plot_spectrum(M2, ax2132);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.2.1 The spectral radius for row-stochastic matrices\n",
    "### Theorem 2.8 (Geršgorin Disks Theorem)\n",
    "\n",
    "Similar to before, the Geršgorin Disks are now visualized for a row-stochastic matrix and another matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax2211, ax2212) = plt.subplots(1,2, figsize=(custom_figsize_square[0]*2, custom_figsize_square[1]))\n",
    "lib.plot_gersgorin_disks(M1, ax2211)\n",
    "lib.plot_gersgorin_disks(M2, ax2212)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.3.3 Applications to matrix powers and averaging systems\n",
    "\n",
    "### Theorem 2.13 (Powers of non-negative matrices with a simple and strictly dominant eigenvalue)\n",
    "\n",
    "Here is an example for Theorem 2.13, which shows, how the powers of primitive, row-stochastic matrices converges to rank 1. This is also done for the wireless sensor network example.\n",
    "\n",
    "#### Example 2.14 Wireless sensor network\n",
    "In the book it is shown, that the wireless sensor network matrix is primitive. Here, the eigenvectors and eigenvalues are printed again and compared for the semi convergence result $\\lim_{k \\to \\infty} A^k = \\mathbb{1}_n w^T$ to demonstrate Theorem 2.13 for a row-stochastic matrix. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Left eigenvectors of A:\")\n",
    "lib.matprint(eigw_cor)\n",
    "print(\"\\n\")\n",
    "\n",
    "print(\"Eigenvalues (left) of A:\")\n",
    "lib.matprint(lambdas2_cor)\n",
    "print(\"\\n\")\n",
    "\n",
    "print(\"Normalizing dominant eigenvector:\")\n",
    "dom_eigv = eigw_cor[:, 0] / sum(eigw_cor[:, 0])\n",
    "lib.matprint(dom_eigv)\n",
    "print(\"\\n\")\n",
    "\n",
    "print(\"Convergence result of A:\")\n",
    "lib.matprint(np.linalg.matrix_power(A, 50))\n",
    "print(\"\\n\")\n",
    "\n",
    "print(\"equals 1n*w^T\")\n",
    "lib.matprint(np.ones((4,1))@dom_eigv[:, None].T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below is a randomly generated example to show, that primitive, row-stochastic matrices always converge to rank 1. *Note: The code is not robust for semisimple eigenvalue of 1*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Creating a new random primitive (positiv), row stochastic matrix here\n",
    "n_M11=5\n",
    "M11 = np.random.uniform(0, 1,(n_M11,n_M11))\n",
    "M11 = M11 / M11.sum(axis=1, keepdims=1) # Row-stochastic\n",
    "print(\"Random primitive row-stochastic matrix M:\")\n",
    "lib.matprint(M11)\n",
    "print(\"\\n\")\n",
    "\n",
    "print(\"Left eigenvectors of M:\")\n",
    "l_M, m_eigv = spla.eig(M11.T)\n",
    "m_eigv = correct_close_to_zero(m_eigv)\n",
    "l_M = correct_close_to_zero(l_M)\n",
    "lib.matprint(m_eigv)\n",
    "print(\"\\n\")\n",
    "\n",
    "print(\"Eigenvalues (left) of M:\")\n",
    "lib.matprint(l_M)\n",
    "print(\"\\n\")\n",
    "\n",
    "# Here we check the position with numerical unprecision of the eigenvalue 1\n",
    "print(\"Normalizing dominant eigenvector:\")\n",
    "idx_dom = np.where(abs(l_M - 1) < 0.005)[0][0]\n",
    "dom_eigv_M = m_eigv[:, 0] / sum(m_eigv[:, 0])\n",
    "lib.matprint(dom_eigv_M)\n",
    "print(\"\\n\")\n",
    "\n",
    "print(\"Convergence result of M:\")\n",
    "lib.matprint(np.linalg.matrix_power(M11, 500))\n",
    "print(\"\\n\")\n",
    "\n",
    "print(\"equals 1n*w^T\")\n",
    "lib.matprint(np.ones((n_M11,1))@dom_eigv_M[:, None].T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises 2.18\n",
    "\n",
    "This section is similar to exercise 1.4, however, here we actually visualize the graph and its node values. Additionally, we show that the values converge to the initial values multiplied b the dominant left eigenvector as presented in Theorem 2.13 for row stochastic matrices.\n",
    "\n",
    "First, we define the graphs and their adjacency Matrix A and simulate the results. Then, for each graph a cell can be executed for the interactive visualization. A plot of the states is available in the Jupyter Notebook script for Chapter 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define x_0\n",
    "xinitial = np.array([1., -1., 1., -1., 1.])\n",
    "\n",
    "# Defining the 3 different systems.\n",
    "# Complete graph\n",
    "A_complete = np.ones((5,5)) / 5\n",
    "\n",
    "# Cycle graph\n",
    "A_cycle = np.array([\n",
    "    [1/3, 1/3, 0, 0, 1/3],\n",
    "    [1/3, 1/3, 1/3, 0, 0],\n",
    "    [0, 1/3, 1/3, 1/3, 0],\n",
    "    [0, 0, 1/3, 1/3, 1/3],\n",
    "    [1/3, 0, 0, 1/3,  1/3] ]  )\n",
    "\n",
    "\n",
    "# Star topology. center = node 1\n",
    "A_star = np.array([\n",
    "    [1/5, 1/5, 1/5, 1/5, 1/5],\n",
    "    [1/2, 1/2, 0, 0, 0], \n",
    "    [1/2, 0, 1/2, 0, 0],\n",
    "    [1/2, 0, 0, 1/2, 0],\n",
    "    [1/2, 0, 0, 0, 1/2]   ])\n",
    "\n",
    "\n",
    "# Defining simulation time\n",
    "ts = 15\n",
    "\n",
    "# Defining graphs for plotting later\n",
    "n = 5\n",
    "G_star = nx.star_graph(n-1)\n",
    "pos_star = {0:[0.5,0.8], 1:[0.2,0.6],2:[.4,.2],3:[.6,.2],4:[.8,.6]}\n",
    "G_cycle = nx.cycle_graph(n)\n",
    "pos_cycle = {0:[0.5,0.8], 1:[0.35,0.6],2:[.4,.3],3:[.6,.3],4:[.65,.6]}\n",
    "G_complete = nx.complete_graph(n)\n",
    "pos_complete = pos_cycle.copy()\n",
    "\n",
    "# Simulating and saving each network\n",
    "states_complete = lib.simulate_network(A_complete,xinitial, ts)\n",
    "states_star = lib.simulate_network(A_star,xinitial, ts)\n",
    "states_cycle = lib.simulate_network(A_cycle,xinitial, ts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Complete graph**\n",
    "\n",
    "Showing complete graph interactive simulation and Theorem 2.13"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax2181 = plt.subplots(figsize=custom_figsize)\n",
    "\n",
    "# If this cell is executed twice we are making sure in the following, that the previous widget instances are all closed\n",
    "try:\n",
    "    [c.close() for c in widget2181.children]  # Note: close_all() does also affect plot, thus list compr.\n",
    "except NameError:  # Only want to except not defined variable error\n",
    "    pass\n",
    "\n",
    "widget2181 = lib.interactive_network_plot(G_complete, states_complete, pos_complete, ts, fig, ax2181)\n",
    "\n",
    "display(widget2181)\n",
    "\n",
    "# Verifying the results\n",
    "eigval, eigvec = np.linalg.eig(A_complete.transpose())\n",
    "idx_dom = np.argmax(eigval)\n",
    "dom_eigvec = eigvec[0:5,idx_dom]/eigvec[0:5,idx_dom].sum()\n",
    "print(\"Showing Theorem 2.13 for the complete graph\")\n",
    "print(\"Dominant eigenvector: \\n\", dom_eigvec)\n",
    "print(\"Final values : \\n\", xinitial@dom_eigvec*np.ones(5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Star graph**\n",
    "\n",
    "Showing star graph interactive simulation and Theorem 2.13"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax2182 = plt.subplots(figsize=custom_figsize)\n",
    "\n",
    "# If this cell is executed twice we are making sure in the following, that the previous widget instances are all closed\n",
    "try:\n",
    "    [c.close() for c in widget2182.children]  # Note: close_all() does also affect plot, thus list compr.\n",
    "except NameError:  # Only want to except not defined variable error\n",
    "    pass\n",
    "\n",
    "widget2182 = lib.interactive_network_plot(G_star, states_star, pos_star, ts, fig, ax2182)\n",
    "\n",
    "display(widget2182)\n",
    "\n",
    "# Verifying the results\n",
    "eigval, eigvec = np.linalg.eig(A_star.transpose() )\n",
    "idx_dom = np.argmax(eigval)\n",
    "dom_eigvec = eigvec[0:5,idx_dom]/eigvec[0:5,idx_dom].sum()\n",
    "print(\"Showing Theorem 2.13 for the star graph\")\n",
    "print(\"Dominant eigenvector: \\n\", dom_eigvec)\n",
    "print(\"Final values : \\n\", xinitial@dom_eigvec*np.ones(5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Cycle graph**\n",
    "\n",
    "Showing cycle graph interactive simulation and Theorem 2.13"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax2183 = plt.subplots(figsize=custom_figsize)\n",
    "\n",
    "# If this cell is executed twice we are making sure in the following, that the previous widget instances are all closed\n",
    "try:\n",
    "    [c.close() for c in widget2183.children]  # Note: close_all() does also affect plot, thus list compr.\n",
    "except NameError:  # Only want to except not defined variable error\n",
    "    pass\n",
    "\n",
    "widget2183 = lib.interactive_network_plot(G_cycle, states_cycle, pos_cycle, ts, fig, ax2183)\n",
    "\n",
    "display(widget2183)\n",
    "\n",
    "# Verifying the results\n",
    "eigval, eigvec = np.linalg.eig(A_cycle.transpose())\n",
    "idx_dom = np.argmax(eigval)\n",
    "dom_eigvec = eigvec[0:5,idx_dom]/eigvec[0:5,idx_dom].sum()\n",
    "print(\"Showing Theorem 2.13 for the cycle graph\")\n",
    "print(\"Dominant eigenvector: \\n\", dom_eigvec)\n",
    "print(\"Final values : \\n\", xinitial@dom_eigvec*np.ones(5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises 2.19\n",
    "\n",
    "This exercise is about $n$ robots moving on the line trying to gather at a common location (i.e., reach rendezvous), where each robot heads for the centroid of its neighbors. The visualization of the algorithm deals with the discrete part of the task, where on can explore values of the sampling period $T$ for the Euler discretization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setup - change these parameters if wanted\n",
    "n_robots = 8\n",
    "# Number of timesteps and sampling period T (If T is too small, need very high n_dt)\n",
    "n_dt = 25\n",
    "T = 0.3  # Play around with this value, something interesting is happening around (2(n_robots-1)/n_robots)\n",
    "#T = 2*(n_robots-1)/n_robots\n",
    "\n",
    "# Set up initial position matrix and further saving variables\n",
    "current_positions = 2*np.random.random((n_robots,1))-1\n",
    "new_position = current_positions.copy()\n",
    "all_positions = np.zeros((n_dt, n_robots, 1))\n",
    "all_positions[0] = current_positions.copy()\n",
    "\n",
    "for tt in range(1, n_dt):\n",
    "    for index, own_pos in np.ndenumerate(current_positions):\n",
    "        new_position[index] = own_pos + T*(1/(n_robots-1)*(np.sum(current_positions)-own_pos) - own_pos)\n",
    "    all_positions[tt] = new_position.copy()\n",
    "    current_positions = new_position.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig, ax219 = plt.subplots(figsize=custom_figsize)\n",
    "# Set colors of robots for tracking\n",
    "all_colors = np.random.rand(n_robots,3)\n",
    "\n",
    "\n",
    "def plot_robot_pos(ax, pos):\n",
    "    # Set xlim, ylim and aspect ratio\n",
    "    ax219.set_xlim(-1.5, 1.5)\n",
    "    ax219.set_ylim(-0.5, 0.5)\n",
    "    ax219.set_aspect('equal')\n",
    "    # Add horizontal line\n",
    "    ax219.axhline(y=0.0, color='k', linestyle='-')\n",
    "    for i in range(0, pos.shape[0]):\n",
    "        # Add a robot as circle\n",
    "        bug = mpl.patches.Circle((pos[i], 0), radius=0.06, ec='black', color=all_colors[i])\n",
    "        ax.add_patch(bug)\n",
    "\n",
    "def interactive_robots(timestep):\n",
    "    ax219.clear()\n",
    "    plot_robot_pos(ax219, all_positions[timestep['new'], :]) # Take the new value received from the slider dict \n",
    "    return None\n",
    "\n",
    "# Plot initial configuration\n",
    "plot_robot_pos(ax219, all_positions[0, :])\n",
    "\n",
    "# Widget\n",
    "# If this cell is executed twice we are making sure in the following, that the previous widget instances are all closed\n",
    "try:\n",
    "    [c.close() for c in widget219.children]  # Note: close_all() does also affect plot, thus list compr.\n",
    "except NameError:  # Only want to except not defined variable error\n",
    "    pass\n",
    "\n",
    "widget219 = lib.create_widgets_play_slider(fnc=interactive_robots, minv=0, maxv=n_dt-1, step=1, play_speed=500)\n",
    "\n",
    "display(widget219)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can even compute the convergence solution in advance. Please refer and solve the given Excersice first to try to understand the following calculations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A_robot = 1/(n_robots-1) * np.ones((n_robots, n_robots)) - n_robots/(n_robots-1)*np.identity(n_robots)\n",
    "eigval, eigvec = np.linalg.eig(A_robot)\n",
    "idx = np.argmin(abs(eigval))\n",
    "z_eigvec = eigvec[:,idx]/np.sqrt(np.sum(eigvec[:,idx]**2))\n",
    "\n",
    "final_values = z_eigvec[None, :] @ all_positions[0] @ z_eigvec[None, :]\n",
    "print(\"Final values :\", final_values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.11"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": true,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}