lectures/lecture11/lecture-11.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Lecture 11. Intro to iterative methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## The main topics for this part\n",
    "\n",
    "Concept of **iterative methods** for linear systems: \n",
    "- Richardson iteration and its convergence\n",
    "- Chebyshev iteration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Iterative methods\n",
    "\n",
    "- If we want to achieve $\\mathcal{O}(N)$ complexity of solving sparse linear systems, then direct solvers are not appropriate.\n",
    "\n",
    "- If we want to solve partial eigenproblem, the full eigendecomposition is too costly.\n",
    "\n",
    "- For both problems we will use iterative, Krylov subspace solvers, which treat the matrix as a **black-box** linear operator."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Matrix as a black box\n",
    "\n",
    "- We have now an absolutely different view on a matrix: matrix is now a **linear operator**, that acts on a vector,  \n",
    "and this action can be computed in $\\mathcal{O}(N)$ operations.\n",
    "\n",
    "- **This is the only information** we know about the matrix: the <font color='red'> matrix-by-vector product (matvec) </font>\n",
    "\n",
    "- Can we solve linear systems using only matvecs?\n",
    "\n",
    "- Of course, we can multiply by the colums of the identity matrix, and recover the full matrix, but it is not what we need."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Richardson iteration\n",
    "The simplest idea is the **\"simple iteration method\"** or **Richardson iteration**.  \n",
    "\n",
    "\n",
    "  $$Ax = f,$$\n",
    "  $$\\tau  (Ax - f) = 0,$$\n",
    "   $$x - \\tau (Ax - f) = x,$$\n",
    "   $$x_{k+1} = x_k - \\tau (Ax_k - f),$$\n",
    "   \n",
    "   where $\\tau$ is the **iteration parameter**, which can be always chosen such that the method **converges**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Connection to ODEs\n",
    "\n",
    "- The Richardson iteration has a deep connection to the Ordinary Differential Equations (ODE).\n",
    "\n",
    "\n",
    "- Consider a time-dependent problem ($A=A^*>0$)\n",
    "\n",
    "$$\\frac{dy}{dt} + A y = f, \\quad y(0) = y_0.$$\n",
    "\n",
    "- Then $y(t) \\rightarrow A^{-1} f$ as $t \\rightarrow \\infty$, and the **Euler scheme** reads\n",
    "\n",
    "$$\\frac{y_{k+1} - y_k}{\\tau} = -A y_k + f.$$\n",
    "\n",
    "which leads to the Richardson iteration\n",
    "\n",
    "$$\n",
    "    y_{k+1} = y_k - \\tau(Ay_k -f)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Convergence of the Richardson method\n",
    "\n",
    "- Let $x_*$ be the solution; introduce an error $e_k = x_{k} - x_*$, then  \n",
    "\n",
    "$$\n",
    "     e_{k+1} = (I - \\tau A) e_k,\n",
    "$$\n",
    "\n",
    "therefore if $\\Vert I - \\tau A \\Vert < 1$ in any norm, the iteration converges. \n",
    "\n",
    "- For symmetric positive definite case it is always possible to select $\\tau$ such that the method converges.\n",
    "\n",
    "- What about the non-symmetric case? Below demo will be presented..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Optimal parameter choice\n",
    "\n",
    "- The choise of $\\tau$ that minimizes $\\|I - \\tau A\\|_2$ for $A = A^* > 0$ is (prove it!)\n",
    "\n",
    "$$\n",
    "  \\tau_\\mathrm{opt} = \\frac{2}{\\lambda_{\\min} + \\lambda_{\\max}}.\n",
    "$$\n",
    "\n",
    "where $\\lambda_{\\min}$ is the minimal eigenvalue, and $\\lambda_{\\max}$ is the maximal eigenvalue of the matrix $A$.\n",
    "\n",
    "- So, to find optimal parameter, we need to know the **bounds of the spectrum** of the matrix $A$, and we can compute it by using **power method**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Condition number and convergence speed\n",
    "\n",
    "Even with the optimal parameter choice, the error at the next step satisfies\n",
    "\n",
    "$$\\|e_{k+1}\\|_2 \\leq q \\|e_k\\|_2 , \\quad\\rightarrow \\quad \\|e_k\\|_2 \\leq q^{k} \\|e_0\\|_2,$$\n",
    "\n",
    "where \n",
    "\n",
    "$$\n",
    "   q = \\frac{\\lambda_{\\max} - \\lambda_{\\min}}{\\lambda_{\\max} + \\lambda_{\\min}} = \\frac{\\mathrm{cond}(A) - 1}{\\mathrm{cond}(A)+1},\n",
    "$$\n",
    "\n",
    "$$\\mathrm{cond}(A) = \\frac{\\lambda_{\\max}}{\\lambda_{\\min}} \\quad \\text{for} \\quad A=A^*>0$$\n",
    "\n",
    "is the condition number of $A$.\n",
    "\n",
    "Let us do some demo..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Maximum eigenvalue = 3.6825070656623633, minimum eigenvalue = 0.08101405277100539\n",
      "Condition number = 45.45516413147915\n",
      "[0.9186479  0.94484203 0.95170654 0.95466793 0.95595562 0.95651542\n",
      " 0.9567591  0.95686534 0.95691172 0.95693198 0.95694084 0.95694472\n",
      " 0.95694642 0.95694716 0.95694748 0.95694763 0.95694769 0.95694771\n",
      " 0.95694773 0.95694773 0.95694773 0.95694774 0.95694774 0.95694774\n",
      " 0.95694774 0.95694774 0.95694774 0.95694774 0.95694774 0.95694774\n",
      " 0.95694774 0.95694774 0.95694774 0.95694774 0.95694774 0.95694774\n",
      " 0.95694774 0.95694774 0.95694774 0.95694774 0.95694774 0.95694774\n",
      " 0.95694774 0.95694774 0.95694774 0.95694774 0.95694774 0.95694774\n",
      " 0.95694774 0.95694774 0.95694774 0.95694774 0.95694774 0.95694774\n",
      " 0.95694774 0.95694774 0.95694774 0.95694774 0.95694774 0.95694774\n",
      " 0.95694774 0.95694774 0.95694774 0.95694774 0.95694774 0.95694774\n",
      " 0.95694774 0.95694774 0.95694774 0.95694774 0.95694774 0.95694774\n",
      " 0.95694774 0.95694774 0.95694774 0.95694774 0.95694774 0.95694774\n",
      " 0.95694774 0.95694774 0.95694774 0.95694774 0.95694774 0.95694774\n",
      " 0.95694774 0.95694774 0.95694774 0.95694774 0.95694774 0.95694774\n",
      " 0.95694774 0.95694774 0.95694774 0.95694774 0.95694774 0.95694774\n",
      " 0.95694774 0.95694774 0.95694774]\n",
      "Theoretical factor: 0.956947735792311\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaMAAAEcCAYAAABqCdtUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl4FFXWBvD3JOwICWFTQIUkoOwQggsIUUnALSoQQMcZdUYhuHyKqAQcncVxhKCDu5jA6My4jLIpMqKQoLK5kYRVBEMCoqKAQALIDuf7o6qhabqT7qQ61V39/p6nn5DuStW5ldCnb91b94iqgoiIyE5RdgdARETEZERERLZjMiIiItsxGRERke2YjIiIyHZMRkREZDsmIyIish2TERER2Y7JiIiIbFfL7gDCRbNmzbRt27Z2h0FEFFYKCwt/UdXmlW3HZOSntm3boqCgwO4wiIjCioh85892vExHRES2q3YyEpHGVgRCRESRq8rJSETaicgmAGUisktEBru9NlREii2JkIiIHK86PaNJADIBNAEwAsDNIjIEAFR1NoCE6odHRESRoDrJqEBVF6lquarmq+pwAE1F5ErzdRZKIiIiv1QnGZV5PqGq0wCIiPSsxn6JiCjCVKtnJCIjRaRYRNq6nlTVReY/y6sTGBERRY4qJyNVXQlgBoDhqrrFy2u9qheaM5Tu3I9JH24Ay7sTEflWrand5njRSh+vba7Ovp3i4w078MriEkxdXGJ3KEREIYs3vQbZHZe1Q3r3VnhqwUZ8smGH3eEQEYUkr8sBicgrAJJR/RlxAqBEVUdUcz9hS0QweWg3lO7cj/v+uxLv3dsXCc3PsjssIqKQIhzLqJiIpANIT0xMHFlcXPX7eH/YcwDXv7gcsfVr462Rl+DsmHrWBUlEFKJEpFBVkyvbjpfpKqGq81R1VExMTLX206ZJA7zy217YvvcQrn9xGVZ9f8bMeCKiiMVkVIMuaheHOXf3RZ1aURie8zneW/mj3SEREYUEX2NGkwC0s+gYpao6waJ9hb0Lzm6E9++9DHe9UYgx76zCii278ei1nVC/TrTdoRER2cZXPaOJAOIsOsZui/bjGHEN6+D1Oy7GPxZuRM6SUny1eTeev7knOp7DBdCJKDJxAoOfkpOTNRjF9ZYW78TYGatRfvAoHkzrgDv7xSM6Siw/DhGRHTiBIUz0a98cH93fD5d3aI6JH27A8JzPsfmXX+0Oi4ioRgWUjESksYjMEJHj5mOFiDwYrOAiRdOz6iLnd73w7IgeKN6+D1c/twTTl5bi+An2WokoMgTaM5oMYAWA4QDGA9gM4CmzuF53q4OLJCKCG3u2Rt7YFPRNaIYnPvgGGa98hk079tkdGhFR0AWajApV9SlVnW1+Ha6qUTCK7M1yX72bqqZl43qYflsynh3RA5t/+RXXPLcML35cjKPHT9gdGhFR0PhMRiIyVUQGi4j7FC+v141UdRaA3jCSElXTyV7SAylI69QSTy/8Fte/uBxrf2BVDiJypop6RpkAZgPYY9YsmgogVkQmettYVcsAcGlqCzVvVBcv3ZKEnN/1wq79h3Hjy8sxcf43OHjkuN2hERFZqqJklG9egusNIBdAAowxoyxXchKRO0WkBwCYl+hYwygIBnU+G3ljUzA8uQ1ylpTiqueW4LNNv9gdFhGRZXzeZyQi7bzVJDJLiqfCSFKpAGJhXL4rAjDMs9CeUwTrPqNAfVbyCx6ZsxZbdh3AiORz8cg1HRHToLbdYREReeXvfUbVvulVRGJU1fGDGaGSjADg0NHjeCb/W0xfuhlNGtTB4zd0xtVdzoYIb5YlotBSYze9RkIiCjX1akdjwtUdMfeevmjZuC7ufrMIo14vxM/lh+wOjYioSrgCQxjr0joGc+/piwlXX4ilxTuRNmUx3vzyO5zgzbJEFGaYjMJcregoZKYkYMGY/ujaJgZ/fHcdbsr9AiU799sdGhGR35iMKiEi6SKSW14e2lcjz2/aEG/eeTEmZ3TDxu37cPWzS/Hix8U4cow3yxJR6OOq3X4KpQkMldmx7xD+Om89PljzEy5o2QiThnZFz/Oa2B0WEUUgrtodwVo0qoeXfpOE6bcmY++hoxgy9TP8dd7X+PXwMbtDIyLyqsrJSERiROQhK4Mha6V2aomFD/TH7y45H68t34KBzyzBJxt32B0WEdEZqtMzioNx4yuFsEb1auPxG7pg1uhLUb9ONH7/2gqMeXsldu0/bHdoREQn8TJdhEhuG4cP7rsM9w9ojw/W/oTUKYsxp+gHcMyQiEIBk1EEqVsrGg+kdcAH9/VD22YNMXbGatz66lf4fvcBu0MjogjHZBSBOrRshFmj++Cv13dG0Xd7MPAZVpYlInsxGUWo6CjBbX3aYuHYFFya0BRPfPANhry8HOu37bU7NCKKQExGEa51bH3887ZkvHBzT/yw5yCuf3EZnlqwAYeOsmYSEdUcJiOCiCC9eyvkj03BjT1b46VPSnDNc0vxRekuu0MjogjBZEQnNWlYB08P64437rgYR0+cwE25X2DCnLUoP3jU7tCIyOGYjOgMl7VvhoVjUjCqfzzeWbEVaVMW46N1P9kdFhE5GJMReVW/TjQeuaYj5t5zGZqdVRej3yhC5usF2L6XNZOIyHrVTUYsLepwXdvEYO69fZF11YX4dONOpP6DNZOIyHpVTkaquhlAloWxUIiqHR2Fuy43aiZ1aW3WTJrGmklEZJ1q9YzMhEQRom2zhnhr5MXIHtoVG37ai6ufM2omHT3OmklEVD0cM6KAiAhG9D4P+Q+mIK1jSzy98Fukv7AMq74vszs0IgpjTEZUJS0a1cNLtyQh93e9sOfAEQx5eTken7eeNZOIqEqYjKhaBnY+G3ljU/Cbi8/Dq8s3Y+AzS/ApayYRUYCYjKjaGterjSdu7IqZoy9FvdpRuP21FXjgnVXY/esRu0MjojDBZESW6d02DvPv74f7BrTH/9ZsQ+qUxXhv5Y+smURElbIkGYnIDBEZ4vb9ABFpbMW+KbzUrRWNsWbNpPObNsCYd1bh9tdWsGYSEVXIqp5RjqrOcX2jqosAJFu0b0uJSLyIjBORVPNrrN0xOZGrZtJf0jthxZbdGPTsEry6bDNrJhGRV1Ylo1IRaQsAIrJQRFYASLNo31abqaqTVTUfQC6AaXYH5FTRUYLb+7ZD3tgUXNwuDo//bz2GTP0MG35mzSQiOp1U9Xq+iBQDyIfx5v6xiPQAEA9jiaB8VS23Lswzjp0BoMxMKJ6vJcHolZWa8ZS6tjNfy1bVNLft96hqk8qOmZycrAUFBVY1IeKoKt5fvQ2Pz1uP8oNHMTolAfdemYh6taPtDo2IgkhEClW10itltapxjNkAVgAYLiI5ADYDaAdgWJATUSqM3swwL6/F48xkM1NESlXVlZw8787cLSJJqloUrJjJuFn2hh6t0b99czzxwTd48ZNNmL/uJ0wa0g0XtYuzOzwisll11qYbr6qzVXW0qrYHkAlgMoBHRGSTiEy1LEqcHOvJgZFQdvvYLBNAjsdzOQCyzX/zXc9mTRrWwT+Gd8frd1yEo8dPYHjO53jk3bXYe4g1k4gimWVTu1V1s6pOU9XhqpoIIzFZRlVLVTVTVXMr2CwDgGcPp8B8HjCSmOeEBSYoG/Rr3xwLxvTHyH7t8PZXW5H6j8X4aN3PdodFRDYJOBn5O2W7phdRNWfFndFrUtUy8/V4GONIZyQfXqKzR4M6tfDHazvhvXv6oulZdTH6jUKMfr2QNZOIIlBAyUhEZgDY45o55/b8lRbGVFVxwKnk40W8Z9IxE9QZkyCoZnVrE4v37+2LcVddgI837kDqlMX471dbWTOJKIIE2jPKAzBcVbd4PF8oIg/ZfKOrv/cLjXTdZwTj8t1IXxuKyCgRKRCRgp07d1oSJHlXOzoKd1+eiAVj+qPTOY0xYc5a3DztC5SyZhJRRAg0GcXCSEinUdVyVX0awHBLogoiVS1y3WdkfvVZ+0BVc1U1WVWTmzdvXpNhRqx2zRrivyMvwaQhXbH+p7246rmleOmTTayZRORwASUjVX0KwEwRmSgiVwQppmrhigrhLypKcNNF52HR2BSkdmyBpxZsRPoLy7CaNZOIHCvQMaNXYNzUmgZgkYgcF5EVIjJVRCYC6BWMIP1Uan49bYKCW3IqBYWVFo3r4eVbeiHHrJk0+OXl+Nv/1uPAEdZMInKaQG96LVHV0a5vzHEX10MBDLAwtoCoapmIlML71O0y86bXgIlIOoD0xMTE6oZIVTSo89m4NKEpsj/cgH8u24wFX/+Mvw/uipQOvHRK5BSBjhmddp3EHHcZby71MAr2jxnl48wFWpNQjRlzqjpPVUfFxMRUKzCqnsb1auPvg7tiRualqFMrCre9+hXGsmYSkWMEmozyReROH6/tqW4wAYiD99lzWThzmaBM83lygIvaxWH+ff3wf1cm4v3VrJlE5BQBL5QqIjEw3vBnqOpe87l2AEpglJK4y/IocXLsZwKMG1szYIwB5QPIU9VZbtslARgBY928eABF3hZUDRQXSg09G37ei/Gz12LV92VI6dAcfx/cBW2aNLA7LCJy4+9CqVVetdvLAQfAWCG7RldeqClMRqHp+AnFfz7fgqcWbAQAPDTwAtzWpy2io8TewIgIgP/JyMq16RY5NRFR6IqOEvy+bzssfKA/LjJrJg1lzSSisGNZMgIAq1fqDgUiki4iueXlQauKQRZo06QBXru9N567qQe27j6A655fhn8s3IhDR4/bHRoR+aFal+nM5X/iYUwocNUSampRbCGFl+nCx+5fj+CJ/63HnJU/Ir55Q9ZMIrJRtceMzMqtriSTAGP2mnvi8ZzNJgBUVR1ZupPJKPws+XYnHnl3LX7YcxC3XHwesq6+EI3r1bY7LKKIYkUyWohTN7O6RoNLzUcZjISUA6NkQxmMSq8z/DloOGIyCk8HjhzDlIXf4tXlm9G8UV08fkMXDOp8tt1hEUWMak9gUNWBMG5inQagiapGqWqiqg5U1eEACsxierPNyQuuqdZEIaNBnVp49LpOePfuvohrWBeZrxfi7jcLsYM1k4hCSoUTGFR1lrn8zzDzst1pL3vZfryVwRFZpfu5p2om5X+zAwOmLMbbX23lzbJEIcKv2XSqOh1AeYgU0atRnE3nHK6aSR/d3w+dzmmM8WbNpM2//Gp3aEQRz++p3aq6WVU/FpEe5ooLEYFr0zlPfPOzTtZM+nrbXgx6dglrJhHZLOD7jFR1lXlza5HNlV2Jqsy9ZtKAC42aSde/uBxrfmDNJCI7VPmmV1WdBqCpiFzpSkoiMsSyyIhqQIvG9TD1t0bNpN2/HsaNLy3H3z9gzSSimlatFRhcl+4ApJmF96ZZExZRzRrU+WzkjU3BTRedh2lLN2PQs0uw5NuddodFFDEsWQ5IVWcDyIX3sg5EYaFxvdp40qyZVDs6Cre++hXGzliFPayZRBR0Vi6UWgRgkVX7I7KLq2bSvVck4v1V2zBgymLMXcWaSUTBZOlCqTAK2TkKp3ZHpnq1o/HQoAvwv/suw3lxDXD/26vw+3+twA97DtgdGpEjWVbPyOm4HFDk8qyZ9PCgC3DrpayZROSPGq9nRORUrppJC8b0R++2cfjrvPXIeOUzbPx5n92hETmG156ROTMuGV6W/Al0/wBKVHVENfdjO/aMCABUFXNXbcPj/1uPfYeO4q6UBNxzZSLq1nLkYvVE1VbjZcedjsmI3LnXTEpo3hCThnZD77asmUTkiZfpiIIormEdTBnRA//+w0U4dPQEhr3yOR59by32HTpqd2hEYYnJiKgaUjo0x8IH+uOOy9rhrS+3Im3KEuSt3253WERhx9eY0SQAVi2GWqqqEyzal214mY4qs/r7MmTNXoMNP+/DtV3PwZ+v74QWjerZHRaRrao1ZiQiMTDKi1tht6qG7U06IpIOID0xMXFkcXGx3eFQiDt6/ARyl5TiuUXFqFcrCn+8tiOGJ58LEU4Dp8jECQwWY8+IAlG6cz8mzFmLLzfvxqXxTfHkkK5o16yh3WER1ThOYCCykatm0sQhXbFuWzmuenYJXv6UNZOIfGEyIgqSqCjBzRedh/yxKbj8guaY/BFrJhH5wmREFGQtG9dDzu+S8cpvk7Brv1Ez6cn53+DgkeN2h0YUMpiMiGrIVV3OQd7YFIzofR5yl5Ri4LOLsaz4F7vDIgoJTEZENSimfm1MHNIVb4+6BLWiovDbf36Jh2auZs0kinhMRkQ2uCS+KT68vx/uuSIB7638EWnPLMb7q7exZhJFLCYjIpvUqx2NhwddiHn/dxlax9bHff9diTv+XYAfyw7aHRpRjWMyIrJZx3MaY87dffHYdZ3weckuDJyyGP/+bAuOn2AviSIHk1ElWOmVakJ0lOCOy9ph4QP90attHP78/tcY9spn+HY7ayZRZLBsBQYR6QGjBlKBqq6yZKchhCswUE1RVby36kc8Pm899h8+hrsuT8Q9VySwZhKFJTtWYHgEQC6AWSJSLCJDLNw3UcQQEQzu2Qb5Y1NwXbdWeH5RMa59fhkKtuy2OzSioLEyGWUB6KWqiaraHsAiC/dNFHGanlUXz4zogX/9vjcOHjmODNZMIgezLBmp6mZVXen2PQdZiCxw+QUtsPCB/vhD33Z4kzWTyKE4gYEoDDSsWwt/Su+EOXf1QUz92hj5nwLc/WYhduw7ZHdoRJao0gQGc7KC13pHqvpxdYMKRZzAQKHiyLETyF1Sguc/3oR6taLw6LWdMCy5DWsmUUgKSj0jERkKYIbrWy+bqKo6csoPkxGFmpKd+zFh9lp8tWU3+iQ0xZODu6ItayZRiAlWMtpk/jMLQKm3bdzHjZyEyYhC0YkTirdXfI+J87/BkeMn8EBaB9x5WTvUiuYVeAoNwUpGJwBkqOqc6gQXjpiMKJRt33sIf5q7Dgu+3o7OrRoje2g3dGkdY3dYREG7zygfPsaKiMg+7jWTduw7jBteWo6JrJlEYaRWgNtnApghIjEw7iM6o2Slqm6xIC4iqoKrupyDSxOaYdKH3yBnSSk+XPcznhzcFZe1b2Z3aEQVCvQy3QAAMwH47P9zAgNRaPi8ZBceeXctNv/yKzJ6tcGj13ZEbIM6dodFESZYY0YFAOJhLPtT4m0bVZ3m9w7DgIikA0hPTEwcWVxcbHc4RAE5dPQ4Xvi4GDmLSxHboDb+nN4Z13U7h9PAqcYEcwLDKFWdXp3gwhF7RhTO1m/bi/Fz1mDND+W48sIW+NuNXdA6tr7dYVEECNYEhiIAXK2RKMx0atUY797dF49e25E1kygkBZqMRgLIFpErRKRxMAIiouCIjhLc2S8eCx/oj6Tzm7BmEoWUQJPRIgAJMKZ47xGR4x6PY9aHSERWOjeuAf7zh4swZXh3lP7yK659fimm5H2Lw8c4DZzsE+jU7mkA2K8nCnMigiFJbdC/Q3P87X/r8fyiYsxf+xOyh3ZFr/N5KyHVvEAnMDRW1b1BjCdkcQIDOdknG3fg0XfXYVv5Qfz24vMx7qoL0KhebbvDIgcI1gSGQhF5sIoxEVGIusKsmfT7Pu3wxpffIW3KEuSzZhLVoECT0ccAbgpGIERkL8+aSXf+pwD3vFWEnfsO2x0aRYCAkpGqZgIQEflIRM4PUkxEZKOe5zXBvP+7DA8N7IC8r7cjdcpizCj4HlWpfUbkr4CSkYjshrECw0AApZxNR+RMdWpF4d4r22P+/f1wQctGGDdrDX77zy/x3a5f7Q6NHIqz6YjIp8QWZ+HtUZfgvyu2YtL8DRj07BI8kNoBd7BmElmsqmXHrwQwDEYvqRTATKeWG3fhbDqKdD+XH8Jjc9chbz1rJpH/gjWbDiIyFcZNryMANDW/5onI2wFHSURh4+yYesj9XS9MvYU1k8h6gY4ZjYRR0yhNVeNUNVlV42AkpOGc9k3kbCKCq7ueg/wHUjCsVxvkLCnFoGeXYPmmX+wOjcJcoD2jTADjVHWR+5OqOgvAeACjrQqMiEJXTIPamDS0G/478hJECXDL9C/x8MzVKDtwxO7QKEwFmoySYKzc7U0JjDEkIooQlyY0xUdj+uOuyxMwZ+WPSJ2yGPNWb+M0cApYoMkoH0CGj9dGwHeiIiKHqlc7GllXXYj37+2Lc2Lq4//+uxIj/1OAn8oP2h0ahZFAk9F4AKNF5GUR6SEijc2vUwEMBZBlfYj2EpF0EcktLy+3OxSikNa5VQzevbsPHr22I5Zt+gVpU5bg9c+34ARrJpEfAp7aLSIZMMqOu+Z0uuoXZzqt5Lg7Tu0m8t/WXQfwx/fWYmnxL+h1fhNMGtIV7Vs2sjssskFQyo57HCAVQDsY9xkVqKqjuw5MRkSBUVXMKfoRf/tgPX49fAz3XJGIuy5PQN1a0XaHRjUo6Mko0jAZEVXNL/sP4/F56/H+6m1o3+IsTBraDb3Ob2J3WFRDgnbTq7nzxr4eVdkfETlXs7Pq4vmbe+K123vj18PHkPHKZ/jz3HXYf5hLWdIpgd702lNEigHsqeBBRHSGKy5sgYVjU3DbpW3xny++Q9qUxVj0DWsmkSHQSq8FMO4lyoVxX9EZnDqJgZfpiKxTtHUPxs9eg2+378d13c7Bn9M7o3mjunaHRUEQlDEjETkBYJSqTq9OcOGIyYjIWkeOncAri0vw4sebUL9ONB69tiMyerWBiFT+wxQ2gjVmlA9gd9VCIiI6pU6tKNw3oD3m338Z2rc4Cw+zZlJECzQZZQHIFpHBnKxARFZIbNEIMzIvxd9u7ILV35dj0LNLkLukBMeOn7A7NKpBgV6m2w0gFr4L7KmqBlqwLyzwMh1R8P1UfhCPvfc18r/Zji6tG2PSENZMCnf+XqZjpVciChnnxNTHtFt7Yf7an/Hn97/GDS8tx8h+8RiT2h71avNmWSfjTa9+Ys+IqGaVHziKJ+d/g3cKvsf5TRtg4uCu6JPYzO6wKEBBvemViCjYYhrURnZGN7w18mIIgN9M/xLjZq1G+YGjdodGQcBkREQhrU9Cs5M1k2YX/YgBUxbjgzU/sWaSwzAZEVHIc6+ZdHZMXdzzVhFG/qeQNZMchMmIiMJG51YxeO/uvvjjNR2xbNNO1kxyECYjIgortaKjMLJ/PBaOSUGPc2Px2NyvMTznc2zasc/u0KgamIyIKCyd17QBXr/jIjw9rDuKd+zHNc8tw3P5xThyjDfLhiMmIyIKWyKCjF5tkD82BYO6nI1n8r/FdS8sReF3LCAQbpiMiCjsNW9UFy/c3BOv3p6M/YeMmkl/ef9r1kwKIz6TkYicEJHjAT74myci21x5YcuTNZP+/fkWDJyyGB9vYM2kcFDRckBPgUv/EFGYOatuLfzl+s5I794K42evwR/+VYD07q3wp+s6sWZSCIvY5YBEJE9V0/zdnssBEYWfI8dOYOqnJXjxk2I0rFsLj17bCUOTWrNmUg3ickA+iEiqiIwCkGp3LEQUXHVqReH+1PaYf18/JDQ/Cw/NXI3f/fMrbN11wO7QyEOVekYiMgRG+fEzqOrTVQpEJANAmarme3ktCUAygFLzuKXetgvweKqqfn88Ys+IKLydOKF448vvMPmjjTh+QvHgwA64vU9b1IqOuM/kNSpYZcfbASgA0MR8SgG4v6Grqga8zruIpAKYCWCYZ5IRkXgAOe6X1ERkJoAsVS0N9Fhu+2AyIopA28oO4k9z1yH/mx3o2joGk4Z2RedWrJkULMG6TJcDYDOAXqoaBaAcQBKAOBglyTMDDDJeRHJg9HZ8lTPPNI/rGUd2IMciIgKAVrH1Me3WZLz4m574qfwgrn9xOSZ9uAGHjh63O7SIFmjP6ASAVFX92Py+AMCdqrpKRAYAmKSqvasUiEgJgEwvPaMSAGnuvSARiQWwx9WzMceAEirYfZ6X/bJnRBThyg4cwZPzv8GMgh+MmklDuqJPAmsmWSlYlV7LALj3Z0thjOWsglGOPCnA/VXITDpn9JpUtUxEICLxqlqqqrlWHpeIIkNsgzqYnNEdN/RojUfeXYvfTPsSI5LPxSPXdERMg9p2hxdRAr1Mlw/gJrfvZwDIMic0ZMNIVlaKA4zk4+N1r5MoiIgC0TexGT66vz8yU+Ixq+gHDJiyGPPXsmZSTQo0GWXBGDMCAKjqLABbAMyCkThGWhaZIdbi/UFEkkRknPnvbHPyhK9tR4lIgYgU7Ny50+pQiCiE1K8TjQlXd8Tce4yaSXe/yZpJNcmSm15FpJ2qbq58ywr3ccaYkTmlu9Db2I6IKIyxpGpN8fYXx4yIIsex4yfw6vLNmJL3LWpFRSHr6gtxy0XnISqKN8sGqkZveq1uIqqMOXZERFQjakVHYVT/BCwY0x/dz43BY++tY82kIAsoGYlI48oeFsfnmkEX5xFHrMfrRESWO79pQ7xxx8V4KqPbyZpJzy9izaRgCLRnVAZgTyUPy5gTF0px5thRHIzVGoKejEQkXURyy8vLg30oIgpBIoJhyecif2wKBnZuiSl5Rs2koq2smWSlQJPRaB+P6TBWYgjoplc/5cOYPu4uyXw+6FR1nqqOionhHdpEkax5o7p48TdJmH5rMvYdOoahU1kzyUoB3WdUwf0800RkEyq+8bQycfA+ey4LxlJB7sfORHASHxFRhVI7tcTF8XF4esFG/PvzLchbvx1P3NgFV1zYwu7QwpplJSREpCeAGaraPoCfiQUwAcb9QhkwLsnlw1gxYZbbdkkARgBYYW5bVFOz6Fw4m46IPBV+twfjZ69B8Y79uL57K/wpvROancWaSe6CslBqJQd8GMZyQAEvlBoOmIyIyJvDx45j6qcleOmTTWhYtxYeu7YThrBm0klBmdpdUSlyGCswFFU1YCKicFS3VjTGpHbA/Pv6Ib5ZQzw4czVuffUrfL+bNZMCEehCqdnwXYq8BMZlOkdNOxORdADpiYmJI4uLi+0Oh4hCmKtmUvaHG3BCgbFpHfD7vpFdM6nGL9M5HS/TEZG/tpUdxGPvrcOiDTvQrU0MJg3phk6trL4NMzyw7DgRkU1axdbH9NuMmknbyg4i/cVlmPwRayZVxOfU7qqupqCqe6seDhGRM4gIruvWCpclNsMTH3yDlz8twYfrfsbEIV1xSXxTu8MLORX1jPxZbSEUSLXrAAAQ+0lEQVSoKzAQEYW72AZ18PSw7njjjotx/ITiptwvMH72GpQfOGp3aCGlopteR+PMyQqxMGbNZeH02kWJAB4GMM7S6IiIHOKy9s2wYEx/PJv/LaYv24xFG3bg8es746ouZ3MaOAKfTTcVRkmH6V5eexhAnKpOsDA+23E2HRFZbd2P5ciavQZfb9uLgZ1a4vEbuuDsmHp2hxUUwZrAMAIeJcDdFMFYRcFRuDYdEVmtS+sYzL2nLyZcfSGWFO9E2pTFeOOL73DiROTObq7KbLpRPp4fBo9SD0RE5F2t6Chkphg1k7q2icGj763DiNzPsWnHfrtDs0WgyWgSgIEi8pGIDBaRHiJypYi8AyNJ+VpIlYiIvDi/aUO8eefFmJzRDd9u349rnluKFyKwZlLAN72KyDgA42FMZlAYpSPKAExU1acsjzBE8KZXIgq2HfsO4a/z1uODNT/hgpaNMGloV/Q8r4ndYVVL0FdgMFfSbgegVFVXVmknYYTJiIhqSv767Xhs7jr8vPcQbu/TFg8NvAAN6wZU8Sdk+JuMqtw6VS0CF0YlIrKcq2bSUws24l+fbcHCr7fjicFdcMUFzq2Z5HPMSEQWisgxj+d8rtptPhxX8pBlx4nIDo3q1cbjN3TBrNGXon6daPz+tRUY8/ZK7Np/2O7QgsLnZTpzbChZVYe7PVfRqt0AAFUdb2mEIYKX6YjILoePHcfLn5Tg5U834ay6tfDYdZ0wuGd41Eziqt0WYzIiIrt9u30fxs9eg6KtZejXvhmeHNwV58Y1sDusCgWruF5bEenh8dxDIvKOiNwRaJBEROS/Di0bYeboPvjr9Z1R9N0eDHxmCaYvLcVxB9wsG+h9RjkAMl3fiMhCGGvVJQDIFZEnLYyNiIg8REcJbuvTFnljU3BpQlM88cE3GPLycqzfFt4FEwJNRmkA8gBARHoCSAWQZXbBJsBYhYGIiIKsVWx9/PO2ZLxwc0/8sOcgrg/zmklVWQ6o1PyaCmMyg2vVhRIA8VYERURElRMRpHdvhfyxKbixZ2u8/GkJrnluKb4o3WV3aAELNBmVAhhuFt4bAaDIrZhePE4lKiIiqiFNGp6qmXTMrJk0Yc5alB8Mn5pJgSaj8eZjD4AkGHWNXDIBzLIorpDB+4yIKFy4aiaN6h+Pd1ZsRdqUxfho3c92h+WXqqxN1w5GIspX1XK350cCmOH+nJNwajcRhZO1Pxg1k9b/tBeDOhs1k1o2rvmaScGqZwQY40TJAN4xL9dBRGJgJKjwn19IROQAXdvEYO69fZF11YX4dONOpE5ZjLe+3BqyNZMCvc9oAIyJCmnmIw4AzN5QAowSE0REFAJqR0fhrsuNmkldWsXgkXfX4qZpX6BkZ+jVTAq0Z5QNYJbZ5fJchyIXxqQGIiIKIW2bNcRbIy9G9tCu2PDTXlz93FK89MkmHD0eOjWTAk1GSQDe8fHaHhg1joiIKMSICEb0Pg/5D6YgrWNLPLVgI9JfWIZV35fZHRqAwJNREYzxIm8yAeRXLxwiIgqmFo3q4aVbkjDt1mSUHTiKwS8vx+Pz1uPXw/YWXQi0ntEkADPEWCpWAcSYa9Vlw7gJtpfF8RERURCkdWqJS+LjMPmjjXh1+WYs+Ppn/H1wF1xuU82kgHpGqjoLwHAYvSCB0VMqAtAbwHBVXWV5hEREFBSN6tXG327sgpmjL0W92lG43ayZtPvXIzUeS3XKjveEueoCy44TEYW3w8eO46VPSjDVrJn0p/ROuLFH9Wsm2VLPSER6OK13JCLpANITExNHFhcX2x0OEVFQfbt9H7Jmr8HKrWXo36E5/n5jl2rVTLLkplcRaWzWKtolIsW+ahaZ280AUFjFeEOWqs5T1VExMTF2h0JEFHQdWjbCrNF98Jf0TijYshsDn1mC91dvC/pxKxszKoJxc+siGGNEuSLyoPsGIjIRxrTuVBjr1hERURiLjhLc3rcd8samoG9iU8Q3axj0Y/qcTSciD5v/bOdab05ECgFMFhHXDa7ZMJLUeFV9KtjBEhFRzWkdWx/Tb+tdI8eqaGr3CACTPBY+zQKwEEAZjCQ0GcBEpy6OSkRENaOiZJQEwHP62Arzay6M3hCTEBERVVtlY0anrRPhlnxymIiIiMgqVSkhQUREZKnKlgMaKiKePSCFUXrc67xxVZ1uSWRERBQxKktGvmbI+ZrCrQCYjIiIKCAVJSMuekpERDXCZzKKhPXmiIgoNFi6Np2TichOAN9V8cebAfjFwnDCRSS2OxLbDERmuyOxzUDg7T5fVZtXthGTUQ0QkQJ/Fgp0mkhsdyS2GYjMdkdim4HgtZtTu4mIyHZMRkREZDsmo5qRa3cANonEdkdim4HIbHckthkIUrs5ZkRERLZjz4iIiGzHZERERLarbDkgqgYRSQKQDKAUQDyAUlXNtzcqa4lIBoy2JZhfc1R1lsc2jj0PIhILIFtVMz2ed2SbzXaNALALQFMYv+9Sj9cd024RSYVRTgcw2luiqrke24R9m83/x2Xe4vanfZacA1XlIwgP8xeS5/HcTADxdsdmYRszACS5fR8LoATAqEg5DwByAMyMhN+9+fvO8Wy/U9sNIwmlejkHjvr7BpAKYI9nW/1tn1XngJfpgicTxhuVuxwYpdqdIl5Vi1zfqGoZjPa5t9ux50FE4gHEeXnJcW02e4DT1K0HaH6aTnXbzGntzlSPT/dq9PrT3LdBmLZZROJFJAdGMtntYzN/2mfNObA7Kzv1AaOHEO/xXKxxyu2Pz4L2xQIoBBDr8Xw8jNXb451+HgCMMh+ePSPHtdl8Y8n28rz7J2RHtdv8+z7j0z3cegFOabPZDm89o0rbZ9U5YM8oCMxPkWd82lCj5+D6RB3WzLbEmw+vnHwezLGEGV6ed2qbRwFY4fmkmuNFDm13PoA8czwEwMne4Ezz305s80n+tM/Kc8BkFBxxwKlfiBdh/UfqoqpN1O0ynSkVxkBoKZx9HmJ9tMupbY4FUCoio0Qkw/XV7XXHtVtVswCUASgUkXGu9uqpCQyOa7MHf9pn2TngbLrgiLU7ABtlApho/tuR50FEMtRjxqAbx7XZ7dNtvNsbMUQkW0TizOcc124AUNVeIpIH4zJlEYABbi87ss1u/GmfZeeAPSOyjIiMArBbVSfbHUuwmJclfH0KdCpfbzjvIAwG6qvD/JvOgjFpIR5GLyncezwhickoiMw3rohg/gfNVNU0L6856TwMVz/un3BYm0s9vgIAzEu0se5vzk5qtznTLF9Vi8zfeTsY56DQYzvHtNkbf9pnxTlgMgoO13/a06b9uv3CSuE82Tj9EgbgsPNgDmQXVLKZo9oMnDYeUNG4gKPa7Ypb3W7oVdUy88NWqTmBxVFt9sKf9ll2DjhmFASqWiYipTjz8kYcTg3uO4b5CTLLcxDTgechGUCCiIxwey4JQLyIZANYoaqzHNZml5N31nt7zaG/6xIfr+XAnMDisDafxt/2WXUO2DMKnnwYf9DuksznHcO8pp6tpy8Jk+p26cYx50FVc1U1y/0BIA9Akfm9a1KDY9rsJgenlsUBcLKn6P6G46R2l8JY4sqbWBiTGQBntdkbf9pnzTmw+2Yrpz5g/MF6LpGRhzBaJsSPNmbAuP8kye2RitOXiHH0eYBxedLzplfHtdlsU4mXNmU4td0wErDnckDxcLv51ylthrEcUIaX5yttn1XngPWMgshtUckVMP6IXQOhYc+8JrzHx8ulqprgtq3jzoPZ88sCMBzGf8ZcGEm4yHzdyW0ugdFrmOnZJqe12+z5J8BYGBYweoLeFkoNuzab/4cnwIg5A0ZvMB9GYpnltl2l7bPiHDAZERGR7ThmREREtmMyIiIi2zEZERGR7ZiMiIjIdkxGRERkOyYjIiKyHZMRERHZjsmIiIhsx2REljPXplPzMcrL6/G+XgtyXNkiEvJ3eYtIjojsqewcmdVH82oytkCFQ4zuzL/NPPPc59gdTyThqt0UbFkwlsohP4jITBhr/I00n6po1ePeMNYCDGXhEONJaiz6mmZ+aAmbJOoE7BlRMOXCKK+QYXcgYSQDRjmOWeajyNeGqjpMVcX9ObP3V+OVSH0d11uMoc5cZw1wzsrbYYHJiIJpJoxP9hPsDiQcuL2ZV6msubnw5Th4lHoINruOG0SpMGs02R1IJGEyomDLBpDk9mmTKNSlgb2iGsdkREFlLrdfBiMp+eRtcoGIjDIHkmM9ttsjIkkiUmi+Xugq6Oc2+Fzi6/Kg+bN55n5KzBLSnttkeO7fMwbz3znmPjwrXVZ0XNc+R7m9lo1TlUVd21R4edP9nJljTa6SHjPNn9/jsX3AbTLP6Uy3CRV57pfjKjuuj9+rz/PgGY+XbT0L/GWYsaq3+KooFeZ4kXkOxpmPCv+GqXqYjKgmTASQasGbhEssjEuAEwFkwqifMtPtMczcbqaPn58Go3Balvm95xtshvmz+ea+Sj23MbcrhFHhMqeySzrmPgvNfaUBeAdAjtuMrYnm8zDjSnCvKeOHkW4/nwmjBk8vC9qU6bH/OLMdfh3Xkx/nwZ3r95yD03/Prn25vp9lHtPVrko/GFQQ38nxIvMDxnBVnWy2i2OfwWR3hUE+nPeA8clSYVbJhPHmoDArwMJ4U1EAo9x+Jtv4czxtP6PM7WI9t4NbBU6358a5PZdhPpdUyTHi3WMzn9sDt2qe5nOFrufcjpfn6xx4OSd73I/h4zy5Yjmj4qaPfZ7Wnop+3qo2wRgXOu0YlRzXM8ZKz4NHPBlenov1+B3HVhRzgH+742D0UGM9/j5jrTwOH2c+2DOioFPjE3YugFH+XM7yU4Hbv12Xt9yv87umRMdVEpurumUycPKTcSyAcW6XfhSnyqq7y4QfzMthsTA+4bsfOx/GJcxh3n7OKla2SU/N7gu4l1vF8+D+Oy3x8doi8zKaFeOSrl7ecHWr6KqqZcoJDUHFZEQ1xXW93ZIbXX28MVT1zaIUp95cXV97wbg04/447c3STGT+cO3T2/buxw6WarXJHJeZaY7Z+Co1H0gcfp+HihKA+ZoreWQDKJTq32CbCvOyn7cxKgoe3vRKNUJVS0VkFoxp3qE2Uykep94gT75RBpBsKuPaTzwAz/uG4nF6Ly8Yqtwm8809HsabvSvOQt8/4Vcclp0Hs1fVy+xxjwKQLSLj1BjnCYhb4pmoqmUiMg7GmFRCoPuiwLFnRDVpIozLNCP83N7nQLhVzDexVJgJ0rwM5XX2XzUuMbreZE+7BGYO5rsG6a2w29uTVW2TOUEgFcZNuLnmfrwlM6/H9SJo58G8jDYZRpLrXcXdeN5fVArzMq85q44TGIKIPSOqMapaJCJFMAaJPZUAxlpmMBJDKiy6pOfJnI6chVOf+MtgJEqXkTCmKbtmcsXCeAMthZ/jRO7MT9mZMGaNxcGYQeY69izz0321mccBgEzXV1V1XcYKuE1mb7YMRm8DMM5TlpftKjqu53aWnQdzSniW2R7XZb4k83tXMs2G0dPxuZKFG2/3F7mS72ljSGQ99oyopp3xZmaaAeNTbTaMT8hNUYU3fj+9A+M+kjwYn+p7uY9NqDGluheMN+w8GFPBSyuIvVLmG1kaTk1PzoTR47B68sJkGIk8G26XwqrRppEwegczcep34/dxPVl8HmbAGN/JdNvXZLekEQtjxl2yn/uLg9vkCvOclZofkPzt/VEViWrIL2JMRFQlIlICYJifPSOyEXtGRORoTEThgcmIiBzJnJxR5UurVLN4mY6IiGzHnhEREdmOyYiIiGzHZERERLZjMiIiItsxGRERke2YjIiIyHZMRkREZLv/B1ATi/H/4oO6AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.rc(\"text\", usetex=True)\n",
    "import scipy as sp\n",
    "import scipy.sparse\n",
    "import scipy.sparse.linalg as spla\n",
    "import scipy\n",
    "from scipy.sparse import csc_matrix\n",
    "n = 10\n",
    "ex = np.ones(n);\n",
    "A = sp.sparse.spdiags(np.vstack((-ex,  2*ex, -ex)), [-1, 0, 1], n, n, 'csr'); \n",
    "rhs = np.ones(n)\n",
    "ev1, vec = spla.eigsh(A, k=2, which='LA')\n",
    "ev2, vec = spla.eigsh(A, k=2, which='SA')\n",
    "lam_max = ev1[0]\n",
    "lam_min = ev2[0]\n",
    "\n",
    "tau_opt = 2.0/(lam_max + lam_min)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "plt.close(fig)\n",
    "\n",
    "niters = 100\n",
    "x = np.zeros(n)\n",
    "res_richardson = []\n",
    "for i in range(niters):\n",
    "    rr = A.dot(x) - rhs\n",
    "    x = x - tau_opt * rr\n",
    "    res_richardson.append(np.linalg.norm(rr))\n",
    "#Convergence of an ordinary Richardson (with optimal parameter)\n",
    "plt.semilogy(res_richardson)\n",
    "plt.xlabel(\"Number of iterations, $k$\", fontsize=20)\n",
    "plt.ylabel(\"Residual norm, $\\|Ax_k - b\\|_2$\", fontsize=20)\n",
    "plt.xticks(fontsize=20)\n",
    "plt.yticks(fontsize=20)\n",
    "print(\"Maximum eigenvalue = {}, minimum eigenvalue = {}\".format(lam_max, lam_min))\n",
    "cond_number = lam_max.real / lam_min.real\n",
    "print(\"Condition number = {}\".format(cond_number))\n",
    "print(np.array(res_richardson)[1:] / np.array(res_richardson)[:-1])\n",
    "print(\"Theoretical factor: {}\".format((cond_number - 1) / (cond_number + 1)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "- Thus, for **ill-conditioned** matrices the error of the simple iteration method decays very slowly.<br>\n",
    "- This is another reason why **condition number** is so important:\n",
    "\n",
    "    * Besides  the bound on the error in the solution, it also gives an estimate of the number of iterations for the iterative methods.\n",
    "\n",
    "- Main questions for the iterative method is how to make the matrix **better conditioned**.\n",
    "- The answer is <font color='red'> use preconditioners </font>. Preconditioners will be discussed in further lectures."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Consider non-hermitian matrix $A$\n",
    "\n",
    "Possible cases of Richardson iteration behaviour:\n",
    "- convergence\n",
    "- divergence\n",
    "- almost stable trajectory\n",
    "\n",
    "**Q:** how can we identify our case **before** running iterative method?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Spectrum of the matrix = [0.5+1.32287566j 0.5-1.32287566j]\n"
     ]
    }
   ],
   "source": [
    "# B = np.random.randn(2, 2)\n",
    "B = np.array([[1, 2], [-1, 0]])\n",
    "# B = np.array([[0, 1], [-1, 0]])\n",
    "x_true = np.zeros(2)\n",
    "f = B.dot(x_true)\n",
    "eigvals = np.linalg.eigvals(B)\n",
    "print(\"Spectrum of the matrix = {}\".format(eigvals))\n",
    "\n",
    "# Run Richardson iteration\n",
    "x = np.array([0, -1])\n",
    "tau = 1e-2\n",
    "conv_x = [x]\n",
    "r = B.dot(x) - f\n",
    "conv_r = [np.linalg.norm(r)]\n",
    "num_iter = 1000\n",
    "for i in range(num_iter):\n",
    "    x = x - tau * r\n",
    "    conv_x.append(x)\n",
    "    r = B.dot(x) - f\n",
    "    conv_r.append(np.linalg.norm(r))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0,0.5,'Residual norm')"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAETCAYAAAAcboCHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl4VOXd//H3nT0QshK2sIaAiCwhCXWpaxtstUUtRqK29rGtLF2sWwtqn19ta1sLtNqqVUG7PK1tZXFp7eKCW+tOwiKgsoVdlEAWIAESwv3745zIkCZkmCVnZvJ5XddcyZycOfM9M5PznXs31lpEREQCFed1ACIiEt2USEREJChKJCIiEhQlEhERCYoSiYiIBEWJREREgqJEIiIiQVEiERGRoCiRiIhIUBK8DqAr9O7d2w4dOtTrMEREokplZeUea21uZ/t1i0QydOhQKioqvA5DRCSqGGO2+rOfqrZERCQoSiQiIhIUJRIREQmKEomIiARFiURERIKiRCIiIkFRIhERkaAokZxAY9MRfvj0WmobmrwORUQkYsV0IjHGTDbGLKivrw/o8VXVDfzpzW3csHAlLUe1tr2ISHtiOpFYa5+21k7PyMgI6PFj8jL4wSWn8e/11fxq6foQRyciEhtiOpGEwlWfGERZ8UDufXEjL77/kdfhiIhEHCWSThhj+PFlYxjdP50bH1vJ9ppGr0MSEYkoSiR+SEmM56EvFWOBbz+2giMtR70OSUQkYiiR+GlwTg9+8oWxrNhWx70vbPA6HBGRiKFEchIuGT+AsuKB3P/SRt7eXON1OCIiEUGJ5CT94JLTGJTdgxsfW0F9Y7PX4YiIeE6J5CSlJSfwqysnsHv/YW5/ajXWanyJiHRvSiQBKByUyU2TRvKPd3bx1MqdXocjIuIpJZIAzTxvOBOHZvH9p9ayo1ZdgkWk+1IiCVB8nOHuqYVY4OZFqzSFioh0W0okQRiU3YMfXnIab2+uYcG/q7wOR0TEE0okQZpSlMfnxvbn7ufXsWZnYJNDiohEMyWSIBlj+MkXxpDdM4kbF67kUHOL1yGJiHQpJZIQyOyRxM+vGM/G3Qf42b/e9zocEZEupUQSIueMyOWrnxzG71/fwivrq70OR0SkyyiRhNCsz57CyL5pfGfxKmq0qqKIdBNKJCGUkhjPL8snUN/YzO1PaNS7iHQPSiQhNnpAOt/5zEieWfshSyp3eB2OiEjYKZGEwXVn53NGfjY/+Ntatu3VqHcRiW1Rm0iMMZnGmFJjTJkxJtPreHzFxRl+MbWQuDjDzYtWaiEsEYlpEZtIjDFFbe6XuYljlrvpNmvtUmApML3LA+xEXmYqP75sDBVba/n1S5u8DkdEJGwiMpEYY0qBxT73iwDcxFHn3s90t9UBw72IszOXFuYxZUIe9yxdzz/e2eV1OCIiYRGRicRNGL6TV5UDde7vVUApTkLJdKu1IvYr/0+njKV4SBY3LVpJxRatqigisSciE0k7MgHfq3AOMB8noZQCC7wIyh8pifE8/OUS8jJT+Z/fvs1bVXu9DklEJKQSvA4gUNbaKo4vtRzHGDMdt+1k8ODBXRVWu7J7JvHY9DO4+uE3+Z/fvc29V07gwtP6BX3cPQcOs/6j/Wze08C2vY3UH2zmwOEjAPRMSqBXSgKDsnswrHdPxg/MJKNHYtDPKSLSVrQkkjog2/09E+j0a721dgFuSaWkpMTzkYF901NYOONMvvK7ZUz/YyXXnT2MmyaNpGeyf2+BtZYtext5Y9NeKrbUULmtlq0+XYuTEuLISE0kLTkBAzQ2tVB/sJmD7iSSxsCoful8elQfPjeuP6P69cIYE45TFZFuJloSyUKgxP09H6enVtTpnZbM4plncuff3+WRVzfzt1Uf8LWzh/GFCXn0SU85bt+jRy1baxpZvrWW1zft5fVNe9hVf8g9ThLFQ7L40ulDGD0gnWG9e9IvPYW4uOMTg7WW6v2H2bj7ABVba3l90x4eeHkj97+0kTF56Xzt7GF8buwAkhKipYZTRCKRicRpPIwxZcDDwDRr7RJ323Scqqx8t7Tht5KSEltRURH6QINQubWWOf96n7fdBvi8zFT6Z6QQH2eoP9jM9ppGGpqc0kRWj0TOGt6bM4fncNbwHIb17hlwaWLPgcP8c/Uu/u/1LWyqbqB/Rgo3TxrJlKKBxMephCIixxhjKq21JZ3uF4mJJFSMMZOByQUFBdM2bNjgdTjt2rh7Py+8t5s1H+yjpuEwzUcsmT0SGZCZyuj+6YwdmMEpfXv9V2kjWEePWl7ZUM0vl25g1fY6RvXrxZ2XjWHi0OzOHywi3YISiY9ILJFECmst/1i9i7v++T4f1B/kK2cN47ufOYXUpHivQxMRj/mbSFQ53s0ZY/j8uAE8d9O5XHPGEH772mY+d+9/WPfhfq9DE5EooUQiAPRMTuBHl47hz9edzv7DR7js16/x15U7vQ5LRKKAEokc56yC3vzj+rMZm5fBDY+tZM4z72tdFRE5oZhOJMaYycaYBfX19V6HElX6pKfwp2mnc/Xpg3nw5U3cvGgVTUc0g7GItC+mE4m19mlr7fSMjAyvQ4k6ifFx/OSyMXznwpE8uWIn1/2hgkPu4EYREV8xnUgkOMYYvvWpEcy5fCz/2VDN9D9WKpmIyH9RIpFOlU8czJwp45RMRKRdSiTil6kTB32cTL79lxW0HFUDvIg4YjqRqLE9tKZOHMQdnx/Nc+9+xA/+tla9uUQEiPFEosb20Lv2k8OYcV4+f3xzKw+8HLHriYlIF4qW2X8lgsz+zCg+qj/EvGfXMSSnB58fN8DrkETEQzFdIpHwiIszzC0bT8mQLL67+B3e/WCf1yGJiIeUSCQgSQlxPPClItJTE5jxaAW1DU1ehyQiHonpRKLG9vDq0yuFh75UzEf1h7n+Lys40qLR7yLdUUwnEjW2h9+EwVn8+LIxvLpxD/e+EJlrvohIeAXV2G6MKeTYWurHsda+GMyxJXpMnTiIt7fUcN9LGzljeA5nDe/tdUgi0oUCSiTGmMuBRa1329nFAloZqRv54SWnsXxbLTc+tpJ/3XAOOWnJXockIl0k0KqtOcBmYCpQ3M6t0xW1JLb0TE7gvqsmUNfYzC2LV3FUI99Fuo1Aq7bygTJr7ROhDEai22kDMvje507ljr+t5bevbea6c/K9DklEukCgJZKldNA2It3bl88cQumpfZj77DrWf6TlekW6g0ATyQxghjHmFmNMoTFmaNtb6EIMnLr/dj1jDHdNGUdacgI3L1pJs7oEi8S8QBNJPjAcmAtUApvauXlO3X+9kdsrmZ9cNoY1O/dx34sbvQ5HRMIs0DaSOe7PeURI0pDIctHY/kyZkMevX9rIp0f1YfygTK9DEpEwCTSRFAHTrbWPhDIYiS13XHIab1Tt5eZFK/nHt88hJVE9wkViUaBVW8uBmlAGIrEnIzWRuWXj2FTdwNxn1nkdjoiESaCJZBowxxhzgTEmPZQBSWw5Z0QuXz5zCL99bTOvb9rjdTgiEgaBJpIXcBrblwK1xpiWNrcjoQtRot1tF51Kfu+efGfRKvYdavY6HBEJsUDbSB7GmQZFpFOpSfHcXV7I5Q++zg/+upa7ywu9DklEQijQRPITa61WMxK/FQ7K5JsXFHDvCxsoHd2Xi8f29zokEQmRQKu2Ko0xt4Q0kjDQgMTIcv2nChg3MIPbn1zN7n2HvA5HREIk0ETyInBlKAMJBw1IjCyJ8XHcPbWQg00tzHr8HaxV7ahILAgokVhrZwDGGPOMMWZIiGOSGFbQJ43bLz6Vl9dV86e3tnkdjoiEQECJxBhTgzNNyoVAlXptycm45owhnDOiN3f+/V3W7FS1o0i0U68t6XJxcYZflhfy+fte5et/quTpb51NZo8kr8MSkQAFWrU121p7K/AckAFMcH8+Z6291f2bSIdy0pJ54ItFfFh/iJsWrtRCWCJRLNDGdowxD+IMSCwHctyfzxtjHgtRbBLjJgzO4vuTT+OlddX86O/vqvFdJEoFumb7NJw1SSZZa1/w2V4GLDLGLLPW/iJEMUoM+9Lpg9m6p4FHXt1Mbq9kvnlBgdchichJCmZhq1m+SQTAWrsEuBWYGWxg0j0YY7j94lO5rHAA855dx6+WbghbyaTlqKXpyFEONbfQoqo0kZAJZhr5WR38bRNOjy4Rv8TFGX5+xXji4+K4Z+l6dtQ28qNLx5CaFNi089ZaNlU3ULGlhmVbann/w318tO8wexsO05qj4gxk90yiX0YKo/qlM2ZAOmeP6M3w3DSMMSE8O5HYF2giWQqU4QxMbKscZ5p5zxljJgOTCwpUXRLpEuLjmFc2jrzMFO57aSMrttfx/c+P5tyRuZ0+tunIUdZ8UP9x4qjYUkNtozM5ZE7PJMbkZTBuYAa5ackku2uiHGpuYW9DE9trGnl5XTVLKncAkJeZyqWFA7i8eCDDc9PCd8IiMcQEUo1gjCkCKoCHgAVAFU4pZAYwHaftpL0k44mSkhJbUVHhdRjip/9sqOb2J1ezveYgY/MymDy+P+MHZtIvI4U4Y6hpaGJrTSOrd9Sxakc9q7bXcfiIszb8sN49KRmSxcSh2ZQMzWJY755+lTC21zTy7w3VLH33I/69YQ8tRy3nn5LLty4ooGRodrhPWSQiGWMqrbUlne4XaH2027C+AKfbL0Drf+sMa+3DAR00TJRIos/hIy0sXLadRRXbWbOz/flBkxLiGN0/neIhWUwcmkXxkGxyeyUH/dy79x9i0bLt/Pa1LdQ0NHF2QW/umDyaEX17BX1skWgS9kTi80SlwDCcUkmFtTbihiorkUS3j/Yd4t1d+6g50ETLUUtOWhL9M1IZ0TeNxPiAe7B36mBTC396ayv3vrCBhqYWrj1rKLdcOJIeSYHWCItEly5LJNFAiUSCUdPQxLxn1/HYsm0MzenJPeWFFA7K9DoskbDzN5EE9XXOGJPe0S2Y44pEkuyeSdw1ZSx/vu4MDje3cPmDr/PAyxs1gFLEFeikjROMMRuA2hPcRGLKmcNz+NeN5/LZ0/ox95l1XP+XFRxsavE6LBHPBTNpYw4wD2fciEi3kJGayP1XT+C0V9KZ9+w6Nu9p4HfXTqRPeorXoYl4JpgBidOttY+EMhiRaGCM4RvnFzCqXy++9ecVlD30Bo9+7XQG5/TwOjQRTwTaRrIUqAllICLR5lOj+vLodadTf7CZsodeZ/1H+70OScQTgSaS2cAcY8wX1LAu3VnR4CwWzzwTgKsffotN1Qc8jkik6wWaSF4AhgNLgFqtkCjd2ci+vfjztDOw1vLFh99ie02j1yGJdCmtkCgSAgV90nj0utO5csGbXPXwmzz+9bPoqwZ46SY0IFEkhN7ZUcdVC95kaO+eLJpxJj2TNQpeoleXDEgUkeONG5jJ/VcX8d6ufXz7Lyu07ol0C0okIiF2wag+/PDSMbzw/m7u/Pu7XocjEnYxXe7WeiTilWvOGPLxEsKjB6QztWSQ1yGJhE1Ml0istU9ba6dnZGR0vrNIiN160Sg+WZDD/z61hjU7I25SbJGQielEIuKlhPg47r1yAr17JjHz0UpqG5q8DkkkLJRIRMIoJy2ZB75UzO59h7lh4UqOqvFdYpASiUiYFQ7K5I5LRvPv9dU8/J8qr8MRCTklEpEucPUnBnPRmH7Me3Yd7+yo8zockZDqNJEYY462MwVKZzdNkSLiwxjDXVPGktsrmRseW0nDYf2LSOzwp/vvPDQdikjQMnskcU95IVc9/CY/+Nta5l0x3uuQREKi00RirZ3dFYGIdAdn5OfwzfMLuP+ljZx3Si6fHzfA65BEgqY2EpEudkPpCAoHZfK9J9ewq/6g1+GIBC2oke3GmClAfnt/s9b+PJhji8SqxPg47ikv5OJf/YfvLF7FH796OnFxxuuwRAIWUCIxxgwDKoAsd5MFfP8TLKBEItKBYb178v3Jo7ntidX89rXNXHdOu9/HRKJCoFVb84HNQLG1Ng6ox1nHPRtnGd4ZoQlPJHZdOXEQpaf2Ze4z63j/w31ehyMSsEATSSkwy1q7wr1fBWCtrQPmoEQi0iljDHMuH0t6aiI3PraSQ80tXockEpBAE0kd4DsTYhXQuvhJJk7pREQ6kZOWzLyycbz/4X5+/uw6r8MRCUigiWQpcKXP/UXAbLfxfQ5OohERP1wwqg/XnDGER17dzGsb93gdjshJCzSRzMZpIwHAWrsE2AIswWknmRZ0ZCLdyO0Xn0p+bk9uWbSKukbNEizRJaBEYq3dbK29tc22ScBwa222tfaJkEQn0k2kJsXzq/IJ7DlwmO89uQZrNZmERI+QDki01m7ufC8Rac/YgRncNGkk/1i9iydX7PQ6HBG/BZRIjDHpnd1CHahIdzDzvOFMHJrF9/+6lu01jV6HI+KXYHpt1XZyCztjTJExprQrnkukK8THGe6eWgjALYtW0aKFsCQKBDpFykzanxG4BKehfbq/BzLGFFlrl/vcL8NJVEXW2rmdPLwcWObvc4lEg0HZPfjhJadxy+JVzP/3Jr5xfoHXIYmcUECJxFq7oIM/PWyM2QgM9+c4bmlifuv+xpgi9/hLjTH5bZNMOxbSwVxfItFsSlEeL76/m7ufW8+5I3IZk5fR+YNEPBKO2X+XAmX+7GitXYo7Kt5VzrExKFU4I+gxxpS1uWWGMmCRSGOM4SdfGEPvtGS+/qdKqvcf9jokkQ6FI5GUEngpIROo8bmfA844lTa31mRTCkxUYpFYlNkjiYeuKaZ6/2Gu+79lNDZpVUWJTIH22upw+V2cke0nqo4KGWvtXGvtbJ/EIhJTCgdlcu+VE3hnZz3T/lChZCIRKdDG9hMtv7sJZ8qUQNThjIwHp3SyN8DjYIyZjtvoP3jw4EAPI+K5C0/rx8/LxvPdJav48m/e5jf/M5GMHokhOfah5hY27j5A9f7DVB84THPLUQCSE+LJ7ZXMgIwUhvXuSUK81sCTjgXa2B6u5XcXcmzyx3yc9paAuB0CFgCUlJSoD6VEtcuLB5KaFM8Nj63gc/f9h/uvLqJw0MnX6NY0NFGxpYbKrbVUbK1l9Y56mtzk0ZGkhDhO7Z/OOQW9uWBULhMGZWkhLjmO8XIqBrer78PANHe+rtaSRBWQf4LeYSelpKTEVlRUhOJQIp5avq2W6/+8gg/qDzK1eBAzzssnPzet3X2ttWze00DF1loqt9SybGsNVdUNACTGG8bmZTBxaDbjB2XSLyOF3LRkkhOcksfB5haq9x9me20j736wjxXb6lixvY6Wo5aBWalcUTyIq04fRJ9eKV127tL1jDGV1tqSTvfrLJEEOkrdWhsxK/UokUgsqT/YzH0vbOD/3thCc4tldP90xg3MILdXMgbY29DE5j0NrN5Zz/5DTptKRmoixUOyKBmaRcmQbMYNzCAlMf7knrexmZfW7WZx5XZe27iX5IQ4rvrEYGaeN5x+GUoosSiUieQoHbeHdMhae3Kf0jAwxkwGJhcUFEzbsGGD1+GIhNRH+w7x1IqdvLK+mvd27aPuYDMAWT2SGJiVypi8DMbmZVA8JIuC3LSQVkdt3tPAgy9v5InlO0mIN1z/qRFcd84wkhM8/7eXEAplIpnOfyeSTJzeWbM5fu2RAuC7OKsnRsya7SqRSHdgreWodaZZ6Srbaxr56T/f419rPmRY757MLRvHxKHZnT9QokLIEkkHB38QqLTWPtLO374LZFtrbzvpA4eJEolIeP17fTX/+9QadtQ28q0LCrj+0yNIVE+vqOdvIgn0nS7n+IGDvpbj58h2EYkN547M5Z83nMMXJgzk3hc3ctWCNzUavxsJ5itDRxMzXsGxsSCeMsZMNsYsqK+v9zoUkZiXlpzAL6aO51dXFrLmg3ouvf9V1uzU/153EGgi+RlwoTHmGWPMF4wxhcaYTxljFuIkmJB02w2WtfZpa+30jAxNeCfSVS4tzGPJzLOwwBUPvcFL7+/2OiQJs0CX2p0L3Ap8AngcqMQZPDgJmB1J7SMi0vXG5GXw1299kuF9ejLtDxU8veoDr0OSMAq4asud5yobZyT6VKDYXa99XsiiE5Go1adXCn+edgZFQ7L49mMrWLhsm9chSZgE3a3CWrvcWvu4tXZFKAISkdiRnpLIH776Cc4dkcutT6zmKa1FH5M6TSTGmOeMMUfabOtw9l/3FhFTlKqxXcR7KYnxzL+mmDOG5XDL4lU8s+ZDr0OSEPNnQOIsoMRaO9Vn2xw6Ge1urb01JBGGgMaRiHiv4fARrvnNW6zeWc/vrv0EZ4/o7XVI0omwDkiMNkokIpGh/mAz5fPfYGfdQZ74+lmM6NvL65DkBMI6INEYM9QYU9hm23eMMQuNMV8L5JgiEvsyUhP5zbUTSUmM5yu/X8aeAxq0GAsCbWyfD8xovWOMeQ5n7q3hwAJjzE9DEJuIxKC8zFQe+XIJew4cZtofKjjU3OJ1SBKkQBPJJOB5AGPMBJy102e7RaDbcEa3i4i0a/ygTH5ZXsiKbXX88Ol3vQ5HghRM998q92cpTsN762j2TTirG3pOvbZEItdnx/TnG+cP5y9vb2NxxXavw5EgBJpIqoCp7qJX5cByn4Ws8jmWZDylKVJEItstF57CJwty+N+n1rD2A33hi1aBJpJb3VstUISzLkmrGcCSIOMSkW4gPs5w75UTyO6ZxMxHK6l3F+eS6BLoXFtLcBrWpwJZ1toXff48B2dSRxGRTuWkJfPrLxaxq+4Q33tyNd1hSEKsCaaNxOLMs7WwdV13Y0wGTglFnwQR8VvR4CxumjSSv7+zi8eXaxqVaBPoOJJP4zSqT3Jv2QDW2nqckopKJCJyUmaeN5zTh2Vzx1/XsGVPg9fhyEkItEQyB1jidvdtu0D0ApwGeBERv8XHGe4pLyQhPo4bHltBc8tRr0MSPwWaSIqAhR38rRbIDPC4IaXuvyLRZUBmKj+bMpZVO+r55dL1Xocjfgo0kSzHaR9pzwycRa48p+6/ItHnorH9mVoykAdf3kTFlhqvwxE/BLPU7q3GmLtwGtYz3OV2nwUu5/juwCIiJ+X7k09jYFYPblq0kv2H1CU40gXT/XcqTunD4JRQlgMTganW2pUhi1BEup205ATuKR/PztqD/EhTqES8YJbaXeIutVvM8UvtPh6y6ESk2yoeks03zi9gceUOnlmzy+tw5ARCsdTuirZL7badYl5EJBA3lI5gbF4Gtz2xmt37DnkdjnTAr0RijEl31xrZa4zZ0NGaI+5+i4DKkEYpIt1SYnwc95QXcrC5hVmPv6NR7xHK3xLJcpyBhy/gtIksMMbc4ruD2/BeizMbcMQssysi0a2gTxq3X3wqL6+r5tE3t3odjrSj00RijPmu++swa+1Ua20BsBKYa4zpZYy5zhizF6fh/Va3nWReGGP2m8aRiMSGa84Ywnkjc/nJP99jU/UBr8ORNvwpkZQDP3OnP2k1G6dkUoczkv1hnEQTEQmklcaRiMQGYwzzysaRmhjPTQtXatR7hPEnkRQBFW22LXN/LsCZ/ffWNolGRCSk+qSncNeUsbyzo557X9jgdTjiw982kjrfOz5JY74SiIh0lc+O6U9Z8UB+/dJGKrdq1HukCLr7r4hIV7pj8mgGZKZy08JVHDh8xOtwBEjwc7/LjTFtSx4WZ7nddufcstY+ElRkIiLt6JWSyD3lhZTPf4M7n36XOWXjvA6p2/M3kXTUiN5RN18LKJGISFhMHJrNzPOG88DLmzg9P5spRQO9Dqlb8yeRFIc9ChGRk3Rj6UhWbKtj9uPvkJeZyun5OV6H1G2Z7jBStKSkxFZUtO14JiLRrr6xmSkPvsaeA00smXkmI/r28jqkmGKMqXQXMDwhNbaLSNTK6JHI7679BInxcVy54E3e27XP65C6JZVIRCTqVVUf4IuPvEVjUwu/vrqIs0f0Dup4DYePsGp7HRt2H2Dr3ka21zZy4NARGptbAEhPSSA9NZFhOT0p6JNG0eAsBmWnYkzblcejm78lkphOJMaYycDkgoKCaRs2aACTSCzbXtPIV3+/jI3VB5h53nC+dUEBPZP960+098Bh3qyq4Y2qPVRurWPdh/s46l4aUxPjGZzdg4zURFKS4rHWsv/QEeoam9hee5AWd8e8zFRKT+3DJYV5FA3OjImkokTiQyUSke6hsekId/x1LYsrd9CnVzJfPH0IF4/tR35uGvFxzoXdWsuO2oOs/aCeZVtqeX3T3o+rxHomxVM0JIsJg7MoGpzJ6AHp5KYld5gUmo4cpWrPAd7eXMOrG/bwyvpqDh85SkGfNKadM4zLJuSRnBDfZecfakokPpRIRLqXyq01/HLpBl7duAdrISUxjuweSRy1UNPQRJM7V1dSQhwlQ7L4ZEFvzhyew7i8DBLiA2863n+omWfWfMjvXtvCu7v20adXMrdcOJKy4kEfJ7JookTiQ4lEpHvaWXeQN9wSR/3BZgyQ3TOJwTk9OG1ABqP69SIlMfQlBmstr23cy93Pr2P5tjpG90/nzsvGUDwkK+TPFU5KJD6USETEC9Zann5nFz/753t8uO8Q088dzk2TRkRNdZe6/4qIeMwYwyXjB/DczedRPnEQD72yiUvvf43Nexq8Di2klEhERMIsLTmBu6aM43fXTuSjfYe45L5XWfruR16HFTJKJCIiXeSCUX14+vqzGdK7B9f9oYIHX94UE+vQK5GIiHShgVk9WDLzLCaPH8CcZ97nR39/l6NHozuZ+Dv7r4iIhEhKYjy/Ki8kNy2Z3762mZqGJn5xxfiguh57SYlERMQDcXGG//f5U+ndK4m5z6wD4O6phVE53kSJRETEI8YYvnF+AQBzn1lHQlwc88rGERdlyUSJRETEY984v4AjLZa7n19PSmIcP75sTFTN1aVEIiISAb796RE0NrXw0CubyMtK/bikEg2USEREIsSsz5zCrvqDzH1mHXmZqVxamOd1SH6Jzi4CfjLGTDbGLKivr/c6FBGRTsXFGeaWjeOM/Gy+s3gVb1Xt9Tokv8R0IrHWPm2tnZ6RkeF1KCIifklOiGf+NSUMyu7BN/+8nF31B70OqVMxnUhERKJRRmoiC64p5mBTCzMfXc7hIy1eh3RCSiQiIhGooE8vfjG1kFXb67jjr2u9Dudp5UmNAAALLklEQVSElEhERCLUZ8f045sXDOexZdtZtGy71+F0SIlERCSC3TzpFM4ansMdf1vLxt37vQ6nXUokIiIRLD7OcE95IalJ8Vz/l5Ucao689hIlEhGRCNc3PYWfXzGO93bt42f/et/rcP6LEomISBT41Ki+fOWTQ/n961siblEsJRIRkShx60WjGN0/ndmPv8OeA4e9DudjSiQiIlEiOSGee8oL2X/oCN97cnXErK6oRCIiEkVO6deLWy4cybNrP+LJFTu9DgdQIhERiTrXnZPPxKFZ3PG3tXxQ5/0UKkokIiJRJj7O8PMrxtNy1DL78Xc8r+JSIhERiUJDcnpy+8Wn8p8Ne3j0rW2exqJEIiISpb54+mDOHZnLT//xHlv2NHgWhxKJiEiUMsYw9/JxJMYbbl60kiMtRz2JQ4lERCSK9ctI4c7LxrB8Wx0PvbLJkxiUSEREotylhXlMHj+AXy7dwOodXb8irBKJiEgMuPPS0+idlsyNC1dwsKlrJ3aM2kRijMk3xhQZY2YZY/K9jkdExEuZPZKYd8U4NlU3MOeZrp3Y0fNEYowpanO/zBhTaoyZ1clDi4AqYClQFq74RESixTkjcrn2LGdix1fWV3fZ83qaSIwxpcBin/tFANbapUBd2yTjy1q7xFpbB5QCS8Idq4hINLj1olGM7JvGTQtXsqu+a0a9e5pI3IRR5bOpHKhzf6/CSRKtpRTfW6a7vRRYaq31PYaISLeVkhjPA18s5nBzC9/403KajoS/S7DnVVttZAI1Pvdz4OPSh++tzk0is4EZxhhVbYmIuAr6pDG3bDwrttXx03++F/bnSwj7M4SJW5pZ6nUcIiKR6HPj+vPOznz69koJ+3NFWiKpA7Ld3zOBvYEeyBgzHZgOMHjw4OAjExGJMrdddGqXPE+kVW0tBFq78uYTRInDWrvAWltirS3Jzc0NSXAiIvLfvO61VQaUtLZxWGuXu9tLgbrW+yIiErk8rdqy1i6hTddda+0Cj8IREZEARFrVVkgZYyYbYxbU13f93DMiIt1FTCcSa+3T1trpGRkZXociIhKzYjqRiIhI+CmRiIhIUGI6kaiNREQk/Iy11usYws4YUw1sDeIQvYE9IQonGnS38wWdc3ehcz45Q6y1nQ7E6xaJJFjGmAprbYnXcXSV7na+oHPuLnTO4RHTVVsiIhJ+SiQiIhIUJRL/dLfR9t3tfEHn3F3onMNAbSQiIhIUlUikWzLGzPL5vcwYU9rZNpFI1HZJcn8/z6H8jCuRnECsX0yMMdPd2xyfbTF/UXVnl57k/l4EHy+UVmeMKWpvm2fBhoB7TmW+K4nG+vvscy7T29kWM+fsfpYX+9z36/Mc6s+4EkkHYu1i0pbPevcLgHz3nynmL6rtKMdZUA2gCijtYFs0u82daTu/qy4sXnJjr3LPpSqWz7n1HH02+ft5DulnXImkY7F2MWkrn2PnVOXej/mLqjGmyP3na5UJ1Pjcz+lgW1RySyHLAKy1c901fmL+fQZaS9n53eicwf/Pc0g/40okHYuZi0l73BUkW3tzFAEVxPhF1ZXd+S4xZSKQ434Db62+ien32U0cVcaYWo6dU0yfs9eUSLo5tzi/vDusRtlOaQScb6StySUT2NvBtmi212f10bLOdo52xphMnPfwLuBhY0x+Jw+JJf5+nkP6Gfd0hcQIF2sXk46UWmtnu793dM6x8jrkuxeVbCDbTaILgdbpI/KB1kTT3rZotJdjdeh1OCWUWH+fpwN3WWvrjDFVQBmxf86tTubzHLLPuEokHVuI8wJD9F9M2mWMmW6tnev+Xkr75xwzr4O1donb6AzOhQOfb+qlQJ21dnl72zwJODSWcOz9y8RpL4np99mX+37XEaPn7JYwS1pLmv5+nkP9GdeAxBNwuw5W4TTYxdSIWJ9ugzU438qusNYube+cY/l16A7c968GmNha+oz199ltD6oCsk90frF0zl5SIhERkaCoaktERIKiRCIiIkFRIhERkaAokYiISFCUSEREJChKJCIiEhQlEjkhd1Zg696mt/P3/I7+Fua45hhjIr7vujFmvjGmtrPXyBgzyxjzfFfGdrIiMUb38/e8+/rO9zqe7kpTpMjJmE33XKo0IMaYxTgTYk5zN1WdYPeJRP7ssxEXo7W2CpjkfqmIqCTXnahEIv5qXbck5if9C6EyYHbr1CwnmobCWnuFtdb4bnNLXV0+4WBHz9tejJHAZx2RqJviJFYokYi/FuN8o77N60Cigc+FuO6EO3b8+ExgFk6Jpst49bxBKsVZyCqg11qCp0QiJ2MOUGSicCU5iWmTUGnEU0ok4jd3Urs6jq0+1672GsLdteGt+43Xd79ad9GlSvfvlW4Dv28j6qaOqtTcxz7vHmeTOxll233K2h6/bQzu7/PdY2S2PcYJnrf1mL5rg88BNrl3W/c5YZWg72vmtq3Uun9a7D6+ts3+J31O7mu62Kfx/3nfKqzOnreD97XD16FtPO3sG6ovJKW47SPuec5ybyf8nEroKJHIyboLKA1h3X0mTrXZXcAMnOm8F/vcrnD3W9zB4x8G5uN0BADnwu17cSxzH7vUPVZV233c/Spx1meY31kViXvMSvdYk3CmI5/v02voLnc7blzDfaav98c0n8fPAIYDxSE4pxltjp/tnodfz9uWH6+Dr9b3eT7Hv89B8W0fcb8ATHWXRhiO00YlXcFaq5tuHd5wvu1ZnAWwwLkgWJyLEzgXBAtM93nMHOejddxxprv7Zbbdr/XYbbbN8tlW5m4r6uQ58n1jc7fVAnPa7FfZus3n+Z4/idek1vc5OnidWmMp8/OYx53PiR4fqnPCaQc57jk6ed62MXb6OrSJp6ydbZknitGP120WTukvs81nMDPYY+vm/00lEjkp1vlmuwCY7k8VkJ8qfH5vrRLyrfNu7TZ7wvXWrdMVdCnuym/ut9VMYJY5NhbG4lxA21arzMAPbhVSJs43a9/nXopT7XdFe48LlVCekz3Wi+ykS5cBvg6+7+mmdv4eiNYS1FTrs56ItbbOqvG9y2gciQRiDk4JYzrOCnxB6eAfPtCLQBXHLyEKTvVM2+PVtInhRGM8fLUes739qwjgonySgjontzqq3D1OMLGe9OsQpgt7KTAXKHar8qbZ6F7RMiopkchJs9ZWGWOW4HQFjrTeMvkcu7h9fJE7iUTRmdbj5ANtL1j5HF+6CoeAz8kdlZ6P80WgNc7Kjh/hVxxevQ6+7SOt67PPwml3GR7u55bjqWpLAnUXTtVGuZ/7d9hoGypuVVspbnJzv5m228ssiGq51gvkcdVG7jf91gblUKhpb2Og5+Q2xJfiDJBc4B6nvUTU7vO2o6tehxNpO36kCrf60+29pcb2LqISiQTEWrvcGLMcp7GzrU3w8brZS3H+4cMyF5fbZXU2x75p1+EkuVbTcLqytvYYysS5+FXhZ7uIL/eb7wyc3knZOD2VWp97idtGEDT3eQBmtP601ra2B5z0ObmlyDpgjnu8Oo71dPP3edvuF9LXwU12c3BKGP5UT7U3fqQ1OR7XZiLhpRKJBOO/LkSuRTjVHXNwvpnmEMBF208LccYQPI/zbbrYty7eOt1ui3Euts/jdBeuOkHsnXIvUJM41oV1Bs43/VA3tM/FScJz8Kk+CuKcpuF8Y1/MsffG7+dtKwyvQyZOD72SznZ0ZePT2O++LlXuFxh/S1YSAsbaiJ9AVUS6CWPMJuAKNZhHF5VIRCSiKIlEHyUSEYkIbmeBgKscxTuq2hIRkaCoRCIiIkFRIhERkaAokYiISFCUSEREJChKJCIiEhQlEhERCYoSiYiIBOX/A/e6UxqWuVQuAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.semilogy(conv_r)\n",
    "plt.xlabel(\"Number of iteration, $k$\", fontsize=20)\n",
    "plt.ylabel(\"Residual norm\", fontsize=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'$x_0 = (0, -1)$')"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaUAAAEuCAYAAAAupQ21AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3V9sHed95vHnJ1tNlWy6R0ycAhGgdY4auIDhxKIZxEVvUpiMdxs0UBzKAta53JC9TLK1GTVGYi+cdSh3s7lbUN7L+EKmkiXaNZCEKjY3QRyIlJQIAWqkZGID7KJVQxFb2EKtOO9ezDvSaHj+zMyZOTPvnO8HICi+Z+acd0TyPHz/zPuac04AADTBgborAABAjFACADQGoQQAaAxCCWgJM+vWXYe8QqwzqkUoAS1gZtOSpuuuRwHTvu6AJMmYfQeEzcw6kpadc4up8mlJM5K2JXUlbTvnLtRQRZnZvKS9Xq9vZiuSlpxze+OvGZqGlhIQvtOSlpMFvlts2Tl31jl3wTl3VtJiHd1lZjYr6cUBhywrugaAUAJaYNo5t50qW5S0kipbUSq8qmRmXd8K6kra7XecrztdeJBE9x0QNN8tNuVbQsnyLUlzybDy3XzXnXM25mrG9Vns131oZguSlL4OTB5aSkDYTkl6OVngw2df6yQes2nojLeXJZ2suxKo3911VwBoEv+GPSupI+mYc27RzJ6StKeo5dG0N87pHhMEpqTbIdRDV9Hkh8Zwzu2Z2Uzd9UD9CCXA8y2M2bgLyczWzWxd0V/ws5LmzaxTdJaYma0qCrusFnuMFWWR5zWARiGUgNseT41pTEm64EPovJkdG2XactmtLB+ibZpGvT1K6KMdCCXgtpdTX09LWoq/SE0a6Eqal3TJH3e2hjfTgbPaynqDzzkGtTvCa+4quqZLBc9HCxBKgJd8M/X31mjAzaarzrmH/LEbiu7DGfd4U783/zg8p5LH+JZV8vGh4vudctTpoqQzOY5Po5U04QgloLc59fmL3a+UcKuF4gfpZ4c9YQVjSnHL4g6+Pts9XmtK0aoKmUPJHzuusB3Y8sNkIJQAz8y6iTfsWUkbicc6kmZ8y6mr/X/R75rZtHOub9dT2WNKPnym+jx8QdESQ8n6TPvypppiPAncpwToVnfdll+FIF5dIPkGuZDoyusXBHXo1+pZ0v4WzqISY2Rm1jGzrfjG1YpNaXgrsVHT1FEPWkpAZFvSWUUtJEl6RNJy4g07OStvV727xuqw0auF5ltRS2a2rGicJ14LL/3GP6Woq7L0lRR86/K0f+2Oov/POUnrzrnzqWNnJZ0ruw4ID8sMATn5ltSL8UQHX3bdOXe4prosplcIz/kcC3Uv7+PXyOsVmpgwdN8BOaVbJX6GWi1jNb4uTVw2KK8pAgkSLSWgEN9CmVW99ynFdZmV1C3S2vGB2q1rnyVfhwXVuNcTmoUxJaAA30KJW0y1vpk65y6Y2XLBm2Vn6+y68+NOHQIJMVpKQEs0YWwoLzObT096wGQjlAAAjcFEBwBAYxBKAIDGYKJDTu9///vdvffeW3c1ACAom5ub/+ycu2fYcYRSTvfee682NjaGHwgAuMXMXs9yHN13AIDGIJQAAI1BKAEAGoNQAgA0BqEEAGgMQgkA0BiEEgCgMbhPCZVbu7yjZ//m57r+1s19j5kkJ+mASb/tsQxjr8fjsrvM9I5zOtI5pCcfvU8njh+p7BoAjAcLsuY0MzPjJvXm2UHh0hSDwqvXYwQaMB5mtumcmxl6HKGUT9tDKYTgqcKgAOscOqhnPn0/wQWMgFCqSFtCaVLDZxTp4KL7EMiOUKpIiKG0dnlHL3z/Ne3s3eg7doPREVpAf4RSRUIIJVpBzcQEDUwyQqkiTQ2lUINolNl3bUFYYRIQShVpUiiFEkRVzXRLdksOmmHX67EQMEMQbUIoVaQJobR2eUdL3/mZ/vU3vx3r67btTbJXqDU5wGhRIWSEUkXqDKWn167qpVffqLzriinQt4UQXG37YwHtRChVpI5Qenrtqr796huVPDcBNJomhhYhhSYilCoyzlAqO4wIoPEbNu41DoQUmoBQqsg4QqnMMCKImqnOsCKkUAdCqSJVhtLa5R39xepP9ZsR3pkIobDVcaMzIYVxmJhQMrNpSTOStiV1JW075y5kOG/eH3/Mf15xzp0fdl5VoTRq6+hzDx/VcyceKLFGaII6WlSEFKowEaFkZnGYzCXKViUtOee2B5w3ryi8LvmvO5I2JS07584Oes0qQumJF3+sH23t5j7v4AHphZMP8qYxgcbZoqL1jTJMSigtS7qYbOGY2aykRefcyQHnPeWcO5MqW1AUcDboNcsOpblv/lC/+Kc3c51DGCFtHCFFCwqjmJRQ2pI0l2wV+VbP9X7h4h//W0mPOOf2EuVdSVuSjg1qZZUZSnlbSIQRsiKk0DStD6U4fCQdToaLf8xpQLiY2XVFoXQpUTbWUMozhnRA0jdPEUYoruqQIqAwzCSEUlfSVq8WkQ+luSwTHhLnLCgaUzo86LgyQmnt8o6+cO5KpmP/+NiUXvr8H430ekAaIYVxm4RQmpa0WWIobUo6lx5r8o8tSFqQpKNHjz70+uuvF6+4pO6XX1GWVeuYUYdxqTqkmCwBQilHKPnQOZmcxdfPqC2lrN12BBLqVNUK9LSgJlfWUDowjspUyY8tjXJ+V9FsvaGBVAYCCSE4cfyILn/1k/rWqQd1pHNIUrSn1ajiP4F39m7oC+eu6MFnf6C1yzujPzFaI+SWUjzR4Viv2Xfp8gHPsyrp8+nJEv2M0lLKMpbEGBKarNfNvGVsukj3Xvu1vvtOujUl/GSPWXSbwyYs+GNXFE1uGBpesVFC6cN/+YpuDhhMuvuA6e//658Wem6gLmV29dG9116T0n13QdESQ0nTvnygxGy7ZCtr1odaJQYFkiT91cmPVvXSQGXK7Oqjew+ht5Q6klZTywytKxoj2k4cc8cSQn6ZoSlJySbPlKJW1+Kg1xylpXTvl1/p+9i77j6g1577D4WeF2iiMltQdO+FL2tL6e5xVKYqzrk9M1uKlxtStLBqr+64qfgfcZD1ecrM3XhlW/7sR+p6aaASJ44f0YnjR0qZbr5346a+cO6KvnjuCt17LRd0S6kOVbWUfvWNTxWtEhAUWlCTaVLGlAAEJjkGdfjdB0d6rrgFxdhTexBKAGpR5gQJwqk96L7LqaruO26YBSJldO/Rrdc8E3GfUh2qCiWJcSUgqYxw4r6n5mBMKUB0OwC3ldG9x31P4aGllNMoLaX7v/o9vfn2O30fPyBpm9YSMBDde2GipdRAX//M4DGj3ypaRRxAf2XM3mNiRHMRSmOU5a+yrLvRApOOcGonQmnMPvfw0aHHfPzr62OoCdAOZYXTF89doaeiAQilMXvuxAO6e8ho7T/+y9ua++YPx1MhoCVGnRjhFPVU0GqqFxMdchp151kp275KEnsrAWUoOjGCyRDlYqJDg504fkR/fGxq6HE/2trVEy/+eAw1AtqraPce4031oKWUUxktpdjHv76uf/yXt4ce9/vv/R395Ctj2a0daD1aTvVgRYeKlBlKktQ9/UqmpfwPmPTNxx9szS/EEy/+WD/a2i10LiGNMjy9dlUvvfpG7q3cCadiCKWKlB1KWceXYm0YZxq23FJehBSKGuVGXMIpn4kJJTObVrQl+raiTf62nXNZtkMvdF7ZoSRFf7HluT/pwx94j9a/9IlS6zAuZQdSGgGFIoqGk0l6gsWUM5mIUDKzrqSV1Hboq5KWeuw+O/J5UjWhJOUPprsPmP7q5EeD+iut6kBKY+V15MV4U3UmJZSWJV10zp1PlM1KWnTOnSz7PKm6UJLyd+VJ0rvuPqDlz34kiF+GcYeSJN1t0t8/z3qCyIdwKt+khNKWpLlk68bMOpKuO+f63jZX9Dyp2lCKZZ2VlxTCWFMdoRQL4f8HzUM4laf19yn5EOlKumMKl3Nuzz/eLfO8cfrJV+Yy3ceU9KOtXd375VdYJqWPH23t6g9O1xeKCNMo9zixbFExwYaSpCnpdpj00C9ccp9nZgtmtmFmG9euXctd0SJe+vwf6VunHsz9Dfr2q2/o3i+/ouP/hRv+0n7j6m2tIVxxOH3u4aPKunIRyxYVE2z3nZ89t9mru83MnKLuuX2z6YqeFxtH913a3Dd/qF/805uFzm1St9WogZDcmbdIF2e/5wLyoEuvmNZ3302S9S99Ind3Xizu1rv/q9+r/a+1MoPgJ1+Z06++8Sn96huf0u+/93dyn89K7CiKLr1qBR9KfoxobOfVJe7OO1jwO/bm2+/oC+eu6EM1jztV0UKJA+rDH3hP5nNGaWUBUrFwoktvuJC77zqSrks61msWXbp81PNidXTfpa1d3tGTq1d087ejPc97fucuff0zD9TSnfCRr31P/+9f+28N30vWQMvaTfitU+1Ztgn1K7Js0SR16U3SlPCTzrlLibKuojGjw2WfJzUjlGJrl3f0pXNXNGI21X5X+oe+/MrQX+S8LawsoWeSfsnYEkpUZLyp7t+/cZmUUFpRFCRnE2Xzkk4NuXm20HlSs0Iplnc1iGEOv/ugvvZn4f/1lmUyBBMeUIUi4dT2VtOkhFJH0mpquaB1RSszbCeO2ZS0HIdQlvP6aWIoxcrq1ks6YNJ//HiYf8VlCWtCCVXK26XX5lbTRISSdGuK9ylJFxXdY3QpOaXbB9AvFa1rdzbref00OZRiVYRTLJSQyjplnFBC1Wg1RSYmlMYthFCKrV3e0TN//XPt3ci/LH8eTeruy3NP1++96y797Nl/X3GNgEjecDp08C49/1g9E5GqQChVJKRQSqqy9dTLOFtURW8uppWEOuTp0jv87oO6/NVPVl6ncSCUKhJqKMXG1XrKKuu09LInc/zuXaa/+/qflvZ8QB55Wk1t6cojlCoSeiglrV3e0env/kw3xtV8ahBaSWiCrK2mNkyAIJQq0qZQSpqkgCKQ0CST0moilCrS1lBKa2NI0WWHJnvw2R9k6lYPtdVEKFVkUkIpLfSQYmt0NF30O3ZVN25mW34rtFZT1lC6exyVQfhOHD9yxw9/0yZM9PLhD7xH61/6RN3VADKJf7+yduXt3bip09+9ese5bUBLKadJbSllVXeLqkn7RwFFtXHaON13FSGURjdqcBE8mAR5JkCE0D1NKFWEUAIwTllbTU0fY2LnWQBogedOPKD/nmEjwbbsbEsoAUDDxbvcdg4NDqZ4Z9uQg4lQAoBAPPPp+2UZjnvp1TeC3W6dUAKAQJw4fkRPPHx0aDA5RVPLQ0QoAUBAso4xXX/rph589gfBtZgIJQAITDzG9LmHjw48LsTJD8Gu6OB3jp2RtK1o59jtLDvHmtm8P/6Y/7zinDtfZV0BoArxvUmDtnVxicebfi+TFGgomVlX0rJzbi5Rtmpm28657QHnzSsKr/P+646kTTObSm6VDgCheO7EA/rfP/2/Q5f8CiWYQu2+W5S0kipbkbQ85Lyuc+5S/IVzbs+fk34uAAjGM5++X4cO3jX0uBBm5YUaSvOSLqXKNnx5T75VdMp/TrrgH++WWkMAGJMTx4/o+cceGDr5IYRZecGFkg+VrqTdZLlv9fQNF/9413/kfc0FM9sws41r167lrzQAVCzr5Ifrb91s9MSH4EJJ0pR0O4R66Bs6zrnDye47b1bS3qCxKOfcWefcjHNu5p577sldYQAYl+dOPDA0mJq86kOIoZTufhvVoqTnS35OAKhNyMEUYiiVxswWJO06587UXRcAKNNzJx4YulZeEyc+1BpKZtbN8dFJnTtSi8mPPS0mp5UDQJsMWyuviRMfartPKb7XKMcpFyWdUXSzrBSNLd0aV0qEVN+xoZRlSY/keH0ACMqJ40e08fruwJtr44kPTbl/qbZQ8hMLThY4b8/MtrV/bGlKQyYsxMxsRdLSgMkSANAKWVZ9aNKNtaGOKV1QtMRQ0rQvH8iPIy0nw8vMZrlPCUBbZZn40JTxpVBDaUn7W1mLvlxS1J1nZls+hOKy+ObajplN+49ZSSeztLAAIFTDJj40ZXwpyLXvfBfekpktKxpritfCSwfLVPwPP+a02ucpCSQArffMp+/XF89dkevz+PW3bmrt8o5OHD8y1nolBRlKkuRvgk3fCJt8fE/S4dTXWTZtBIBWyjLx4Sv/62qtoRRq9x0AoIBh40tvvv1OrTfVEkoAMGGGzbKrc9IDoQQAE2jYpIcXvv/a+CqTQCgBwAR65tP3D3x8Z+/GmGpyJ0IJACbQieNHht67VMfYEqEEABOqiWNLhULJzP6TmT1mZr9XdoUAAONzpHOo72N1jC0VbSnNSDov6bqZ/cLM/oeZfYaQAoCwPPnofQNv4Bz32FKhUHLO/blz7oCkj0k6K+mYpO/ozpD6aIn1BABU4MTxI3piyNjSOLvwRhpTcs5dcs694Jz7ZCKkLsefzeyimb23jIoCAKoxbGxpnGvilTrRwYfU45LO+ZBalXSJbj0AaLZBY0vX37o5tnoUnejwoO+i+wszu7fHIU6S/DbjpySdLlxDAEDlnnz0voGPj6sLr2hL6S8VLW7655K2fDfd835G3mOKuu8k3Vo4lYVQAaDBThw/okMH+0fCuLrwiobSup/s8AeSPqyom25O0v9U1Cp6XpJ8UN0r6dejVxUAUKXnH/tI38fG1YVXNJQu+MB50Dm37Zw745ybcc5NOec+5py74o97SFFgbZVT3dv8Bn0LftfYBb9ZX97n6Pit0QFg4tW5ZUWs0H5KzrlfSjptZo8MOe6ThWo1hN+6fNk5N5coWzWz7Zw7yC4rsREgAKC/cWwAOOqU8L8tqyI5LUpKt3BWFIVMJj7YCCQASBi0evg4xpVCXftuXvt3nd3w5VnNSlovrUYA0AKDVg8fx7hScKFkZh1JXUm7yXK/3XncAhr2HLOSXs7xmgtmtmFmG9euXctZYwAIR93jSsGFknyXWxxCPQwNJUmdAefv45w76ydyzNxzzz1ZTwMA5BRiKHVGOdnM5p1z58uqDAC0zaAbS6u+iTbEUCrMd/1lbiEBwCRyAx6rerJDoSnhZcky/pOwm+xyM7NcXXDe4865sznPAYCJcqRzqO+WFVVPdqgtlOJ7jXKcclHSGUnxfUhTSrR6fCtIicfTrzetaIYeAGCAJx+9T184d2X4gRWoLZT8Ta4nC5y3Z2bb2j+2NCVpb8DNszOSjpnZqUTZtKSumS1LushYEwBEM/AmLpRGdEFRyCTvVZr25T316rYzs6ckfcw5t1R6DQEAuYU60WFJ+1tZi75c0q117bbMbGHA87yvisoBAIoJsqXku/CW4m43RfcmLffouuu5jJAfz1qS9LikeFHWFb/NBgCgJkGGknRrn6a+IeJn5h3u89i2opbVYjW1A4BwDboXqerN8ULtvgMAVGTQvUiD7mEqA6EEALjDoHuRjnQOVfrahBIAILMnH72v0ucnlAAAmTV6kz8AQLtUveDqMIQSAOCWF77/Wq2vTygBAG7ptxCrNHir9LIQSgCAWwbdhzRoq/SyEEoAAEnReNKg+5DGsVU6oQQAkFT/eJJEKAEAvLrHkyRCCQAg6em1qwMfH8d4kkQoAcDEW7u8o2+/+sbAY8YxniQRSgAw8QYtwCpVv95dEqEEABNu0AKspurXu0sKdj8lM5tWtCX6tqJN/radc323Q+9x7ilJv1a0++xKjw0CAaD1ho0lPfHw0bF13UmBhpLfOXbZOTeXKFs1s+1h4WJm85LmnHOLibIVseEfgAnz9NrVgWNJhw4e0HMnHhhjjQINJUUBspIqW5G0LOlkv5PMrCPpRefc4UTZvKTZKioJAE21dnlHLw2Z3PD8Yx8ZU21uC3VMaV77t0Lf8OWDnJZ0NlngnDsvaa734QDQTs/+zc8Hrt7QOXRwrN12seBCybd2upJ2k+XOuT3/eHfA6QuSLqYLM3T5LZjZhpltXLt2LX+lAaBBnl67OnRyw7juS0oLsftuSrodQj10FU1+6KUjadvMFhSF2pSkXd9a6ss5d1a+hTUzM1P1FvUAUJlh40jS+Cc3JIUYSp0iJyVaUF0fMnH5splNJcsAoI2yBNLnHj469skNScF1342gX5idUzRBAgBaK8vEhs6hg7UGklRzS2nI+E/abrLLzsw6A7rwetlOfZYkOecumVnHzLrcqwSgrYZNbKhzHCmptlCK7zXKccpFSWd0O1SmJN0RUv6fPYPFObdnZkqekzJoLAoAgjVsYoNU7zhSUm2h5Fslfe8pGnDenplta3933JSkvSGtnVurP/R5DABaJYRxpKRQx5QuKFpiKGnalw+y4o+7xS85NCzMACA4oQWSFG4oLWl/K2vRl0uKuvPMbMtP/46d1f7lhJYlfb6SWgJATUKZ2JAW4pTwuAtvycyWFY01xWvhpVs7Uz3Om/Nr3W1JOubPy7SQKwCEYO3yjv7zyz8NYmJDWpChJEWz5rR/qaHk43uSDvco3xaLrwJoqafXruqlV98YGEhScyY2pAUbSgCAO2UZQ5KaN46URCgBQODWLu/omb/+ufZuDJ72bYpaSE0NJIlQAoCgrV3e0envXtWNm+8MPO4uM/23xz/ayC67JEIJAAIVT2h4xw0eQTIpiECSCCUACFLWCQ1Scyc19EIoAUBAso4fSWGMIaURSgAQgDxhJEmH331QX/uz+4NpIcUIJQBouDxddaFMaOiHUAKAhsrbOgppQkM/hBIANFDWqd6xePwo5ECSCCUAaJysU71joY4f9UIoAUBDTMpkhkEIJQCoWZGxo9CmemdFKAFATfKGkdTO1lFSsKHkd4ydUWKL8yz7IpnZrG7vPvs+SVvOubOVVRQAUoqEUehTvbMKMpTMLN7Uby5Rtmpm24O2NfdBJufcmUTZvJktEEwAqlYkjCTp0MG79PxjD7Q+kKRAQ0nRJn0rqbIVRVubp7dJv+M859wdG/w5586b2aqirdIBoFRrl3f0wvdf087eDZmU6QbYpLZ316WFGkrz2h9KG5LWh5w3Y2bdHq2pTmk1AwD1bhXlCaRJC6NYcKFkZh1FY0i7yXLn3J6ZqU/oxC5IWjezk347dZnZvKTVIa+5IGlBko4ePTrqJQBosaJddLFJDaNYcKEkaUqKQqjP411Fkx/2cc4t+YkOm2a2FB83bDzJP35WkmZmZvK2vgG03KhddBJhFAsxlEbqanPOPWRm64rGny5JeqSUWgGYOKN20UmEUVqIoTQS3xW3pKjFtaqo1TQ3aNYeAMRGbRXF5xzpHNKTj95HGKXUGkp+andWu8kuOzPrDOjC6/d6K4qmkm/7rz8kH0ySDud5LgCThVbReNQWSvG9RjlOuSjpjG6PF01JuiOk/D97tnjix5MtIh9qc2a2aWazWW6+BTA5yhgrkgijPGoLJR8Og+4p6nfenplta//Y0pSkvQHdcDOStvo8ttLj+QBMoH5BlCeQ6KIrLtQxpQuKQuZSomzal/czKAQ7Q84F0GJlBFGMVtFoQg2lJUVjQcmp3Iv+Q9Kt7rpNRWNIZ51z2/4+pju66Xw34vuY6ABMnjLGiWgVlSvIUPJdeEtmtqxorCleCy8dLFOp8xbNbMGf92tfvOecW6q+1gCaoKxxIolWURWCDCVJ8isyXBrw+J56zKhj4VVg8pTVPUerqHrBhhIApMXh8w97N/RvDx2UmXT9rZsEUUAIJQDBSofQm2//RjffiWKHcaIwEUoAgtKvK67oAqhpjBPVi1AC0HhlTtlOo1XULIQSgEZIdsV9sHNIf/KH9+j//N01gmjCEEoAatOvBbSzd0PffvWNW8cRRJODUAJQqV4z4vbeurlvYkIZwXPwgOnf/O7dt54/fq0PEkTBIJQAlCrrjLiyJibQAmoXQglAblnHf8oKnjSCqL0IJQBDDWr9VDH+0wtBNBkIJQB9Wz7DuuDKlgyeZB0YE5ochBIwAfKETrrlU2UIJScmEDyQCCWgFZoaOhIz4pAPoQQ0TDpgnnz0PklqbOj0wvgPijLnqhqWbKeZmRm3sbFRdzUQiEEtmF6Bkw4YKWppyHRHWd3S3W6M/2AYM9t0zs0MPS70UDKzeUUb9WXaztzMphVtpb6taHPA7aznSoRSG+VtmWQ95k/+8B59Z3NHN26+0/e1mxg4vTD2g1FNRCiZ2ayibdFPZgkWv/X5inNuLlG2Kmkp63bohFJ1soRD2WW9giNLUGQ5ZtRdTceJlg+q1upQ8uGyJGnTf17MGErLki46584nymb9+SezvDahVI21yzs6/d2rQ8Oh7LKQgmMUhA7qljWUgpzo4Fs1i5JkZks5Tp2XtJIq25C0XlLVUNAL339tXzfXzd/uj4uyy0IIpKwtN0IHbRBkKBVhZh1FY0i7yXLn3J6Zycy6/brwzGxB0oIkHT16tPK6TqJ/2LtRdxUqMawl1q9Flx6/kYaPcRE6aIOJCSVJU1IUQn0e7yqa/LCPc+6spLNS1H1XSe0m3Ac7h7RTUzClg6OsMaVDB+/SZx86sm/yw7AJE/0ChtDBJJikUOrUXQH09+Sj99UyptQrOMqafZen9ULgAJFJCiU0WPymPO7Zd4OCI0tQECZAuWqdfedn0WW126vrzcy2lGH2nb8/aVPS4fTzmJmTNJdlBh+z7wAgv8bPvvOBtJzjlIuSzozwkvF40ZSkW6HkJ0AkHwcA1KS2UPIz3TLdG1TS6+2Z2bb2jy1NKVoRglACgJodqLsCY3ZB0RJDSdO+HABQszaE0pR6zKwzs46Zbfl7jGJL2t86W/TlAICaBTn7zo8DnVZ0b1FH0rKZzUlaTy4hJH9vUsx34S3Fyw3585fpugOAZggylPzsuYGtG3/M4R7llyRdqqhqAIARtKH7DgDQEkGuEl4nM7sm6fW669HH+yX9c92VKEFbrkNqz7W05TokrqUu/845d8+wgwilFjGzjSw3pzVdW65Das+1tOU6JK6l6ei+AwA0BqEEAGgMQqldztZdgZK05Tqk9lxLW65D4loajTElAEBj0FICADQGoQQAaAxCCQDQGEEuM4Qw+Y0WZxTtXdWVtJ1lY8Wi51VphGuZ98cf859XUus1jl0Z/79+Pcpl59xiBVXMWofC1+HPPSXp15Lep+j7UtuamCP8fM2DeXLfAAAFh0lEQVQq2vlAiq5jyzkX1mQI5xwfgX4o+uFbkDQbf8543rykpyStSFqXND+GunYVLZibLFuV1K3ivIZey7yk6cTXHUlbkhZCu5Yez7MiaTXE6/Dfl5X09YR2Lf79YDZVNl/nz1eh66+7AnwU/MYF9saoaJfh+VTZ7LA3sqLnNfRanupRthD9bRjWtaSO7/qfvRC/Jx1J11Nl84paGKFdS88grfP7Uuj6664AHwW/cYG9Mfrg66bKOsNet+h5TbsW//impE6qvCvJZfmLvinX0ufnZ6HmUCr687WsqNsxXV5nS7zotWz2qrdSf7w2/YOJDuGa1/4tODZ8eU++3/+U/5x0wT/eLbWGd75uV9JustxF24v0fd2i51WpaJ38413/0Qhl/P/6MYyXK6lgRiNex4KivdXu4GoaTxrxWi5IWvfjUfHzzStqxQaDUApQgG+MU4nX76VffYqeV6XCdXLOHXbRfl5Js5L2anoTLOP/tzPg/HEZ5To6krbNbMHM5uPPpdcwu1F+vpYk7UnaNLOn4utwgU10YPZdmLL84PZ8k3PO7dv4UNW/Me7brr7i86pUdp0WJT1f8nNmNdK1mNm8q3nmoFfoOhJ/vHWTb9xmtmxmUzW9mY/0PXHOPWRm64q6JS9JeqSUWo0RLaUwtemNcWKZ2YKkXefcmbrrkpdvrdfdQhpVv9+jc4re1IPjf6aWJM0p+uN0s45u7lEQShNunG+MPcayKj2vSqPWyb9RLDrn5kqq0ih1KXItj7ua7xVLK3Ad26nPkiTfxdqp8828yPfEzFYkXXDOXfLfmw8purbNsutXJbrvapbzB3832WVnZiP15yfeGB8q+hwZxb/0U0r8dZ34xevXbVj0vCqVVadl1d+1Uuha/ED6RrVVy6XQdTjn9sxM6t/i69sNXqGi35OOdOcEDf/eMGdmm2Y227Q/IvohlGrkQyFPN8FFSWcU2Buj/+Xf1v7ukikNGMsqel6VyqiT/4t2qe4JAiNcy4ykY2Z2KlE2LalrZsuSLo5zrGnE78mtFRP6PDZWI35Ptvo8ttLj+RqLUKqR/wE7WeC8EN8YLyj6xUnOPpv25VWcV6XCdfLdpcvJ75GfVr1d0wy83NfSawKAmT0l6WN+Blgdin5PVtLH+ZZgXTMipWLXMui9pDPk3Gap+0YpPop9KPplWkiVzSvDDYyK7s1I35w3my4rub4d7V+BYj35muqxukSW82r4vy96LfP+/3468TGrepe0KXQtPZ5nOcvPXtOuIy7rcV7lS29VcC0r2r/MUFc9bg5u8geb/AXKd9WtusRAuZ8Kuuj8X3j+mE1FP5Rnfdm8ohZVckxgStJJV/FimolFLy8q+mWJB2ST1/RLRS24s1nPq0Pea/FfX+/zdNvOuWNV17mfot8X/1hX0WyvxxW9UZ5VFLLp+7EqN8LPV3wNW4oWyl0N7ecrUb6g6Bp+7Yv20t+zpiOUAtamN0YAkAglAECDcJ8SAKAxCCUAQGMQSgCAxiCUAACNQSgBABqDUAIANAahBABoDEIJANAYhBIAoDEIJQBAYxBKAIDGIJQAAI3BJn9AwPy2C7OKto045pxb9Bvu7Umac87l3kQSqBOrhAOB8luRPJ7YK2vdP3RSUVCtSjrsat52HciD7jsgXI+nNnCbUrSn1p5z7ryilhOBhKDQUgICZWadZOiYmVPUZVfrrqnAKAgloAXMbFbSunPO6q4LMAq674B2mJN0qe5KAKMilIBA+Zl3sVlJG4nHOr71BASFUAIC5ANny8y6Zjbti5OTGhYYW0KIGFMCAuRbSUuSNn3Ry5KWk18z8w4hIpQAAI1B9x0AoDEIJQBAYxBKAIDGIJQAAI1BKAEAGoNQAgA0BqEEAGgMQgkA0BiEEgCgMf4/XLiKOx7Gf9YAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter([x[0] for x in conv_x], [x[1] for x in conv_x])\n",
    "plt.xlabel(\"$x$\", fontsize=20)\n",
    "plt.ylabel(\"$y$\", fontsize=20)\n",
    "plt.xticks(fontsize=20)\n",
    "plt.yticks(fontsize=20)\n",
    "plt.title(\"$x_0 = (0, -1)$\", fontsize=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Better iterative methods\n",
    "\n",
    "But before preconditioners, we can use **better iterative methods**. \n",
    "\n",
    "There is a whole **zoo** of iterative methods, but we need to know just few of them."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Attempt 1: The steepest descent method\n",
    "\n",
    "- Suppose we **change** $\\tau$ every step, i.e. \n",
    "\n",
    "$$\n",
    "   x_{k+1} = x_k - \\tau_k (A x_k - f).\n",
    "$$\n",
    "\n",
    "- A possible choice of $\\tau_k$ is such that it minimizes norm of the current residual\n",
    "\n",
    "$$ \\tau_k = \\arg\\min_{\\tau} \\|A(x_k - \\tau_k (A x_k - f)) - f\\|_2^2.$$\n",
    "\n",
    "- This problem can be solved abnalytically (derive this solution!)\n",
    "\n",
    "$$ \\tau_k = \\frac{r_k^{\\top}r_k}{r_k^{\\top}Ar_k}, \\quad r_k = Ax_k - f $$\n",
    "\n",
    "- This method is called **the steepest descent**.\n",
    "\n",
    "- However, it still converges similarly to the Richardson iteration."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Attempt 2: Chebyshev iteration\n",
    "\n",
    "Another way to find $\\tau_k$ is to consider\n",
    "\n",
    "$$e_{k+1} = (I - \\tau_k A) e_k = (I - \\tau_k A) (I - \\tau_{k-1} A)  e_{k-1} = \\ldots = p(A) e_0, $$\n",
    "\n",
    "where $p(A)$ is a **matrix polynomial** (simplest matrix function)  \n",
    "\n",
    "$$\n",
    "   p(A) = (I - \\tau_k A) \\ldots (I - \\tau_0 A),\n",
    "$$\n",
    "\n",
    "and $p(0) = 1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Optimal choice of time steps\n",
    "The error is written as \n",
    "\n",
    "$$e_{k+1} = p(A) e_0, $$\n",
    "\n",
    "and hence\n",
    "\n",
    "$$\\|e_{k+1}\\| \\leq \\|p(A)\\| \\|e_0\\|, $$\n",
    "\n",
    "where $p(0) = 1$ and $p(A)$ is a **matrix polynomial**. \n",
    "\n",
    "To get better **error reduction**, we need to minimize\n",
    "\n",
    "$$\\Vert p(A) \\Vert$$ \n",
    "\n",
    "over all possible polynomials $p(x)$ of degree $k+1$ such that $p(0)=1$. We will use $\\|\\cdot\\|_2$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "##  Polynomials least deviating from zeros\n",
    "\n",
    "Important special case: $A = A^* > 0$.\n",
    "\n",
    "Then $A = U \\Lambda U^*$, \n",
    "\n",
    "and \n",
    "\n",
    "\n",
    "\n",
    "$$\\Vert p(A) \\Vert_2 = \\Vert U p(\\Lambda) U^* \\Vert_2 = \\Vert p(\\Lambda) \\Vert_2 = \\max_i |p(\\lambda_i)| \\overset{!}{\\leq}\n",
    "\\max_{\\lambda_\\min \\leq \\lambda {\\leq} \\lambda_\\max} |p(\\lambda)|.$$\n",
    "\n",
    "The latter inequality is the only approximation. Here we make a <font color='red'> crucial assumption </font> that we do not want to benefit from distribution of spectra between $\\lambda_\\min$ and $\\lambda_\\max$.\n",
    "\n",
    "Thus, we need to find a polynomial such that $p(0) = 1$, that has the least possible deviation from $0$ on $[\\lambda_\\min, \\lambda_\\max]$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "##  Polynomials least deviating from zeros (2)\n",
    "\n",
    "We can do the affine transformation of the interval $[\\lambda_\\min, \\lambda_\\max]$ to the interval $[-1, 1]$:\n",
    "\n",
    "$$\n",
    "\\xi = \\frac{{\\lambda_\\max + \\lambda_\\min - (\\lambda_\\min-\\lambda_\\max)x}}{2}, \\quad x\\in [-1, 1].\n",
    "$$\n",
    "\n",
    "The problem is then reduced to the problem of finding the **polynomial least deviating from zero** on an interval $[-1, 1]$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Exact solution: Chebyshev polynomials\n",
    "\n",
    "The exact solution to this problem is given by the famous **Chebyshev polynomials** of the form\n",
    "\n",
    "$$T_n(x) =  \\cos (n \\arccos x)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## What do you need to know about Chebyshev polynomials\n",
    "\n",
    "1. This is a polynomial! \n",
    "\n",
    "2. We can express $T_n$ from $T_{n-1}$ and $T_{n-2}$: \n",
    "\n",
    "$$T_n(x) = 2x T_{n-1}(x) - T_{n-2}(x), \\quad T_0(x)=1, \\quad T_1(x)=x$$\n",
    "\n",
    "3. $|T_n(x)| \\leq 1$ on $x \\in [-1, 1]$.\n",
    "\n",
    "4. It has $(n+1)$ **alternation points**, where the maximal absolute value is achieved (this is the sufficient and necessary condition for the **optimality**) (Chebyshev alternance theorem, no proof here).\n",
    "\n",
    "5. The **roots** are just \n",
    "\n",
    "$$n \\arccos x_k = \\frac{\\pi}{2} + \\pi k, \\quad \\rightarrow\\quad x_k = \\cos \\frac{\\pi(2k + 1)}{2n}, \\; k = 0, \\ldots,n-1$$\n",
    "\n",
    "We can plot them..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Interval $x\\\\in[-1.1, 1.1]$')"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEICAYAAAC55kg0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJztvXmUHNd93/u9vW8z07MCgx0DcCdFcTCiRIlaLINRZL3IfjZI2lIkO1YEKfJJ4vjF1HGW915iJwpoS8dSbFmkbCeKpfhRpB07tjYSlC1RqwWClCkBIMEZAsQ6e0/39F5V9/1Rdat7Zqq6q7qW7ur5fc7BwUz1VtN9u373t31/jHMOgiAIYnsS6vYJEARBEN2DjABBEMQ2howAQRDENoaMAEEQxDaGjABBEMQ2howAQRDENoaMAEEQxDaGjIABjLFVxli2i6+fZYzN+vyaxxljx31+zaOMsYf8fM1+hdZssOiltU9GoEOCuvhakOWcP2rnAYyxaTvvw+b7c85P2nk9whm0Zq2tWafr2sp9emntkxHonA91+wSsIhYgY+whty4EjLGjAH4DgKXdp937E55Aa7bNGvRiXff62o90+wR6Gc29fhrAKQBHATzBOf8oY+wRANOMsWcBnNSOPQT1SzYH4H7tKZ4GcBLAlPb7I2IHwBib5ZwfYow9rt2+AuB+znnO5b/hEQCP290xtYNzfpIxNgWLC9vu/YnOoDVrjpU16MW67vW1T55Ae6YBnOCcHwJwHAA45x8CMMc5P6J9maYBvE67z4cAnGh67DLn/H4Aj0H7omn3P6091/2c8yMAHoG6W2iLFk88yhg7of3/uMn9TgD4aC+5noQv0JolLEOeQHtynPM57ecVxljWYOfzIIApxthTBo99GAA4509oC1zc/zFAdxU/BHVnNYc2aDu9U5zzHGPshPaFNn2c27s0IhDQmiUsQ0agPSsW7/cxzvkT4hdt4W9+7GltR3W0aTf2Uai7rSlY2FWJL4jmXp7Sjpl9oaYMKhBy7dxsLQZr5ro+Sl/SnofW7EZozbaAjEDn5JpK8h6D6k4/Aeius9EifwzqDkp80WYAnNZ2SDNWXrTpNY8BeFa8Huf8tMHdfyB2dXZwOxZL9Ay0ZoktUE6gc05CXdAntMX8OGPsWS3xZvjl0HZdx6HGUgHgiwCOaY+5z+LrHgfwALQvrOaaG+6qOOcPa9UV0xaf2zLa694H4D7tZ3H8WaPXM7s/4Su0Zs3XbNbCfSyvayvP2SswGirT/2iLdwaqu2zoWjPGHupkB2b2eia7PKP7uva6RP/g95q1ek5W17XF5+uJtU/hoG2AtnBdW7ytYIwdpcoOwil+rlkr9PO6pnAQIci50ZRj54vSi64xEShcWbNWcNsA9NLap3AQQRDENoY8AYIgiG1Mz+cExsbG+IEDB7p9GkSf8uyzzy5xzsf9fl1a14SX2FnXPW8EDhw4gFOnTnX7NIg+hTF2sRuvS+ua8BI765rCQQRBENsYMgIEQRDbGDICBEEQ2xgyAgRBENsYV4xAK50PxtgxtmmeptExgiAIwn8cGwGt881sQMQ0oHfb5Zg6Mm7LMafnQBAEQXSGYyOgXczNtMEfBCB0vOegjrszOkYQBEF0Aa9zApuHVIyaHHPMYqGKtVLdjaeyxVqpjq/+6DoUJVjyG5Ks4AcXVnDuer7bp2KKonD8xXNXUKxKHT/H575zAf/7h1ddPCuC6D7X1yr4+JMvYnZx3fFz9WRimDF2nDF2ijF2anFx0dJjfvm//wD//i9/5PGZbeTM1Tze/vG/xYc//yy+eOqSr6/tlIf+7O9x/2e+i3d96ls4P1/o9ukY8rcvLeBXH3sev/Wlsx0/x+e/dxFfeeGai2dFEN3nSq6M//r1l/HqSsnxc3ltBHIARrSfswCWTY5tgHP+KOd8hnM+Mz7evvO5Ksk4cy2PuSXnVtEOn/vOBVQlBXfsHsLHn3rJ0Y7VT85ey+N/PXcF9x/Zg1Q0jBNfPdftUzLkz05fAQD86d+9iucvdTYdUOEcoRBz87QIoutU6zIAIBkNO34uT4zAphF2U9rPU1AnGxkdc8TsQhGywnEtV3H6VJaRZAVPnrmOt988gX/3rluwWKji5Nl5317fCY9+cw6ZWAT/9l234ENvncLJswuYc8GtdJO1ch1PnZnHz03vQYgBXz+30NHzKBwIMTICRH9RkVQjkOgFI8AYOwZgRvtf8DSgD4YQFUQ5zvlpo2NOz+ElLZyxXKyhollIr/n+KytYLdXxztt3YubACAbiEXxvzup87+6hKBzfeGkRR2/dgWwqhnffuRsA8Mz5pS6f2UZ+8MoKapKCB2b2YN9ICrMLnRkpWeEIkw0g+oxyTQHgjifgWEBOm0H6xKZjR5p+3jIWzu2h0OeuN2La8/kK9o+m3Xx6Q755fhHRMMNbbxpHOMQwc2AY339lS2Sr5zh7PY+VYg33Hh4DAOwbTWH/aArPnF/EL77xQHdPrgkR2rtp5wAOjWfwsgMjQOEgot8Qm91E1HkwpycTw3Z5sanC5apPIaFz1wo4PDGAVEy1o6+fGsXcYhGLhaovr98p39J2/G/SjAAA3Ht4DN+bW0FdVrp1WluYWyxiJB1DNhXD4YkMXlkqQurg/DjnCFM4iOgzyvUeCgf1AucX1nH77kEAwPV82ZfXPHc9j1t2Dui/v/6gmus+daG3Q0LfmV3GDRMZ7BxK6MfuPTyG9aqEH11Z6+KZbWRuqYipMdWjOzSRQU1WcGnV/mcrc045AaLvqJAR2Ei+XMdtk0MA/PEEVoo1zOeruGVyUD92y+QgGNsYmupFfnx1Da/dm91w7I496nt35lrv9AzMLRYxNa4agcMTGQDoKC8gK6BwENF3VCXVK6ZwkEalrmAkE8NQMopra957Aue0i+XNkw1PIBENY+9wCi/3WJVNMwuFCpbWaxuMFwDsziYxkIjgbI8YgXyljqX1KqbG1Yu/MAKdvLcK5wj3xSoniAblmowQA2IuLO7Afz1khaMmK0hEwpgcSuD6mveewFltt3/zzo0X0xsmMnh5vneNwNlr6nlvNgKMMdyyc1C/vdvMLRYBQA8HDSaiGEhEOvpsFcoJEH1IpS4jEQ2DubC2A28ERGwsGQthYjDhS2L2/HwBI+kYxgfiG447SWD6gdjp37rJCACqV/Pi9UJPyF9cXFaNwMGxRpXXeCbe0WcrK9yVLwpB9BJlzQi4QeCNQHOWfCARQcGHrt3Lq2XsG0ltOe4kgekHZ6/lsWsogaFUdMttt0wOYr0q4XIPnPt8Xt3xNyevxzJxLK7bNwKKwhGmnADRZ1Tqiis9AkAfGIHmLPlAPIL1ivdG4NJqCXuGk1uOi9h1r2rxvHi9gJsNvABArccH0BOCcvP5KlKxMDLxRhvL+EAcSx0YAZmTESD6j4okI+5CUhjoMyOQiUc81++RFY6ruTL2DG/1BPQqFi2m3UsoCscrS0UcGjdupBPx94vLzgWpnDKfr2DHYGJDGGcsE8NSB+Egko0g+pFKTSZPQFCpN9qn0/EIijUZsodx7YVCBXWZG3oCg4kohpJRXMl1/0K6mflCBVVJMe2mzqZiGExEcHGl+wZsIV/FxKZ8y1gmjnxFQlWyJwuiKBzkCBD9RkWinIBOual9eiChhg+KNe+8AREz32uQEwCAPcPJnoirb+bCkmqYDrSQ1Dgwlu4NT6CgegLNjGlGYXm9Zuu5KBxE9CPlmuxKjwDQD0ag1pBUFTFkL/MClzT9biNPAFBr7q/0oBEQFTcHxoyNFwDsG0l13Qhwzk09AQC28gKcc3AKBxF9CCWGm9iQE9A8gXUP8wJil787a2wE9gyncHm1DM67X2rZzCvLRcTCIUwOGZ83oHoJV3LlrmoIFaoSynV5qyeQiQGwZwREVJA8AaLfUBPDZAQAbCwRFZ5AwUNP4PJqCeMDcdN43O7hJMp1GatdGHXZiotLJewdSba8IO4fTUFWeFc9mQWtPHRi0NgTsNMrIHJDZAOIfoMSw01URWK4qaTQywqha2sV7BpKmN4uwkSXV7sfW2/mwnJxQ/OVESJpfGG5e8nh+bx6kd/sCYjGvCUbOQFF88ZIO4joNyqSQjkBge4JREK+hIMW8lVMDLY3Ar2UF+Cc49JKyTSZLdg/qt5+yYW5pZ2yUFA9gc1GQHh6nXgCJBtB9BvlmoxEhDwBAM2yEf4khhcKlS1Jy2b2ZNULaS9VCK2V6yjWZNM8hmAsE0ckxHDVB/0lMxY0T2CzJAcADKejyJU68ATICBB9BOccFUlGMtYjRoAxdowxdpQx9pDBbdOMMc4Ym9X+PaIdP6H9f9zp6zc8gTAG4qocglfSETVJwWqpjokBc09gMBlBJh7BlVzvGAFxLu2MQDjEsGMwgWtdPPeVUg2xSAhpgwU+lIwib8PAK1p+m8JBRD9RkxVw7s4sAcDheEnG2DQAcM5PMsamGGPTm2YGj3DOWdN9c9rx49pM4g85eX1ANQKxSAihEEM6rr4pXnkCQrtmc9KyGcYYJgbjelijFxAzFnablLU2szub7KonsFqsYSQVMxR9G0xEkS9bT7jLXISDXDs9naYNzCHO+Ue1Y8egrvFpzvnD7r8qQQCVmpgl0BuewINoXNjnABxtvpFzfrLp1xnO+Zz28wc554c23d4R1bqCRET9MyLhEJLRMNar3lTm6JUrLcJBALBzMKEnOHuBq9rOflcbTwAAJrMJX2YymLFSrGM4HTO8bTARRb5iwwiInIDLngBj7CiAk9qs7CnNE9Y3RABy4neCcJuK5N58YcC5EcgCaJ6nOGp0J+1L88WmQ1NmISTt/scZY6cYY6cWFxdbnkC5tjE2lo5HsF61Jy1glQUtKdkqHASoSU2hhNkLXM2VEYuEMGpycW1mciiJ62uVrklKrxSrpuc5mIwgX7bu5YleDQ+kpKfQ2PDMab+33BARhFs0N8i6gV+J4fs45+ILAs75w9qOaVQzEBvgnD/KOZ/hnM+Mj4+3fOKKtLFediAR8aw6SDcCLcJB4vaFfLVnGsau5MrYnU1auhjuyiZQl3lHip1usFpq7QmsdRIOctkT0Nbno9qv0wBOwcKGyM7mhiDMaHgCvWEEcgBGtJ+zAJZN7qe7xtoX4Zj26zLUXVTHqBoajTcjE49g3UbIwA6L+QoYQ9sd9Y6BBGqy0jMNY1dyZezKtvZeBLu0juJu5QVWijWMGMw7ANTEcLkuoyZZ62j2ukRUC/mc3pQHM8XO5oYgzBCimb0SDnoMjYv4FICTAMAY0yeZM8Y2X+RPifsBOKT93jFq08QmI+ChJzCajiPSZq6nGIbSKyGhq7myfnFvx6RmLLpRIVSXFayVW3gCSa36y6KRF46Yh9VBR0VSGNY3RAThCBEO6glPQOyAtJBOrmlH9PSmu85teswDmjcwa3UXZUZlk5peJhHxTDZiPt+6R0CwQwsX9YIRqMsKFgpVTFpICgMNT6AbJa45zXNqlRMAYLlM1EvZCMbYcVEBpK1/ww0RQbiN2+EgRyWigOriGhw70vTzHDaVgho9plPKdVkXFwO89QSW1muGTUybEYnjhR6oEFosVME5MNlC6qKZbCqKWCSk5z/8ZFVrBGuVEwBgOS/gVU5Au+ifYIx9FOru/37O+WnG2IzBhoggXEV4AimXmsUcG4FuU9k0cDkRDesxM7dZKdZwgzY9rBUTPeQJiHPY0SaZLWCMYWIgrpfD+omYFTCSah0OstoroCjedAxrRQ3DBsdd29wQhBlik5uOuXP5DrxsRLm+sTooFQuj7NFQmeViFSMWyizjkTBG0jFc7wkjYK2stZmJgXhXPYGRjPF7PCSMgMWcAElJE/1ISRiBOBkBAGqmPNHkFiWjYZTrsuvlmaWahEpdMb1AbWY809lgdLcxE2RrxcRAoitGYKXYxhNICE+g+zkBgugWRZfDQX1gBDaq6SVjYShc1ddwE3GBstJwBQCjmZjtUYhecH2tgkiIWT5vQA1ndSOUtaq9x1nTcJC687GaEyABOaIfKdUkhEMM8UhvlIh2HTUn0PgzRGhIJE/cQt+lpq3F1kfSMSwXu28E5rVRjXbKJHcMJlCoSLpCq1/kynWkY2HETBZ3MhpGJMQsh4O8ko0giG5SrMpIx8KudcIH2ggoCoek8A0XDSEhUXb5ArasGwFrO+qxHgoHtZp/YISogPK7uilfruvJXyMYYxhMWheRo6EyRD9SrEqu5QOAgBsBSdvpRcM+eALrNsNB6RgKFclyd6tXzOcrliuDBKIXwm8l1Hylrsf9zbAjJ03hIKIfKdVk1/IBQOCNgHqBbXb3Rbmo256AHg6ymBge1WbirnQ5JDSfr9pKCgNNfQ4+J4fzZUmP+5uRiUcsewIiLUSTxYh+oliT9AFabhBoI1CX1Z1epMkICAvptiewXKwhGmYYsPjmj2rGopshoUpdxlq5bt8IDIpwUO95Aul4GCWLJcB6dVCgVzlBbKRYlZByqUcACLgRkI3CQR7lBFaKVQybDDsxQnQxdzM53GpUYytGUjFEQsx/T6DSOicAiI5wa58t5zRjmOg/ilVZH6DlBoE2ApK8NRzkZXWQ1aQw0KgiWu6iJ7BU1IxAxp4RCIUYRjMx372YfFnCYKL1Dicdj6BoURZEpsQw0YeUauQJ6NR1T8D7nMBysaaHeKwg7tvNXoElbSc/ZtMIAKoR8/PcFYWjYMETSMcj9sNB5AkQfUSxJlN1kECWRR1448/wKiewWqxZ7hEAgIF4BLFwSN+Nd4Ml7SI+NmDdeAnGMv72ORRrEhSO9jmBWNiyQKDikYAcQXSTYlVCmqqDVOpadVCzJ5D0yBPIlesYNhl2YgRjakhlpYuegAhF2QljCdRmN/8MmJD/blcdlI5HUKkreiiwFQpVBxF9hqJwtUSUPAEV4e5HQt4mhmWFY61cR7ZNqGIzoz7vpjeztF7FYCKCeMT+rmE0HffVgIku4HaegCiNK1rw9GR9xrDDkyOIHkFc1zKUGFapGySG45EQGHM3HFSo1MG5uaaNGaPpeHcTw+s1jNmsDBKMZmIo1mTfpCOEKJyVnAAAS8lhhWQjiD6jqOXDeioxzBg7xhg7yhh7yOT2E9r/x60+xiqSvDUxzBhTlURdNAJi4lXWRjgIULuGl7qZGF6vYsxGHqMZ0RntlycjGsAGLFQHAbCUHPZqqAxBdIuiVh7dMyWi2qBtMWQjJ37fxHHG2Cy0EZMWH2MJyWSnJ+Sk3ULo3Ns2Ahk1ru62rLVVltarHSWFgUbHs1+ejNVwkEiIWekVEPMEqDqI6BeEB9xLnsCDUAdsA+pF/qjBfT7IOT+kXfStPsYSIjkY3TT4PeGyEciVhSdgMxyUiaNSV1ByuVLJKkvrtY7KQ4FGMtlvT8CLcBA5AkS/IK4lvSQbkQWw0vT7qMF9pjaFfto+hjF2nDF2ijF2anFx0fTFJWWrbAQgpou5d+FdE+Egu4nhdPd6BWqSgrVyHaMdhoPGfO5zEKJw7cJBYvFbKRMlKWmi32jkBHokHGQFzvnDmhcwqg3htvKYRznnM5zzmfHxcdP76UYgvCkcFHPZEyi1HnZihtiF+1lqKRDCdZ2Gg3RPwK9wULmOVCy8xavbjB1PQCYVUaLPKLo8WhJwbgRyAEa0n7MAlptv1Hb0x7RflwFMtXuMHUQ4KBIyCAe56AmslkS82t4b382uYSH50Gk4KBOPIBYJ+aaCakU8DmgkxKyUiHJKDBN9RiMn0DuewGNQL+zQ/j8JAIyxrHbslDgG4JD2u+FjOsHUE4iGXS1tXCvXMZiIINJml7qZRlzdf0+gYQQ68wQYYxjzsbrJiow0AKRjNjwBrVmMPAGiX2g0VdoLTbfCkRHgnJ8GAC3MkxO/A3i66fYHNG9glnN+usVjbCPJW5vFANVKupmMzZVqtkNBAPR4fDfKRHXJiA49AUCdnbDikwGz6gmkYmEwZjMcFOhuGIJosFauI8SAjIvVQY6fiXP+qMGxI21u33KsE8RQGSNPwN0S0brt8lBAzU2kY+GuhoNGHRiB0XTct+qgQkWyJHnNGEM6FrGUGNabxcgTIPqEfLmOgUTUVWXcQO+RJIOhMgCQiLkbDsqV6x15AoB6Ee5GOGh5vYpENORIaGo0HfOxOqhuOeeSjodRstQnQDkBor/IV6yFTe0QbCOgewKbwkEuJ4bXSjXb5aGC0Yx/F9JmRI+A1SE4RvjZ7NZuyHwz6XgE61Y6hhWhHURGgOgP1sp1DLmYDwACbwQ02YjNHcOxMEp12bWLl+oJdGgE0v4PZwHUcJCTUBCgzhTwo9mNc458RWrbIyBIx6wNlvHaE9jc7e6WHApBmJEvW8ud2SHYRkA2/pInomFwDlSl9nLD7eCcI+/A+mZTMV17yE8WC1WMd1gZJBAlrl6XiZZqMmSFW17c6XjYVnWQFzkBrbDh8abfXZNDIQgzyBPYhFAR3RwOEjMF3MgLFGsyFN6+k9WMkXRM1x7yk+Vi55IRAlFe6rUno+sGWQ0HxSKWvBPFw+og7WI/13TINTkUgjDDahWdHQJtBGSD8ZJAo5HCjTCGrmnT4RufTUVRlRTXJ521QlE4VmyOwzRCTFLz2hPQZaQtvscJix3hir/jJa1IqBCEI6z209gh0EbAVEXUxcEydnepmxnWqor89AZWSzXICnfsCfilfdR4j60tbquJf11KmhLDRB9QkxSU6zKFg5oR4aCogWwE4M5gmYJFYTMzxEhKP42AqO13bAQy/iiJ2vW2rGpD6Z6APyWibeVQrAojEoQRTjekZgTaCMgKB2Nbv+Ru5gSch4PUC6mfyeGlgmgUcxYOSsUiSEbDnovI2V3cSYsd4Qr3VUa6rRyKVWFEgjBiTbsWkSfQRF3mW7wAwOWcQADDQW55AoCa2PYvJ2DN20pGw6hJip4TMkPm3Mvy0GMAZoRAoptyKARhhNMNqRnuZhh8RlaULZIRQFM4yAVPQBdschwO8s8TEAZnuMMu52bGMjEs+RQOGrAaDmr6fFsN11AU7llSmHP+BIAnNh1zRQ6FIIxYK9vLnVkl8J6A0U5PJIbdDAdZvUBtRg8H+aTBAzQSucMdNrg1M5yOYdXjcy9UJSSjYcQi1paj8PTa5XxkxTtPgCD8RgxeonBQE5KiGA4hsXqRsEK+IiERDVm+QG0mFlH1e/z2BIaSUdvS10YMp7zvc1AlI6zvbhIWcz5qToCMANEfeBUOCrQRkBW+RTwOaIQL3MgJFFxozlC7hv3zBFaKNX2WgVOyqajnSW27DTBiyHa7z1fhnOYLE33DmsU53HYJtBGoy8ZGwM2cQL5sXdPGjOF0FCsBNQLDqRjWqxJqLkhwmGH3PU7G1GXb7vOlcBDRTyyv15COhfXrm1s4NgLtRLO02ujjjLETTcdOiNucvLYkK4Yhj3gkBMZcyglUrKtbmqGGVPwLB60Ua64khYFGXiFX9s6I2X2Pk1HhCbTWD1I8rA4iCL9ZLlYxZmHmhl0cGYF2ollaudxJrWpiqmnQ/HHG2Cw2aq/YRjIJBzHGkIq6M13MDdW+4S6Eg0ZdCwd53+dg9z22mvhXOCcZaaJvWFqvuva9bsapJ9BONGuq6dgcGs00H+ScH9KMR8dIMjcsEQWsd5W2o2BD4tiM4VTU8wobAeccq6Uahl0MBwHw9PztDsqwmvORFU6SEUTfsLzuXBTSCKdGoKVoltYhKWqnp6EOmgc0r6BNCKlte73qCRj/CYloGBWXmsWchoOyqRjyFQmS7F1cXVCoSqjL3EVPwNs+ByHVbS8xbLVElKaKEf2DGzNCjPAlMayFifQh85zzhzUvYLQpRKRjtb1eMmkWA9ybM5wvSy6Eg9THi+y+l4gdu2ueQFqEg7zxBMp1GZLCbRlaqyWinHMaMk/0BbKmDOx0RogRTr8ibUWzNI5yzj8K6Lv8Y9rxZTRCRLaRTKqDAHW36DQnUKnLqMmKC9VBQjrCeyMgJCPc8gS87ni2KyMNWJcFkbl3HcME4SerpRoUjp70BAxFsxhjWXEHxthxzvnD2s9HoYaERC7gEBohIttIitIyHOTUE3BLta+RXPU+L+C2JyA6eb06d7sy0oD1EmDKCRD9glAB6LmcQAvRrKebjp9gjM0yxlabHvOA5g3MOhHaapcYdloialfYzAw/9YOE2NuISyWijDE1se2REShU7HdBhkMM8UiobU5A4dwvGWmC8BQx3c+pMrARjpWIjESzOOdHtP9PAhi28phOkBSOlIk0QioWxpVVZ0agkwuUEX4qiepGwMXF4mWfg25obXpbVqq/FIUGyhD9gTACPecJdBs1HGT8JU+40CcgBJucqvaJChs/wkErpRpiYVWvyC1U6QiPw0E2va2khc9X5uq8CYIIOkt6OKj3EsNdpVViOBl1IxzkjieQiUcQCTF/wkHrqmSEm01S3noCneVdrHkC1DFM9AfL61VEQsx1BVEg6EZA4YYqooA7JaL6LAGHbzxjDNmU95LMAFxtFBN4ee75Dsd3Ji30gXg5VIYg/OT6WgU7BhOedMAH2wjIiumXPKXtFDlvPX2qFSJU4bREFABG0t4lV5tZdlEyQjCciiJXrjt6L83Il+uIR0KIR+yFr6yUAJOUNNEvXF4tY89w0pPnDrQRqLeoDkrEwuAcqDpQv8yX64iEmC5T4ISsTyJyq0X3PYHhVAyywvVdu5t02pFtpQRYnSzW6ZkRRO9webWEPcMpT5470EbAbJ4A0DSC0EFyuFCRMJiMuuKCDXuYXG3GTfE4gZeJbbUj276nlYqFabIYsS2oSQqu5yvkCRihykaYl4gCQMlBXiBfqbsSCgL8kZOuywryFck1GWnBiIcdz516AlZyPtQxTPQD19cqUDjICBghKRzRFiWigDNPwA0ZaYGYLuZFXF0gcg4jaXcrCLIe9jl0+h4nLeQEOCWGiT7g8moJACgcZIQkc4RNZCOSFkXGWlGwKXHciuFUFHWZo+iCsqkZeqNY2t2GkmEvw0FayM0uyWik7WcrK+QJEMHn8moZAHkChtRlBVGTxLDVObStyFfqGIi7s6v2Q5d/RdcNctcTaJy7B+Ggcr2jnEAyFmpb/SVzkGwEEXgur5YQDjFMDiU8ef5AGwFZaaUdZG0ObSvyZfc8gUZy1bu8wIpa3Z/uAAAgAElEQVSuIOquJ6Amx933BDjnevLdLqlYBLLCUWsxo0FROEyWB0EEhkurZewcTJjmP50SWCPAOYektAoHqRfvcps5tK3IV9zLCTTkpL3zBFY98gTCWqei24nhqqSgJisdvcf6TIFaCyNAOQGiDzh3vYAbdmQ8e/7AGgFJUcMAZolhMYe2U09AkhWUajIG3DICupKol+GguvZa7uuLqNVN7p57QzKisxJRACjVzY28rNCMYSLY1CQFLy8UcOvkoGevEVgjIGtGIGyaExDVQZ01ixVcEo8TZH3JCVQxmIiYSmk4QRWRc9cTaHRkd1YiCrSu/lI4zRMggs35hQLqMsetu7wzAu5c4bpAXYsFR1sMlQGAUofhoLxLMtKCbNL7mQIrpbpe0+82w6kY5vMVV59zzcG8BiuDZfq5WexKrozra2XcvHMQ6bh3X+NyTcaL8wUsr1cRDjHsHUlh73AKsUhv7R/LNRlXciVcyalrNBkNY89wEpND3ujtWKVSl3HmWl7v5L91clBfu1Y4czUPAJ56AoE1AsITaDVjGOi8RLTQobCZGZFwCIOJiKddw15IRgiyqShevF5w9TmdTG6zMmye92l10H/79iv42JfPoSYryKai+M2fvh3/6M5drj1/TVLwF89dwZdeuIbvzi2jtkl6JRpmuOfQGN55+0686zWTrm2U7DK7uI7//fxVfPP8In54KQfFoFBsLBPH3QeH8e47d+FtN03YugA74aX5Aj7zjVl8+YVrqNQb799APIJjM3vwL3/yBj060Ioz1/JIxcLYP5r27FwdX+G0CWE5ANNijGS729s9xgp1WQsHmXzJY5EQIiHWcYlopxLHrRhOe9s1vFKsYVfWmzIyT3MCHTaLAW08Ae6vdpAb67odZ6/l8Zt/fQZvuXEcv3D3PnzmG7P453/6HFZLNbz/ngOOnltWOP7n372KT//Ny7i2VsHBsTTe94b9uPvgCHYMJiDJCl5dKeHM1TyePDOP3/jzF/Cfv3QW73n9Pnzg3oOYGPRm7TXDOceTZ+bxx996Bd9/ZQUhBrxmTxa/8hOHcXgig13ZJEKMoVST8MpSEc9fyuGbLy3iyy9cRzYVxfvvOYBfeuMBzzzmYlXCb3/tRXzuuxeQiobxf961B2+9cRw7hxKYz1fwtR9dx+e+cwF//ffX8On3TuN1B0ZaPt/zl3K4ZXLQU4/WkRFgjE0D6gQxxtgUY2y6eVyk0e3iNrPHWEXUh7d6c5zISbsdDgKEiJyHnkCphts8ih0Op6Io1WRUJdm24qcZTob2JKPth837OWO43XfBLf7zl89iIBHF7z74WmRTMfzETRP4yBdO4//+yx9DUTh+6U0HO3relxcK+PUn/h7PvZrD6w4M48TPvQZvvmFsSyhl5sAIfnYa+LfvugU/vLyGP3xmDp99Zg6f++4F/OI9B/Dhtx7yxBsVF/9PnjyPM9fy2DuSxEP/8CYcO7IHEwPGxufNN4zj/feoRR7fmV3G5793EZ96+jw++805/OIbD+BDb5ly9VyfOb+I3/jzF3B5tYz337Mf/+rojVue/x237cQv33sQ//xPn8N7P/t9fOZ903j7zTsMn29pvYrnL+XwL3/yBtfO0QinnsCDAJ7Sfp4DcBTA6Ta3j7Z5jCVkzQi06ghNWhAZM8OtqWLNDKei+sBot+GcY7lY82yHI1zXXKmOHYMuGQEXPIFW4T5F8XXGcLvvgmOurZXxzPkl/F/33ah/HrFICJ9+7zR+5X+exv/7V2egcOCX77VuCOqygke+MYtPPf0y0vEwPvnzr8W779zVNo7OGMNr92bxe++ZxsXlIj558jwefWYOX/j+q/jAvQfxT9980JXKOkVRL/6felq9+B8cS+MTD9yJd9+5y3LdfCQcwltuHMdbbhzHywsF/N7XX8Yj35zFF753ER9480F84F5n57pWruM/fekMvnjqMqbG03j8w/e03OHfvnsIf/GRN+Ef/9H38aE/eRZ/8N4jOHrrVkPw9XML4Bw4eouxkXALp9mdLICVpt9HLdze7jFgjB1njJ1ijJ1aXFw0fGG9OqidEejUEyh3XrlihhchFUGpJqMmKZ7lBLyYk5yv1BGLhDqK0+oloi2rg3ydMezKum7F02cXAADvvGPnhuOxSAi//55pvOO2HfiPf30Gf/jMnKXn+9GVNfz0730bv/PkS/gHt+3AU7/2Vvz0a3fbTqTuH03jEw++Fk/+6lvw5hvG8Mmnz+PND/8NPvON2Y43YYrC8eUXruGnPvUMPvz5Z1Guy/jEA3fiqX/1Fvzs9J6OG6cOTwzgd3/+LnztV9+Ce28Yw++eVM/1D/521nYRCefqOd73iW/gz05fwT972yF8+V+8uW2IBwCGUlF8/p++HrdODuKffeFZPPnj61vu89SZeUwOJTzz7gU9mRjWBtE/CgAzMzOGugCKlmtptdNLRp15AoypiRy38KLMUtDQDfLKCGjVTS5KR6gy0p0ZWSslojLnMCke6wpW1nUrTp6dx4HRFA6Nb20cikVC+L33TONf/Olz+K0vncV8voKH/uHNhuXCy+tV/Nevv4w/+d5FjKZjeOR9R/CO23ZuuZ9dbtgxgD/4x0fwwuU1fPypF/FfvnIOf/jMK/iVnziEnzuyx9Jnna/U8ZfPXcH/+O5FnF9Yx9R4Gr/74Gvxf7xm0tWO2RubzvUTT72IE189hz/6lnqu73n9vpYhT845Tr+6iv/ylXP4wYVV3Do5iD/+pdfh9t1Dts5hKBnF//jA6/H+P/47fOQLp/H7753WP4czV/N4+uw8PvjmKc+rm5xe4XIAhNnLAli2eHurx1hC0cNB5vdx4gkUKnVk4hFXwwnDqRjWqxJqkuJ6iZ2uIOpBoxjQHA5yzxMoVOodh9uslIgq/grItfsuOKIqyfjO7DLe94b9pheFaDiET/3CXfiPf3UGn33mFfzNi4v48FsP4d7DY0jFw5hdWMdXf3wdX/jeqyjXZfz86/bioXfcjKGUu9U9d+wZwn//J3fjBxdW8NtfexH/4a/O4GNfOYe33zSBNx0exR17stg1lEBCC9deXi3hR1fy+OZLi/j27BIqdQW37RrEp37hLrzrjklPk6J37BnCf/snd+PUhRX8zpPquX7y6fM4essOvP3mCRwaz2D3cBKyzHFptYTvzS3jiWcv49z1AsYyMXzsZ+/A/Uc690yGklH8yQfuxvv/6O/wK184jV9/x034iZsn8G/+1wsYTEbxkbcddvkv3opTI/AYgBnt5ykAJwGAMZblnOfMbjc5ZgvZamK44+qgznepZuhqnOWaaTKrU5Z1yQiPPIG0+30O+Urn73E8EkKItfYEpBZDhzzAbK27wvn5ddQkBUf2D7e8XzQcwm/+zO14203j+E9fPot//fgPN9weDjG847Yd+LX7bsLhCe+kCADgdQdG8NjxN+D5Szn85fNX8eUXruGrBmEPwf7RFB6Y2Yv7j+zFHXvs7aqdMnNgBP/f8XvwndklPH7qMr724+t44tnLhve9Y/cQfutnbsfP3LUbGRciBYMJ1RD82hd/iI995Rw+9pVziEVC+J3773TdQBvh6C/gnJ9mjM0wxo4CyDVVQzwN4IjZ7SaPsYWitE8Mp2LhjsMvbg6UEYgLdK5Ud90IrOricQHKCZQ7GygDqInJdtVfqsCgP/GgFt8FV7DbNPST2k72uUs5/PjKGqqSgl3ZJO4+OIKxjLsCg61gjOGufcO4a98w/p9/dCuurlXwwuUclos1lGsy4pEQ9oykcHg8g70j3ujl2+GNh8bwxkNjqEkKzl3P48JyCddyZU3FM4nX7Bny5DwHElE8+r4jOHVxFS/NF/CmQ2M4MOZdb0Azjq9yWpxz87EjbW7fcswuVjyBRDTsoFms8wuUGeJCuuKBdMSKx55AIhpGIhpyNRyUr9QdaaQnY5GWieG6rPjpCbiyrs04cy2PdCyMfTYuQIwxTO8bxvS+1t6DXzDGsDubxO6sN7r4bhKLhPCaPVm8Zk/Wt9dkjOF1B0YsJZbdpIfSZvbQE8MtvuMpR9VBnc2+bYWXs3pXSzWEQ8z1c27G7RGZ+bLkqPoqGQu1NPKtpMaDxpmredwyOdiXHdBEdwmuEbDSJxBtP4LQDDdlpAWNkIr7FUIrxRqGUzFPKwnEiEy3yDtIDANAKhoxLetrJzUeJBSF48y1vKciYsT2JbDfEL1PoFViOBZx1CfgVTjIi16BlWLNs3yAYDjl3kyBSl3ta3BiaBOxMMp1Y5XYdlLjQeLqWhnrVQk37yQjQLhPYI2A7gm0qQ6qSYpuMCw/t8KxXpVcTwwnY2pc3Qs56dVi3fVhMptxs9nNiXicIBkNoWLi6TUEBgO7xHXEjNn9o91PnBL9R2C/IZbCQR2OmCzWJCjcXd0ggdtxdcFKyTvJCEE2FXXNgOUdyEgLUrGI6VAZITXuZ2LYK7weNE5sbwJrBMRo2dayEWLEpD0j4IVukMDtuLpA5AS8ZCQdw1q5rpfnOsEdT8C8D0SSW0uNB4nLqyUwBkwOkREg3CfARkCEg8zvY0VawIiCg4lX7XAzri6QFY5cyfucQDYVg8IbF3AnOBGPE7QSCBQ5gX7xBHYMJHpukAvRHwR2VelS0m2axQD74aBGqMIDI5COuZ4TyJfrULh3PQKCxpxk50ZADO0ZcuBttWoWk7Qa4v7ICZQoFER4RmC/IbLFxDBgf8SkkwHo7VA9AXeNwErJW/E4gZvVTW7Ma0jFzEuA9XBQn3gCZAQIrwiuEbAgG2Fl+pQRhar7A2UEwyn34uoCvVvY45yA8DTc8GR0b8tBTiARDaMqKYbvpR4OCnhOQJIVXFurYM8wVQYR3hBYI6BYFJAD7M8ZFhcot0tEAXfj6gKvZaQFboaD8pU6YuEQ4g7i3PpgGWnr5yvp1UGBXeIAgOv5CmSFkydAeEZgvyGKpeqg9oNHjPBioIxgRKvld1M/aNUnI+CmnHS+rAr0OelwbpX415vFAu4JXF+rAAAmA6C3QwSTwBoBkRNodQ3puDqoKiEZDXtSjZH1QDrCr5zAYCKCcIi5lBOQHHdktwr3iZxA0GUjltarAICxjLefLbF9Cew3RLEkG9FpdZD7MtKCYQ+Gs6ys15CKhTsa02gHxhiGU1GsuDBdLF+uOxa7a2Xk63p1ULA9gUVtJrWf8s/E9iKwRsCKlLReImq7Wcx93SCBm3F1wUrJ+0YxgVvNbm68x8kW08XkPukTWNY8Aa+9PGL7ElgjIApCWlUHJSKd5QQKFfdlpAVuVtgIVoveS0YI3CpxVT0Bl8JBRp5AnySGl9aryKaihrOCCcINHK8sxtgxxthRxthDJrcf1/6daDp2QtzW6es2JouZ3ycUYkhEW2vOG6GGg7zxBAbiEURciqsLVkp1zxvFBMOpWMfT2ppRcwIOw0Etwn1ynySGl9drFAoiPMWREWCMTQMA5/wkgJz4ven2owBOahOXprTfAeA4Y2wWwFynr21FShpo3VVqhhtJSzMYY6oQm5vhoGLVc8kIwXAq5kplkxuGtmV1kGxtffQ6S+tVSgoTnuLUE3gQQE77eQ7A0U23TzUdm9N+B4APcs4PacajI2QLKqJAZ4NlChXnSctWuC0it1qs+5cTSEeRK9V12Y5OqNRlVCUFQx7mBEQ4KOhhlOX1GkbJEyA8xOk3JAtgpen30eYbOeePNs1dnQZwSvt5ykII6RRj7NTi4qLhC3MLiWFAExmz4Qlwzh2PPWyHWmHjjhGoSjLWq5Lef+A1w6kYarLS8cQ2oNGH4dQItNKGkvukY3hxvYpxMgKEh/iyTdLCRKc556cBgHP+sOYFjDaFiHQ04zHDOZ8ZHx83fE5ZnzFswQjYuGBVJQU1WfFEN0jgVlwdgP48fuUERjSPw4kRy2lGQMxc7pREq8RwH1QHVSUZhYrkW6iP2J60vdKZJG/nRB4AwIh2LAtg2eRpjnLOP9r0fCuc8ye0+0+ZPKYlDQG51vdLRSO2jIAbEsftGE7F8PylXPs7WmBZqyP360IhLty5Uh17R9rc2YQ1lzyB1jmB4FcHic92bIA8AcI72hqBpnCOEY8BmNF+ngJwEgAYY1nOeU77+Tjn/GHt56NQQ0IiIXwIwCOdnLjeLNbGE0jEwvqF3QpioIxXzWLAxri608HwosrIr5yAXuLqIKchvJds0tk5R8MhRELMuGPYYuFAL+O3gSe2J462SSK8o13cc+J3AE83HT/BGJtljK02PeYBxtgxALNNj7GFFQE5QJ1Da8sTcGHiVTtEXL3oIK4u8Es8TtBoduvcCLjlCQDmOR9RHRTkxPBSUW0Uo8Qw4SWOt7tGngLn/Ij2/0kAw1YeYxeR+Gu3k07FIrYSw41wkHeegIirrxZryMSdvY7uCfhYIgo4a3YTlVGuGAGTEZNyH8hG5F3KnRBEKwK7TbLqCSRslojm9YlX3n3xmuPqTller4ExIOvh+TYj3hcnfQ75ch2MuRNyM/ME6n0wVGbNh/wUQQTWCFgZNA+oZYR2Oob1L56HF1URullxoVdgtVTDUDLq2xjFSDiEwUTEUZ/DmiYZ0WoqnFXMPIF+GC/p5YQ7ghAE9huiWKwOEh3DVpub/KgO0o2AFvN1wkqxpoeX/GIkHcOKA08gV667FuIwzQn4VCJq0CXfUkbFDvmKhEQ0hHjEW3VYYnsTXCNgYbwkoF4kZIWjJlyHNqyV64hHQp7KMotEn6j+cMJKseZbPkDgtON5rVx3Ldxm6gn4EA7SCh8eb/q9pYyKXdZKzkX2CKIdgTUCupS0BdkIAKjUrBmBvIsXKDMGExFEwwzLLnQNd0NbxqmSaK7kshHoUomodrFv1r9qJ6Nii3zF+7VIEIE1AronYEE2AgBKdcnS87q5SzWDMYaRdEzXindCN1Qmh1MxrDoYLOOmoTUvEVUQCTHHfRg2aSmjYpe1sndzLQhCEFwjwK3t8uwOlvHrizeSjjvWD5JkBSsl/wXGhtMxx30CbnoCFZMZw71WHmpFE6uZvMdChgQBuNAn0C1kzlvOEhAkWihNGpGv1DExkHByapYYy8Sw5DAnsFKqgXNgvAvhoFJNRlWSbSctOeeuJ4ZLJs1iUYeSEW0kU4xoKaOi9cc8CgAzMzNtKxXyZQmHxzPWT5ggOiCwRkBReNukMNDwBKz2CqyV67hhYsDRuVlhJB3DxeWSo+dY7tL82aw+J7mOHYP2jECxJkNWuPeJYUVB2KEn0EFTo6GMSqdQOIjwg8CGg2SFWwoHpbWO3PWqxZyAi0nLVoym445zAkva4/0WGGuUuNr3ZERVkVPdIEEyFkZVUvQckaAuc8/F4zTpkxnt/1YyKrZRFI4CJYYJHwiuJ8Dbl4cC0GUZStX2noCicBSq3s0XbmY0E0OxJqNSlzsuRxVGwG+BMSdGwO1mPL36S5KRijU+N1lRPO8R0JRwn9h0zLEkCgCs1yQonLqFCe8JrCegWMwJiHBQ0YInUKhK4NzbbmGBuHA7KRNdKnRHaliUpC514MmsuayHkzQJ90ly7yWG7UDdwoRfBNYIWA0HZWyEg9yaeGUFUdGz4iA5vFSsIhYJYcChCJ1dRtOdN7utldx9jxMmMwUkhQdaQTRf9l7DiiCAIBsBbi8nYMUT8EM3SCBCKksOpCOWCjWMpWN+18JjKBlFOMQchYPcurgJT2+zPpSkKIGeJUDicYRfBNYIWB3IEg2HEIuEsF6zbgT82H2JkIoT6Yil9WpXpk6FQgzDqRiWOzBgbo2WFJgNm1cTw8E1An7MtSAIIMBGQFZ4W8kIQSYesZQY9jMc5IaI3HKx6nt5qKDTPoe1ch3RMNMv3k4Rz7M5JyAHPhxEngDhD46/Je1UExljJ7T/j1t9jBVkxbouTCoWthUO8sMIZOIRxCIhZ55Aoda10YOjmc5kL9RuYfdCWPqw+S2eQLDDQWK9puOkIEp4iyMjYFE18ThjbBaa0JZbSosK521lpAWZeMRSYtjPnABjDKPpWMfVQZxz1RPo0hDy0XS8o3NX+zDcS2TrOQGD6qBogKuDxOjRtM9Jf2L74dQTsKKa+EHO+aGmVvu2j7GisaKWiFr7kqfjERQt5gTCIYZ0zJ/dV6e7aUCtHqnLvGvhIPXcOwsHuelpmeUErFaP9SqlmoRwiCEeCW5IiwgGTleYFdXEqU2hn7aP4Zw/yjmf4ZzPjI+PG76wnZxAOh7BupWcgNah6Ve1zagDEblF0S3ss26QYCwTx3pVsjW1DQBy5ZouO+EGZjmBuqIEOidQrMpIxcK+V34R2w/PvyWc84c1L2BUa6d3BTUcZDUxHEbJUjhI8rUuezTduYicLhnRLU+gw2Y31z0BkxJRWQl2dVCxKuk9LgThJW1XWRslxZaqidpjV7T2+mWoolotH2MVRWk/UEaQikUsJ4b9lO4dzcQ69gS6bQT06qb1GnZnk5Yf57Y2k1mzWF3mCHusHeQlpZqs5zsIwkvaXvHaaKEYqiYyxrKc8xyAU2hMXjoE4BHtmGOlRZlzWPWUrSaG8z6rNo6k4yjXZZRq0gbdGyuIePxol8JBouPZTrObrHDkK+56W9FwCNEw25ITkGQl4IlhiZLChC842iq1UE18uun2BzSVxVnO+Wm3lBYVG4m/dDyMYq39sHk/Rks2M+qgYWxpvYoQU6d8dYNOmt0KFXcbxQSJaNiwTyDIieFiVULa5saAIDrB8Soz8hQ450fa3O5YadGqbASgJoZlhaMqKS0VO/3Wb2+Oq+8dSdl67NJ6FSPpeNcudMITsFPdlHNZN0iQioW35AT6ITG8K0uNYoT3BPZbYlVKGoC+o2qVF+Cc+zJfuBldRK6DruGl9VrXKoMAIB0LIx4J2UoMr2izBIZdbnAzGjYvB1w2opMQIUF0QnCNgGJNShpoFpEzL2cs12VILk68soLwBDqpEFpa755kBKA2u41l4rbkpIViqttdzgmD6WL1HpwxbIf1qkzdwoQvBNYI2In5ZrQvU6vksJ+SEQJxEe9El181At3zBAD7DWOiEmrEbU8gttUTkGTF88liXlKqUU6A8IfAfktkmx3DAFp2DXdDujcZC2MgHsFC3p4R4JxjsVDFeJckIwQjaXslriJ0JOYRuEUqttUTkALsCSgKV0tEqTqI8IHAGgFuMzEMtPEEPEpatmN8MI7Fgj0jkK9IqNQV7BhMeHRW1rA7J3mlWEUyGtYbvNzCKCcgBTgnIP4Wv+RLiO1NYI2ArFj3BMTkrfWKuRHIV7ozyWliII75fMXWYxa0+3fbCIxlYlgq1tqW3gqWizXXQ0GAlhMwGCoTCWh1UENBlDwBwnuC+S0BIHNYlo0QZZ9iUIcR3cgJAMDEQAILNj2B6z1iBEYzMdQkxVIjHqDmBLwwAqbhoIB6Ag0FUfIECO8JrBFQFA6rIV9xYRdzW41Y69Jg74mBOBYKFcu7aQCY13IIOwa7mxOwO2vYKyOwORwkyQo4R2D7BIQnQCWihB8E81sCVUDOak4gHgkhFg7pF3oj1ko1MAYM+DzJaWIwjkpdQcHibhqAHj6aGOi+JwDA8pjJ5XVvhuAkNnkCFUkBANeml/mNHg4iI0D4QGCNgKxYmzEMqDXtg8lIy3DQSqmG4VTM9w5cEdKxUyG0kK9gMBFxPcFql0aJa/c9gaqkQFFUb0oYhEQ0mMu7ROEgwkeC+S2B5gnY0FofTET1ua1GrBRrGHZZ08YKosxzoWA9OTyfr2LnUHe9AKBhBKxUN5VrMsp1GSMe9DakNo2YFBISrSRCehlRykyJYcIPAmsE7AqEDSajegWQEV7tUtshQjp2PIH5QqXrSWFArQ5iDJYS2yJk5EU4aPN0saAbgZLW2U5S0oQfBNYIcBvVQYBqBFrlBFaL9e4YgUH7nsBCvtr1fAAARMIhjGXieslqK1aL6ns/4nKjGLB1poAwBkHNCYhqKxoqQ/hBYI2A2jFs/f6DiQgKLYyAVzXs7RiIR5CMhvWKn3YoCsdCoaIbj26zYzCul6y2YnFdvY8X8w8ymzrCK3U1MRxUT0A3YuQJED4QXCNgY8YwIMJBxkaAc45VLTHsN4wxTGYTuLZWtnT/pfUq6jLHLhvTvLxkx0DCkgFb0Mta3fdgMomNzYCNi2gwl3e5JiPEgFhAS1yJYBFYf1NRrM8YBkRiWALnW6uK8hUJssK74gkAwK6hJK7mrIWDLudUY7E72/1wEABMDCbw/KVc2/uJvMG4B8qnwhMQZbZ+5QSaRq8e4px/VDt2DOoI1WnO+cOdPG+lLiMRpSHzhD843mowxo4xxo4yxh4yuG2aMcYZY7Pav0e04ye0/43mF1tC4dZnDANqw1hNVvRQQTNeqVtaZXIogas5a56AuF/PeAKDcSwXa6hJW9/XZhYKFQynoohF3N/dZjbJgvhhBLTJeCe1AUlT2ndgGgDE/G3xu13KmhEgCD9w9I20sOhHOOeMc34IwP0ATmjHjzPGZtGYP2wbmXPYUQoWncBGISFhBNwedmKVyWwSi+vVthdSoGEE7Ax395KdWnhnsY2QnJfJbD0ctMkT8DgxPAXgqPbznPb7g1C9AHHsaPMDGGPHGWOnGGOnFhcXTZ+4UlcCm9QmgofTbVnLRa8ZB8EM51xc9D/IOT+06XYdK18WxYaAHNCQiDbqFVgVnkCX5vXuGkqAc1gSkruaq2AgEfG9s9kMEeNvd+4LhapnyWw9MawZgUazmHcXUs75o01jUqcBnAKQBbDSdLdRg8fMcM5nxsfHTZ+7IsmIB7TRjQgeTlday0Uv0FznLzYdmjILIQHWvix2ZgwDrUXkuh0OEqGda2vtjcCVXLlnvACgqcS1nRHIVzybfyDkFQp6Ytg/2QjN+z3NOT/t1nNWajJ5AoRv+JUYvq951y8SZoyx+xhjR808glbY9QSEiJxRr4CYfds9I6Dupq1UCF1ZLfdMPgBoeALXWxgwzjkW170LB4VCDJl4ZEs4KO4w/2CSs5rbtF6PiqQwVK94RMYSxukAAAuDSURBVPs5C2C5k9etSJQTIPyjrRFo80Wwuuj1XIH2fCuc8ye0+0/ZOmMNO4PmAbVPADBWEl0t1hCLhLrWoTk5pF7UrVQIXV0rY3p/1utTssxIKoZYJISrLYzAaqmOuswx4eEktEw8siExHI+EbFWPGdEU7jGEMXa8aUNzFMBjAGa0m6cA2N7cAGo4i3oECL9oawTafBEMFz1jLMs5z2k/b77In0IjIXwIwCN2TligykZYv39Wi/evlraKnS0WqhjPxLtWkpeORzCYiLStECpWJeRKdd1o9AKhEMOebBJXVs3PXXRDe9nglo6HN3gCXl9EtYv+CcbYR6FuhO7nnJ9mjM1ot+U6DRFV6krXvFJi++EoHNRi0T8N4EjTXec2PeY4Y2wFwGynXxS1Osj6RTubjCIcYoZD3VUtnu524O4dSeHSaqnlfcTt+0ZSfpySZXYPJ3G5xbmLRjEvpS4yiajeJ1Cuy0hEvDUCmic8bHC8pfdgBTUxTJ4A4Q+OcwJGi55zfqTp5zkAH2r3mA5e11afQCjEMJqOYamw1ROYz1dxeDzj9JQccWA0jTPX8i3vc2GpqN+3l9gznMSTV83PXeQLdnooejcQj2BdS/pX6kqgwymUGCb8JLB1aHZmDAvGMnFjTyDffU9g32gKl1dLkGTzXoELy+pue/9Yb3kCe4ZTWC7Wtox4FFxeLSHEgEkPu5ybE8NlLScQVCqSEthZCETwCORK45yriWGbib+xga1GoFyTUahImOiyNPOB0RTqMm9ZJnpxuYjRdEzveegVRMnqlZxxSOjSahmTQ0lPxz1mEhEUqw0p6SB7AuWa9+EsghAE0ghoA6RshYMAVf9+8xQskbTstj7/vhE1xHNx2Ty2/spSEQfGeisUBKjhIAC4bJIcvrxawu5hb5PZmXgEBT0cFNyLKOccFSnYRowIFgE1AqoVsFsBOJ6JY7FQ3TDUvVeGth/QQjwXloum97m4XML+0d4KBQFqOAhoZQTK2Dvs7XkPJNRwEOcc5QB7AjVZAefBlcEmgkcgjYCsuQK2w0GZOGqysmHCmJA76LYnsGMggXgkhIsmRqBSl3FtrdJzSWEAmBiIIxpmhtVNNUnB9XxF9xa8Ih2PQOFqPiDI2juVWrBnIRDBI5BGQHgCdofCC9mC5ryAMAJeNjJZIRRi2DeSwitLxuEg4SH0YjgoFGI4MJrG3OJWA3ZtrQzO4bkRaFYSLdeCq71TkYTuUTDPnwgegVxpwhOwnxPQjEDTTNzFQhWxSEiXlegmN+4YwEvzBcPbzmrlozftGPDzlCxzeCKD2YX1LcdFiGivx70NA4nGTIGqFNwSS138LqA5DSJ4BNIIiMSw3QbfsQG1C7M5OXxdKw/thQEet0wO4NWVkp7gbObstQJi4RCmxnvPEwBUI3BxpYSqtLFM9NKK6tl4LXqnD5bRPIGghlOEJxDUnAYRPIJpBJTOwkG6J9AUDrqwVMT+kd64sN68cxAADL2Bs9fyuGFHxtMySyccnshAVjgubApnvTS/jmQ07LkREDILS4UqKlKAcwL6fOTe/JyJ/iOQK03uMCcwnIohFg7hqqbWyTnH3GKxZ3bXt+xSjcCZa8ZG4JbJQb9PyTKHJ9SO65c3hYTOXsvjpp0DjsXc2iGUVV9dKUFWeGAvon7MQiCIZgL5TWmUiNq7sIRDDAfGUnrsemm9hkJVwsEeSbbuGkpgMBHR4/+ChUIFS+u1njYCh8YzYGyjEeCc49z1PG6Z9D6PMZaJIxJieEWT1gjqRbSRGA7m+RPBI5hGQFNWsGsEAC2BqVWxzC2qF6ypLusGCRhjuGVyED++srbh+AuX1d9v7WEjkIiGsW8ktcGALRSqWC3V9TCXl4RDDDsGE3jxuupFiRxB0KjUfBmNSRA6gTQCjXCQ/cceGs/g4nIRVUnGnLZrnOoRTwAA7j44gheurG0YfvPtl5cRj4Rw177emSNgxJH9w/jBhRW9GU8YBL88mJ1DCZx+dRVAIzwVNMgTIPwmkEZAJIY79QQUDlxYKuGVpSJikVBPTeq69/AYFA58d7Yxn+dbLy/i7oMjPX9heMPBUSwXa3pI6MeasuhNO/0pa50cSkDS1saNPr2m21BimPCbQK60TpvFANUTANTY9csL6zg4mu7oebzirn3DSMfCeOb8IgB1Nu9L8+t40+GxLp9Ze14/pQ6Z+/4r6tjpr59bwK2Tg771YEwOqV3fam6l+30fnVCmcBDhM64YAW3YttltxzYPlTc6ZgfZgScgKoGee3UV3355CXcfHGnzCH+JRUJ4w9Qovn5uATVJwZdeuAYAePMNvW8E9o2ksGMwjmfOL+L6WgXPXlzFO2/f6dvri4lrfnkeXkDhIMJvHBsBbarY4ya3TQP6FKYcY2za6Jjd19SrgzrYwadiEbx2bxZ//O1XUJUU/Oz0btvP4TXvu2c/rq1V8Mg3ZvHpv53F6w+O9HRSWMAYw8/ctRtPnpnHw187BwB45x1+GgHVE7jJh0S0V4jEcJDnIRDBwvFK0y7mcyY3Pwh1GD20+xw1OWYLMXfFrmyE4D+8+zZwqAnh1+7tvWTrW28cx5sOj+LjT72ExUIVv/6Om3qio9kKH3nbYWSTUfz56St40+FRHJ7wb1cupCn8KEn1CjFQJiifNxF8vK6jywJYafp91OTYBhhjxwEcB4B9+/ZtedJ0PIyfumNnx5Oq7tybxW8fuxM7BxM9+WVjjOF37r8TT52Zx+GJDGYO9FbIqhVDySh+/z3TuLBcwgMze3x97dt2DeKz75/B224a9/V13eTGHQP4qTsmu30axDaiJ4uptRnEjwLAzMwM33z7nuEUPv3eI1seZ4djR/y9QNllciiJ999zoNun0RFvPDyGNx72/3UZY7jv1h3+v7CLHDuyp+fXJtFftDUC2q58M3NaGKgdOQBiG5sFIOoejY4RBEEQPtPWCGi7clswxrKc8xyAxwDMaIenAAjDYXSMIAiC8Bk3qoOOAZjR/hc8DQCc89PafY4CyHHOTxsdc3oOBEEQRGc4zglwzp8A8MSmY0eaft7iSXTiXRAEQRDuQ8XIBEEQ2xgyAgRBENsYMgIEQRDbGDICBEEQ2xgmtN97FcbYIoCLJjePAVjy8XTM6JXzAHrnXHrlPIDW57Kfc+57i3GbdW1Gt9/Tbr4+/e32sLyue94ItIIxdopzPtP+ntvjPIDeOZdeOQ+gt87FCd3+O7r5+vS3e/faFA4iCILYxpARIAiC2MYE3Qj0StNZr5wH0Dvn0ivnAfTWuTih239HN1+f/naPCHROgCAIgnBG0D0BgiAIwgGBNQJ25xoT24d+XRv9+ncRG/H7cw6kEbA719jjc2n5oTDGTmj/G81l8PJ1fbsodOs9MHmtnlkbbtLtv8vPz7iba7vba7kbn3MgjUAHc409weKHcpwxNgvz83X9df282HXrPTCjV9aG23Tz7/LzM+7m2u6FtdyNzzmQRqANbWcYu4iVD+WDnPNDFiexufW6fl7suvUedIKfa8NPvP67/PyMu7m2e30te/I596MR8BMrH8qUB65ru9f182LXrfeA8A8/P+Nuru1tuZZ7ctC8R3ONuwLn/GEAYIzdxxg72gO7Yd/pofegp9ZGM0Ff8z30GXtKl/9OTz7nnjQCHs017og2X86WH4r22BVt+tqydj5u0G4x+HlR6NZ7YBmv1oabdHvN99A67+ba7sm17PX6DWQ4yO5cYyevxTl/1OCfePMfQ2Mh6B8KYyyrHTuFxgd1SPvdDdq9ruHtHtGt98AQP9eGn3j9d/XQOu/m2u76Wu7G+qWOYYdou4M5AFNiN8cYe1bMWRa7B+32h31+3Q23e0W33gPCP/z8jLu5trfjWiYjQBAEsY0JZDiIIAiCcAcyAgRBENsYMgIEQRDbGDICBEEQ2xgyAgRBENsYMgIEQRDbGDICBEEQ25j/H4TAu7oOiRMgAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "x1 = np.linspace(-1, 1, 128)\n",
    "x2 = np.linspace(-1.1, 1.1, 128)\n",
    "p = np.polynomial.Chebyshev((0, 0, 0, 0, 0, 0, 0, 0, 0, 1), (-1, 1)) #These are Chebyshev series, a proto of \"chebfun system\" in MATLAB\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2)\n",
    "ax1.plot(x1, p(x1))\n",
    "ax1.set_title('Interval $x\\in[-1, 1]$')\n",
    "ax2.plot(x2, p(x2))\n",
    "ax2.set_title('Interval $x\\in[-1.1, 1.1]$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Convergence of the Chebyshev-accelerated Richardson iteration\n",
    "\n",
    "Note that $p(x) = (1-\\tau_n x)\\dots (1-\\tau_0 x)$, hence roots of $p(x)$ are $1/\\tau_i$ and that we additionally need to map back from $[-1,1]$ to $[\\lambda_\\min, \\lambda_\\max]$.\n",
    "This results into \n",
    "\n",
    "$$\\tau_i = \\frac{2}{\\lambda_\\max + \\lambda_\\min - (\\lambda_\\max - \\lambda_\\min)x_i}, \\quad x_i = \\cos \\frac{\\pi(2i + 1)}{2n}\\quad i=0,\\dots,n-1$$\n",
    "\n",
    "The convergence (we only give the result without the proof) is now given by\n",
    "\n",
    "$$\n",
    "   e_{k+1} \\leq C q^k e_0, \\quad q = \\frac{\\sqrt{\\mathrm{cond}(A)}-1}{\\sqrt{\\mathrm{cond}(A)}+1},\n",
    "$$\n",
    "\n",
    "which is better than in the Richardson iteration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaMAAAEcCAYAAABqCdtUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8W9d54P3fo42SLVIgtXhXJMhJ7CRObJBu0yWbBaZZOpnEIaWk+2caEeky77xNEyJKO5/2nbaRSSfpMtOkINPpTDt9G4lM+nYyadOSchYnbRqRdOo2jmObkCwvsbWQICXborbn/eNcgBcgQALkBS8BPt/PBx8I915cnAtCeHDOec45oqoYY4wxYVoTdgGMMcYYC0bGGGNCZ8HIGGNM6CwYGWOMCZ0FI2OMMaGzYGSMMSZ0FoyMMcaEzoKRMcaY0FkwMsYYE7p1YRegVmzbtk137doVdjGMMaamjI6OnlHV7QsdZ8GoTLt27WJkZCTsYhhjTE0RkSfLOc6a6YwxxoRuycFIRJqCKIgxxpjVa9HBSER2i8gTQEZEzorIe3z73isijwdSQmOMMXVvKTWj+4AE0AzsB94vIvcCqOrngT1LL54xxpjVYCnBaERVj6rqlKoOq+o+YKuI3OPtt4WSjDHGlGUpwShTuEFV+wERkbuWcF5jjDGrzJJqRiJyQEQeF5Fd2Y2qetT759RSCmaMMWb1WPQ4I1V9SETSuOa6E0X2tS61cMaYOvTwEZg5B3f/YtglMSvIkga9quoU8FCJfceXcu5aNjMzw8TEBOfOnePKlSthF8fUobVr19LY2EhLSwsNDQ1hF6cy3/lLuPiCBSOTx2ZgCNjMzAwnT56kubmZXbt2sX79ekQk7GKZOqKqXLp0ienpaU6ePMnOnTtrKyBNPQ3X3xF2KcwKUzQYicifAG0sPSNOgHFV3b/E89SMiYkJmpub2bZtW9hFMXVKRNiwYUPuMzYxMcENN9wQcqnKpOqC0SvfHnZJzApTNBip6geXuyD14ty5c9iEqma5NDU1ceLEidoJRi+ehcsXoOnmsEtiVhibmy5gV65cYf369WEXw6wS69evr61+yamn3P0WC0YmnwWjKrA+IrNcau6zNvW0u7dgZAqU6jO6D9gd0GukVfVgQOcyxtSyqWfc/ZZbwi2HWXFKZdMdAloCeo2JgM5jjKl1U0/Buk1wTVBfL6ZelEpgmMJmUDDGBG3qaddEV2vNi6bqrM/IGLN8ssHImAIVBSMRaRKRIyJyxbsdE5Ffr1bhTO0ZGxtDRObcWltb6e3tnXN8b28v7e3tZZ8/mUyG0mnf19eHiJDJzJkf2FTCgpEpodKaUS9wDNgHfBQ4DtzvLa73uqALZ2pXT08P4+PjjI+PMzo6SiKRIJVKsWfPnrwv9GPHjjE8PBxiSc2yuTwD55+zYGSKqnQ6oFFvmYg8ItIBDIpIe+GkqSuNiMRws0sARIBBVU2HWKS6FIlEiEajucexWIx9+/axe/duDhw4wMDAAEDu3qwC08+6ewtGpoiSNSMR+YyIvEdEmnybi04PpKqDwN24lV9Xuriq9nm3XiAZdoFWi0gkQjwet5rQamVjjMw85mumSwCfBya9NYs+A0RE5FCxg1U1A4xXoYxziEiHiMRL7IuJSJeIxLP3BYckRCSyDMU0RYyNjdHSYmm9q1IuGNkYIzPXfM10w6r6Vq9Zay/Qjlfz8ZrlhoFR3HpG3/EW2GsFPlvNAnvBpR/oLLIvCvSoartv24CIpH1NcT3AcRHJ1oiWrWb0/3zxuzzy7PRyvdyivOrGJn7r37060HNmMhlGRkZIJpOk02lSqVRuXzKZpLe3F9XZSvfY2BjJZJKRkREA2tra6OnpIRaLzTlmeHiYWCxGf39/3v50Op3bn8lkiMfjpFKpvKbDZDJJX18fk5OTJBIJhoeHGR0dJRKJ0Nvby6FD7nfXvn372LNnz5zrGhwczF0TMOc1CsuYSCTo6uqa8/pHjx6d91rqxrQXjJpuDLccZkVaqGaEqo6p6v2q+lZVXYMLOH3AVlxCw5iIXAEGcF/0VSEiURFJAVFKD6RNAKmCbSl/uVS1DzeoN+Hd7Gd6FSQSiVwmXXNzM+3t7YyNjdHT05P3hVxoeHiY1la3LmN/fz89PT1EIhEOHz6cd1xnZ2cuKSKdTtPZmf/bJBvw+vv7GRoaYmJiInfeQq2trYyMjJBIJHKBKJlM0tbWRn9/P62trSST+b9Zsq/Z0dHB6OgoAwMDRKPRXHLG4OAgra2tRKNRhoaG2L9/P4lEgkQivyU7k8kseC11Y+ppuHY7rN8UdknMSqSqS7oBW5Z6jkW85jiu76fY9mjBtoi7zNzjbt+/u3BLXCz4mq2trVqORx55pKzj6tXo6KgC2t3draOjozo6OqpdXV0K6OTk5Jzju7u71fv7qKpqNBrVeDxe8vzZ4wcGBuZsK3b+wnIVe17h60UikTnbenp68l5jYGBg3teMRCLa1dWVt21oaEgBHRoaWtK1FKqZz9xf3KuaelPYpTDLDNd6tuB37JIHvaqbrSF0Xj/QnFqTur6sbM0qDoz59vXhsgDrsE0kXHv27CEWixGLxejpcRXTI0eOzPucdDqda15bSDw+2xVYrAmtULbZK9uk5udvNhwbGyOTycypwUQi+d2M2dffu3cvvb29jI3lPla5psHCc8TjcSKRyJwMwkqvpWbZGCMzj3qagaEFZoNPEdlANSfwqOrY3MNNUCKRCN3d3QsGmWygKCfBoTA4FDM4OEhnZyetra00NzeXPM7fj5Qtg39bqdcfGhoCXN9Pa2trbvDufOeIRqNzAmI511LzsovqWfKCKaGegtGC/6O9oJP2suy6RKQbOLzQ88zSHTzoJm6fLyBlv7wnJpY+t257ezvJZJL29nb6+/s5evRoWc+rpAzxeJzR0VEmJyfp6elheHiY3t7e3DmK1cLS6fSCga4uXcjAxfPQdFPYJTErVD0Fo7Ko6qD6xhnNVyvyAtaIiIycPn16OYtZdyKRCAcPHqS3t7fklDrRaJRoNJpr1vOrZBqedDrN8PBwLlkiFouVHQBisRiRSCSv6Q5gfLz0qIVszS8Wi3Hs2DHa2tyY6sJzDA4O5hIWVh0bY2QWUOkMDCueiETmaaqriNen1AfQ1tZWdMCvKV93dzeHDh3Km4GhUCqVor29nc7OTvbv38/ExETu2Gyz2EKi0SiRSCRXC4tEIkUDXCkHDx4kmUySSCTo7OxkbGxszrx6fX199PT0kEgkck1vY2NjuYy8VCpFIpFgYmKC/fv35/rCOjo68vqIVg0bY2QWUE81o2ybSF6Hg2+Aq035swIcPHiQwcHBvA5/v2zTV7YGkU3tLqxlLKS/v5+JiQk6OztJJpMV1Ua6u7vp7u7myJEjJBIJzp49OyeY7du3j46ODlKpFJ2dnaRSKbq7u3Np611dXQwNDeVStVOpFD09Pat3+iOrGZkFiOrifvCLyBbggKp+ItgilfXa40BCVYeLbO/0N715A2FHVbV0D3YZ2traNDsIcz7f+973uP3225fyUsZUpCY+c0O/Bd/6NPzG87Cmnn4Dm4WIyKiqti103FI+FS24+ehWkmFmJ0HNinnbjTFhmXraJS9YIDIl1Oono4Xi2XNJ5k4TlMAmQzUmXDbGyCygZhIYvL6fg7jxQhGgR0TagSF1s4ajqhkRSYpID27dpexcddZfZEyYpp6C3W8KuxRmBauZYORlyC1Yw/H6i2wQqzErxZVLcO4HELFMOlNarTbTGWNqxfQzoFetmc7My4KRMaa6bIyRKYMFI2NMdWWecveRneGWw6xoFoyMMdU15QUjm5fOzMOCkTGmuqaegmt3wPqNYZfErGAWjIwx1ZV5yjLpzIKWGowkkFIYY+rX1FOWvGAWtOhgpKrHsZkNjDHzyS2qZ2ndZn5Lqhl5AckYY4p74QxcvmCZdGZB1mdkjKmeqZPu3prpzAIsGJmqGRsbo7Ozk+bmZkSE5uZm2tvbGR6enUQ9mUwisvSux6DOE/Zr1J3sGCNrpjMLsGBkqqKvr4/W1lbGxsY4ePAgAwMDHDx4kImJiYoXyjM1LDv7gmXTmQXUzESppnYMDw+TSCSIx+Nzlgrv7u4mkwlkVXhTC6aegg2NsLHYii/GzLKakQlcIpEgGo3OCURZkYh9Ma0amadcE501b5oFBBKMROSIiNzre7xXRJqCOLepLYODg6TTaZJJy/o3uJqRNdGZMgTVTJdS1aPZB6p6VETuAR4I6PyBEpEO/+Ps4nxV93cfhef+dVleatGuvwPeft+in37s2DEA9u3bV9HzxsbGSCaTDA8PE4vF6O/vJxaL5R0zODjIoUOHGBsbIxaL0dPTQzweL3qekZERWlpaSKVSxONxkskkvb29TE5O5tXMent7SSaTue2Dg4Mkk0nSabceYzweJ5VKEY1GAylrIpGgr68PVc07vrAcdWPqKbj57rBLYWpAUM10aRHZBSAi/yAix4D2gM4dKBHphlwAGsatHmsCkv0Sr/QLtbOzk0QiQSqVIp1O09mZv3r84OAgnZ2dxONxBgYGiEajtLe3514v68CBAyQSCXp6egByxyQSCcAlVvgdPnyYeDxOJBLJvW5HRwejo6O51yns41pKWbPlGBzM//1z+PBhYrFYfQWimfPw0qTVjEx5VHVRN+Bx4DPAPd7jO4F7gfcCWxZ73jJfuwOIl9gXA7qAePbety8CTC7mNVtbW7UcjzzySFnH1at4PK7uY1We7u5uBXRgYGDOtsnJydy2SCSi3d3dec+NxWK5bdnn+I2PjyugXV1dueOj0Whu/+TkpAKaSqVUVXVgYGDO61ajrJFIROPx+JxyZstRqRX7mXv+EdXfalJ9eGDhY03dAka0jO/YpdSMPo+rWewTkceBXqAHGFfVqSWcd14iEgf6S+yLAj2q2qeqw6raByS87QBtuFpch4jERaTbt88EINtcVWnGnL+5bc+ePXn7xsbGyGQy9Pb2IiK529jYGGNjpVeYj0ajxONxRkZGADh48CDpdDr3nGwtqaurK68Me/fupbe3t+S5l1rWrq6uvLFW2VpSpU2bK54tqmcqsJS56T6qqp9X1Q+q6suBBC4gfUxEnhCRzwRWSlygEZEUEAUmShyWAAoHsaRwQRLvuTFgWFWHgT6geMqXWZS773b9A0eOHKnoefM1T2Wb4kZHRxkfH8+7DQwMzHveaDSae35HRweRSITDhw8Drmmso2O2+zASieQyAJPJJK2trbS3z21tXmpZC5vq/E2FdSXjzb5gzXSmDIGldqvqcVXtV9V9qnorLjAFRlXTqprwajuldACFP2dHvO0AaSCtqhnvnBkgarWj4HR0dBCNRnN9NkHwJw9Eo9G820Jf4Ol0Ou/5+/bto6+vj0wmw9jYWC4wZMXjcUZHR5mcnKSnp4fh4WF6e8v/KJdT1uzjVCpVshx1YeopWLMONl8XdklMDag4GJWbsq3LPImqiEQoUmvKBh4v4KSLPNVGYAasVMd+1nxNa8VkO/aLpYvP1xyYyWQYHh7Oa1ZLJpNkMhn6+vqIRCJzsvGyIpEI3d3dxGKxXIZgkGVNJBIMDw/nmuv8NbS6MfW0W911zdqwS2JqQEXBSESOAJPZzDnf9nsCLNNitcBs8CkiqqppIOMFrmwAS3vbTUCyWWSDg4O0trbS19fH4OAgvb29tLe3s3fv3orP2d/fz/DwMJ2dnQwPDzM4OEh7e/ucL/3Ozk7S6TTDw8Ps3buXSCTCwYOzCZPRaJRYLMahQ4fm9NH09fWxZ88eent7c+UdGxsr2lS31LJmg8+hQ4fqMxCBt6iezdZtylROlkP2BhwA3ltk+xbgw0BTJedb7A0YpyCbDtcXpCWO1+zxuNpTDy7TrgcXpEq9TheumW9k586dZWWOrNjMphCMj49rR0eHRiIRBTQSiWhHR4eOj4/njimWBZdKpYpmtY2Ojuay9SKRiHZ1deWO8We6RaNRBTQej+e9VuH5R0dH87ZPTk5qd3d37vnRaDQvKy6osmbFYrE52XmLsWI/c5+8XfULibBLYUJGmdl04o4tj4h8BDfAdbrE/g+o6mfLPuEiicg4kFCXhJDdFgNGVXXOvCMiokC7//hKtbW1aTYraz7f+973uP322xf7MmYZ9Pb2kkqlGB8fD7sogViRn7nLF+F3d8CbkvAWG8q3monIqKq2LXRcRc10qno/MCAih0TkLYsuXRVlm+CMKSWVStVnwsBKMv00oNZMZ8pW0XRAIvIngOBmV0h6NY4xXFNWBjeotOo1oxKy/T4t+JISfMHJ+oVWsez4omPHjpFOp3Nji0yVZNcxsrRuU6ZK56YbV9UPZh94A1CzNwUq75kOiKpmRCSNC4h+LUBGLUlhVctkMhw4cACAgYGB+hvTs9LkxhhZzciUp9JglJep5vXBDAOIyF3APsKrGeGVpY38sUYxb7tZxWKxGJOTk2EXY/XInARZ41K7jSlDpeOMhkXkAyX2Lef/9Bbm1oAAkkDh4JaEt90Ys1wyJ10gWrs+7JKYGlFpAsNxXALDB/yDX0VkN65PpjXg8uWISEREekRkABeIekQk5V8OQt0Yo6R3XIc3Q3ePNdEZs8wyJ62JzlSk4vWM1E2C+tmCbcdFpJ0qJglkA00Zx40xd0qgZaWqiK1saZZBJUMzllXmJOx+Q9ilMDUkqMX1UN/ieqvZ2rVruXTpEhs2bAi7KGYVuHTpEmvXrrDpdq5cgnPP2mzdpiKBTZQKEPRM3bWosbGR6emiY4KNCdz09DSNjY1hFyPf9DOgV62ZzlRkScFIRJpE5E4RucdLbKizBVkq19LSwuTkJGfOnOHixYsrtxnF1CxV5eLFi5w5c4bJyUlaWlrCLlI+S+s2i1CymU5E7sRlrUWBPbikgahvW2E2m+DGGq1qDQ0N7Ny5k4mJCU6cOMGVK1fCLpKpQ2vXrqWxsZGdO3fS0NAQdnHyWTAyizBfn1Evs4NZs73xae92FBeQUrglGzLAcaCyFdXqVENDAzfccAM33HBD2EUxZvnZGCOzCCWDkaq+1UubjgNJLVhKXET+RFX7C7bZ4FJjVrvMSWi8EdZZEo8p37x9Rqo66E3/0+k12+XtLnL8R4MsnDGmBtkYI7MIZSUweMtCTK2QRfSMMStZ5qRNkGoqVvY4I2/2heNeDWlqoeONMavQlcsw/azVjEzFFjMDw3cARGRMRJpKLbRnjFmFpp8BvWLByFRs0eOMvOSFrd4YoyYAEbk3sJIZY2qPpXWbRVrSoFdVPa6qDwDt3sJ7/Qs9xxhTxywYmUUKZDogVf080EfxZR2MMatF5iQg0HRz2CUxNSawuem82bJtslRjVrPMSWi8wcYYmYoFOlEqbiE7Y8xqZWOMzCIFGoy89O+aISKpsMtgTF2ZsmBkFifomlHNEJE40BV2OYypG1cuwdQzFozMohQdZ+RlxrWx9Fm4BRhX1f1LPE/+Sd2ceRlVnTMXnojEcGVP4yZzTRceJyIRZid4NcYEYeppN8aoeVfYJTE1qGgw8uajW5G8Gk0/0FlkXxToUdV237YBEUmrqn9J9LiqDtrS4MYEKPOku29+WbjlMDWpZprpRCTq9fFEcbWaYhK4ZS38UkCP7zxxwGYXNyZok14wilgwMpWrmWCkqmlVTahq3zyHdQBjBdtGvO255jlVteY5Y4I2eQLWrLN1jMyilOozug/YHdBrpFX1YEDnKskLNHNqTaqaEZFsE14MaBGRNm93RES6gOGCZjxjTKUyT8KWm2FtxVNeGlNyotRDuOXFg1CqSS1oLeCCT4n9UVUd9G8QkdQCNS1jTLkmn7QmOrNopRIYpqi9ZSLKnorIq0V1ef/uBgaL1Yy8WlMXwM6dlq5qzLwmT8Bt7wi7FKZG1UyfUZBUNaOqvaoq3n3RJjpV7VPVNlVt2759+3IX05jaMXMeXjxjad1m0eouGHm1HmPMcsrN1m3NdGZx6ikYZWs3eX1dvuBkCQrGVEtujNGuUIthalfdBCMvcSHN3L6jFtxsDRaMjKmWyRPu3oKRWaS6CUaeYdxUQH4xbJCrMdU1+SSsvxau2Rp2SUyNqtVg1ELx7Lkkc6cJSnjbjTHVknnSTQNkU2yZRaqZ0Wle389B3MDWCNAjIu3AUHb8kDfANSkiPcAx79gea6IzpsomT1jyglmSmglGXp/QgjUcb8XZwimBjDHVouqa6Xa/KeySmBpWq810xpiV4sWzcOkFm63bLElgwUhE7hSRD4jInUGd0xhTA7KZdNZMZ5YgyJrRx4A+YFBEHheRewM8tzFmpbK0bhOAIPuMksAhVX0IQES2BHhuY8xKlR3wasuNmyUILBip6vGCx7U20aoxZjEmT8A126Bhc9glMTXMEhiMMUsz+aQlL5glW1TNyEtSKLrekao+sKQSGWNqS+ZJuKk17FKYGldRMBKR9wJHsg+LHKLA2qUWyhhTI65chqmn4dWWr2SWptKaUQ9wHJesYLMaGLPaTT0FVy9Dy+6wS2JqXKXBKAp0qOoXqlEYY0yNmfTylpotGJmlqTSBYZgSfUXGmFVowgtGVjMyS1RpzSgBHPHGEB0FMoUHqOqJAMpljKkFk8dhbQM03hh2SUyNW0wz3R6gd55jLIHBmNVi4rhL615jo0TM0iwmgQHgfmA84LIYY2rN5AnrLzKBqDQYxYAuVf1sNQpjjKkhqq5mtOvHwy6JqQOVBqMxYKIaBVkuIhLDLU0eAe4Gkrb4njGL8MJpb+kIqxmZpas0GB3AJTBMAqOqOl2FMlWNt1psm6r2eY/jwBCuH8wYUwnLpDMBqjQYHcXVKIYBZO5696qqVV89VkQ6gIyqDhfZl635pHEJF2nfcVHcgN0+7/EIEBWRiLeSrDGmXDbGyASo0sDRj5vyJzRebaYf6CyyLwr0qGq7b9uAiKRVNa2qYyLS7ntKGy6oWSAyplITxwGxSVJNICoNRr8XVtOcF2iSwCil+60SQKpgWwqXBdgJUNA/lMA1PRpjKjV5HJpugnUNYZfE1IFKBweMisivV6UkC/BqNolsf08JHbgkC78Rb3seEekCDqvqYIDFNGb1mDhu/UUmMJUGoweA91WjIEvlJSdEKag1ZZvgvJpV9tg4ri/JApExizV53JYaN4GpKBipagIQEfmyiKy0huIWmA0+RUQhl+AwkU1q8JIhjDGVmDnnUrutZmQCUul6RtlaRwxIh5VNV0JkoQO82tFRIOIrexooWkPymvK6AHbu3BlMKY2pB5Mn3L1l0pmA1Fw23VJ4yQvNFRzfh5cG3tbWVrPXbUzgbIyRCVhFwUhVkwAicg8uOy2Kq1kMrJTlxm3MkDHLwMYYmYBVPNWuiHwGN+h1P7DVux8Skc8FXLZKZVO289Zb8hIb/PuNMUs1cRw2NcOmBVvHjSlLRcFIRA7gxua0q2qLqrapagsuIO0LK+0bcokLaeb2HbXgBrZaMDImKJPHrVZkAlVpzSgBdKvqUf9GL0X6o8AHgyrYIg3jZlXwi3nbjTFBmTgOLdGFjzOmTJUGoxhzB5VmjeOlTy+DFopnzyWZO01QwttujAnClUsw9bQlL5hAVZpNN4ybzaBYssJ+SgeqJfP6fg7iAl4E6PHmmRvKDl5V1YyIJEWkBzjmHdtjTXTGBChzEvSKNdOZQFUajD4KjHhjdPqYnRk7AbwXaC/91KXx+oQWrOGo6hhVDIrGrHpnvUWet9rKKyY4laZ2j4nIPlwgSnibs6NHEyslvdsYU0UTXjBqsWBkglPxbAlek9igN7/bblztaERVp4IunDFmBTo7Dg1NcO22sEti6siip+4ptrCdMWYVmBh3mXRzpwMzZtEqHvQKICJNpW5BF7CuTD8LXzkEly+GXRJjFu/suPUXmcBVOuj1LhF5HJic52ZKGfkz+Np98M0/DLskxizO5Ysw9ZT1F5nALWai1K3A/bhxRaYS6a+4+6/fD69+D2y7NdzyGFOpyROgV61mZAJXaTCKAV2q+tlqFKauvTQJz4xC6y/Ad/8a/s//DT//RWt3N7XFMulMlVTaZzRMwUqqpkzHH3S/KF/7Pmj/L3DiQfjOX4ZdKmMqM+GNH7eakQlYpcEoiZv54D2WrFCh8QdgQyPc3AZ3/Rzs/FH4+9+YHUBoTC04Ow4bI3BNy8LHGlOBSoPRUWAPbmXUSRG5UnC7HHwR60T6K7D7DbB2PaxZA+/6I1izFv60HZ4eDbt0xpRnwjLpTHWsqpVeQzORdh2/r/+V2W3bXg6/OAT/6174nz8JHX8Gr3xbaEU0pixn07Dz9WGXwtShRa30aio07mXR7XlL/vate1xA+stO+Nz74Y0fgTd8GNZtWP4yGrOQSxe8tO6fCrskpg4tatCrqVD6K7DlFthaJJV78w74hS/BHfvgaz3Qfw/84OHlL6MxC5k8Aag105mqsGBUbVcuQ/rrEH1z6TTuhs1wbwre/zl44RT0vwW+/DGXDm7MSmFp3aaKLBhV2fOP/hPMTPHiLW9a+OBXvh1++VvwuvfDtz4NfxSDb/e7xcyMCVtu6Qhb4dUEb9ETpdYqEYniFggcww3i7fPWSqqKE9/+IttV+NEjV2gcfoDbr29iZ8s13NJyDbe0bOKGLZu4YctGtmxaj4i4lNl//9/gh7rg7z8Gf/thF5jefBBe816XgWdMGCbGYVMLbGoOuySmDonq6kqOE5FRVW31/h0B+lW1cKnyOdra2nRkZKTi15ucOMPJf/tHvnn1VXz32Wm+/9w5np58kQuXruYdt2n9Wq5rauC6po1cv2Uj1zdtZEdjA3e8+E+85tE/4prJR7m6/TbWvPkg3P4ulx5uzHL6Hz8Jly/AB2zCflM+7zu3baHjarJmJCIdQKbYMhYiEgPamF2FNp09ztuXm0HCW6Y8Xs2yNrdso/mN7+J1vm2qyunzMzw9+RLPTV3g2cxL/GDqAs9NX+D5qQuMPjnJqXMzXLx8FYgg/CbvWPNtPvT8AHsGfp7ja17GFyM/zZM72tm2ZRM7GjdyXVND3v2mDVaDMgGbSMOuN4RdClOnai4YecGjH5hTm/Ga4HpUtd23bUBE0qqaDU6FTXITIhLzlitfFiLCjsaN7GjcWPIYVWXqpUs8N32BU9MzPD99J38//TO4PJ18AAAfr0lEQVTsePJL/Nizf8b/NfFxjk/+Of/10rvpv/x6rhZ0/zVuXMd1Xu1qR6OrcW1vbGBH00auy943NXDNhpr7CJgwXHwRpp+xTDpTNTXzTeQFmiQwSun58RJAqmBbCujBBa+amcNERIhcs4HINRu47Xr/nlfC1f8E3/1rdn/9E3zq9H/j/u1/x6k7f5XHr3s7p164wqlz2QB2gVPnZhjJq2nl29ywjh1NswHLH7hywaxpI5sbauajYqohOyddiyUvmOoo+Q0jIlepfLYFVdWqfGt5NZsEgIiUGnzbwdxgNAIMef+eACIF+2smQOWsWQt3dMCr74VHv8jar93PDV/5NW5o/kN4w6/Dj79vzsDZbE3r1LmZvEDl7l3weuhkhuenLzBTJGhdu2EtO3xBytWu5gauzQ3rXCKGqS9nH3f3214ebjlM3ZovcNxPDU394yUjRCmoNXn9QtmaVZoiwWc5m+gCtWYNvOrfu4SG7/8dfL0X/vd/hK/1wo//Gtz1M7CuAcivab3iusaSp1RVpi9c5pQXrE6du8DzvuB1enqGh592QaswCQNmEzF2NG7MBasdebUsV9NqtKBVW8484e6LDdw2JgAlg1ENTv3TAi74lNgfVdVh/xegF6BqPzVIBG57hxun9MSwm8nhSx+CBz/pglLs53JBaeFTCVs2rWfLpvW8fIGgdW7GC1rTM75a1uy/v/vsNA88eooXL16Z8/yN69fkJVzs8CVgzPZ1baRpkwWtFeHs49B0E2y4NuySmDpVTx0Bhc1vpRwQkW5mxxkdKHWgiHQBXQA7d+5ccgGrTgRe3g63xt0URF/tceOUHvwUvOFDcNfPwvrSSROVvZTQtHE9TRvXc+uO0kEL4NyF2ebBbJNgtsZ16twFvveDab722AznZ+ZO+t6wbo2rYfkC1o6CzMHrmhpmx2mZ6jjzuNWKTFUtapyRiNyLaxKbQ1U/sdRClfH640DCn9rtpW2PquqcbyQRUaC9WCp4uRY7zihUqnD8a/DV++DkP0HjjS4oVVBTWk4vzFzO1aqen77AaV8tyx/IzhUJWhvWrcnLHMwmXuQeewEtco0FrYqpwn074bX74J2fDLs0psZUZZyRiOzGJQRkh2Ar4P+frUDVg9F8RCRSzRkVaoqImxNv95vg+Nfhq4dcTekbv+8SHe762RU1Q/i1DevY3bCO3dvmbwp66eKVvJrV89MzuT6u56cv8Pip83zziTNMXygStNau8VLc5wlcjQ00X7OBNWssaAFw/hTMTMO2V4RdElPHKm2mSwHHgbiqPiQiE8A9wAngiHcLi5d7Sgu+sUReYoN//+ojAtE3we43QvqrLih96UPwjT+AN33EzYW3dn3YpSzbpg1rednWa3nZ1vKCVmEG4anpCzx/7gLp0y/wrfQEUy/Nnftv3RopGJflayr09WltvXYVBK1sJp0105kqqjQYxfECkfc4DbmMtR7gPuCzAZavbF4Z0hRP3c54qeGrm4hbUyn6ZnjiKHzld1323Td+H978MW/uu/qZZqjcoHXh0hWvSfACz03N5AJYtqnwybMv8u0TE2ReLB60tm1ucP1XviDlHs/2cW29toG1tRq0zlhat6m+SoNRBtjie5zGTb3zHVwQiAVUrsUaxpXHn6odox4y5oIkAi+Pw617XUr4V34PvvAB+Man4C2/Abe9s/RyF3Vo4/q13sS118x73GzQmpmT+n7q3AxPTbzIyIkJJosErbVrhG2bN+QC1fYimYPXNTWwdfMKDFpnn4B1m6Dp5rBLYupYpcFoGHgf8Nfe4yPAIa+5roe5U+1USwvFs+eSwADQ59uW8G6mUDYl/BVvg+9+Ab7ycTj803BTG8R/yzXrmZxyg9bFy1c5fX428eK0r3/r1LkZnslc4DtPZThz/uKc564R2JqtaTXO9mcVzj24bfMG1q1dplrsmcfdNEB1VGs2K09F2XReAkNCVT/q2zYE7MUFog+o6hcCLyW5vp+DuCy+DlytbBgYUtVB33ExYD9wzDt2bClZdFk1mU1XqSuX4V/+X5d9N/0MRN8C8d+GG+8Mu2R16eLlq5w5nz9G6/S0N8j43Oz4rbMvzFD431QEtl6bTcKYDVTbfYkY1zU1sG1zA+uXGrT+6C64/rWw738u7TxmVSo3my6QJSREZLeqHl/yiVawVRGMsi5dgGOfhQc/4VabfU0H3POb0LI77JKtSpevXOXM+Yten9bswOLCpsKz52e4WjRobcg1C/qbBLf7ZsnYvrmBDeuKBK3LF+H3rndDAu75zeW5YFNXlnUJiXoPRKvO+o3wo78KsZ+Fb/4h/NOn4ZG/gbs/AG/8CFy7NewSrirr1q5xa1xt2chr5+m2uXzlKmdfuFhy7sFT52Z45NlpzhQJWgAt126YbRb0MghfLs/ybr3CcW5k/eSLbG9soGGdLU9igldpM13TQseo6vSSSrRCraqaUaHpH7h08If+AjY0wht+DX74lwKbzcEsrytXlbMvzOTNhDGb9j47uPj0+Rn28m36Nvw+75r5HR5Wt3xE8zXrZ5ck8dW4soOLdzS6fRvXW9AyVWqmK2cmb1Wty0/gqg5GWae+B8O/DY99GbbsdEkOr3nvqsq8W02uXFVeeuATbP7G7/L1947xgwvrXa3Ll0F4ykt/v1ykqrVl0/qicw8WzkloQau+VauZ7oMUD0ZtuDneuio8n6klO26HnzoM6a/BP/wGfP4X4Z//BH7iENxyd9ilMwFbu0bYfP44bL6ON95RelG9q1eViRcv5gLV6YK5B5+fniE9fp7T52e4dGXu10fTxnVzpm6y1YtXn4qCkar2ldjVLyJPALYM5GoQfRN0fQ3+5a/g6O/An8bhjn0u827LTWGXzgTpzOMLTgO0xhv4u21zA6+idEv+1atK5qVLeTNh+O+fn77At49PcPrcDBevzF2epDG3EGT+8iS2enF9CPKvNowbd3QwwHOalWrNWrde0qve7WZw+Mf/Ct/7osu6+tH/COs3hV1CE4Szj7u/cQDWrBFart1Ay7UbuP2G0sepKpkXL+Wlt89OnOtqWmMnJ3l+ep7Vixsb5tSwCgOZrV68sgT514hTYiZvU8caNsPe/+xmAh/6z242h4f+An7i43DbT1p/Ui174axL7V/maYBEhOZrN9B87QZuu770caVWL/YvU/KdpzKcOld8Icjs6sX5E+XO7eOy1YuXR6Wzds+XwCDA6JJLZGpT88tg35+72cH/LgmHf8YNmn3H/TanWa3KTZC6Mv9+la5efLqgH6tw9eJT0zO8dGnuQpDXbFibG5+VX7vKnyWjaaMFraWotGY031Lk44Q7a7dZCXa/ERIPwsifwgO/B5/+EfiRX3Hjkxo2h106U4nTj7r77bW9dIR/9eL5FoJUVc7PXJ6dumk6f+7B3OrF05WtXly4ppatXlxcpQkMtbYUuQnD2nXwwwl49XtcKvg3/wD+dQDedh/c/u+s6a5WnH7MTZC6pQZWOQ6AiNC4cT2NG9dz6475fzi5oHWh6OrFz0/Pv3pxdiHI65rmBi5/bWu1LQRpPXimejbvgHd/GmI/D1/6dTjys25J9HfcDy3WvbjinX7UNbHaBKlzbG5Yx+btm9mzff6glV292K2hNTsuK9tE+P3nzvHgY2dKrl68fXP+Wlp5GYRe4Gquk6BVMhiVM9tCMfU6A4NZgp0/DF1fhWP9runuj18Pb/ww/Nh/WpHLnxvPmcdg54+EXYqaVu7qxS9evJzLHCxcvfjUuQs8cfo8/zhefPXi9WslN+uFf7Xi65o2st0LZNc1rfzVi+erGWVYYLaFEmxkmplr7Tp4/S+5NOG/P+iy7h4+DO/8lBu3ZFaWmfMw9RRs//mwS7IqXLNhHbu2rWPXAkGr5OrF51yN6/iZF/jn44tYvdh3H9bqxfMFo2KzLURw6xYlyV+76FbgI0B3oKUz9afpBuj8H26M0pc+DH/+Ljdg9q2/C43XhV06k3XmMXe//bZwy2HyLGb1Yn8tK5uc8eTZFzlWYiFI/+rF2Znd333nTfzQ7pZqXZZ73VI7is22ICKfAbpUdc7S4iJyBrDpnE15bo3DL/+TGzD7jd+Hx/7ejVdq+w9uQK0J1+nvu3sLRjWp3IUgZy7nr178fEEG4dOTLzJ2cpI7b4mEF4xK2A8Mldg3BvwJK3wGBm/xvTZcLe9uIKmq6XBLtUqt3wRv+ZirGX3pQ/C3H3YDZt/xSZvrLmynH4U166HZ1rCqZw3r1nJz8zXc3Dx/0FoOi0mTKTUZaiduOfAVy1sttk1V+1S1F0hROria5bLtVvi5v4GO/w7nT7m57v7mV9y/TTjOPAZbb3V9fcYsg0qD0X3AW0XkyyLyHhG5U0TuEZHDuCBVaiLVBYlIh4jES+yLiUiXiMSz94t8mSiuvytrBIh6QcqEScQtR/Grx9zcdv/yOfijmGvCu3Qh7NKtPqcfrfnBrqa2VDrotdfLZ/8o8FZcgoPgkhmSqnr/YgrhBZd+XO2qcF8U6FHVdt+2ARFJV9q8pqpjItLu29QGZFQ1U+o5Zpk1NLpkhtgvuLnuhn8bRv47vCkJr32f/VJfDpcuwOQJuGPOf0djqqbiZjpV7VXVFtwX+T6gVVVbFhOIRCQqIilcjWWixGEJXHOaXwqX1VexggCWwK3DZFaabbfC+//KNd9ds9U12/3xD8G/HIYrczOATIDOPgF6dcGlI4wJUkUrvVaTiIwDCVUdLrK93R9EvGa1SVUV73EX86+lNFTkvF3AhKoOllM+W+k1RKrw/b+Fr3wcnv83aLwR7v5FaP0FuHZb2KWrP//2eRj8D/DBb8L1rwm7NKbGLXmlVxH5B+AeVV3n27bQsuPqP36pvKAzp9akqhkRQUSiqpqeZ9G/UueNA+nCAGVWKBG47Z3wirfDE0Pwrc/AA78DX+uBPXvhNffCK9/umvjM0p3+Psgal8BgzDKZL3AMkz+wFeaftbsaWsAFnxL7o0BF/UZeaveEqo55jzvKrR2ZkK1ZA6/4CXc79SiM/Tl896/hsb+DtRvg5h+C3W+AXT8ON7zOgtNinX4UmnfB+o1hl8SsIvMNeu0tsm25Z+0ONMvNS4Y4CkR8EwumgaLByGvK6wLYuXN1zFxcM3bcBm/7uEt2ePrbbpXZEw/CV+8jl1fTEnXNTNteAS173C/95l2uaa8OJpasmtOPwbZXhl0Ks8pUurjeLiCiqt/xbfswbvDoP6jqnwZauoB5/U7NFRzfh5eu3tbWtjI610y+NWtg5+vdDdzKpCf/GZ77V3juYfjBwy5QqW+lz3WbIHILbLkFttwETTd79zfBlpvd/YbwBwGG4spll8Dwip8IuyRmlam0fyeFq0n8EuT6lfYCDwF9IrJHVT8WbBFd35GlX5uybGqGV77N3bIuX3SpyhPjkDnp3Z6Eqadd0HqhyODaTS0uMG25xd1HboHITu/2Mvc69Vi7mjwOVy/BdqsZmeVVaTBqBzoAROQuIA50q+onRKQblyYdZDDK9ge14Ou/8g1StWl8zMLWbXADOEsN4rx0Ac49C9PPwtQzMP20C1RTT7sv5+Nfh4vn8p/T0OSWWm/eDS27XZNgyx7Yugcab6jdQPX8d929BSOzzBaT+ZYNAHFc43w2k20cl1AQGC9rLs3cvqMW3GBVC0Zm6dZv9IJJiY+vKlzIQOap2VrV5Al3O/0oPPZluHJx9vgNm11Q2vYK96W+/TZ3a4mu/Elgnxl1ySDXWUq3WV6VBqM0sM8LEPuBMd9iehVntpVpGDfAdsy3LeZtN6b6RFyz3KZmuOG1c/dfveJqURPjcHbc9bmceQxOfsstt561bqMLSte9Bq6/w53r+jtWVtbfM2OuTLbooVlmlQajjwJHcPO7Ca52lJUABoo9qUwtFM+eS3rn9Y8lSng3Y8K3Zq3XZPcy2HNP/r6Zcy4wnfoePP8InPquq0l95395B4hb2vumVrgx5mYrv+6OcKY9unoFnn0I7vrp5X9ts+pVOjfdoIjswauZqOqUb3cPLlCVzev7OYirVUWAHm/uuKHs2B+vqS4pIj3AMe/YHmuiMzWhodEFmptaZ7epwrnnXLbfs99xAeCJo/Avf+X2r7/GHf+yH3Njpm6+e3nG/Jx+FC69kF9WY5bJYn5+Ka7Z7ICI7FPVaRHZggtQhys6kcuQW3DskjdAdWyh44ypCSJuxdumG2ZTqFVdU9/T33ap6U99C77eC1+7D9Y2wM4fdrWuPXtdM9+axaz+soBnRt39TQvO3GJM4CodZ7QX+AdcKvdduKa1aVWd8mpM9wG/HHgpjal3Il76+C1uKQ2AlzKu3+nEg5D+qpvBfPi3XbbeK94Gr3wH7H5jcLWmZ0Zh45bSiRzGVFGlNaMeYFBV93vz1Pn14cYhWTAyJgibIvljps4955rzHvsyPHwERv/MpZjf9k549b0QfbNLY1+sp0ddv1U1al3GLKDSYBQDPl5i3yQBT99jjPFpvN4lF9z1025s1IkH4bv/Hzz6RdffdM02eN374K6fgR23V3buiy/AqUfgDR+qTtmNWUClP4HGcP1FxSSwdGtjlsf6jfDydnj3H8OHH4f3fw5e9qPwzyn49Ovhs+1uKYhy1376wcOgVyx5wYSm0prRfcARcbOMKrBFRO7ENd/FAfskG7Pc1jW4JTRe+XZ44Qw8fBiOfdatSdR4I/xwwq3/NN94pmzywo2x5SmzMQUqqhl56db7cLUgwdWUxnATpe7zT6BqjAnBtdvgR34FfnUU3n/YjWEa/i34g9fCg590456KeWbEzcPXeN3yltcYz2KWHR/0lh1vJX/Z8c8HXjpjzOKsWeMSH37+f8MHHoCb2+Dof4E/fJ2rNV25nH/8M6PWRGdCtei0GVV9SFU/r6oPZbd5TXbGmJXk5lb46QE48ABsvx2+9OuQegM8+iU3Oez5U27OPQtGJkTzBiMRaRKRwyJyVkQeF5FfnOe4I8BoVUppjFm6m1rhF/4P7PsLlz33uZ+CT93ubtn9xoRkoQSGMdzA1mFcWneft7bQJ7MHiMghoBuYws1dZ4xZqUTgVe+Cl78VnvwmTKTd7ONXLrpph4wJSclgJCIf8f65OzsHnYiMAr0i0oebtbsHl8jwUVW9v9qFNcYEZP1GuHUvbm1MY8I3X81oP3BfwWSoSdx0QBlcEOoFDhUcY4wxxlRkvmAUA0YKth3z7vtwtSELQsYYY5ZsoWy6jP+BL/ikLBAZY4wJis2IaIwxJnQLZdO9V0QKa0CKW3q86Bx1qvrZQEpmjDFm1VgoGJXKkCuVwq2ABSNjjDEVmS8Y2Qg4Y4wxy6JkMPJP82OMMcZUk6hq2GWoCSJyGnhykU/fBpwJsDi1YjVe92q8Zlid170arxkqv+6Xqer2hQ6yYLQMRGREVUstSli3VuN1r8ZrhtV53avxmqF6122p3cYYY0JnwcgYY0zoLBgtj76wCxCS1Xjdq/GaYXVe92q8ZqjSdVufkTHGmNBZzcgYY0zoLBgZY4wJ3ULTAZklEJEY0AakgSiQVtXhcEsVLBHpwF3bHu8+paqDBcfU7fsgIhGgR1UTBdvr8pq969oPnAW24v7e6YL9dXPdIhLHLacD7nrHVbWv4Jiav2bv/3GmWLnLub5A3gNVtVsVbt4fZKhg2wAQDbtsAV5jBxDzPY4A40DXankfgBQwsBr+9t7fO1V4/fV63bggFC/yHtTV5xuIA5OF11ru9QX1HlgzXfUkcF9UfincUu31IqqqY9kHqprBXZ//uuv2fRCRKNBSZFfdXbNXA+xXXw3Q+zUd9x1Wb9ed0IJf9+pq/e3+Y6jRaxaRqIikcMFkosRh5VxfMO9B2FG5Xm+4GkK0YFvEveXhly+A64sAo0CkYHsUN3t7tN7fB6DLuxXWjOrumr0vlp4i2/2/kOvqur3P95xf9/hqAfVyzd51FKsZLXh9Qb0HVjOqAu9X5JxfG+pqDtlf1DXNu5aodyuqnt8Hry/hSJHt9XrNXcCxwo3q9RfV6XUPA0NefwiQqw0OeP+ux2vOKef6gnwPLBhVRwvM/kGKqOkPaZaqNquvmc4Tx3WEpqnv9yFS4rrq9ZojQFpEukSkI3vv2193162qSSADjIpId/Z6dTaBoe6uuUA51xfYe2DZdNURCbsAIUoAh7x/1+X7ICIdWpAx6FN31+z7dRv1fREjIj0i0uJtq7vrBlDVVhEZwjVTjgF7fbvr8pp9yrm+wN4DqxmZwIhIFzChqr1hl6VavGaJUr8C61WpL5zD1EBH/VJ4n+kkLmkhiqsl1XqNZ0WyYFRF3hfXquD9B02oanuRffX0PuzTMsZP1Nk1pwvuAfCaaCP+L+d6um4v02xYVce8v/lu3HswWnBc3VxzMeVcXxDvgQWj6sj+p81L+/X9wdLUnx7ymzCgzt4HryN7ZIHD6uqaIa8/YL5+gbq67my51TegV1Uz3o+ttJfAUlfXXEQ51xfYe2B9RlWgqhkRSTO3eaOF2c79uuH9gkwWdmLW4fvQBuwRkf2+bTEgKiI9wDFVHayza87Kjawvtq9O/9bjJfal8BJY6uya85R7fUG9B1Yzqp5h3AfaL+Ztrxtem3qP5k8JE/c13dTN+6Cqfaqa9N+AIWDMe5xNaqiba/ZJMTstDpCrKfq/cOrputO4Ka6KieCSGaC+rrmYcq4vmPcg7MFW9XrDfWALp8gYooamCSnjGjtw409ivluc/Cli6vp9wDVPFg56rbtr9q5pvMg1ddTrdeMCcOF0QFF8g3/r5Zpx0wF1FNm+4PUF9R7YekZV5JtU8hjuQ5ztCK15XpvwZIndaVXd4zu27t4Hr+aXBPbh/jP24YLwmLe/nq95HFdrGCi8pnq7bq/mvwc3MSy4mmCxiVJr7pq9/8MHcWXuwNUGh3GBZdB33ILXF8R7YMHIGGNM6KzPyBhjTOgsGBljjAmdBSNjjDGhs2BkjDEmdBaMjDHGhM6CkTHGmNBZMDLGGBM6C0bGGGNCZ8HIBM6bm069W1eR/dFS+6pcrh4RWfGjvEUkJSKTC71H3uqjQ8tZtkrVQhn9vM/mkPfep8Iuz2pis3abakvipsoxZRCRAdwcfwe8TfPNenw3bi7AlawWypijbtLXdu9HS80E0XpgNSNTTX245RU6wi5IDenALccx6N3GSh2oqp2qKv5tXu1v2VciLfW6xcq40nnzrEH9zLxdEywYmWoawP2yPxh2QWqB78t8UcuaexNfdlOw1EO1hfW6VRTHW6Mp7IKsJhaMTLX1ADHfr01jVrp2rFa07CwYmaryptvP4IJSScWSC0Sky+tIjhQcNykiMREZ9faPZhf083U+j5dqHvSeO+SdZ9xbQrrwmI7C8xeWwft3yjtH4UqX871u9pxdvn09zK4smj1m3uZN/3vm9TVll/QY8J4/WXB8xdfkvacDvoSKIX9z3EKvW+LvWvJ9KCxPkWMLF/jr8Mqqxcq3SHG8/iLvPej2bvN+hs3SWDAyy+EQEA/gSyIrgmsCPAQkcOunDPhund5xAyWe349bOC3pPS78gu3wnjvsnStdeIx33ChuhcvUQk063jlHvXO1A4eBlC9j65C3Ha9ce/xrypThgO/5CdwaPK0BXFOi4Pwt3nWU9bqFyngf/LJ/5xT5f+fsubKPB73XzF7Xgj8M5ilfrr/I+4GxT1V7veuyvs9qCnuFQbvV3w33y1LxVsnEfTko3gqwuC8VBbp8z+lxH8e883R5x0UKj8O3AqdvW7dvW4e3LbbAa0T9ZfO2TeJbzdPbNprd5nu9oVLvQZH3ZNL/GiXep2xZ5qy4WeKcedcz3/ODuiZcv1DeayzwuoVlXPB9KChPR5FtkYK/cWS+Mlf42e3G1VAjBZ/PSJCvY7e5N6sZmapT9wu7D+gqpzmrTCO+f2ebt/zt/NmU6JYFypZd3bINcr+MI0C3r+lHmV1W3S9BGbzmsAjuF77/tYdxTZidxZ4XlCCvSWez+yqu5S7yffD/TcdL7DvqNaMF0S+ZreXtU9+KrqqaUUtoqCoLRma5ZNvbAxnoWuKLYbFfFmlmv1yz9624phn/Le/L0gtk5cies9jx/teuliVdk9cvM+D12ZRaar6ScpT9PswXALx92eDRA4zK0gfYxvGa/Yr1UZnqsUGvZlmoalpEBnFp3istUynK7Bdk7ouygmCzkOx5okDhuKEo+bW8alj0NXlf7lHcl322nKOln1FWOQJ7H7xaVatX4+4CekSkW10/T0V8geeQqmZEpBvXJ7Wn0nOZylnNyCynQ7hmmv1lHl+yIzwo3pdYHC9Aes1QRbP/ltDEmP2SzWsC8zrzs530QZgotnGx1+QlCMRxg3D7vPMUC2ZFX7eIqr0PXjNaLy7I3b3I0xSOL0rjNfN6WXWWwFBFVjMyy0ZVx0RkDNdJXGgc3FxmuMAQJ6AmvUJeOnKS2V/8GVygzDqAS1POZnJFcF+gacrsJ/LzfmUncFljLbgMsuxrD3q/7pfMex2ARPZeVbPNWBVfk1ebzeBqG+Dep2SR4+Z73cLjAnsfvJTwpHc92Wa+mPc4G0x7cDWdkjNZ+BQbX5QNvnl9SCZ4VjMyy23Ol5nnCO5XbQ/uF/JWFvHFX6bDuHEkQ7hf9a3+vgl1KdWtuC/sIVwqeHqesi/I+yJrZzY9OYGrcQSdvNCLC+Q9+JrClnBNB3C1gwFm/zZlv26hgN+HI7j+nYTvXL2+oBHBZdy1lXm+FnzJFd57lvZ+IJVb+zOLJKorfhJjY4xZFBEZBzrLrBmZEFnNyBhT1ywQ1QYLRsaYuuQlZyy6adUsL2umM8YYEzqrGRljjAmdBSNjjDGhs2BkjDEmdBaMjDHGhM6CkTHGmNBZMDLGGBM6C0bGGGNC9/8Di2aZaow3Ru4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "niters = 64\n",
    "roots = [np.cos((np.pi * (2 * i + 1)) / (2 * niters)) for i in range(niters)]\n",
    "taus = [(lam_max + lam_min - (lam_min - lam_max) * r) / 2 for r in roots]\n",
    "x = np.zeros(n)\n",
    "r = A.dot(x) - rhs\n",
    "res_cheb = [np.linalg.norm(r)]\n",
    "\n",
    "# Implementation may be non-optimal if number of iterations is not power of two\n",
    "def good_shuffle(idx):\n",
    "    if len(idx) == 1:\n",
    "        return idx\n",
    "    else:\n",
    "        new_len = int(np.ceil((len(idx) / 2)))\n",
    "        new_idx = good_shuffle(idx[:new_len])\n",
    "        res_perm = []\n",
    "        perm_count = 0\n",
    "        for i in new_idx:\n",
    "            res_perm.append(i)\n",
    "            perm_count += 1\n",
    "            if perm_count == len(idx):\n",
    "                break\n",
    "            res_perm.append(len(idx) + 1 - i)\n",
    "            perm_count += 1\n",
    "            if perm_count == len(idx):\n",
    "                break\n",
    "        return res_perm\n",
    "\n",
    "# good_perm = good_shuffle([i for i in range(1, niters+1)])\n",
    "# good_perm = [i for i in range(niters, 0, -1)]\n",
    "good_perm = [i for i in range(niters)]\n",
    "# good_perm = np.random.permutation([i for i in range(1, niters+1)])\n",
    "\n",
    "for i in range(niters):\n",
    "    x = x - 1.0/taus[good_perm[i] - 1] * r\n",
    "    r = A.dot(x) - rhs\n",
    "    res_cheb.append(np.linalg.norm(r))\n",
    "    \n",
    "plt.semilogy(res_richardson, label=\"Richardson\")\n",
    "plt.semilogy(res_cheb, label=\"Chebyshev\")\n",
    "plt.legend(fontsize=20)\n",
    "plt.xlabel(\"Number of iterations, $k$\", fontsize=20)\n",
    "plt.ylabel(\"Residual norm, $\\|Ax_k - b\\|_2$\", fontsize=20)\n",
    "plt.xticks(fontsize=20)\n",
    "_ = plt.yticks(fontsize=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "### What happened with great Chebyshev iterations?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- Permutation of roots of Chebyshev polynomial has crucial effect on convergence\n",
    "- On the optimal permutation you can read in paper (V. Lebedev, S. Finogenov 1971)  ([ru](http://www.mathnet.ru/links/42ab1a270ad179d43359ef8fab53422c/zvmmf6870.pdf), [en](https://www.sciencedirect.com/science/article/abs/pii/0041555371901698?via%3Dihub))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### [Chebfun](http://www.chebfun.org/) project\n",
    "\n",
    "- Opensource project for numerical computing ([Python](https://github.com/chebpy/chebpy) and [Matlab](https://github.com/chebfun/chebfun) interfaces)\n",
    "- It is based on numerical algorithms working with piecewise polynomial interpolants and Chebyshev polynomials\n",
    "- This project was initiated by Nick Trefethen and his student Zachary Battles in 2002, see [paper on chebfun project in SISC](http://www.chebfun.org/publications/chebfun_paper.pdf)\n",
    "- Chebfun toolbox focuses mostly on the following problems\n",
    "    - Approximation\n",
    "    - Quadrature\n",
    "    - ODE\n",
    "    - PDE\n",
    "    - Rootfinding\n",
    "    - 1D global optimization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Beyond Chebyshev\n",
    "\n",
    "- We have made an important assumption about the spectrum: it is contained within an interval over the real line (and we need to know the bounds)\n",
    "\n",
    "- If the spectrum is contained within **two intervals**, and we know the bounds, we can also put the optimization problem for the **optimal polynomial**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Spectrum of the matrix contained in multiple segments\n",
    "\n",
    "- For the case of **two segments** the best polynomial is given by **[Zolotarev polynomials](https://www.encyclopediaofmath.org/index.php/Zolotarev_polynomials)** (expressed in terms of elliptic functions). Original paper was published in 1877, see details [here](http://gpntb.dlibrary.org/ru/nodes/4245-zolotarev-e-prilozhenie-ellipticheskih-funktsiy-k-voprosam-o-funktsiyah-naimenee-i-naibolee-otklonyayuschihsya-ot-nulya-spb-1877) \n",
    "\n",
    "- For the case of **more than two segments** the best polynomial can be expressed in terms of **hyperelliptic functions**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## How can we make it better\n",
    "\n",
    "- The implementation of the Chebyshev acceleration requires the knowledge of the spectrum.\n",
    "\n",
    "- It only stores the **previous vector** $x_k$ and computes the new correction vector\n",
    "\n",
    "$$r_k = A x_k - f.$$\n",
    "\n",
    "- It belongs to the class of **two-term** iterative methods, i.e. it approximates $x_{k+1}$ using 2 vectors: $x_k$ and $r_k$.\n",
    "\n",
    "- It appears that if we **store more vectors**, then we can go without the spectrum estimation (and better convergence in practice)!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Crucial point: Krylov subspace\n",
    "\n",
    "The Chebyshev method produces the approximation of the form\n",
    "\n",
    "$$x_{k+1} = x_0 + p(A) r_0,$$\n",
    "\n",
    "i.e. it lies in the **Krylov subspace** of the matrix  which is defined as\n",
    "\n",
    "$$\n",
    "   \\mathcal{K}_k(A, r_0) = \\mathrm{Span}(r_0, Ar_0, A^2 r_0, \\ldots, A^{k-1}r_0 )\n",
    "$$\n",
    "\n",
    "The most natural approach then is to find the vector in this **linear subspace** that minimizes \n",
    "certain **norm of the error**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Idea of Krylov methods\n",
    "\n",
    "The idea is to minimize given functional:\n",
    "- Energy norm of error for systems with hermitian positive-definite matrices (CG method).\n",
    "- Residual norm for systems with general matrices (MINRES and GMRES methods).\n",
    "- Rayleigh quotient for eigenvalue problems (Lanczos method).\n",
    "\n",
    "To make methods practical one has to \n",
    "1. Orthogonalize vectors $A^i r_0$ of the Krylov subspace for stability (Lanczos process).\n",
    "2. Derive recurrent formulas to decrease complexity.\n",
    "\n",
    "We will consider these methods in details on the next lecture."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Take home message\n",
    "\n",
    "- Main idea of iterative methods\n",
    "- Richardson iteration: hermitian and non-hermitian case\n",
    "- Chebyshev acceleration\n",
    "- Definition of Krylov subspace"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Questions?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "from IPython.core.display import HTML\n",
    "def css_styling():\n",
    "    styles = open(\"./styles/custom.css\", \"r\").read()\n",
    "    return HTML(styles)\n",
    "css_styling()"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}