seminars/seminar9/seminar9.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Семинар 9\n",
    "\n",
    "##  Сдвиги для ускорения методов вычисления собственных значений"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Степенной метод: напоминание\n",
    "\n",
    "- Задача поиска максимального по модулю собственного значения\n",
    "- Сходимость линейная, коэффициент $\\frac{|\\lambda_2|}{|\\lambda_1|}$\n",
    "- Может не сходиться. В каком случае?\n",
    "- Сложность пропорциональна сложности умножения матрицы на вектор"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Сдвиги для ускорения сходимости степенного метода\n",
    "\n",
    "- Степенной метод \n",
    "\n",
    "\\begin{align*}\n",
    "& x_{k+1} = A x_k\\\\\n",
    "& x_{k+1} := \\frac{x_{k+1}}{\\|x_{k+1}\\|}\n",
    "\\end{align*}\n",
    "\n",
    "- Метод обратной итерации со сдвигами\n",
    "\n",
    "\\begin{align*}\n",
    "& x_{k+1} = (A - \\sigma I)^{-1} x_k\\\\\n",
    "& x_{k+1} := \\frac{x_{k+1}}{\\|x_{k+1}\\|}\n",
    "\\end{align*}\n",
    "\n",
    "- Куда сойдётся?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import scipy.linalg as splin\n",
    "\n",
    "def shift_inverse_iteration(A, x0, num_iter, shift=0.):\n",
    "    x = x0.copy()\n",
    "    conv = [x]\n",
    "    n = A.shape[0]\n",
    "    for i in range(num_iter):\n",
    "        x = np.linalg.solve(A - shift * np.eye(n), x)\n",
    "        x = x / np.linalg.norm(x)\n",
    "        conv.append(x)\n",
    "    return x, conv\n",
    "\n",
    "def shift_inverse_iteration_fast(A, x0, num_iter, shift=0.):\n",
    "    x = x0.copy()\n",
    "    conv = [x]\n",
    "    n = A.shape[0]\n",
    "    lu_factor = splin.lu_factor(A - shift * np.eye(n))\n",
    "    for i in range(num_iter):\n",
    "        x = splin.lu_solve(lu_factor, x)\n",
    "        x = x / np.linalg.norm(x)\n",
    "        conv.append(x)\n",
    "    return x, conv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.01811231 0.04519876 0.04519876 ... 7.95480124 7.95480124 7.98188769]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARAAAAD8CAYAAAC/+/tYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAYhUlEQVR4nO3dfXBV9Z3H8fdXAoj40ICRpkDX6FAVt4o1E7F2qlWoWhlhttORsrbpalfHsbbpbuvItDNuZ8cdttvdssrSig/duD4wjHXXjLPaAj4w1goNbVbBiDyklQCaiE9UWwX57h/3BE/CzU3OOffhnHs/r5nMvfeXc3K/CfHj+Z3vyfmZuyMiEscRlS5ARLJLASIisSlARCQ2BYiIxKYAEZHYFCAiEltqA8TMLjGzLWa2zcxuqnAt083sCTPrNrPNZvatYHySma02s63BY31on8VB7VvM7OIK1DzGzH5nZo+ktVYz+4iZPWhmLwY/23PTWGfw3t8O/u03mdkDZnZkWmo1s7vNrM/MNoXGItdmZmeb2fPB5241Mxvxzd09dR/AGGA7cBIwDvg/YGYF62kEPhU8PwZ4CZgJ/BC4KRi/Cfjn4PnMoObxQFPwvYwpc81/B9wPPBK8Tl2tQDvw9eD5OOAjKa1zKtADTAherwK+lpZagc8CnwI2hcYi1wZsAM4FDHgUuHTE9y7nL3WEH8i5wC9CrxcDiytdV6ieh4G5wBagMRhrBLbkqxf4BXBuGeubBqwFLgwFSKpqBY4N/qO0IeOpqjN4r6nATmASUAc8Anw+TbUCJw4JkEi1Bdu8GBr/MnD7SO+b1inMwD/YgN5grOLM7ETgLGA9MMXd9wAEjycEm1W6/qXAjcDB0Fjaaj0J6Ad+Fky17jSziSmsE3ffBfwIeBnYA7zl7r9MY60hUWubGjwfOl5QWgMk39yr4tfcm9nRwM+BNnd/u9CmecbKUr+ZzQP63H3jaHfJM1aOWuvIHXb/xN3PAt4hd6g9nEr+TOuB+eQO+T8GTDSzKwvtkmes4r+/geFqi1VzWgOkF5geej0N2F2hWgAws7HkwuM+d38oGH7VzBqDzzcCfcF4Jes/D7jczH4PrAQuNLN7U1hrL9Dr7uuD1w+SC5S01QkwB+hx93533w88BHw6pbUOiFpbb/B86HhBaQ2Q3wAzzKzJzMYBC4GOShUTnI2+C+h2938LfaoDaA2et5I7NzIwvtDMxptZEzCD3AmqknP3xe4+zd1PJPdze9zdr0xbre7+CrDTzE4Jhi4CXkhbnYGXgdlmdlTwu3AR0J3SWgdEqi2Y5uwzs9nB9/jV0D7DK8dJqJgnhb5ArtuxHfhehWv5DLnDueeAruDjC8BkcicrtwaPk0L7fC+ofQujOJtdorov4MOTqKmrFZgFdAY/1/8B6tNYZ/DePwBeBDYB/0Wui5GKWoEHyJ2b2U/uSOLqOLUBzcH3tx1YxpAT3Pk+LNhRRCSytE5hRCQDFCAiEpsCRERiU4CISGxlD5A0/ZGciCRT1gAxszHAfwCXkvujni+b2cwR9rmmHLUllZU6ITu1ZqVOyE6txa6z3EcgLcA2d9/h7u+Tu1Jy/gj7ZOIfhuzUCdmpNSt1QnZqzXSApOGPjESkSMp6IZmZfQm42N2/Hrz+CtDi7jcM2e4agqQ8YsKxZ9cddwJ1RxgnNxzNuLpc5nXveZsDB3O1f3zSURw3YWzZvo98+vv7aWhoqGgNo5WVWrNSJ2Sn1nx1bty48TV3j1V8XVGqGr1R/ZGRu68AVgCMb5zhja1LATjy6LFs+P7nAZjzr0+yrf8dAD4A/mHRLOadoYMZkajM7A9x9y33FCbRH8n1/XE/ly59EoB/+qtPUh8cdRwEbri/i0ee21X0gkVkeGUNEHc/AHyD3F2QuoFV7r650D6nTDlm0I0Kul95h0uXPklL02Q6bvgM9UflQsRRiIiUW9mvA3H3/3X3T7j7ye5+y0jbj6s7gtsWzRo01v3KO5xzy2oAbv/K2RwZnBdRiIiUVyauRJ13xlSWDQmRV/e9zxeXP01L02TuubplUIh8UyEiUhapD5DNu9+m/Zkdh0IkPJ3p++N+FixbdyhEdE5EpLxSHyAH3bm5o5uv3b2eM6fVs+7Gz1EXqrqrdx/n3LKaxuMm6JyISJmlPkAGPPnSa3z/v59n+qSjuP9vZw/63MB0Zvqkow47J6LpjEjppD5APnbckYeeP7X1NZaueZGWpsl5pzMD3RlNZ0TKI/UBMvno8fzg8tMOvV66ZjtL17zIvDOmsu7Gz6nFK1JBqQ8QgNZPn0TbnJMPvV66Zjtfuzu3GoBavCKVk4kAAWibc+qgEBk4JzJci/eKn/5KLV6REstMgEAuRMLTmYFzIvlavLvffl/nRERKLFMBArnpzNBzIuEW79jQdzQwnVGLV6Q0MhcgcPg5kXCL9z61eEXKJpMBAsNPZ9TiFSmfzAYI5J/OqMUrUj6ZDhDI3+Jtf2YH0ycdlbfF23r3s3mnM20ru9jQs7eMlYtkX+YDBA5v8d7c0T2oOxP21Et7ueL2w1u8Bw7Cojue1ZGISARVESBw+DmR8HRmaIis73lz0DmRCWM/DBFNZ0RGr2oCBIa/YvXMafWHhUi4xdt+VYuuWBWJoaoCBKJfsaqbEonEV3UBAtGuWFWLVyS+qgwQUItXpByqNkCgeC1eTWdE8qvqAIFkLV5NZ0QKq/oAgfgtXk1nRAqriQCBeC1e0E2JRAqpmQABtXhFiq2mAgTU4hUpppoLEIje4tVNiUTyq8kAgWgtXt2USCS/mg0QGLnFO3Q6oxavyGA1HSBQuMU7dDqjFq/IYDUfIKB1Z0TiUoAE1OIViU4BEqIWr0g0CpAhRlp3Ri1ekQ/FDhAzm25mT5hZt5ltNrNvBeOTzGy1mW0NHutD+yw2s21mtsXMLi7GN1AKhdadUYtX5ENJjkAOAH/v7qcBs4HrzWwmcBOw1t1nAGuD1wSfWwicDlwCLDezMUmKL6Wo05kFy9ZpOiM1J3aAuPsed/9t8Hwf0A1MBeYD7cFm7cCC4Pl8YKW7v+fuPcA2oCXu+5fDSFes1oV+el29+9TilZpTlHMgZnYicBawHpji7nsgFzLACcFmU4Gdod16g7F8X+8aM+s0s87+/v5ilBhboRbv/UOW0VSLV2pN4gAxs6OBnwNt7v52oU3zjHm+Dd19hbs3u3tzQ0ND0hITG67FO7CMZphavFJLEgWImY0lFx73uftDwfCrZtYYfL4R6AvGe4Hpod2nAbuTvH85qcUrcrgkXRgD7gK63f3fQp/qAFqD563Aw6HxhWY23syagBnAhrjvXwlq8YoMluQI5DzgK8CFZtYVfHwBWALMNbOtwNzgNe6+GVgFvAA8Blzv7h8kqr4CorZ4r/jpr9Tilapl7nlPQ6RGc3Ozd3Z2VrqMw7Q/s4ObO7oPvW6bczJtc07lked2ccP9XYNO7pz20Yk82nYBG3r2cu09G3njT/uB3Emh2xbNYt4Zec8li5SFmW109+Y4++pK1JhGavGODf1kte6MVCsFSAKFbkp0X54Wr9adkWqjAElouJsS5WvxPvXSXnVnpKooQIogyrozms5INVGAFInWnZFapAApoqg3JfrrO36tK1Yl0xQgRRblitWX3/izzolIpilASmCkFu8RoRTRFauSZQqQEinU4r31y7MGbaubEklWKUBKKOq6M5rOSNYoQEosyrozavFK1ihAykDrzki1UoCUiVq8Uo0UIGWkFq9UGwVImY10UyK1eCVLFCAVUOimRGrxSpYoQCpE91iVaqAAqaCRrlhVi1fSTgFSYWrxSpYpQFJALV7JKgVISqjFK1mkAEkRtXglaxQgKaMWr2SJAiSF1OKVrFCApJRavJIFCpAUi9riXbBsXd7pjEJESkUBknJRWrxdvftYsGydWrxSNgqQDBjpnMiY0LZdvft0TkTKRgGSEYVavNuXXHbYWrxq8Uo5KEAyZLgWL8CPF84atK1avFIOCpCMUYtX0kQBkkFq8UpaKEAyqtC6M/lavNfd+xtNZ6ToEgeImY0xs9+Z2SPB60lmttrMtgaP9aFtF5vZNjPbYmYXJ33vWjfSujNhj27q03RGiq4YRyDfArpDr28C1rr7DGBt8BozmwksBE4HLgGWm9kYJJFC684MDRFNZ6TYEgWImU0DLgPuDA3PB9qD5+3AgtD4Snd/z917gG1AS5L3l5zhrlg9c1o9q66dPWhb3ZRIiinpEchS4EZyR8IDprj7HoDg8YRgfCqwM7RdbzAmRTDcFastTZPz3pToi8uf1hWrkljsADGzeUCfu28c7S55xnyYr32NmXWaWWd/f3/cEmuOWrxSbkmOQM4DLjez3wMrgQvN7F7gVTNrBAge+4Lte4Hpof2nAbvzfWF3X+Huze7e3NDQkKDE2hO1xasrViWJ2AHi7ovdfZq7n0ju5Ojj7n4l0AG0Bpu1Ag8HzzuAhWY23syagBnAhtiVy7CitHgH7rGqFq/EUYrrQJYAc81sKzA3eI27bwZWAS8AjwHXu/sHJXh/YeQWr+6xKsVg7nlPQ6RGc3Ozd3Z2VrqMzGp/Zgc3d3zYZW+bczJtc05l5+vvcv6/PMHB0D//aR+dyKNtF7Dz9Xe5fNnTvPHufiB38uq2RbOYd4bOeVcjM9vo7s1x9tWVqFWu0E2Jht5jVS1eiUoBUgOirjujFq+MlgKkRqjFK6WgAKkhI607oxavRKUAqTGF1p1Ri1eiUoDUIC2jKcWiAKlRI12xOnQZTf0Vr+SjAKlhavFKUgqQGqcWryShABG1eCU2BYgAavFKPAoQOUQtXolKASKDqMUrUShA5DBq8cpoKUAkr0I3JcrX4m29+9m805m2lV1s6NlbxsqlnBQgMqwo68489dLeQdOZgRA5cBAW3fGsjkSqlAJECoq77sw9V7cwYeyHIaLpTHVSgMiICq07ky9EFixbR0vTZNqvatEVq1VOASKjEuWK1a7efYdCRFesVjcFiIzaSC3e8DqlXb371OKtAQoQiaRQi3f7kssYG/qNUou3+ilAJLLhWrwAP144a9C2hVq8ms5knwJEYknS4tV0pnooQCS2uC1eTWeqhwJEEona4j3/h48DuilRtVCASGJRWrx/eP1PuilRFVGASFGM1OIN/6LppkTVQwEiRVOoxfuUbkpUlRQgUlSF/oo3302Jvrj8abV4M0wBIkU3UotX91itHgoQKYlCLd5891hVizebFCBSMoXWnRk6nVGLN5sUIFJSavFWNwWIlJxavNUrUYCY2UfM7EEze9HMus3sXDObZGarzWxr8Fgf2n6xmW0zsy1mdnHy8iUrCq07oxZvdiU9Avl34DF3PxU4E+gGbgLWuvsMYG3wGjObCSwETgcuAZab2Zi8X1WqUtR1Z9TiTb/YAWJmxwKfBe4CcPf33f1NYD7QHmzWDiwIns8HVrr7e+7eA2wDWuK+v2STltGsLkmOQE4C+oGfmdnvzOxOM5sITHH3PQDB4wnB9lOBnaH9e4Oxw5jZNWbWaWad/f39CUqUNBpp3Rm1eLMjSYDUAZ8CfuLuZwHvEExXhmF5xjzfhu6+wt2b3b25oaEhQYmSVmrxVockAdIL9Lr7+uD1g+QC5VUzawQIHvtC208P7T8N2J3g/SXj1OLNvtgB4u6vADvN7JRg6CLgBaADaA3GWoGHg+cdwEIzG29mTcAMYEPc95fqoBZvtpl73lnE6HY2mwXcCYwDdgB/Qy6UVgEfB14GvuTurwfbfw+4CjgAtLn7oyO9R3Nzs3d2dsauUbKh/Zkd3NzRfej1BZ84nn9c8EkAPvvDJwbNdaccM44HrzsPgMuXPc0b7+4HcnPk2xbNYt4ZeU+tyTDMbKO7N8faN0mAlIMCpHYsXfMiS9dsP/T6/BnH0371OTzy3C6+cX/XoG1POHosG77/eTb07OWrd23gzwcOArn/e92qEIkkSYDoSlRJDbV4s0cBIqmiFm+2KEAkdaLclKj7lXe47t7f6IrVClGASCpFWXfm0U19ms5UiAJEUqvQTYlWXTtb05kUUIBIqg13xWrjcRPyTmfOuWU1oCtWy0UBIqkX5YrVgb/i1RWr5aEAkUxQizedFCCSGVFbvAuWrdM5kRJTgEimRGnxdvXuY8GydWrxlpACRDJnpBZv+DZ3Xb37NJ0pIQWIZFKhFu/2JZcxNvSbrRZv6ShAJLOGa/HufP1dfrxw1qBt1eItDQWIZJpavJWlAJHMU4u3chQgUhUKrTuTr8V7/g8f17ozRaAAkaoRZd2ZgXusqsWbjAJEqorusVpeChCpOoWuWM23jKZavPEpQKQqRV13Ri3eeBQgUrXU4i09BYhUNbV4S0sBIlVPLd7SUYBITVCLtzQUIFIz1OItPgWI1BS1eItLASI1Ry3e4lGASE1Si7c4FCBSs9TiTU4BIjVNLd5kFCBS89TijU8BIoJavHEpQEQCavFGlyhAzOzbZrbZzDaZ2QNmdqSZTTKz1Wa2NXisD22/2My2mdkWM7s4efkixRVl3ZnuV96h9e5na3o6EztAzGwq8E2g2d3/EhgDLARuAta6+wxgbfAaM5sZfP504BJguZmNyfe1RSpppHVnwp56aW9NT2eSTmHqgAlmVgccBewG5gPtwefbgQXB8/nASnd/z917gG1AS8L3FymJQuvODA2RWp7OxA4Qd98F/Ah4GdgDvOXuvwSmuPueYJs9wAnBLlOBnaEv0RuMiaTScFesnjmtPm+InLdkDVBbV6wmmcLUkzuqaAI+Bkw0sysL7ZJnzIf52teYWaeZdfb398ctUSSxKFes7nrzvZq7YjXJFGYO0OPu/e6+H3gI+DTwqpk1AgSPfcH2vcD00P7TyE15DuPuK9y92d2bGxoaEpQoklyhFu+qa2cPWkaz1lq8SQLkZWC2mR1lZgZcBHQDHUBrsE0r8HDwvANYaGbjzawJmAFsSPD+ImUzXIu3pWkyj3/n8BbvObesrokrVpOcA1kPPAj8Fng++ForgCXAXDPbCswNXuPum4FVwAvAY8D17v5BoupFyihKi3fgD/CqvcVr7nlPQ6RGc3Ozd3Z2VroMkUOWrnmRpWu2H3rdNudk2uacyiPP7eKG+7sGndg77aMTebTtAjb07OXaezbyxp/2A7kTgrctmsW8MyrfRzCzje7eHGdfXYkqElGhFm++P8BbsGxd1bZ4FSAiMUS5KVFX776qvSmRAkQkJt2USAEikshILd660HymGlu8ChCRhIa7KVHjcRN44rvV3eJVgIgUQZSbElVTi1cBIlIktXiPVQWISBEVuilRvhbvFbf/KtMtXgWISJFFafGu73kz0y1eBYhICdRKi1cBIlIitXBORAEiUkJR153JWotXASJSYtXc4lWAiJRBtU5nFCAiZVKNLV4FiEgZRbkp0fqeN4ddd6ZtZRcbevaWsfL8FCAiZZZk3ZmBEDlwEBbd8WzFj0QUICIVEHfdmXuubmHC2A9DpNLTGQWISIVEXXdmoMXbflVLaq5YVYCIVFDWr1hVgIhUWJZbvAoQkRTIaotXASKSEsVq8ZZzOqMAEUmRJC3eSkxnFCAiKRO3xVuJ6YwCRCSF4rR4ofw3JVKAiKRUFlq8ChCRFEt7i1cBIpJyUVu85y1ZU7abEilARDIgSot315vvle2mRAoQkYwo1OJdde1sxob+ay7XdEYBIpIhw7V4W5om8/h3Dp/OlLrFqwARyZgo686UusWrABHJoGK2eI+YcEx93DpGDBAzu9vM+sxsU2hskpmtNrOtwWN96HOLzWybmW0xs4tD42eb2fPB5241Mxv6XiIyelFbvAuWrct7TqTu2Cknxa1hNEcg/wlcMmTsJmCtu88A1gavMbOZwELg9GCf5WY2JtjnJ8A1wIzgY+jXFJGIoqw709W7L++6MyT4X/mIAeLu64DXhwzPB9qD5+3AgtD4Snd/z917gG1Ai5k1Ase6+6/d3YF7QvuISAJJ1p05+fiJ+MEP9sd977jnQKa4+x6A4PGEYHwqsDO0XW8wNjV4PnQ8LzO7xsw6zayzv78/ZokitaPQdGbVtbOpCx1lhFu8a79zAfv7ep6L+77FPoma72DIC4zn5e4r3L3Z3ZsbGhqKVpxINRvuitWWpsk88d38Ld6k4gbIq8G0hOCxLxjvBaaHtpsG7A7Gp+UZF5EiitXiHTN2XNz3ixsgHUBr8LwVeDg0vtDMxptZE7mTpRuCac4+M5sddF++GtpHRIooaot37KSpp8R9r9G0cR8Afg2cYma9ZnY1sASYa2ZbgbnBa9x9M7AKeAF4DLje3T8IvtR1wJ3kTqxuBx6NW7SIFJbvnEj7MzvytnhtTF3sIxDLNUXSy8z2AVsqXccoHA+8VukiRikrtWalTkhprUdMrG+oO2byxwH84MEP9vft6AVeY8zYceOO/4tPYnDgrT4+ePetWM3cuqJWWxpb3L250kWMxMw6s1AnZKfWrNQJ2am12HXqUnYRiU0BIiKxZSFAVlS6gFHKSp2QnVqzUidkp9ai1pn6k6gikl5ZOAIRkZRSgIhIbAoQEYlNASIisSlARCS2/wcSZfowveuVtAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import scipy.sparse as spsp\n",
    "import scipy.sparse.linalg as spsplin\n",
    "n = 32\n",
    "A = np.random.randn(n, n)\n",
    "A = A.T @ A + np.eye(n)\n",
    "\n",
    "L1d = spsp.spdiags((-1*np.ones(n), 2 * np.ones(n), -1*np.ones(n)), [-1, 0, 1], n, n)\n",
    "L2d = spsp.kronsum(L1d, L1d)\n",
    "plt.spy(L2d, markersize=1)\n",
    "A = L2d.toarray()\n",
    "# A = A + A.T\n",
    "# Q, _ = np.linalg.qr(A)\n",
    "# d = (np.random.rand(n) + 10) * 10\n",
    "# d[0] = 1\n",
    "# d[1] = 2\n",
    "# Lam = np.diag(d)\n",
    "# print(Lam.shape)\n",
    "# A = Q @ Lam @ Q.T\n",
    "true_spectrum = np.linalg.eigvalsh(A)\n",
    "print(np.sort(true_spectrum))\n",
    "x0 = np.random.randn(A.shape[0])\n",
    "num_iter = 40"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.01811230970766158\n",
      "0.01811230970766158\n"
     ]
    }
   ],
   "source": [
    "smallest_eigv, conv_sm = shift_inverse_iteration(A, x0, num_iter)\n",
    "print(smallest_eigv @ A @ smallest_eigv)\n",
    "smallest_eigv, conv_sm = shift_inverse_iteration_fast(A, x0, num_iter)\n",
    "print(smallest_eigv @ A @ smallest_eigv)\n",
    "\n",
    "# %timeit shift_inverse_iteration(A, x0, num_iter)\n",
    "# %timeit shift_inverse_iteration_fast(A, x0, num_iter)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fec26312810>]"
      ]
     },
     "execution_count": 97,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAD4CAYAAAAKA1qZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAjiElEQVR4nO3deXRV9b338fc3M4QwJQEhAxlIUERlyIMi8xBFxQtttQVrtZWW2ooI3Furvfd56n267qO1CooT4qx1uNR6q3WOTAFBZRCQQUIIhARkHgNk/j1/cLRpJAjk5OyTnM9rraxk/3Ky92ftBflk7/07e5tzDhERCU1hXgcQERHvqAREREKYSkBEJISpBEREQphKQEQkhEV4HeBsJCQkuLS0NK9jiIg0KytXrtznnEs81feaVQmkpaWxYsUKr2OIiDQrZlbc0Pd0OkhEJISpBEREQphKQEQkhKkERERCmEpARCSEqQREREKYSkBEJISFRAnsPlLOjLwCCveUeR1FRCSohEQJnKisYda8zazaftDrKCIiQSUkSiC5Qysiw42t+455HUVEJKiERAlEhIfRLT6Wor06HSQiUldIlABAekIsRXt1JCAiUlfIlEBGYizF+49TU6tnKouIfC1kSiAzoQ2VNbXsOHjC6ygiIkEjZEogPTEWgC37dF1ARORrnpeAmcWa2UozG9OU28lIOFkCW3VdQETkG34vATN71sz2mNm6euOjzWyTmRWa2V11vvVbYK6/c9TXMTaKdq0iKdKRgIjIN5riSOB5YHTdATMLBx4DrgJ6AhPMrKeZjQI2ALubIMc/MTPNEBIRqcfvj5d0zuWbWVq94f5AoXOuCMDMXgPGAm2AWE4Wwwkze9c5V+vvTF/LSIxl2Zb9TbV6EZFmJ1DPGE4CSuoslwKXOucmA5jZT4F9pyoAM5sETAJITU1tVIjMxDa8sWoHxyuraR3VrB6vLCLSJAJ1YdhOMfbNhH3n3PPOubdP9YPOuTnOuRznXE5iYmKjQqT7Lg7rlJCIyEmBKoFSIKXOcjKwM0Db/kaGb5qo7iEkInJSoEpgOZBlZulmFgWMB94K0La/kRYfi5mOBEREvtYUU0RfBZYBPcys1MwmOueqgcnAB8BGYK5zbr2/t/1dYiLD6dqulaaJioj4NMXsoAkNjL8LvOvv7Z2tjMRYnQ4SEfHx/B3DgZbhe6+Ac7qRnIhI6JVAYhvKKqrZe7TC6ygiIp4LwRLwTRPVKSERkdArAb1XQETkH0KuBLq2a0VMZJgeNSkiQgiWQFiYkRavGUIiIhCCJQAnrwvomoCISKiWQEIbth84TmV1k92wVESkWQjNEkiMpabWUXLwuNdRREQ8FZIloBlCIiInhWQJZCS0AdAMIREJeSFZAu1aRxIfG6UZQiIS8kKyBMA3Q0ing0QkxIVuCSS00S2lRSTkhWwJpCfGsq+sksMnqryOIiLimZAtgYwEPWpSRCR0SyBRM4REREK2BFI7tiY8zHQkICIhLWRLICoijJQOrTRDSERCWsiWAJw8JbRFp4NEJIR5WgJmlmFmz5jZ615sPz0hlm37j1Fbq+cNi0ho8nsJmNmzZrbHzNbVGx9tZpvMrNDM7gJwzhU55yb6O8OZykiMpbyqlq+OlHsVQUTEU01xJPA8MLrugJmFA48BVwE9gQlm1rMJtn1WdA8hEQl1fi8B51w+cKDecH+g0PeXfyXwGjD2TNZnZpPMbIWZrdi7d69fs3790HnNEBKRUBWoawJJQEmd5VIgyczizWw20MfM7j7VDzrn5jjncpxzOYmJiX4N1SkumtiocM0QEpGQFRGg7dgpxpxzbj9wa4AyfIuZaYaQiIS0QB0JlAIpdZaTgZ0B2vZppSfoofMiEroCVQLLgSwzSzezKGA88FaAtn1aGYmx7Dh0gvKqGq+jiIgEXFNMEX0VWAb0MLNSM5vonKsGJgMfABuBuc659f7e9rlIT4jFOdi2X0cDIhJ6/H5NwDk3oYHxd4F3/b29xsr03Uhu695jnH9eW4/TiIgEVkjfNgLqPHRe1wVEJASFfAnERkfQuW20ZgiJSEgK+RKAk+8c1gwhEQlFKgGge6c2bNp1lP1lFV5HEREJKJUA8JMB3aiqqeU//77B6ygiIgGlEgCyO8cxeXgWb63ZyUcbdnsdR0QkYFQCPr8alkmPznH8x9/WcaS8yus4IiIBoRLwiYoI44/XXcyeo+Xc996XXscREQkIlUAdvVPaM3FQOq98up1lW/Z7HUdEpMmpBOqZntuD1I6tueuNtZyo1P2ERKRlUwnU0yoqnPt+cBHF+48z86MCr+OIiDQplcApXJ6ZwIT+qTy9uIg1JYe8jiMi0mRUAg24++rzSYyL5rd/XUtlda3XcUREmoRKoAFtYyL5r3EX8eWuozyxcIvXcUREmoRK4DRG9ezMtZd05dEFm3VaSERaJJXAd7jn2p60axXFuMc/5tcvr2TdjsNeRxIR8RuVwHeIbxPN+1MH86uhmSwu2MeYR5bws+c+Y2XxAa+jiYg0mjnnvM5wxnJyctyKFSs82/7hE1W8tGwbz368jQPHKrk0vSOTR3RnUPcEzMyzXCIip2NmK51zOaf8nkrg7B2vrObVz0qYk7+F3UcquCipHWN7d+XKC88jpWNrr+OJiPyToC0BMxsHXAN0Ah5zzn14utcHSwl8raK6hjdW7eCFpdv4ctdRAHoltWX0hecxutd5dO8U53FCEZEmKgEzexYYA+xxzvWqMz4aeBgIB552zt13BuvqADzgnJt4utcFWwnUtW3fMT5Yv4v31+/i8+2HAMhMjOXKC89jXJ8ksjurEETEG01VAkOAMuDFr0vAzMKBAiAXKAWWAxM4WQj31lvFLc65Pb6fexB42Tm36nTbDOYSqGvX4XLyNpwshE+KDlBT6xh1QWcmj+hO75T2XscTkRDTZKeDzCwNeLtOCQwA7nHOXelbvhvAOVe/AL7+eQPuA/Kccx818JpJwCSA1NTUfsXFxeec1wsHjlXy4rJtPPfxNg6fqGJQ9wR+PTyTARnxupgsIgFxuhLw9xTRJKCkznKpb6whtwOjgOvM7NZTvcA5N8c5l+Ocy0lMTPRf0gDpGBvF1FHZfHzXCH539fls2n2UG576lB88sZT5X+6mOV2YF5GWJ8LP6zvVn7YN/pZzzs0CZvk5Q1BqEx3BpCGZ3DQgjb+sLGX2wi3c8vwKenZpyx2jsriiZ2cdGYhIwPn7SKAUSKmznAzs9PM2mrWYyHB+clk3Fv5mGA9cfwnlVTX88qWVjH3sYxZu2qMjAxEJKH+XwHIgy8zSzSwKGA+85edttAiR4WFc1y+ZD6cN4U/XXcyBY5X89Lnl/PDJZXxSpKeaiUhgnHMJmNmrwDKgh5mVmtlE51w1MBn4ANgIzHXOrfdP1JYpIjyM63NSmP+vw/jDuF5sP3Cc8XM+4canP+Xz7Qe9jiciLZzeMRxkyqtq+PMnxTyxcAv7j1Uy6oJOTMvN5sKu7byOJiLNVNC+Y/hshUIJfO1YRTXPfbyVOflFHCmv5pqLujAtN0vvQhaRs6YSaMYOn6ji6cVFPLtkKyeqahjXJ4mpI7NJjdc9ikTkzKgEWoD9ZRU8mV/EC0u3UVPruD4nhSkju9OlXSuvo4lIkFMJtCB7jpTz2IJCXvlsO2bGDf1T+fXwTDrFxXgdTUSClEqgBSo9eJxH5hXy+qpSIsONmwek8cuhmXSMjfI6mogEGZVAC7Zt3zFmzdvM31bvoFVkOD8bmM4vBmfQrnWk19FEJEioBEJA4Z6jzPxoM++s/Yq4mAh+PiiDWwalERejMhAJdSqBELLxqyPMzCvgww27ad86kklDMrh5QBqx0f6+TZSINBcqgRC0tvQQM/IKWLhpL/GxUfxqWCY3XtaNmMhwr6OJSICpBELYyuKDzMwrYEnhPhLjorltWCYTLk0lOkJlIBIqVALCJ0X7mZFXwGdbD9ClXQyTR3Tn+n4pREX4+x6CIhJsVAICgHOOpVv28+CHm1i1/RDJHVoxZUQW3++bRES4ykCkpVIJyD9xzrGwYC8z8wpYW3qYtPjWTBmZxdjeSYSH6cE2Ii1NIB8vKc2AmTG8RyfevG0gT92UQ6uoCKbPXcMVMxfx9zU7qa1tPn8YiEjjqARCmJmR27Mz79w+iCd+3JfwMOP2Vz/nqocX8/66r/SUM5EQoBIQwsKMqy7qwnt3DOHh8b2pqqnl1j+vYswjS/how26VgUgLpmsC8i3VNbW8uXonD8/bzPYDx7kkuR3TcrMZmp2Ima4ZiDQ3ujAs56SqppY3VpUya14hOw6doF+3DkzPzebyzHiVgUgzohKQRqmsrmXuihIenV/IriPlXJrekem52VyaEe91NBE5A0E7O8jMhpnZYjObbWbDvMwiDYuKCOPGy7qx8DfDuOfanhTtO8aP5nzCjU9/ysrig17HE5FGOOcSMLNnzWyPma2rNz7azDaZWaGZ3fUdq3FAGRADlJ5rFgmMmMhwfjownfzfDOc/rrmAjV8d4QdPLOWnz33G2tJDXscTkXNwzqeDzGwIJ3+Bv+ic6+UbCwcKgFxO/lJfDkwAwoF7663iFmCfc67WzDoDM5xzPz7dNnU6KLgcq6jmxWXFPJm/hUPHqxh1QWem5WZxYdd2XkcTkTpOdzronO8v7JzLN7O0esP9gULnXJFvw68BY51z9wJjTrO6g0D0uWYRb8RGR/juTprK8x9v46nFRVwzazdX9TqPqaOy6XFenNcRReQ7+Psm80lASZ3lUuDShl5sZt8HrgTaA4828JpJwCSA1NRUf+UUP4qLieT2kVncdHkazyzZyrNLtvL++l2Mubgrd4zMonunNl5HFJEG+LsETjVvsMHzTc65N4A3TrdC59wcYA6cPB3UqHTSpNq1imR6bja3DExjTn4Rzy/dxjtrdzKudxJTRmaRlhDrdUQRqcffs4NKgZQ6y8nATj9vQ4Jc+9ZR3Dn6fBbfOZyfD87g3XVfMXLGIu58fQ0lB457HU9E6vB3CSwHssws3cyigPHAW37ehjQT8W2i+d3VF5B/53BuGtCNv63eyfAHFvK7//mCnYdOeB1PRGjcFNFXgWVADzMrNbOJzrlqYDLwAbARmOucW++fqNJcdYqL4ffXXkj+b4YzoX8qf1lRwrA/LeT3b65j95Fyr+OJhDS9Y1gCbsehEzw6fzN/WVFKeJhx42XduHVoJolxmiAm0hR02wgJStv3H2fW/M28saqU6Ihwbrq8G78ckknH2Civo4m0KCoBCWpFe8uYNW8zb67ZSevIcH42MJ1fDM6gXetIr6OJtAgqAWkWNu8+ykPzNvPO2q+Ii45g4uB0bhmUTtsYlYFIY6gEpFn5ctcRZuYV8MH63bRrFcmkIRncfHkabaL9/bYWkdCgEpBmad2Owzz0UQEfbdxDh9aR/HJoJjcN6EbrKJWByNlQCUiztrrkEDPzClhUsJeENlHcOjSTGy/rRkxkuNfRRJoFlYC0CCuLDzAjr4CPC/fTKS6a24Z3Z3z/FKIjVAYip6MSkBblk6L9zPiwgM+2HaBLuxgmj+jO9f1SiIrw9BlJIkFLJSAtjnOOpVv28+CHm1i1/RDJHVoxZUQW3+ubRGS4ykCkLpWAtFjOORYV7GVGXgFrSw+TFt+aKSOzGNs7ifCwU93UViT0BO0zhkUay8wY1qMTb942kKduyqF1VATT564hd+Yi3lqzk9ra5vNHjogXVALSIpgZuT078/btg5h9Y18iw8KY8urnjH44n/e++EplINIAlYC0KGFhxuheXXjvjsE8MqEPNbWOX728ijGPLCFvw26a0+lPkUBQCUiLFBZmXHtJVz6cNpSZP7qE45XV/OLFFYx97GMWbNqjMhDx0YVhCQnVNbW8sWoHs+ZvpvTgCfqmtmd6bg8Gdo/HTBeQpWXT7CARn8rqWl5fWcqj8zez83A5/dM6Mi03mwGZ8V5HE2kyKgGReiqqa/jv5SU8tqCQ3UcquDwznum52eSkdfQ6mojfqQREGlBeVcMrn27n8YVb2FdWwZDsRKbnZtM7pb3X0UT8RiUg8h2OV1bz0rJinswv4sCxSkae34lpudn0SmrndTSRRlMJiJyhsopqXli6jTn5RRw+UcUVPTszLTebC7q09TqayDkL2hIws1TgUWAfUOCcu+90r1cJSKAcKa/i2SVbeWbxVo5WVHPNRV2YOiqLrM5xXkcTOWtNctsIM3vWzPaY2bp646PNbJOZFZrZXd+xmmzgHefcLUDPc80i4m9tYyKZOiqbJb8dweTh3Vm4aQ9XPJTP1Nc+p2hvmdfxRPzmnI8EzGwIUAa86Jzr5RsLBwqAXKAUWA5MAMKBe+ut4hagBngdcMBLzrnnTrdNHQmIVw4cq2ROfhEvLN1GRXUN3++bzJQRWaTGt/Y6msh3arLTQWaWBrxdpwQGAPc45670Ld8N4JyrXwBf//y/AZ855/LN7HXn3HWneM0kYBJAampqv+Li4nPOK9JYe49WMHvRFv78STE1tY7rc5K5bXh3kjuoDCR4BfIuoklASZ3lUt9YQ94HppjZbGDbqV7gnJvjnMtxzuUkJib6LajIuUiMi+Z/j+lJ/p3D+fGlqfx15Q6GP7CQ//jbF+w6XO51PJGz5u8ndp/q/fcNHmo459YB3/rrXyTYdW4bw3+O7cUvh2by6IJCXvushLkrSrmhfyq/Hp5Jp7gYryOKnBF/HwmUAil1lpOBnX7ehkjQ6Nq+Ff/vexex4N+GMa53V176pJgh9y/gv97ZwP6yCq/jiXwnf5fAciDLzNLNLAoYD7zl522IBJ2Ujq25/7pLmDd9KFdf1IVnlmxl8P0L+OP7X3LwWKXX8UQa1Jgpoq8Cy4AeZlZqZhOdc9XAZOADYCMw1zm33j9RRYJfWkIsM37Ymw+nDWXkBZ2ZvWgLg+9fwIMfbuLw8Sqv44l8i94xLNKENu06ysPzCnj3i13ExUTw80EZ3DIojbiYSK+jSQgJ2ncMny2VgDRX63ce5qGPNpO3YTftW0fyi8EZ/PTyNGKj/T03Q+TbVAIiQWJt6SFm5hWwYNNeOsZGcevQDH5yWRqtosK9jiYtmEpAJMisLD7IQx8VsHjzPhLaRPPrYZnccGkqMZEqA/E/lYBIkPps6wFm5G3ik6IDnNc2httGdOeHOclER6gMxH9UAiJBbumWfcz4sIAVxQdJat+KySO6c12/ZCLD/T2LW0KRSkCkGXDOsXjzPh7MK2BNySFSOrZiyogsvtcniQiVgTSCSkCkGXHOsWDTHmbkFbBuxxHSE2K5Y2QW117SlfCwU92ZReT0AnkDORFpJDNjxPmd+fvkQTz5k35ER4Qx9b9Xc+VD+by9die1tc3nDzcJfioBkSBlZlx54Xm8O2Uwj93QF4DJr3zO1bMW8/66XTSno3gJXioBkSAXFmZcc3EXPpg6hIfH96aiupZb/7ySMY8sYd7G3SoDaRSVgEgzER5mjO2dRN60ITxw/SUcLa9m4gsrGPf4UhYV7FUZyDnRhWGRZqqqppbXV5by6PxCdhw6QU63DkzPzeby7gleR5Mgo9lBIi1YZXUt/72ihMfmF7LrSDmXZXRkem4P+qd39DqaBAmVgEgIKK+q4dXPtvP4wi3sPVrB4KwEpo7Kpl+3Dl5HE4+pBERCyInKGl7+tJgnFm5h/7FKhvVIZHpuNhcnt/c6mnhEJSASgo5XVvPC0mKezN/CoeNVjLqgM9Nys7iwazuvo0mAqQREQtjR8iqe/3gbTy0u4kh5NVf1Oo9pudlkd47zOpoEiEpARDh8oopnlmzl2SVbOVZZzZiLu3LHyCy6d2rjdTRpYioBEfnGoeOVzMkv4vml2yivqmFcnySmjMgiLSHW62jSRIKiBMwsA/h3oJ1z7rqGxk5HJSDiP/vLKngyv4gXl22jqsbxg75J3D4ii5SOrb2OJn7W6BvImdmzZrbHzNbVGx9tZpvMrNDM7jrdOpxzRc65id81JiKBEd8mmt9dfQH5dw7npgHd+NvqnQx/YCG/+58v2HnohNfxJEDO9LYRzwOj6w6YWTjwGHAV0BOYYGY9zewiM3u73kcnv6YWEb/pFBfD76+9kEW/GcaE/qn8ZUUJw/60kN+/uY7dR8q9jidNLOJMXuScyzeztHrD/YFC51wRgJm9Box1zt0LjPFXQDObBEwCSE1N9ddqRaSeLu1a8Ydxvfjl0AweW1DIy59u57XlJdx4WTduHZpJYly01xGlCTTmBnJJQEmd5VLf2CmZWbyZzQb6mNndDY3V55yb45zLcc7lJCYmNiKuiJyJ5A6tuff7FzP/X4dx7SVdee7jrQy5fwH3vreRA8cqvY4nfnZGRwINONUjjhq8yuyc2w/c+l1jIhIcUuNb88D1l/DrYZnMmreZOflF/HlZMT8bmM7PB6fTvnWU1xHFDxpzJFAKpNRZTgZ2Ni6OiASbjMQ2PDS+D3nThjD8/E48uqCQwX9cwEMfFXCkvMrreNJIjSmB5UCWmaWbWRQwHnjLP7FEJNh07xTHozf05f2pgxnYPYGHPtrMoPvm8+j8zZRVVHsdT87RmU4RfRVYBvQws1Izm+icqwYmAx8AG4G5zrn1TRdVRILB+ee1ZfZP+vH27YPon96RBz4sYPAf5zN70RaOV6oMmhu9Y1hEGmVNySFm5BWwqGAvCW2iuHVoJjde1o2YyHCvo4lPULxj2B9UAiLBa2XxAWbkFfBx4X46xUVz2/DujO+fQnSEysBrKgERCZhPivYzI6+Az7YeoEu7GCaP6M71/VKIitAjzb2iEhCRgHLO8XHhfh7M28Tn2w+R3KEVU0Zk8f2+SUSEqwwCTSUgIp5wzrGoYC8z8gpYW3qYtPjWTBmZxdjeSYSHneqtRtIUGn0DORGRc2FmDOvRiTdvG8jTN+XQOiqC6XPXcMXMRby1Zie1tc3nj9CWSiUgIk3OzBjVszNv3z6I2Tf2JSIsjCmvfs7oh/N574uvVAYeUgmISMCEhRmje3XhvTsG88iEPtTUOn718iqueWQJeRt205xOT7cUKgERCbiwMOPaS7ry4bShzPzRJZyorOYXL65g7GMfs2DTHpVBAOnCsIh4rrqmljc+38GseZspPXiCPqntmZ6bzaDuCZjpAnJjaXaQiDQLldW1/HVVKY/M28zOw+X0T+vItNxsBmTGex2tWVMJiEizUlFdw9zlJTy6oJDdRyq4PDOe6bnZ5KR19Dpas6QSEJFmqbyqhlc+3c7jC7ewr6yCIdmJTM/NpndKe6+jNSsqARFp1o5XVvPSsmJmL9rCweNVjDy/E9Nys+mV1M7raM2CSkBEWoSyimpeWLqNOflFHD5RxRU9OzMtN5sLurT1OlpQUwmISItypLyK55Zs4+nFRRytqOaai7owdVQWWZ3jvI4WlFQCItIiHT5exdNLinh2yVaOV9XwL5d0ZcrILDIT23gdLaioBESkRTtwrJI5+UW8sHQbFdU1jOuTxB0js+gWH+t1tKCgEhCRkLCvrILZC7fw0ifFVNc6ruubzOQR3Unp2NrraJ5SCYhISNlzpJzHF27hlc+245zjhzkp3Da8O13bt/I6mieC4lbSZpZhZs+Y2ev1xmPNbKWZjQlUFhFp2Tq1jeGef7mQRb8Zxvj/lcrcFSUM+9NCfv/mOnYfKfc6XlA5oxIws2fNbI+Zras3PtrMNplZoZnddbp1OOeKnHMTT/Gt3wJzzzyyiMiZ6dKuFX8Y14sF/zaMH/RL4uVPtzPk/gX8379vYO/RCq/jBYUzOh1kZkOAMuBF51wv31g4UADkAqXAcmACEA7cW28Vtzjn9vh+7nXn3HW+r0cBCUAMsM859/bpcuh0kIg0xvb9x3lk/mbe+HwHkeHGzQPSmDQkg/g20V5Ha1J+uSZgZmnA23VKYABwj3PuSt/y3QDOufoFUH89dUvgv4BYoCdwAviec6623usnAZMAUlNT+xUXF59RXhGRhmzdd4xZ8zbz5uodxESG87OBafxicAbtW0d5Ha1JNNU1gSSgpM5yqW+soRDxZjYb6FOnMP7dOTcVeAV4qn4B+F4zxzmX45zLSUxMbERcEZGT0hNimfmj3nw4bQgjL+jM4wu3MOiPC5iRV8DhE1VexwuoiEb87Klu8t3gYYVzbj9wawPfe74ROUREzkn3TnE8MqEPk4d35+F5Bcyat5nnlmzlZwPTuGVQeos9MqirMUcCpUBKneVkYGfj4oiIBF6P8+J4/Mf9eHfKYAZnJzBrfiED75vP/e9/yYFjlV7Ha1KNKYHlQJaZpZtZFDAeeMs/sUREAq9n17Y8/uN+fDB1CMPP78QTi7Yw6I/zuffdjS12NtGZThF9FVgG9DCzUjOb6JyrBiYDHwAbgbnOufVNF1VEJDB6nBfHozf0JW/aEK7o2ZmnFhcx+P75/OHtDZQePO51PL/SO4ZFRL5D0d4yHluwhb+t3oFzjlEXdObmy9O4PDO+WTwDWbeNEBHxgx2HTvDyJ8W8tryEA8cq6d6pDTcP6Mb3+ibTJrox82yalkpARMSPyqtqeHvtV7ywdBtf7DhMXHQEP+iXzPU5yZx/XlvCw4Lr6EAlICLSBJxzfF5yiBeXbuOdL76iqsYRGxXORcnt6J3SgT6p7emT0p5ObWM8zakSEBFpYnuPVrB4815WlxxidckhNuw8QnXtyd+vXdvF0Du1PSkdW5MQG018mygS2vzjc8fYKCLDm+5+nioBEZEAK6+qYf3OI3y+/SCrSw6xtvQwuw6XU1nzrRsjANAmOoKIcCMizAgz3+ewf3wemJnAH8b1OqcspyuB4L2SISLSjMVEhtOvWwf6devwzZhzjqMV1ewvq2R/WQX7yirYV1bJ/rJKDp2opKbWffvDOaprHUkdmuZZCCoBEZEAMTPaxkTSNiaS9ITgePRlwB4qIyIiwUclICISwlQCIiIhTCUgIhLCVAIiIiFMJSAiEsJUAiIiIUwlICISwprVbSPMbC9Q3IhVJAD7/BTHn5Tr7CjX2VGus9MSc3VzziWe6hvNqgQay8xWNHT/DC8p19lRrrOjXGcn1HLpdJCISAhTCYiIhLBQK4E5XgdogHKdHeU6O8p1dkIqV0hdExARkX8WakcCIiJSh0pARCSEhUQJmNloM9tkZoVmdpfXeb5mZtvM7AszW21mnj0308yeNbM9ZrauzlhHM8szs82+zx1Ot44A5rrHzHb49tlqM7vag1wpZrbAzDaa2Xozu8M37uk+O00uT/eZmcWY2WdmtsaX6z99417vr4Zyef5vzJcj3Mw+N7O3fctNsr9a/DUBMwsHCoBcoBRYDkxwzm3wNBgnSwDIcc55+sYUMxsClAEvOud6+cbuBw445+7zFWcH59xvgyDXPUCZc+6BQGapl6sL0MU5t8rM4oCVwDjgp3i4z06T64d4uM/MzIBY51yZmUUCS4A7gO/j7f5qKNdoPP435ss3HcgB2jrnxjTV/8lQOBLoDxQ654qcc5XAa8BYjzMFFedcPnCg3vBY4AXf1y9w8pdJQDWQy3POua+cc6t8Xx8FNgJJeLzPTpPLU+6kMt9ipO/D4f3+aiiX58wsGbgGeLrOcJPsr1AogSSgpM5yKUHwH8PHAR+a2Uozm+R1mHo6O+e+gpO/XIBOHuepa7KZrfWdLgr4aaq6zCwN6AN8ShDts3q5wON95ju1sRrYA+Q554JifzWQC7z/N/YQcCdQW2esSfZXKJSAnWIsKNoeGOic6wtcBdzmO/0hp/cEkAn0Br4CHvQqiJm1Af4KTHXOHfEqR32nyOX5PnPO1TjnegPJQH8z6xXoDKfSQC5P95eZjQH2OOdWBmJ7oVACpUBKneVkYKdHWf6Jc26n7/Me4H84eeoqWOz2nWP++lzzHo/zAOCc2+37j1sLPIVH+8x3DvmvwMvOuTd8w57vs1PlCpZ95styCFjIyfPunu+vU+UKgv01EPgX3zXD14ARZvZnmmh/hUIJLAeyzCzdzKKA8cBbHmfCzGJ9F+8ws1jgCmDd6X8qoN4CbvZ9fTPwpodZvvH1fwKf7+HBPvNdUHwG2Oicm1HnW57us4Zyeb3PzCzRzNr7vm4FjAK+xPv9dcpcXu8v59zdzrlk51waJ39fzXfO3UhT7S/nXIv/AK7m5AyhLcC/e53HlykDWOP7WO9lLuBVTh72VnHyyGkiEA/MAzb7PncMklwvAV8Aa33/Kbp4kGsQJ08prgVW+z6u9nqfnSaXp/sMuBj43Lf9dcD/8Y17vb8ayuX5v7E6GYcBbzfl/mrxU0RFRKRhoXA6SEREGqASEBEJYSoBEZEQphIQEQlhKgERkRCmEhARCWEqARGREPb/AYhtRZcRYaSwAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.semilogy([np.linalg.norm(A @ x - (x @ A @ x) * x) for x in conv_sm])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.992228493059952\n"
     ]
    }
   ],
   "source": [
    "shift = 7\n",
    "middle_eigv, conv_m = shift_inverse_iteration_fast(A, x0, num_iter, shift)\n",
    "print(middle_eigv @ A @ middle_eigv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fec24a8d410>]"
      ]
     },
     "execution_count": 103,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAbvElEQVR4nO3de3Sc9X3n8fd3ZiTNaC6SdRnJtmzkO9iQkqA6BNpcCBCTpZC0bAPttiXlhO1JyElPe3aTPe3ZpufsOU276Z4mXZpdk7AkbQrLpqFA4gRyITgJNNjmamNsywZi+aKRfJNGd2l++8c8kkayZFseSc94ns/rHB1pfhrp+fo58uf3zO95nu+Ycw4RESl/Ib8LEBGRxaHAFxEJCAW+iEhAKPBFRAJCgS8iEhARvws4l4aGBtfa2up3GSIil5Rdu3Z1O+cap4+XdOC3trayc+dOv8sQEbmkmNnbM41rSUdEJCAU+CIiAaHAFxEJCAW+iEhAKPBFRAJCgS8iEhAKfBGRgCjLwP/+7uM8sP2Q32WIiJSUsgz8Z97I8MBPFfgiIoXKMvBrqivoGRzxuwwRkZJSnoEfq2BwJMfQ6JjfpYiIlIyyDPxUrAKAMwM6yhcRGVeegR/N94TrUeCLiEwoy8Cv0RG+iMhZFPgiIgFR1oHfMzDqcyUiIqWjrANfR/giIpPKMvB1lY6IyNnKMvArwiGqK8MKfBGRAmUZ+JBf1lHgi4hMKuvA13X4IiKTyjbwUzrCFxGZomwDX0s6IiJTlW3gp6Ja0hERKVS2ga8jfBGRqco68PuGxxgdy/ldiohISVjUwDezuJntMrNbF3pbNTGvY+ag2iuIiECRgW9mD5pZxsx2TxvfYmb7zKzdzD5X8K3PAo8Ws80LVVOtu21FRAoVe4T/ELClcMDMwsD9wC3ARuAuM9toZjcCrwOdRW7zgqSiCnwRkUKRYn7YObfdzFqnDW8G2p1zhwDM7BHgdiABxMlPAgNmts05d9YCu5ndC9wLsHLlyouuTQ3URESmKirwZ7EcOFzwuAN4t3PuPgAzuxvoninsAZxzW4GtAG1tbe5ii1Dgi4hMtRCBbzOMTQS3c+6hBdjmWSZ74ivwRURgYa7S6QBWFDxuAY4uwHbOSS2SRUSmWojA3wGsM7NVZlYJ3Ak8sQDbOadoRZjKSEhH+CIinmIvy3wYeB7YYGYdZnaPc24UuA94CtgLPOqc21N8qXOnu21FRCYVe5XOXbOMbwO2FfO754MCX0RkUtm2VgCvJ/6gAl9EBAIQ+DrCFxHJK+vAT0UjCnwREU9ZB35NrIIz/Qp8EREIQOD3Do2Sy130DbsiImWjrAM/FavAOehVi2QRkfIO/In2CrpSR0QkGIGvE7ciImUe+OqnIyIyqawDX0f4IiKTFPgiIgERiMBXx0wRkTIP/OrKMJGQ6QhfRIQyD3wzI6V+OiIiQJkHPqiBmojIuLIPfB3hi4jklX3g53viq7WCiEgwAl9H+CIi5R/46okvIpJX9oE/ftLWObVIFpFgC0Tgj+UcfcNjfpciIuKrQAQ+qL2CiEhgAl8nbkUk6AIT+DrCF5GgK/vAV098EZG8sg98HeGLiOSVfeCntIYvIgIEIPCTVRHMFPgiIosW+Gb2ETN7wMweN7ObF2u7oZCRiqqBmohIUYFvZg+aWcbMdk8b32Jm+8ys3cw+B+Cc+1fn3CeAu4GPFbPduUrF1F5BRKTYI/yHgC2FA2YWBu4HbgE2AneZ2caCp/y59/1Fo574IiJFBr5zbjtwctrwZqDdOXfIOTcMPALcbnl/DXzPOfdiMdudKwW+iMjCrOEvBw4XPO7wxj4N3AjcYWZ/NNsPm9m9ZrbTzHZ2dXXNS0EKfBERiCzA77QZxpxz7svAl8/3w865rcBWgLa2tnlpcak3QRERWZgj/A5gRcHjFuDoAmzngukqHRGRhQn8HcA6M1tlZpXAncATC7CdC5aKVTA8mmNwRC2SRSS4ir0s82HgeWCDmXWY2T3OuVHgPuApYC/wqHNuT/GlXjy1VxARKXIN3zl31yzj24Btxfzu+VQY+E2pqM/ViIj4o+xbK4B64ouIQEACXy2SRUQCEvhawxcRUeCLiARGIAI/Fc2fm1bgi0iQBSLwI+EQiaoIPQO621ZEgisQgQ/qpyMiEpjAT0bVE19Egi0wgV8Tq9B1+CISaIEKfB3hi0iQKfBFRAIiUIHfM6jAF5HgCkzgp2IV9A+PMTKW87sUERFfBCbwdbetiASdAl9EJCAU+CIiARGYwE+pJ76IBFxgAr8mpgZqIhJsgQl8HeGLSNAFJvC1hi8iQReYwK+KhIlWhBT4IhJYgQl8UHsFEQm2wAW+3gRFRIIqUIGfiuoIX0SCK1CBryUdEQkyBb6ISEAEKvBTetcrEQmwQAV+TayC3qFRxnLO71JERBZdoAJ//G7bXr0RiogE0KIFvpnFzezrZvaAmf3uYm23kO62FZEgKyrwzexBM8uY2e5p41vMbJ+ZtZvZ57zh3wS+5Zz7BHBbMdu9WAp8EQmyYo/wHwK2FA6YWRi4H7gF2AjcZWYbgRbgsPe0sSK3e1EU+CISZEUFvnNuO3By2vBmoN05d8g5Nww8AtwOdJAP/XNu18zuNbOdZrazq6urmPLOosAXkSBbiDX85UweyUM+6JcD3wZ+y8y+Ajw52w8757Y659qcc22NjY3zWljK64mv9goiEkSRBfidNsOYc871AR9fgO1dMB3hi0iQLcQRfgewouBxC3B0AbYzZ7GKMBVhU+CLSCAtRODvANaZ2SozqwTuBJ5YgO3MmZmpvYKIBFaxl2U+DDwPbDCzDjO7xzk3CtwHPAXsBR51zu0pvtT5ofYKIhJURa3hO+fummV8G7CtmN+9UFLRCnp0p62IBFCgWiuAOmaKSHAp8EVEAkKBLyISEIEM/J6BEXJqkSwiARPIwM856BvW3bYiEiyBC/zx9gpa1hGRoAlc4Ku9gogEVeACP6XAF5GAClzgjx/h625bEQmawAa+jvBFJGgCF/ha0hGRoApc4CerIjQkKnm144zfpYiILKrABb6ZcdPGJp55I8PgiC9vrSsi4ovABT7AzZua6Rse47mD3X6XIiKyaAIZ+NetqSdZFeGp3Z1+lyIismgCGfhVkTA3XJHmB3s7GR3L+V2OiMiiCGTgA3xoUzMn+4bZ8dYpv0sREVkUgQ38961vpCoS4qk9x/0uRURkUQQ28ONVEX59XSNP7zmOc2qVLCLlL7CBD7DlymaOnhnktSO6Jl9Eyl+gA//GK9KEQ8b3d2tZR0TKX6ADv7a6kmtX12kdX0QCIdCBD7BlUzMHu/poz/T6XYqIyIIKfODftLEZQMs6IlL2Ah/4zTVR3rmylqf26K5bESlvgQ98yN+E9dqRM3Sc6ve7FBGRBaPAJx/4AE/rKF9EypgCH1jVEOfy5iTf19U6IlLGFPiemzc1s/Otk3Rnh/wuRURkQSxq4JvZR8zsATN73MxuXsxtn8+WTc3kHPzwdS3riEh5uuDAN7MHzSxjZrunjW8xs31m1m5mnzvX73DO/atz7hPA3cDHLqriBXLF0iQr6mJa1hGRsjWXI/yHgC2FA2YWBu4HbgE2AneZ2UYzu8rMvjPtI13wo3/u/VzJMDO2bGrmufYT9AzqDc5FpPxccOA757YDJ6cNbwbanXOHnHPDwCPA7c6515xzt077yFjeXwPfc869ONN2zOxeM9tpZju7urou9t91UT60qZnhsRzPvJFZ1O2KiCyGYtfwlwOHCx53eGOz+TRwI3CHmf3RTE9wzm11zrU559oaGxuLLG9u3rVyCY3JKl2eKSJlKVLkz9sMY7M2l3fOfRn4cpHbXDChkHHzxiYee+kIh0/2s6Ku2u+SRETmTbFH+B3AioLHLcDRIn+nr37vPZcRCRkf/Yef89Iv9faHIlI+ig38HcA6M1tlZpXAncATxZfln8ubU3z7k9dTXRnhzq3/xrbXjvldkojIvJjLZZkPA88DG8ysw8zucc6NAvcBTwF7gUedc3sWptTFszad4LFPXseVy2v45Ddf5Cs/Oai3QRSRS56VcpC1tbW5nTt3+rb9wZEx/tO3XuXJV47ysbYV/LePXklFWDcni0hpM7Ndzrm26ePFnrQta9GKMF/62NW01lfz9z9up+N0P//wu9dQE6vwuzQRkTnT4ep5hELGn968gS/++1/hhTdP8ltfeY7DJ9VGWUQuPQr8C3THNS184w/fTVfvELf9z5/x0wOLe1OYiEixFPhz8J419Tz+qetJJ6P8/oMvcP8z7eRypXsORESkkAJ/jlob4jz2qev4jXcs478/tY//+E+71HtHRC4JCvyLUF0Z4Ut3Xs1f/MZGnnkjw21//zP2He/1uywRkXNS4F8kM+Pj16/i4XuvpW94jI/c/3Mef/mI32WJiMxKgV+kX22t47uf/jWuXJ7iM4+8zF8+uYeh0TG/yxIROYsCfx6kU1H++RPXcvd1rfyfn7/FLX/3U13FIyIlR4E/TyrCIT5/2yYe+vivknOO3/vaC3zqmy9y7MyA36WJiAAK/Hn3/g1pvv/H7+VPb1rPD/d28sG/fZb/9exBhkdzfpcmIgGnwF8A0Yown/7gOn74J+/j+rUNfOF7b3DLl7bzXHu336WJSIAp8BfQirpqHvj9Nh68u42RMcfvfPUX3PPQDr61q4Pu7JDf5YlIwKh52iK44fImrlvTwP9+9hDf/MXb/OiNDGbwjpZaPnh5mhsuT7NpWQqzmd5ATERkfqg98iLL5RyvH+vhR3sz/Hhfhlc7TuMcpJNVfGBDmiuWJmlKRUmnojSlqmhMVlEVCftdtohcQmZrj6zA91l3doif7OvimTcybN/fRe/Q6FnPqYtXkk7mwz8VrSAZjZCoipCIRkhGK0h6X6eiFaRi+c81sfzzIurfLxI46odfohoSVdxxTQt3XNNCLuc41T9MZ88Qnb2DZHoGyXhfd/YM0dU7xNHTA2SHRukdHKV/+Pw3eMUrw6Ri+QlgSXUldfHJj/rE5NcNiSrSySpqYhVaWhIpUwr8EhIKGfWJKuoTVWwkdd7nj47l6Bsao2dwhN7BUbJDo5wZGKFnYISewRF6BkbpGRzhzED+41TfMHuP93Cyb5jT/TM3fKuMhEgnq7yP/LJSOhWlORVlaW2UZTUxmmuiRCu0zCRyqVHgX8Ii4RA11SFqquf+DlyjYzlO9Y9wsm+Yk33DdGeHyPQOkenNv6rI9A7S3pXluYPd9AyevczUkKhkaU2MpTVRltXGaFkSo2VJNSvqYqyoqyYV1buCiZQaBX5ARcIhGr3zAuczMDzG8Z5Bjp0e4OiZ8c8DHD09yFsn+nju4Amy0849pKIRVtRV07Ikxsq66vxHfZzL6qpZviSm9wYW8YECX84rVhlmVUOcVQ3xGb/vnOPMwAiHTw7Qcaqfw6f66Tg1wOGT/Rzs6uPZ/V0MjkzeaRwOGctqo1xWF2dlfTWr6uO0NsRZ1VDNirpqXZUkskAU+FI0M6O2upLa6kquaqk56/u5nKMrO8Rb3X28fbKfX57o9z738d1Xj3FmYPJ8QshgWW1sYoJprY+zujHO6oYEy5fECId0QlnkYinwZcGFQkZTKkpTKsq7V9ef9f3T/cO82d3HWyf6eLO7P/91dx+PvXhkymWqlZEQrfXVrGqIs7oxwaqGOGsaE6xpjFNbXbmY/ySRS5ICX3xXW13JO1dW8s6VS6aMO+fozuYngze7sxzq6uNQdx/tmSw/fiPDyNjkPST18cp8+KfHJ4EEqxvjtCyp1qsCEY8CX0qWmU2cWN68qm7K90bHcnScGuBQd5aDmT4OdmU52JXl6T2dnOg7PPG8ykiI1Q1x1qQTrG1MTHxe3RjXpaUSOAp8uSRFwiFaG/Ine2+4fOr3TvUNc6g7S3smy8Gu/CuC1zrOsO21Y4zfWG4GLUtirEsnWetNAmubEqxNJ3RJqZQtBb6UnSXxSq6J13HNZVNfFQyOjPGmtyTUnsnS3pXlYCbLzw50Mzw2eRVRU6qKtenExGSwLp1gXVOSurjOE8ilTYEvgRGtCHPF0hRXLJ16F/NYznH4ZD8HxieCTJb2TC//b+dh+graV9THK1kzPgF4k8C6pgSNiSq1o5BLggJfAi8csonloZs2Nk2MO+c4dmaQA5ksBzp7ac9kOZDJ8uQrR6fcfVwTq/AmgARr00nWN+VfHTSlNBFIaVnUwDezOLAd+Avn3HcWc9sic2VmLKuNsaw2xvvWN06MO+fI9A7lJ4DOXm9CyPK93cc53T95wjgZjXivBvKvBNY1JVmXTrC0JqqJQHxxQYFvZg8CtwIZ59yVBeNbgC8BYeCrzrkvnOdXfRZ49CJrFSkJZpP3FVy/tmFi3DnHib5h9o+/GujMsr+zlx/u7eT/7pycCBJVEdamE6xvSrC+KamJQBbNBfXDN7P3AlngG+OBb2ZhYD9wE9AB7ADuIh/+fzXtV/wh8A6gAYgC3RdyhB+EfvgSDCeyQ/lXAt6rgvFJoTs7PPGcZFWEtU358wPjE8H6pgTNKU0EMjdF9cN3zm03s9Zpw5uBdufcIW8DjwC3O+f+ivyrgekFfACIAxuBATPb5pzLzfC8e4F7AVauXHkh5YmUvPG219dOu9P4pPeKoHAi+NHeDI/u7Jh4zvjSUOEksL4pSTqpcwQyN8Ws4S8HDhc87gDePduTnXN/BmBmd5M/wj8r7L3nbQW2Qv4Iv4j6REpeXbySa1fXnzURnMgOsb8zf7XQfm9p6OnXO3lkx+R/uVQ0ctYkoKuG5FyKCfyZ/qLOG9DOuYeK2KZIINQnqnhPoor3rJk6EXRnh/KvCLxJIH+y+BgPvzDZgK62uoL13oni8UlgfVOShsT5W2FLeSsm8DuAFQWPW4CjxZUjIufSkKiiIVHFdWumnizuyg5NTAL7O/PLQ9MvH62LV04sDa33rhparxvKAqWYwN8BrDOzVcAR4E7gd+alKhG5YGZGOhklnTz7qqFM79CUSWB/Zy+PvXRkyhvWNCQqWZdOsqF58tXA+nTyot5JTUrbhV6W+TDwfqDBzDrIX0f/NTO7D3iK/JU5Dzrn9ixYpSIyJ4WXj/76uqn3ERw7MzhlaWh/JnvWncXpZFV+EkgXviJIkFSvoUvWBV2W6RddlimyeHI5x9EzAxzozLLPezVwoDPfamJgZHIiWFYTnQj/8WWhdekE8SrduF8qirosU0TKXyhktCyppmVJNR+4PD0xnss5Ok4NeK8E8pPAvuO9PH/oBMOjkxfbLa+NTS4LpfMTwdp0glil2lCXCgW+iJxTKGSsrK9mZX01Nxb0GhrLOX55sj8/ERzPLwsd6Ozlpwe6Jt6cxgxW1lVPLAut914R6P0I/KHAF5GLEg7ZxHsPf2hT88T4yFiOt0/0Tdw/ML5E9JN9GUZz+YkgZNBaH2ddU4INE/cSJFnVEKcyEvLrn1T2FPgiMq8qwiHWppOsTSf58FVLJ8aHR3O82d3nXTU0eY7gB6934s0DRLxJpPD+gfVNSVrrq4mENREUS4EvIouiMhJiQ3P+8s9CgyNjHOoqnAiy7D56hm27J9+hrDIcYnVjnHVNSTYUnCxeWaf3LJ4LBb6I+CpaEWbjshQbl019Y5qB4THaM9mJk8X7j/fy4tunePKVyfs7qyIh1jQmJk4Wb/AmguW1MUKaCM6iwBeRkhSrDHNVSw1XtdRMGc8OjU5OBN7J4n87dILHXjoy8ZzqyvDEu5IVniwOegtqBb6IXFISVRGuXlHL1Stqp4yfGRiZ0mxuf2cvz+7v4lu7CjqPei2oNzQFs/OobrwSkbJ2un+Y/d6VQgcKThaf6Jt8L4KaWMXkjWTpBOubk5d0wzndeCUigVRbXcnmVXVsXlU3ZXx659H9nb1899Vj/PPAZOfRKQ3nmr3JoCnJkku04ZwCX0QCabbOo+MN5/Yd9yaDzEwN56rY0JzwbihL5r9uSpIq8T5DCnwREc/5Gs5NLgvl7yp+dOdh+gsazjWnopOvBJpLr89QaVQhIlLCzIxltTGW1cb4wIapfYaOnB6YuH9gfGnoHw+dYGhan6H1Td4k4GOfIQW+iMhFCoWMFXXVrKir5oNXTO0zdHi8z1DBZPDz9hMMj+UnAj/6DCnwRUTmWThktDbEaW2Ic3NBn6HRsRxvneif0lpi/yx9hr7427/Cu1Yumde6FPgiIoskEg6xNp1gbToxa5+hA5297OvspSE+/5eEKvBFRHw2W5+h+ab2cyIiAaHAFxEJCAW+iEhAKPBFRAJCgS8iEhAKfBGRgFDgi4gEhAJfRCQgSvoNUMysC3j7In+8Aeiex3Lmi+qaG9U1N6prbsq1rsucc43TB0s68IthZjtnescXv6muuVFdc6O65iZodWlJR0QkIBT4IiIBUc6Bv9XvAmahuuZGdc2N6pqbQNVVtmv4IiIyVTkf4YuISAEFvohIQJRl4JvZFjPbZ2btZvY5v+sZZ2ZvmdlrZvayme30sY4HzSxjZrsLxurM7AdmdsD7PL/vrXbxdX3ezI54++xlM/uwD3WtMLNnzGyvme0xs894477us3PU5es+M7Oomb1gZq94df2lN+73/pqtLt//xrw6wmb2kpl9x3s87/ur7NbwzSwM7AduAjqAHcBdzrnXfS2MfOADbc45X2/0MLP3AlngG865K72xvwFOOue+4E2SS5xzny2Buj4PZJ1zX1zMWqbVtRRY6px70cySwC7gI8Dd+LjPzlHXb+PjPjMzA+LOuayZVQA/Az4D/Cb+7q/Z6tqCz39jXn1/ArQBKefcrQvxf7Icj/A3A+3OuUPOuWHgEeB2n2sqKc657cDJacO3A1/3vv46+eBYVLPU5Tvn3DHn3Ive173AXmA5Pu+zc9TlK5eX9R5WeB8O//fXbHX5zsxagH8HfLVgeN73VzkG/nLgcMHjDkrgP4HHAU+b2S4zu9fvYqZpcs4dg3yQAGmf6yl0n5m96i35LPpSUyEzawXeCfyCEtpn0+oCn/eZtzzxMpABfuCcK4n9NUtd4P/f2N8B/xnIFYzN+/4qx8C3GcZKYhYHrnfOvQu4BfiUt4Qh5/YVYA1wNXAM+Fu/CjGzBPAvwB8753r8qmO6GeryfZ8558acc1cDLcBmM7tysWuYySx1+bq/zOxWIOOc27XQ2yrHwO8AVhQ8bgGO+lTLFM65o97nDPAY+eWnUtHprQmPrw1nfK4HAOdcp/efNAc8gE/7zFvz/Rfgm865b3vDvu+zmeoqlX3m1XIa+An5dXLf99dMdZXA/roeuM07x/cIcIOZ/RMLsL/KMfB3AOvMbJWZVQJ3Ak/4XBNmFvdOrGFmceBmYPe5f2pRPQH8gff1HwCP+1jLhPE/eM9H8WGfeSf7vgbsdc79j4Jv+brPZqvL731mZo1mVut9HQNuBN7A//01Y11+7y/n3H9xzrU451rJ59WPnXP/gYXYX865svsAPkz+Sp2DwJ/5XY9X02rgFe9jj591AQ+Tf+k6Qv4V0T1APfAj4ID3ua5E6vpH4DXgVe8/wFIf6vo18suCrwIvex8f9nufnaMuX/cZ8A7gJW/7u4H/6o37vb9mq8v3v7GCGt8PfGeh9lfZXZYpIiIzK8clHRERmYECX0QkIBT4IiIBocAXEQkIBb6ISEAo8EVEAkKBLyISEP8fAPw83q/tomIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.semilogy([np.linalg.norm(A @ x - (x @ A @ x) * x) for x in conv_m])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.992228493059952\n"
     ]
    }
   ],
   "source": [
    "shift = 7\n",
    "large_eigv, conv_l = shift_inverse_iteration(A, x0, num_iter, shift)\n",
    "print(large_eigv @ A @ large_eigv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fec0e835110>]"
      ]
     },
     "execution_count": 105,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAbvElEQVR4nO3de3Sc9X3n8fd3ZiTNaC6SdRnJtmzkO9iQkqA6BNpcCBCTpZC0bAPttiXlhO1JyElPe3aTPe3ZpufsOU276Z4mXZpdk7AkbQrLpqFA4gRyITgJNNjmamNsywZi+aKRfJNGd2l++8c8kkayZFseSc94ns/rHB1pfhrp+fo58uf3zO95nu+Ycw4RESl/Ib8LEBGRxaHAFxEJCAW+iEhAKPBFRAJCgS8iEhARvws4l4aGBtfa2up3GSIil5Rdu3Z1O+cap4+XdOC3trayc+dOv8sQEbmkmNnbM41rSUdEJCAU+CIiAaHAFxEJCAW+iEhAKPBFRAJCgS8iEhAKfBGRgCjLwP/+7uM8sP2Q32WIiJSUsgz8Z97I8MBPFfgiIoXKMvBrqivoGRzxuwwRkZJSnoEfq2BwJMfQ6JjfpYiIlIyyDPxUrAKAMwM6yhcRGVeegR/N94TrUeCLiEwoy8Cv0RG+iMhZFPgiIgFR1oHfMzDqcyUiIqWjrANfR/giIpPKMvB1lY6IyNnKMvArwiGqK8MKfBGRAmUZ+JBf1lHgi4hMKuvA13X4IiKTyjbwUzrCFxGZomwDX0s6IiJTlW3gp6Ja0hERKVS2ga8jfBGRqco68PuGxxgdy/ldiohISVjUwDezuJntMrNbF3pbNTGvY+ag2iuIiECRgW9mD5pZxsx2TxvfYmb7zKzdzD5X8K3PAo8Ws80LVVOtu21FRAoVe4T/ELClcMDMwsD9wC3ARuAuM9toZjcCrwOdRW7zgqSiCnwRkUKRYn7YObfdzFqnDW8G2p1zhwDM7BHgdiABxMlPAgNmts05d9YCu5ndC9wLsHLlyouuTQ3URESmKirwZ7EcOFzwuAN4t3PuPgAzuxvoninsAZxzW4GtAG1tbe5ii1Dgi4hMtRCBbzOMTQS3c+6hBdjmWSZ74ivwRURgYa7S6QBWFDxuAY4uwHbOSS2SRUSmWojA3wGsM7NVZlYJ3Ak8sQDbOadoRZjKSEhH+CIinmIvy3wYeB7YYGYdZnaPc24UuA94CtgLPOqc21N8qXOnu21FRCYVe5XOXbOMbwO2FfO754MCX0RkUtm2VgCvJ/6gAl9EBAIQ+DrCFxHJK+vAT0UjCnwREU9ZB35NrIIz/Qp8EREIQOD3Do2Sy130DbsiImWjrAM/FavAOehVi2QRkfIO/In2CrpSR0QkGIGvE7ciImUe+OqnIyIyqawDX0f4IiKTFPgiIgERiMBXx0wRkTIP/OrKMJGQ6QhfRIQyD3wzI6V+OiIiQJkHPqiBmojIuLIPfB3hi4jklX3g53viq7WCiEgwAl9H+CIi5R/46okvIpJX9oE/ftLWObVIFpFgC0Tgj+UcfcNjfpciIuKrQAQ+qL2CiEhgAl8nbkUk6AIT+DrCF5GgK/vAV098EZG8sg98HeGLiOSVfeCntIYvIgIEIPCTVRHMFPgiIosW+Gb2ETN7wMweN7ObF2u7oZCRiqqBmohIUYFvZg+aWcbMdk8b32Jm+8ys3cw+B+Cc+1fn3CeAu4GPFbPduUrF1F5BRKTYI/yHgC2FA2YWBu4HbgE2AneZ2caCp/y59/1Fo574IiJFBr5zbjtwctrwZqDdOXfIOTcMPALcbnl/DXzPOfdiMdudKwW+iMjCrOEvBw4XPO7wxj4N3AjcYWZ/NNsPm9m9ZrbTzHZ2dXXNS0EKfBERiCzA77QZxpxz7svAl8/3w865rcBWgLa2tnlpcak3QRERWZgj/A5gRcHjFuDoAmzngukqHRGRhQn8HcA6M1tlZpXAncATC7CdC5aKVTA8mmNwRC2SRSS4ir0s82HgeWCDmXWY2T3OuVHgPuApYC/wqHNuT/GlXjy1VxARKXIN3zl31yzj24Btxfzu+VQY+E2pqM/ViIj4o+xbK4B64ouIQEACXy2SRUQCEvhawxcRUeCLiARGIAI/Fc2fm1bgi0iQBSLwI+EQiaoIPQO621ZEgisQgQ/qpyMiEpjAT0bVE19Egi0wgV8Tq9B1+CISaIEKfB3hi0iQKfBFRAIiUIHfM6jAF5HgCkzgp2IV9A+PMTKW87sUERFfBCbwdbetiASdAl9EJCAU+CIiARGYwE+pJ76IBFxgAr8mpgZqIhJsgQl8HeGLSNAFJvC1hi8iQReYwK+KhIlWhBT4IhJYgQl8UHsFEQm2wAW+3gRFRIIqUIGfiuoIX0SCK1CBryUdEQkyBb6ISEAEKvBTetcrEQmwQAV+TayC3qFRxnLO71JERBZdoAJ//G7bXr0RiogE0KIFvpnFzezrZvaAmf3uYm23kO62FZEgKyrwzexBM8uY2e5p41vMbJ+ZtZvZ57zh3wS+5Zz7BHBbMdu9WAp8EQmyYo/wHwK2FA6YWRi4H7gF2AjcZWYbgRbgsPe0sSK3e1EU+CISZEUFvnNuO3By2vBmoN05d8g5Nww8AtwOdJAP/XNu18zuNbOdZrazq6urmPLOosAXkSBbiDX85UweyUM+6JcD3wZ+y8y+Ajw52w8757Y659qcc22NjY3zWljK64mv9goiEkSRBfidNsOYc871AR9fgO1dMB3hi0iQLcQRfgewouBxC3B0AbYzZ7GKMBVhU+CLSCAtRODvANaZ2SozqwTuBJ5YgO3MmZmpvYKIBFaxl2U+DDwPbDCzDjO7xzk3CtwHPAXsBR51zu0pvtT5ofYKIhJURa3hO+fummV8G7CtmN+9UFLRCnp0p62IBFCgWiuAOmaKSHAp8EVEAkKBLyISEIEM/J6BEXJqkSwiARPIwM856BvW3bYiEiyBC/zx9gpa1hGRoAlc4Ku9gogEVeACP6XAF5GAClzgjx/h625bEQmawAa+jvBFJGgCF/ha0hGRoApc4CerIjQkKnm144zfpYiILKrABb6ZcdPGJp55I8PgiC9vrSsi4ovABT7AzZua6Rse47mD3X6XIiKyaAIZ+NetqSdZFeGp3Z1+lyIismgCGfhVkTA3XJHmB3s7GR3L+V2OiMiiCGTgA3xoUzMn+4bZ8dYpv0sREVkUgQ38961vpCoS4qk9x/0uRURkUQQ28ONVEX59XSNP7zmOc2qVLCLlL7CBD7DlymaOnhnktSO6Jl9Eyl+gA//GK9KEQ8b3d2tZR0TKX6ADv7a6kmtX12kdX0QCIdCBD7BlUzMHu/poz/T6XYqIyIIKfODftLEZQMs6IlL2Ah/4zTVR3rmylqf26K5bESlvgQ98yN+E9dqRM3Sc6ve7FBGRBaPAJx/4AE/rKF9EypgCH1jVEOfy5iTf19U6IlLGFPiemzc1s/Otk3Rnh/wuRURkQSxq4JvZR8zsATN73MxuXsxtn8+WTc3kHPzwdS3riEh5uuDAN7MHzSxjZrunjW8xs31m1m5mnzvX73DO/atz7hPA3cDHLqriBXLF0iQr6mJa1hGRsjWXI/yHgC2FA2YWBu4HbgE2AneZ2UYzu8rMvjPtI13wo3/u/VzJMDO2bGrmufYT9AzqDc5FpPxccOA757YDJ6cNbwbanXOHnHPDwCPA7c6515xzt077yFjeXwPfc869ONN2zOxeM9tpZju7urou9t91UT60qZnhsRzPvJFZ1O2KiCyGYtfwlwOHCx53eGOz+TRwI3CHmf3RTE9wzm11zrU559oaGxuLLG9u3rVyCY3JKl2eKSJlKVLkz9sMY7M2l3fOfRn4cpHbXDChkHHzxiYee+kIh0/2s6Ku2u+SRETmTbFH+B3AioLHLcDRIn+nr37vPZcRCRkf/Yef89Iv9faHIlI+ig38HcA6M1tlZpXAncATxZfln8ubU3z7k9dTXRnhzq3/xrbXjvldkojIvJjLZZkPA88DG8ysw8zucc6NAvcBTwF7gUedc3sWptTFszad4LFPXseVy2v45Ddf5Cs/Oai3QRSRS56VcpC1tbW5nTt3+rb9wZEx/tO3XuXJV47ysbYV/LePXklFWDcni0hpM7Ndzrm26ePFnrQta9GKMF/62NW01lfz9z9up+N0P//wu9dQE6vwuzQRkTnT4ep5hELGn968gS/++1/hhTdP8ltfeY7DJ9VGWUQuPQr8C3THNS184w/fTVfvELf9z5/x0wOLe1OYiEixFPhz8J419Tz+qetJJ6P8/oMvcP8z7eRypXsORESkkAJ/jlob4jz2qev4jXcs478/tY//+E+71HtHRC4JCvyLUF0Z4Ut3Xs1f/MZGnnkjw21//zP2He/1uywRkXNS4F8kM+Pj16/i4XuvpW94jI/c/3Mef/mI32WJiMxKgV+kX22t47uf/jWuXJ7iM4+8zF8+uYeh0TG/yxIROYsCfx6kU1H++RPXcvd1rfyfn7/FLX/3U13FIyIlR4E/TyrCIT5/2yYe+vivknOO3/vaC3zqmy9y7MyA36WJiAAK/Hn3/g1pvv/H7+VPb1rPD/d28sG/fZb/9exBhkdzfpcmIgGnwF8A0Yown/7gOn74J+/j+rUNfOF7b3DLl7bzXHu336WJSIAp8BfQirpqHvj9Nh68u42RMcfvfPUX3PPQDr61q4Pu7JDf5YlIwKh52iK44fImrlvTwP9+9hDf/MXb/OiNDGbwjpZaPnh5mhsuT7NpWQqzmd5ATERkfqg98iLL5RyvH+vhR3sz/Hhfhlc7TuMcpJNVfGBDmiuWJmlKRUmnojSlqmhMVlEVCftdtohcQmZrj6zA91l3doif7OvimTcybN/fRe/Q6FnPqYtXkk7mwz8VrSAZjZCoipCIRkhGK0h6X6eiFaRi+c81sfzzIurfLxI46odfohoSVdxxTQt3XNNCLuc41T9MZ88Qnb2DZHoGyXhfd/YM0dU7xNHTA2SHRukdHKV/+Pw3eMUrw6Ri+QlgSXUldfHJj/rE5NcNiSrSySpqYhVaWhIpUwr8EhIKGfWJKuoTVWwkdd7nj47l6Bsao2dwhN7BUbJDo5wZGKFnYISewRF6BkbpGRzhzED+41TfMHuP93Cyb5jT/TM3fKuMhEgnq7yP/LJSOhWlORVlaW2UZTUxmmuiRCu0zCRyqVHgX8Ii4RA11SFqquf+DlyjYzlO9Y9wsm+Yk33DdGeHyPQOkenNv6rI9A7S3pXluYPd9AyevczUkKhkaU2MpTVRltXGaFkSo2VJNSvqYqyoqyYV1buCiZQaBX5ARcIhGr3zAuczMDzG8Z5Bjp0e4OiZ8c8DHD09yFsn+nju4Amy0849pKIRVtRV07Ikxsq66vxHfZzL6qpZviSm9wYW8YECX84rVhlmVUOcVQ3xGb/vnOPMwAiHTw7Qcaqfw6f66Tg1wOGT/Rzs6uPZ/V0MjkzeaRwOGctqo1xWF2dlfTWr6uO0NsRZ1VDNirpqXZUkskAU+FI0M6O2upLa6kquaqk56/u5nKMrO8Rb3X28fbKfX57o9z738d1Xj3FmYPJ8QshgWW1sYoJprY+zujHO6oYEy5fECId0QlnkYinwZcGFQkZTKkpTKsq7V9ef9f3T/cO82d3HWyf6eLO7P/91dx+PvXhkymWqlZEQrfXVrGqIs7oxwaqGOGsaE6xpjFNbXbmY/ySRS5ICX3xXW13JO1dW8s6VS6aMO+fozuYngze7sxzq6uNQdx/tmSw/fiPDyNjkPST18cp8+KfHJ4EEqxvjtCyp1qsCEY8CX0qWmU2cWN68qm7K90bHcnScGuBQd5aDmT4OdmU52JXl6T2dnOg7PPG8ykiI1Q1x1qQTrG1MTHxe3RjXpaUSOAp8uSRFwiFaG/Ine2+4fOr3TvUNc6g7S3smy8Gu/CuC1zrOsO21Y4zfWG4GLUtirEsnWetNAmubEqxNJ3RJqZQtBb6UnSXxSq6J13HNZVNfFQyOjPGmtyTUnsnS3pXlYCbLzw50Mzw2eRVRU6qKtenExGSwLp1gXVOSurjOE8ilTYEvgRGtCHPF0hRXLJ16F/NYznH4ZD8HxieCTJb2TC//b+dh+graV9THK1kzPgF4k8C6pgSNiSq1o5BLggJfAi8csonloZs2Nk2MO+c4dmaQA5ksBzp7ac9kOZDJ8uQrR6fcfVwTq/AmgARr00nWN+VfHTSlNBFIaVnUwDezOLAd+Avn3HcWc9sic2VmLKuNsaw2xvvWN06MO+fI9A7lJ4DOXm9CyPK93cc53T95wjgZjXivBvKvBNY1JVmXTrC0JqqJQHxxQYFvZg8CtwIZ59yVBeNbgC8BYeCrzrkvnOdXfRZ49CJrFSkJZpP3FVy/tmFi3DnHib5h9o+/GujMsr+zlx/u7eT/7pycCBJVEdamE6xvSrC+KamJQBbNBfXDN7P3AlngG+OBb2ZhYD9wE9AB7ADuIh/+fzXtV/wh8A6gAYgC3RdyhB+EfvgSDCeyQ/lXAt6rgvFJoTs7PPGcZFWEtU358wPjE8H6pgTNKU0EMjdF9cN3zm03s9Zpw5uBdufcIW8DjwC3O+f+ivyrgekFfACIAxuBATPb5pzLzfC8e4F7AVauXHkh5YmUvPG219dOu9P4pPeKoHAi+NHeDI/u7Jh4zvjSUOEksL4pSTqpcwQyN8Ws4S8HDhc87gDePduTnXN/BmBmd5M/wj8r7L3nbQW2Qv4Iv4j6REpeXbySa1fXnzURnMgOsb8zf7XQfm9p6OnXO3lkx+R/uVQ0ctYkoKuG5FyKCfyZ/qLOG9DOuYeK2KZIINQnqnhPoor3rJk6EXRnh/KvCLxJIH+y+BgPvzDZgK62uoL13oni8UlgfVOShsT5W2FLeSsm8DuAFQWPW4CjxZUjIufSkKiiIVHFdWumnizuyg5NTAL7O/PLQ9MvH62LV04sDa33rhparxvKAqWYwN8BrDOzVcAR4E7gd+alKhG5YGZGOhklnTz7qqFM79CUSWB/Zy+PvXRkyhvWNCQqWZdOsqF58tXA+nTyot5JTUrbhV6W+TDwfqDBzDrIX0f/NTO7D3iK/JU5Dzrn9ixYpSIyJ4WXj/76uqn3ERw7MzhlaWh/JnvWncXpZFV+EkgXviJIkFSvoUvWBV2W6RddlimyeHI5x9EzAxzozLLPezVwoDPfamJgZHIiWFYTnQj/8WWhdekE8SrduF8qirosU0TKXyhktCyppmVJNR+4PD0xnss5Ok4NeK8E8pPAvuO9PH/oBMOjkxfbLa+NTS4LpfMTwdp0glil2lCXCgW+iJxTKGSsrK9mZX01Nxb0GhrLOX55sj8/ERzPLwsd6Ozlpwe6Jt6cxgxW1lVPLAut914R6P0I/KHAF5GLEg7ZxHsPf2hT88T4yFiOt0/0Tdw/ML5E9JN9GUZz+YkgZNBaH2ddU4INE/cSJFnVEKcyEvLrn1T2FPgiMq8qwiHWppOsTSf58FVLJ8aHR3O82d3nXTU0eY7gB6934s0DRLxJpPD+gfVNSVrrq4mENREUS4EvIouiMhJiQ3P+8s9CgyNjHOoqnAiy7D56hm27J9+hrDIcYnVjnHVNSTYUnCxeWaf3LJ4LBb6I+CpaEWbjshQbl019Y5qB4THaM9mJk8X7j/fy4tunePKVyfs7qyIh1jQmJk4Wb/AmguW1MUKaCM6iwBeRkhSrDHNVSw1XtdRMGc8OjU5OBN7J4n87dILHXjoy8ZzqyvDEu5IVniwOegtqBb6IXFISVRGuXlHL1Stqp4yfGRiZ0mxuf2cvz+7v4lu7CjqPei2oNzQFs/OobrwSkbJ2un+Y/d6VQgcKThaf6Jt8L4KaWMXkjWTpBOubk5d0wzndeCUigVRbXcnmVXVsXlU3ZXx659H9nb1899Vj/PPAZOfRKQ3nmr3JoCnJkku04ZwCX0QCabbOo+MN5/Yd9yaDzEwN56rY0JzwbihL5r9uSpIq8T5DCnwREc/5Gs5NLgvl7yp+dOdh+gsazjWnopOvBJpLr89QaVQhIlLCzIxltTGW1cb4wIapfYaOnB6YuH9gfGnoHw+dYGhan6H1Td4k4GOfIQW+iMhFCoWMFXXVrKir5oNXTO0zdHi8z1DBZPDz9hMMj+UnAj/6DCnwRUTmWThktDbEaW2Ic3NBn6HRsRxvneif0lpi/yx9hr7427/Cu1Yumde6FPgiIoskEg6xNp1gbToxa5+hA5297OvspSE+/5eEKvBFRHw2W5+h+ab2cyIiAaHAFxEJCAW+iEhAKPBFRAJCgS8iEhAKfBGRgFDgi4gEhAJfRCQgSvoNUMysC3j7In+8Aeiex3Lmi+qaG9U1N6prbsq1rsucc43TB0s68IthZjtnescXv6muuVFdc6O65iZodWlJR0QkIBT4IiIBUc6Bv9XvAmahuuZGdc2N6pqbQNVVtmv4IiIyVTkf4YuISAEFvohIQJRl4JvZFjPbZ2btZvY5v+sZZ2ZvmdlrZvayme30sY4HzSxjZrsLxurM7AdmdsD7PL/vrXbxdX3ezI54++xlM/uwD3WtMLNnzGyvme0xs894477us3PU5es+M7Oomb1gZq94df2lN+73/pqtLt//xrw6wmb2kpl9x3s87/ur7NbwzSwM7AduAjqAHcBdzrnXfS2MfOADbc45X2/0MLP3AlngG865K72xvwFOOue+4E2SS5xzny2Buj4PZJ1zX1zMWqbVtRRY6px70cySwC7gI8Dd+LjPzlHXb+PjPjMzA+LOuayZVQA/Az4D/Cb+7q/Z6tqCz39jXn1/ArQBKefcrQvxf7Icj/A3A+3OuUPOuWHgEeB2n2sqKc657cDJacO3A1/3vv46+eBYVLPU5Tvn3DHn3Ive173AXmA5Pu+zc9TlK5eX9R5WeB8O//fXbHX5zsxagH8HfLVgeN73VzkG/nLgcMHjDkrgP4HHAU+b2S4zu9fvYqZpcs4dg3yQAGmf6yl0n5m96i35LPpSUyEzawXeCfyCEtpn0+oCn/eZtzzxMpABfuCcK4n9NUtd4P/f2N8B/xnIFYzN+/4qx8C3GcZKYhYHrnfOvQu4BfiUt4Qh5/YVYA1wNXAM+Fu/CjGzBPAvwB8753r8qmO6GeryfZ8558acc1cDLcBmM7tysWuYySx1+bq/zOxWIOOc27XQ2yrHwO8AVhQ8bgGO+lTLFM65o97nDPAY+eWnUtHprQmPrw1nfK4HAOdcp/efNAc8gE/7zFvz/Rfgm865b3vDvu+zmeoqlX3m1XIa+An5dXLf99dMdZXA/roeuM07x/cIcIOZ/RMLsL/KMfB3AOvMbJWZVQJ3Ak/4XBNmFvdOrGFmceBmYPe5f2pRPQH8gff1HwCP+1jLhPE/eM9H8WGfeSf7vgbsdc79j4Jv+brPZqvL731mZo1mVut9HQNuBN7A//01Y11+7y/n3H9xzrU451rJ59WPnXP/gYXYX865svsAPkz+Sp2DwJ/5XY9X02rgFe9jj591AQ+Tf+k6Qv4V0T1APfAj4ID3ua5E6vpH4DXgVe8/wFIf6vo18suCrwIvex8f9nufnaMuX/cZ8A7gJW/7u4H/6o37vb9mq8v3v7GCGt8PfGeh9lfZXZYpIiIzK8clHRERmYECX0QkIBT4IiIBocAXEQkIBb6ISEAo8EVEAkKBLyISEP8fAPw83q/tomIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "res = np.array([np.linalg.norm(A @ x - (x @ A @ x) * x) for x in conv_l])\n",
    "plt.semilogy(res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Метод Релея\n",
    "\n",
    "- Адаптивный сдвиг на каждой итерации\n",
    "\n",
    "$$ \\sigma_k = \\frac{(Ax_k, x_k)}{(x_k, x_k)} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def rayleigh_quotient_iteration(A, x0, num_iter, shift):\n",
    "    x = x0.copy()\n",
    "    x = x / np.linalg.norm(x)\n",
    "    conv = [x]\n",
    "    n = A.shape[0]\n",
    "    for i in range(1):\n",
    "        x = np.linalg.solve(A - shift * np.eye(n), x)\n",
    "        x = x / np.linalg.norm(x)\n",
    "        conv.append(x)\n",
    "        \n",
    "    for i in range(num_iter):\n",
    "        lam = x @ A @ x\n",
    "        x = np.linalg.solve(A - lam * np.eye(n), x)\n",
    "        x = x / np.linalg.norm(x)\n",
    "        conv.append(x)\n",
    "    return x, conv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "eigv_rq, conv_rq = rayleigh_quotient_iteration(A, x0, num_iter, 7.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.992228533552932\n",
      "[0.01811231 0.04519876 0.04519876 ... 7.95480124 7.95480124 7.98188769]\n",
      "[ 0.25665493 -1.23724391  1.60446583 ...  0.24743374 -0.29562927\n",
      "  1.14999884]\n"
     ]
    }
   ],
   "source": [
    "print(eigv_rq @ A @ eigv_rq)\n",
    "print(np.sort(np.linalg.eigvalsh(A)))\n",
    "print(x0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.89789243e+00 1.21314070e-01 1.32531566e-04 8.11727899e-07\n",
      " 3.23598523e-13 2.40727774e-15 2.27862475e-15 3.17035131e-15\n",
      " 1.77829958e-15 2.01790993e-15 1.77494928e-15 1.56757314e-15\n",
      " 1.60813439e-15 1.52801119e-15 1.61379221e-15 1.59361975e-15\n",
      " 1.62471396e-15 1.70191047e-15 2.30138098e-15 1.50819493e-15\n",
      " 1.52370035e-15 1.52846982e-15 1.81663206e-15 1.74373134e-15\n",
      " 1.90911123e-15 1.57494886e-15 3.20245019e-15 1.53502909e-15\n",
      " 1.87027368e-15 1.64262216e-15 1.86028201e-15 1.63206942e-15\n",
      " 2.32600653e-15 3.18078087e-15 1.89426563e-15 1.78238821e-15\n",
      " 1.67673786e-15 1.90365873e-15 1.54222803e-15 1.56701285e-15\n",
      " 1.63213639e-15 1.55712700e-15]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAD4CAYAAAAKA1qZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAo7ElEQVR4nO3deXRc9X338fd3RsuMZEljLba12JZtDHhFgLEhpE8IZbGDwSltn7DkaWhpKT0lSU/6JCFN2zRpKfSELklwQpzAY9IUODSlxBinkDhQm2LAC8bYeMEIg+VNsmTZlqxtZn7PH3ckj+SRrWWk0Xg+r3PuuTO/We5X1+P5zO/e+7vXnHOIiEhm8qW6ABERSR2FgIhIBlMIiIhkMIWAiEgGUwiIiGSwrFQXMBilpaWuuro61WWIiKSVzZs3H3XOlSV6LK1CoLq6mk2bNqW6DBGRtGJmH/b3mDYHiYhkMIWAiEgGUwiIiGSwtNonICKZoauri7q6Otrb21NdSloJBAJUVVWRnZ094NcoBERkzKmrq6OgoIDq6mrMLNXlpAXnHI2NjdTV1TFt2rQBv06bg0RkzGlvb6ekpEQBMAhmRklJyaB7TykNATO7xszWm9mjZnZNKmsRkbFFATB4Q1lnSQ8BM3vczOrNbHuf9sVmttvM9prZ/bFmB7QAAaAu2bX0OHEQXv57OLp3xBYhIpKORqInsBJYHN9gZn5gObAEmA3cbmazgfXOuSXAV4FvjkAtno6T8N//AAffGrFFiMj5xe/3U1NTw9y5c7n55ptpbm7ueWzHjh1ce+21XHjhhcyYMYNvfOMbRKPR1BU7DEkPAefcOqCpT/NCYK9zrtY51wk8DSxzznWvtWNAbqL3M7N7zGyTmW1qaGgYWlFFk715c7+D5kREegkGg2zdupXt27dTXFzM8uXLAWhra+OWW27h/vvvZ8+ePbzzzju8+eabfOc730lxxUMzWvsEKoH9cffrgEozu9XMfgj8K/BIohc651Y45xY45xaUlSU89cW55eRBfhk0fzS014tIRrvqqqs4cOAAAE8++SRXX301N9xwAwB5eXk88sgjfPvb305liUM2WoeIJtpb4ZxzzwLPjkoFoSkKAZE09M3nd/DuwRNJfc/ZFYV84+Y5A3puJBJh7dq13H333YC3Kejyyy/v9ZwZM2bQ1tZGc3MzoVAoqbWOtNHqCdQBk+PuVwEHR2nZHoWAiAxCW1sbNTU1lJSU0NTUxPXXXw94x+MnOgonXa/XPlo9gY3ATDObBhwAbgPuGKVle0JTYNcLEI2CT8MjRNLFQH+xJ1v3PoHjx4+zdOlSli9fzhe+8AXmzJnDunXrej23traW0tLStOsFwMgcIvoUsAG4yMzqzOxu51wYuA94EdgJPOOc25HsZZ9VaApEOqHlyKguVkTSW1FREd/97nd5+OGH6erq4s477+TVV1/lV7/6FeD1GL7whS/wzW+O3AGOI2kkjg663TlX7pzLds5VOecei7Wvcc5d6Jyb4Zx7INnLPafQVG+uTUIiMkiXXnopl1xyCU8//TTBYJBVq1bxwAMPcOGFF1JaWsrVV1/NnXfemeoyhyRzzh0UmuLNmz+CKYtSW4uIjHktLS297j///PM9t+fOncvLL78MwHPPPceXvvQl7rjjDqZOnTqqNSZD5mwc11gBERkBn/70p6mtrU3LAIBMCgGNFRAROUPmhADoMFERkT4UAiIiGSzzQuD4fm+sgIiIZGAIaKyAiEiPDAsBjRUQkYF54IEHmDNnDvPnz6empoY33ngDgOrqao4ePXrG81etWsVDDz0EQENDA4sWLeLSSy9l/fr1fP/73+93OR/72McA2LdvH08++eQI/CVnl2EhEDdWQESkHxs2bGD16tVs2bKFbdu28atf/YrJkyef9TXdp5cGWLt2LRdffDFvvfUWkydPPmsIvPbaa8DQQiASiQzq+YlkVghorICIDMChQ4coLS0lN9e7zElpaSkVFRU9j3/ve9/jsssuY968eezatQuAlStXct9997F161a+8pWvsGbNGmpqavjqV7/K+++/T01NDV/+8pfPWNa4ceMAuP/++1m/fj01NTX88z//M5FIhC9/+ctcccUVzJ8/nx/+8IcAvPLKK3zyk5/kjjvuYN68ecP+WzNnxDBorIBIOvrF/XD4neS+56R5sOShfh++4YYb+Na3vsWFF17Iddddx2c+8xk+8YlP9DxeWlrKli1b+P73v8/DDz/Mj3/8457Hampq+Na3vsWmTZt45JFH2LdvHzt27GDr1q1nLemhhx7i4YcfZvXq1QCsWLGCoqIiNm7cSEdHR69rGLz55pts376dadOmDWMleDKrJwA6TFREzmncuHFs3ryZFStWUFZWxmc+8xlWrlzZ8/itt94KwOWXX86+fftGpIaXXnqJn/zkJ9TU1LBo0SIaGxt57733AFi4cGFSAgAyrScAXggc2pbqKkRkoM7yi30k+f1+rrnmGq655hrmzZvHE088wV133QXQs5nI7/cTDodHZPnOOb73ve9x44039mp/5ZVXyM/PT9pyMrMnoLECInIWu3fv7vnVDbB169YhnxuooKCAkydPDvp5N954Iz/4wQ/o6uoCYM+ePbS2tg6phrPJvJ5A0eTTYwUKy1NdjYiMQS0tLXz+85+nubmZrKwsLrjgAlasWDGk9yopKeHqq69m7ty5LFmypN9rEc+fP5+srCwuueQS7rrrLr74xS+yb98+LrvsMpxzlJWV8dxzzw3jr0rM0umSaAsWLHCbNm0a3pvseQme/F24+5cweWFyChORpNq5cyezZs1KdRlpKdG6M7PNzrkFiZ6fmZuDQDuHRUTIyBDQWAERkW6ZFwI5+ZBXqp6AyBiXTpuqx4qhrLOUhoCZfdrMfmRmPzezG0ZtwRorIDKmBQIBGhsbFQSD4JyjsbGRQCAwqNcl/eggM3scWArUO+fmxrUvBr4D+IEfO+cecs49BzxnZuOBh4GXkl1PQqEpcGT7qCxKRAavqqqKuro6GhoaUl1KWgkEAlRVVQ3qNSNxiOhK4BHgJ90NZuYHlgPXA3XARjNb5Zx7N/aUv4w9PjpCU2D3L7yxAr7M2yImMtZlZ2cnbUSsnF3SvwGdc+uApj7NC4G9zrla51wn8DSwzDz/APzCObcl2bX0KzQFIh3QWj9qixQRGYtG62dwJbA/7n5drO3zwHXA75jZvYleaGb3mNkmM9uUtK6hrisgIgKMXghYgjbnnPuuc+5y59y9zrlHE73QObfCObfAObegrKwsOdVorICICDB6IVAHxF+RoQo4OErLPpPGCoiIAKMXAhuBmWY2zcxygNuAVaO07DNprICICDACIWBmTwEbgIvMrM7M7nbOhYH7gBeBncAzzrkdyV72oGisgIhI8g8Rdc7d3k/7GmBNspc3ZBorICKSgaeN6BaaAs26roCIZLbMDgGNFRCRDJfBIaCxAiIiGRwCGisgIpLBIaCxAiIimRsCGisgIpLBIQAaKyAiGU8hoBAQkQyW4SEwWWMFRCSjZUwItHaEz2wMTY2NFdDVi0QkM2VECGzd38zV//BrXq9t7P2ADhMVkQyXESFw0cQCxuVm8Rf/+Q4d4cjpB3pCQIeJikhmyogQCOb4+dtPz6W2oZVHX6k9/UBR91gB9QREJDNlRAgAfPKiCSydX87yl/dS29DiNeaOg7wShYCIZKyMCQGAv755NrnZPr7+n9txznmNOkxURDJYRoXAhIIA9y+5mA21jTy75YDXqBAQkQyWUSEAcPsVU7hsSoi/e+Fdmlo7vRA4vh+6ewYiIhkk40LA5zP+/tZ5nGwP8+Cand5YgXA7tOi6AiKSeTIuBAAunlTIH/2v6fz75jp2toW8Rm0SEpEMlJEhAPCFa2cyuTjIt99o8xo0VkBEMlBKQ8DMppvZY2b2s9FedjDHz999eh6vN+V7DeoJiEgGGnIImNnjZlZvZtv7tC82s91mttfM7j/bezjnap1zdw+1huH6xIVlXHfJdJpcAScOv5+qMkREUmY4PYGVwOL4BjPzA8uBJcBs4HYzm21m88xsdZ9pwjCWnTR/uXQWB2wCTXs36gghEck4Qw4B59w6oKlP80Jgb+wXfifwNLDMOfeOc25pn2lMHI4zoSBAbdVvUd2xm/C2f091OSIioyrZ+wQqgf1x9+tibQmZWYmZPQpcamZf6+c595jZJjPb1NAwMqd8zrri93k7Oh33X1+H9hMjsgwRkbEo2SFgCdr63cbinGt0zt3rnJvhnHuwn+escM4tcM4tKCsrS1qh8RbNKOOvun6frLYGeCVhGSIi56Vkh0AdMDnufhVwMMnLSLrScbm0T7iEl8d9Ct74IRzefu4XiYicB5IdAhuBmWY2zcxygNuAVUlexoi4cnoJXz9xKy4Yghf+XJecFJGMMJxDRJ8CNgAXmVmdmd3tnAsD9wEvAjuBZ5xzO5JT6si6anoJhzqDfHjpV2D/6/D2U6kuSURkxGUN9YXOudv7aV8DrBlyRSmyaHoJAC/4r+VPJy+CX/41XPwpCI5PcWUiIiMnY08b0Vdxfg4XTypgQ+0xuOkfoa0J1v5tqssSERlRCoE4V04vYdOHTXSWzoGFfwybHocDW1JdlojIiFEIxLlyegntXVHermuGT34Nxk2I7SSOnPO1IiLpSCEQ58rpxZjB6+83QqAIbngADm6BLU+kujQRkRGhEIgTysvh4kmFbKht9Brm/Q5U/wb86m9g65MQ7kxpfSIiyaYQ6OOq6SVs/vAYHeEImMHN34HCSnjuT+Bf5sH6f4K2Y6kuU0QkKRQCfVw1o4SOcJStHzV7DSUz4E9eg8/+B0yYBWu/Cf80B37xVTi2L5WliogMm0Kgj4XVsf0CtXEnSDWDC66D33sO7n0VZt8CGx+D714Kz/wefPS6TkMtImlJIdBHUV42cyoK2VB7NPETJs2D33oU/mwbXP1FqH0FHr8RHv24FwwdJ0e1XhGR4VAIJHDltBK2fNRMe9dZDg0trIDr/ga+tBNu/i6YD174EvzjLO+w0iPvjlq9IiJDpRBI4KoZJXSGo7zVvV/gbHLy4fLPwR+vgz9cC7OWwpZ/hR9cBY8vgW3P6BoFIjJmDfncQeezK6YV4zPYUNvIVTNKBvYiM6ha4E03/j289VPY9Bg8+0fgz4Fpn/AC4qJPeYPQRETGAIVAAoWBbOZWFvF693iBwcorhqu/AFfdB3Vvws7nYddqeP6L8PyfwZQr4eKlXiiMr05m6SIig6IQ6MeV00tY+T/7aO+KEMj2D+1NfD7vC3/KlXDD38GR7bBztRcIL33dm/JKoXSmN5XE5qUXQmgq+PXPIyIjS98y/bhqegkr1tWy5cNjfOyC0uG/oZl3ZNGked55iZo+gD0vQv27cPQ92LUGTsUdkeTLhsJyyC2E3AJvyhl3+nZuAWTlQlbA29yUFYjdz427HfCm7GDsfmyenaeAERFAIdCvBdXj8fuMDbWNyQmBvoqnwZX39m5rOwZH98LRPd508hB0tEDHCWhtgKZa7xDUjhboah3e8n3ZkJMH2fmxeWzKyesdNjnjIDd2PyfWFij0zq2UGzf36RgDkXSkEOhHwXD3CwxFcDxMvsKbziUahUgHhNsh3BE3tfeZt3nzrrZYWzt0tXsh0nkqbn4KOluh/TgcPwCdLacDiAEMhMst9KZgCAKhPvOi2O3xkDfemweLvbkCRCSlFAJncdX0Eh57tZa2zgjBnCHuFxgpPh/4gt6mnpHknBcQHS2xXsgJb2o/4QVGR2zefb/9OLQ3e72Wtmbvdtep/t/f/F5QBMdDXokXDnklXlj0uh+b8ku9QFFwiCSFQuAsrpxezKP//T6bPzzGx2eOwCahdGDmjYXIyYeCiUN7j3CnFwZtx3pPp5ri7jd594/XwaG34VSj19NJWJMvFhqlsXAo9sIhrzRuXtL7flbOkFeByPlMIXAWV1QXx/YLHM3cEEiGrBxvbMRgxkd090BONXmBcKox7vbR022tjdC4F/a/4d130cTvl1vkBUN+2elwyC/tc78sdr9EO84lY6T8k25m+cA64BvOudWpridefm4W86uKep9MTkZHfA8kNHlgr4lGvR5H61EvKHrmseBobfDaju2DA5u8266fU4MEx58OhfiA6HW77PTmKbMk/eEio2vIIWBmjwNLgXrn3Ny49sXAdwA/8GPn3EPneKuvAs8MtY6R1n2oaGtHmPzclGemnI3P520ayisGLjz3888IjYbTQdFzuxHqd0Hrem+TVcLlZvcOiHET4sJiwunHxk3QpikZc4bzrbYSeAT4SXeDmfmB5cD1QB2w0cxW4QXCg31e/wfAfOBdIDCMOkbUx2eW8v1X3mfdngaWzCtPdTmSTIMNjUiXt0mqtQFa672AaK0/MzyOvue1h9sTv0+g6Mxw6AmRCb2DJLdQvQwZUUMOAefcOjOr7tO8ENjrnKsFMLOngWXOuQfxeg29mNkngXxgNtBmZmuc671R18zuAe4BmDJlylDLHbKF1cWUjsth9TuHFAKZzp/t7RwfyA5y57xDblvrvXBoqY8LitjtlgZo2A371vd/tTp/blzPojsgyhLfzisG3xg7ik3GvGRv36gE9sfdrwMW9fdk59zXAczsLuBo3wCIPWcFsAJgwYIFo37lliy/jyVzy/nZ5jpOdYbJy9EmIRkAs9ggu3FQPP3cz490xW2G6hscsenkITi8zXss2pVgmb7YTu7ucIgFRMLbZd7occl4yf5GS9RvPecXt3NuZZLrSKqb5pfzr69/yK931bN0fkWqy5HzkT92mpDCAfQ2nfP2ZbQ09O5VdAdI9+1jG73b/Y0uDxTF7bfo3iyVKDQmeDvotVnqvJTsEKgD4g/lqAIOJnkZo+6K6mLKCnJZ/fYhhYCknlls1PV4KBvAvozO1n6CIn6z1C74YJ0XLolk551jk1R3aEzQYL40k+wQ2AjMNLNpwAHgNuCOJC9j1Pl9xk3zynnqzY9o6QgzTkcJSTrpPtR2IKctD3d6R0q11PfZfxE3P14HB7f0f4itL6v3ju6eHd9lp293z/NKNSYjxYZziOhTwDVAqZnV4R3n/5iZ3Qe8iHdE0OPOuR1JqTTFbppfzsrX9rF25xGW1VSmuhyRkZGV4106tXAAPd5o1DtsNn7fRUt9n01U9d7O79Z6iHQmeBOLjfhOtEkqQa8je8weSJi2hnN00O39tK8B1gy5ojHq8injmVQYYPW2QwoBEfA2+XSPvD4X57zzTPXaDNVn3toAB7Z4886WxO+TW9i7J5God9F9O2ec9mMMgPphA+TzGZ+aV85PX/+Qk+1dFASyU12SSPowi51NtghKLzj38ztP9eldJDhqqmE37Hu1/0F8WcHEO7njexfdwREcn7GBoRAYhJvml/P4/3zAL989wq2XVaW6HJHzV04e5EyF8VPP/dyew2v7bIaK368x4P0YiXoXfTZVnWfnljp//pJRcNmUEJWhIC9sO6QQEBkrBnN4bTTqDczrCYk++zK6A+ToHq8t4Zls4/ZjnK130R0iY3w8hkJgEMyMT82bxMrX9nG8rYuioDYJiaQVny92NtkSmDDr7M/t3o/RswkqUXA0wIHN59iPUdR7kF6i3kX34ynYj6EQGKSl8yv40foPeGnHYX53wQDPbiki6Sd+P0bJjHM/v/NU78F6QxmP8dln4YLfTOqfcS4KgUGaX1VE1fggL7xzSCEgIqfl5EFO9RDGY8T1LMouHukqz6AQGCQz46b55Ty2/gOOtXYyPl+nBRaRQRrMeIwRprHdQ3Dz/ArCUcdL7x5OdSkiIsOiEBiCORWFTC3JY/W2Q6kuRURkWBQCQ2BmLJ1fzmvvN9LY0s/F0EVE0oBCYIhumldBJOp4cceRVJciIjJkCoEhmlVewPTSfF54J+3PlC0iGUwhMETdm4Q2vN9Iw0ltEhKR9KQQGIab5lcQdfBf27WDWETSk0JgGC6aVEB5UYAtHzWnuhQRkSFRCAxT1fggB5vbUl2GiMiQKASGqSIU5OBxhYCIpCeFwDBVhIIcPt5ONOpSXYqIyKApBIapIhSkK+I4qkFjIpKGFALDVBnyLnx9QPsFRCQNpTQEzOw3zOxRM/uxmb2WylqGqiIUBOBgc3uKKxERGbwhh4CZPW5m9Wa2vU/7YjPbbWZ7zez+s72Hc269c+5eYDXwxFBrSaXTIaCegIikn+FcT2Al8Ajwk+4GM/MDy4HrgTpgo5mtAvzAg31e/wfOufrY7TuAPxxGLSlTGMhmXG6WNgeJSFoacgg459aZWXWf5oXAXudcLYCZPQ0sc849CCxN9D5mNgU47pw70c/j9wD3AEyZMmWo5Y6oilBAPQERSUvJ3idQCeyPu18Xazubu4H/19+DzrkVzrkFzrkFZWVlSSgx+TRWQETSVbIvL2kJ2s56AL1z7htJrmHUVYSCbKs7nuoyREQGLdk9gTog/urrVcB5f67lylCQptZO2jojqS5FRGRQkh0CG4GZZjbNzHKA24BVSV7GmFMRGytwSJuERCTNDOcQ0aeADcBFZlZnZnc758LAfcCLwE7gGefcjuSUOnZVFGmsgIikp+EcHXR7P+1rgDVDrigNaayAiKQrnTYiCSYVBTDTqSNEJP0oBJIg2+9jQkGuegIiknYUAkmisQIiko4UAklSEQpqx7CIpB2FQJJUhrzLTDqni8uISPpQCCRJRVGAjnCUptbOVJciIjJgCoEk0XUFRCQdKQSSpDsEdJioiKQThUCSaMCYiKQjhUCSjM/LJpDtUwiISFpRCCSJmWmsgIikHYVAElWGghzQjmERSSMKgSSqKApqc5CIpBWFQBJVhII0nOygI6yLy4hIelAIJFH3xWWOHO9IcSUiIgOjEEiiSo0VEJE0oxBIonKNFRCRNKMQSKLyIm9zkEJARNKFQiCJAtl+SsflaKyAiKQNhUCSVWisgIikkVELATObbmaPmdnP4tryzewJM/uRmd05WrWMJI0VEJF0MqAQMLPHzazezLb3aV9sZrvNbK+Z3X+293DO1Trn7u7TfCvwM+fcHwG3DKryMaoiFOSQLi4jImlioD2BlcDi+AYz8wPLgSXAbOB2M5ttZvPMbHWfaUI/71sF7I/dPi9GWFWEArR2RjjRFk51KSIi55Q1kCc559aZWXWf5oXAXudcLYCZPQ0sc849CCwd4PLr8IJgK/0EkpndA9wDMGXKlAG+berEjxUoystOcTUiImc3nH0ClZz+FQ/eF3plf082sxIzexS41My+Fmt+FvhtM/sB8Hyi1znnVjjnFjjnFpSVlQ2j3NGhsQIikk4G1BPohyVo63dDuHOuEbi3T1sr8PvDqGHM6T51hA4TFZF0MJyeQB0wOe5+FXBweOWkv9L8XHL8Pp06QkTSwnBCYCMw08ymmVkOcBuwKjllpS+fzygPBXTBeRFJCwM9RPQpYANwkZnVmdndzrkwcB/wIrATeMY5t2PkSk0fGisgIulioEcH3d5P+xpgTVIrOg9UhIK89v7RVJchInJOOm3ECKgMBThyop1wJJrqUkREzkohMAIqQkGiDo6c1MVlRGRsUwiMAI0VEJF0oRAYAZUhXVdARNKDQmAElBfpMpMikh4UAiMgPzeLUF62egIiMuYpBEaIN1ZAA8ZEZGxTCIyQipAGjInI2KcQGCGVoYBCQETGPIXACKkIBTnRHuZke1eqSxER6ZdCYIR0jxU4dFz7BURk7FIIjJDusQI6TFRExjKFwAip0KhhEUkDCoERMqEggN9nCgERGdMUAiPE7zMmFeriMiIytikERlBlKKh9AiIypikERlBFKMAhXXBeRMYwhcAIqggFOdTcTntXJNWliIgkpBAYQVfNKCEcdbyyuz7VpYiIJKQQGEFXTS+hdFwuP996MNWliIgkNGohYGbTzewxM/tZXNssM3vUzH5mZn8yWrWMliy/j6Xzy1m7q16njxCRMWlAIWBmj5tZvZlt79O+2Mx2m9leM7v/bO/hnKt1zt3dp22nc+5e4H8DCwZbfDpYVlNBZzjKizuOpLoUEZEzDLQnsBJYHN9gZn5gObAEmA3cbmazzWyema3uM03o743N7BbgVWDtkP6CMa5mcogpxXn8fOuBVJciInKGAYWAc24d0NSneSGwN/YLvxN4GljmnHvHObe0z9TvnlHn3Crn3MeAOxM9bmb3mNkmM9vU0NAwsL9qDDEzbrmkgv/Ze5SGkx2pLkdEpJfh7BOoBPbH3a+LtSVkZiVm9ihwqZl9LdZ2jZl918x+CKxJ9Drn3Arn3ALn3IKysrJhlJs6y2oqiDp4YZt2EIvI2JI1jNdagjbX35Odc43AvX3aXgFeGUYNaWHmxAJmlRfy87cPctfV01JdjohIj+H0BOqAyXH3qwD91O3HspoK3vqomY8aT6W6FBGRHsMJgY3ATDObZmY5wG3AquSUdf65+ZIKAJ7XJiERGUMGeojoU8AG4CIzqzOzu51zYeA+4EVgJ/CMc27HyJWa3ipDQa6oHs9zbx3AuX63momIjKoB7RNwzt3eT/sa+tmhK2e6paaSv3puO7sOn2RWeWGqyxER0WkjRtNN88rJ8plOIyEiY4ZCYBQV5+fwGzNLef7tg0Sj2iQkIqmnEBhly2oqOdDcxuaPjqW6FBERhcBou372RALZPp1GQkTGBIXAKMvPzeK6WRNZ885huiLRVJcjIhlOIZACy2oqaWrt5NW9R1NdiohkOIVACnziwjKKgtmsGsZRQtqxLIN1vvU8z7e/J1WGc+4gGaKcLB+fmjeJn289SFtnhGCO/6zPrz/Zzo6DJ3j34Al2HDzOjoMn2N90igVTi7lpfjlL5k5iQmFglKqXdBKJOv57Tz0/ff0jXt5dzxVTi7nr6mpumD2RLH96/AaMRB37GlvZeegEOw+dYNehk+w8dIKDx9uZW1nIb148ketmTWRuZSFmiU5pJmdj6TR6dcGCBW7Tpk2pLiMpNrzfyO0/ep2pJXkUBrLJzfKRm+0jN8tPIDZvPtXJjoMnqI87BfWU4jzmVBRSGQqy/r2j7D5yEjNYWF3M0vnlLJ5bTllBbr/Lbe+KcKozQkEgi+w0+RJIprbOCO/Vn2TX4ZPsPnyStq4IC6aOZ9H0EipDwVSXlzSNLR08s6mOf3vjQ+qOtVFWkMsNsyey7r0G9je1UV4U4LNXTuX2hVMozs9Jdbk9wpEoe460sK2umbfrjvPuwePsPnKS9i7vV7/fZ8woy2dWeSEVoSAbP2hiy0fHiDqYWJjLtRdP5PrZE/jYjFIC2Wf/cZVJzGyzcy7hhbsUAikSjToe/MVO6o610RGO0t4VoSMcpSMcoaMrSkc4Sl6On9nlhcypLGJORSGzKwopDGT3ep89R07ywrZDrN52kPcbWvEZLJpWwvSyfJpPdXHsVCfHTnXRfKqTY6c6e/4zARQEshifl8P4vGxCveY5jM/PpiiYzfi8HEJ53rwoL5txOVlEnCMSdYSjjkjEEXGOcDRKNApZfiMny0dulo8cv6/fX2bOudjfG6Uj9rebgc8sNnnXYvDF2qLOW040Smzu1dB9uyMcpTMSpSs27wx7U0c4yr7GVnYf9r74P2xspXtLWneNJzvCAFSND7JoWgmLphdz5bQSJhcHh/XLMhyJ0tjayZET7Rw50UH9ydj8RDtdEUdOlpHt9/VMOVk+cvyGmdHSEaalPUxrR5iTHd68JTYvCmZTXhSkvChAeSg2LwpQEQqyv+kUP339Q9a8c5jOSJQrpxfzf66s5oY5E8n2+4hEHb/eVc8Tr+3j1b1HycnyseySCj73sWrmVBRytKWTQ8fbONjcxsHmdg42t3HoeDuNrR0UBLIZn5fN+HzvM1Lc/dnIz6ErEqX5VFfPZ+54WxfHWjtpbuuivStCKC+H4tjnqzjfe11xfg4FgWw+ONrC2/uPs62umR0HT9AR9j6jhYEs5lYWMau8kIsneWfinTlxHLlZvb/cG1s6eHl3A2t3HmHdngZaOyMEsn1cMGEcgSw/gezYD6tsf+y+r6ftjMez/QSyfLR1RXr+luae/z/ePBx1lBcFqRrvTZWhIJXjg1SNz2N8XvZZPzORqKMrcvqz2hU5fT8ccd7/H3/s/0/35PcNu9emEMgAzjn2HGnhhW0HWbP9ME2tnT1f3n2/5PNy/JxoC8cCondINLd29XwpJkNOlo9cv9fLAfNCLvYFPVrMoLokn4smFnDRpAIunuTNp5bkA7Dr8AneqG3ijQ8aefODJo6d8q4HPaEgl3GBLHAQdQ5HbO6ITd7/ne7/QfH/lcLRKE2tnfTddWMGJfm55Gb56IpEY5PrCa5u2X4jPzeLcXFTfm4WeTl+mk91cfiE9wXdkWA9FuRm8duXV3HnoinMnFjQ73p578hJntiwj//YfIC2rgg5fh+dfbaz52b5qAwFKc7PoaUjTFOr96XY93l9BbJ9hILel31utp/jsc/Z8bbE19oOZvuZW1nI/KoQ86uKmF8Vorokb9Ah3BGO8EZtE2t3HmH/sTbauyKxKUp77AdWT1s4SmQA+9bifywV5eWQ5TMONrdRd6yNlj7/VwLZPrJ9vtM/WmKfk0jUnfFZGAy/z1g8dxLL77hsSK9XCMigdEWiHG/zgsH7NXT6dmtnmCyf4ff5YvPeUzgSPf0Lv+fXeITOcJSoo2dTV/zmr+5f5Jj3HyYa+9KN/w/kM8PnM/xm+H1e78Dvs55596+m3Ngv6p5f1lk+KoqC59zv0i0adextaOGN2kbe+qiZjkgUg169EzMwuuee7u8qi7X4fEbZuBwmFAaYWBhgQkEuEwsDlI7L6fdXnXNe7yrq3Bm/dvt7/rFTXRw63sah5nYOnWgnmO1nydxJ5OcOfHff8bYunt1Sx+Hj7VTEehYVoSAVoWDCX7bOOU51RnoC4dipTrL8Fvui9L74+9sUE459to6d6qSp1QuFKcV5zCjLT8k+iq5I9HRIdEXoCHu3A9k+Qnk5FAWz+91s6pzjRFuYuuZT1B1r48AxrwcVcQ5/7PPqi+vNdn9+vc+n9XxOvcnI8vkIR0//v+ns06udMSGf37q0akh/p0JARCSDnS0EMm/PoIiI9FAIiIhkMIWAiEgGUwiIiGQwhYCISAZTCIiIZDCFgIhIBlMIiIhksLQaLGZmDcCHw3iLUkAn8T83raeB0XoaGK2ngRupdTXVOVeW6IG0CoHhMrNN/Y2ak9O0ngZG62lgtJ4GLhXrSpuDREQymEJARCSDZVoIrEh1AWlC62lgtJ4GRutp4EZ9XWXUPgEREekt03oCIiISRyEgIpLBMiIEzGyxme02s71mdn+q6xlLzOxxM6s3s+1xbcVm9kszey82H5/KGlPNzCab2ctmttPMdpjZF2PtWk99mFnAzN40s7dj6+qbsXatqwTMzG9mb5nZ6tj9UV9P530ImJkfWA4sAWYDt5vZ7NRWNaasBBb3absfWOucmwmsjd3PZGHgz51zs4ArgT+NfYa0ns7UAVzrnLsEqAEWm9mVaF3154vAzrj7o76ezvsQABYCe51ztc65TuBpYFmKaxoznHPrgKY+zcuAJ2K3nwA+PZo1jTXOuUPOuS2x2yfx/tNWovV0Budpid3Njk0OraszmFkVcBPw47jmUV9PmRAClcD+uPt1sTbp30Tn3CHwvgCBCSmuZ8wws2rgUuANtJ4Sim3i2ArUA790zmldJfYvwFeAaFzbqK+nTAgBS9Cm42Jl0MxsHPAfwJ85506kup6xyjkXcc7VAFXAQjObm+KSxhwzWwrUO+c2p7qWTAiBOmBy3P0q4GCKakkXR8ysHCA2r09xPSlnZtl4AfBvzrlnY81aT2fhnGsGXsHb56R11dvVwC1mtg9vE/W1ZvZTUrCeMiEENgIzzWyameUAtwGrUlzTWLcK+Fzs9ueAn6ewlpQzMwMeA3Y65/4p7iGtpz7MrMzMQrHbQeA6YBdaV704577mnKtyzlXjfSf92jn3WVKwnjJixLCZfQpv+5sfeNw590BqKxo7zOwp4Bq8U9geAb4BPAc8A0wBPgJ+1znXd+dxxjCzjwPrgXc4vf32L/D2C2g9xTGz+Xg7NP14PzKfcc59y8xK0LpKyMyuAf6vc25pKtZTRoSAiIgklgmbg0REpB8KARGRDKYQEBHJYAoBEZEMphAQEclgCgERkQymEBARyWD/H8PffhmqY1EdAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "res = np.array([np.linalg.norm(A @ x - (x @ A @ x) * x) for x in conv_rq])\n",
    "plt.semilogy(res, label=\"RQ\")\n",
    "plt.semilogy([np.linalg.norm(A @ x - (x @ A @ x) * x) for x in conv_l], label=\"Shift iter\")\n",
    "plt.legend()\n",
    "print(res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Как найти другие собственные значения?\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.         0.01811231 0.04519876 ... 7.95480124 7.95480124 7.98188769]\n"
     ]
    }
   ],
   "source": [
    "B = A - eigv_rq @ (A @ eigv_rq) * np.outer(eigv_rq, eigv_rq)\n",
    "print(np.sort(np.linalg.eigvals(B)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.992228533552932\n"
     ]
    }
   ],
   "source": [
    "eigv_rq_B, conv_rq_B = rayleigh_quotient_iteration(B, x0, num_iter, 7.)\n",
    "print(eigv_rq_B @ B @ eigv_rq_B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.20120895e-14+0.j 1.24344979e-14+0.j 1.81123097e-02+0.j ...\n",
      " 7.95480124e+00+0.j 7.95480124e+00+0.j 7.98188769e+00+0.j]\n"
     ]
    }
   ],
   "source": [
    "C = B - eigv_rq_B @ (A @ eigv_rq_B) * np.outer(eigv_rq_B, eigv_rq_B)\n",
    "print(np.sort(np.linalg.eigvals(C)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.9909438451461705\n"
     ]
    }
   ],
   "source": [
    "eigv_rq_C, conv_rq_C = rayleigh_quotient_iteration(C, x0, num_iter, 7.)\n",
    "print(eigv_rq_C @ C @ eigv_rq_C)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### И так далее для других собственных векторов....\n",
    "\n",
    "- Какая правильная стратегия использования метода Релея?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Матричные функции: введение"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "\n",
    "$$ B = f(A) = \\sum_{k=0}^{\\infty} c_k A^k$$\n",
    "\n",
    "- Это не поэлементная функция\n",
    "- Нужно уметь вычислять ряд\n",
    "- Устойчивость"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Вычисление степени матрицы\n",
    "\n",
    "- Дана матрица $A$\n",
    "- Нужно вычислить $A^k$\n",
    "- Это необходимо при вычисление матричных функций, которые представляются в виде матричного ряда"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Используем спектральное разложение\n",
    "\n",
    "- $A = U\\Lambda U^{-1}$\n",
    "- Тогда $A^k = U \\Lambda^k U^{-1}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Матричная экспонента: знакомство или напоминание\n",
    "\n",
    "$$\\exp(A) = \\sum_{k=0}^{\\infty} \\frac{A^k}{k!} $$\n",
    "\n",
    "- Про матричные функции будет отдельное занятие!\n",
    "- Важно! Чаще нужно не само значение матричной экспоненты, а результат её умножения на некоторый вектор $v$: \n",
    "\n",
    "$$ \\exp(A)v $$\n",
    "\n",
    "- Это позволяет разработать быстрые методы для вычисления произведения, не считая саму матричную экспоненту!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "72.9093498873854\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import scipy.linalg as splin\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rc(\"text\", usetex=True)\n",
    "%matplotlib inline \n",
    "\n",
    "n = 50\n",
    "A = np.random.randn(n, n)\n",
    "A = A + A.T\n",
    "print(np.linalg.norm(A))\n",
    "expA = splin.expm(A)\n",
    "# print(expA.round(4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def qr_algorithm(A, num_iter, eps):\n",
    "    T = A.copy()\n",
    "    U = np.eye(A.shape[0])\n",
    "    conv = [(T, U)]\n",
    "    for i in range(num_iter):\n",
    "        Q, R = np.linalg.qr(T)\n",
    "        T = R @ Q\n",
    "        U = U @ Q\n",
    "        conv.append((T, U))\n",
    "        if np.sum(np.abs(np.tril(T, k=-1))) < eps:\n",
    "            break\n",
    "    return T, U, conv[1:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5.017469844366763e-12\n"
     ]
    }
   ],
   "source": [
    "T, U, conv = qr_algorithm(A, 4000, 1e-6)\n",
    "print(np.linalg.norm(A - U @ T @ U.T))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.Line2D at 0x7fdca3f6ee10>"
      ]
     },
     "execution_count": 64,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPgAAAD7CAYAAABHRVmlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAL6UlEQVR4nO3cvW4b6RXG8ecEqVwxstXZiMEguQCKdxDuHWjjO5DLII19DbtNkJZ3YKyA9DHTpA2lC0gQAQu741rLyu1JoXd2aYpLDofz8X78f8DCnJEMHtA8O+c5GsrcXQDy9KuhCwDQHRocyBgNDmSMBgcyRoNvMLPJ1vGlmc3M7M1QNWEYm//mKb8PemnwFF4gM5tJ+m7jeCJJ7r6QtN5u/iGZ2VX475uNc9G+xqGuWUr1SvoqPI72fVBH5w2eygsU6rvbOPVK0jo8vpM067umXcKbb+Huc0nj0CTRvsah3q9DbRMzm8Rc7w5Rvg/q6uMKnuoLNJJ0v3H8dKA6to3182t4F46jfY3dfeHur8Ph2N1vFXG9ZjYJ/+OpjBTn+6CWX/fwHCMl/ALFJly5KxNJ7yRdKPLXOIziVaOPFG+9Z0MX0CaWbL9srZ//sUeSPg1WyQ5hrL0NV8Toufu3kl6b2WjoWn7Jjqu3FPn74JA+ruBrpfkCvZM0DY/Hkrb/4Yc2c/e34fFakb7GG3n7Vg/j+JXirXdsZmM91HYWao/9fbBXH1fwd3p4YaSIXyAzu5Q0DX9Wb8hqSbSO6UppZlfhiljVF/NrPNOXzXynSOt192t3vw6Ho3Au2vdBHdbHh03M7EphIbSVIXGkjR/n3euhcb5290Wsr3EYyf8UDi+qhVus9eamlwYHMAyWbEDGaHAgYzQ4kDEaHMgYDQ5krLcGDz8WSUZK9aZUq5RWvSnVukujBm/4Ub/UXqiU6k2pVimtelOq9ZGjGzyxj/oBRTv6Rpfwof334e6pmaRJddvkLs+ePfOXL19qtVrp/Pz8xHL7k1K9KdUqpVVvKrXe3Nz84O6PCm3yYZORjvio38uXL/X3f/xL1zcfdXnxXC/OnjR4SgD7mNn3u853smQLv05oaWbL1Wql65uP+ts//6vrm49dPB2AX9DkCr7WgY/6hQ8PzCVpOp365cVzSVL1J4B+NGnwoz8f++Lsif7y1R8aPBWAUxw9orf1+dgP95/11/f/0Yf7z03+OoAaGv1GlzY+v1vlcklc3YGO9PErm3YilwPdG6zByeVA9/iwCZCxaBqcpRvQvmganJthgPYNlsG3sXQD2hdNg7N0A9oXzYi+jUwOnC7aBieTA6eLZkTfRiYHThdtg5PJgdNFO6LvQi4HjpNUg5PLgeNEO6LvQi4HjpNUg5PLgeMkNaIDOE7SDc7SDdgv6QZn6Qbsl1QG38bSDdgv6QZn6Qbsl/SIvo1MDnwpqwYnkwNfSnpE30YmB76UVYOTyYEvZTWi70IuR8myb3ByOUqW1Yi+C7kcJcu+wcnlKFn2IzpQsuIanKUbSlJcg7N0Q0myz+DbWLqhJMU1OEs3lKTWiG5mk63jSzObmdmbbsrqD5kcOTvY4GY2k/TdxvFEktx9IWm93fypIZMjZwdHdHdfmNndxqlXkt6Hx3eSZpJuO6itF2Ry5KxJBh9Jut84ftpOKcMgkyNnxf2YDChJkwZfSzoLj0eSPrVVTCxYvCEXTRr8naRxeDyWtNj+BjO7MrOlmS1Xq9Up9Q2CxRtyUWeLfilpGv6Uu9+G8zNJ6+p4k7vP3X3q7tPz8/O2a+7c5cVz/fmPv2fxhuSZu3f6BNPp1JfLZafPAZTOzG7cfbp9niVbDWRypIoGr4FMjlQVdy96E9wMg1TR4DVwMwxSxYjeELkcKaDBGyKXIwWM6A2Ry5ECGrwhcjlSwIgOZIwGbwlLN8SIBm8JSzfEiAzeEpZuiBEN3hKWbogRI3pHyOSIAQ3eETI5YsCI3hEyOWJAg3eETI4YMKL3iFyOvtHgPSKXo2+M6D0il6NvNHiPyOXoGyM6kDEafEAs3dA1GnxALN3QNTL4gFi6oWs0+IBYuqFrjOgRIZOjbTR4RMjkaBsjekTI5GgbDR4RMjnaxogeOXI5TkGDR45cjlMwokeOXI5T0OCRI5fjFIzoQMZo8MSwdMMxDo7oZnYVHv7O3d+Gc5eS1pIm7v5td+VhW7V0k8TojoP2NriZzSQt3P3OzL4Lx/eS5O4LMxub2cTdb/soFizdcJxDI/pY0iw8vgvHr/Rw9a7OzR7/NXSlWrq9OHsydClIwN4Gd/e5u8/D4UTSUtJI4SoePO2mNNRBJsc+tZZsZjaRdFt3FDezKzNbmtlytVqdVCD240YY7FP35+CzasGmh/H8LDweSfq0/c3hqj+XpOl06qeViH3I5Nin1ha92pSHJds7SdPw5bGkRXfl4RBuhME+e0f00NDfmNn/zOxHSarG9PC1NRv0+JDLUdl7BXf3haTf7Dg/3/HtiAQ/K0eFe9EzRC5HhQbPELkcFe5FBzJGgxeApVu5aPACcDNMucjgBWDpVi4avAAs3crFiF4gMnk5aPACkcnLwYheIDJ5OWjwApHJy8GIDkmPczk5PQ80OCQ9zuXk9DwwokPS41xOTs+DuXf7C1em06kvl8tOnwMonZnduPt0+zwjOpAxGhySWLLligaHJJZsuWLJBkks2XLFkg3IAEs2nIRMniYaHLWQydNEBkctZPI00eCohQ+opIkRHY2Ry+NHg6Mxcnn8GNHRGLk8fjQ4GiOXx48RHcgYDY7WsHSLDw2O1rB0iw8ZHK1h6RYfGhytYekWH0Z0dIZMPryDDW5ms/DfNxvnLsO5N92Wh5SRyYe3d0Q3s5mkr939tZm9NbNJ9TV3X5jZ2Mwm7n7beaVIDpl8eHuv4O6+cPfX4XAcGvmVpHU4dydp1l15SFmVyV+cPRm6lGLVyuBhFK8afSTpfuPLT1uuCRkjl/erVoO7+7eSXpvZqM73m9mVmS3NbLlarU6pD5khl/frUAafSFIYze8kXelhPD8L3zKS9Gn777n7XNJcevidbK1Vi+SRy/t16OfgM0nVAm0k6d+SFpKqX+42DsdALfysvF+HRvS5pLGZXUmSu19XG/OwYV+zQQfitfcK7u5rhVF76/yjc0ATH+4/6/rmoy4vnrNt7wB3smFQLN26xb3oGBRLt27R4BgUS7duMaIjKtwI0y4aHFEhk7eLER1RIZO3iwZHVMjk7WJEBzJGgyN6LN6ao8ERPRZvzZHBET0Wb83R4Igei7fmGNGRHDJ5fTQ4kkMmr48RHckhk9dHgyM5ZPL6GNGRBXL5bjQ4skAu340RHVkgl+9GgyML5PLdGNGBjNHgyBJLtwc0OLLE0u0BGRxZYun2gAZHlli6PWBERxFKzeQ0OIpQaiZnREcRSs3kNDiKUGomZ0RHsUrI5TQ4ilVCLmdER7FKyOU0OIpVQi5nRAcyVrvBzezNxuNLM5ttngNSl+PSrVaDm9lM0lfh8USS3H0haV0dA6nLcenWJIO/kvQ+PL6TNJN021pFwEByXLodvIKb2SRcrSsjSfcbx0/bLgoYQrV0e3H2ZOhSWlNnRD/rvAogQjlk8r0NvuPqLUlr/dz0I0mfdvy9KzNbmtlytVq1USfQuxwy+aEMPjazsR4a+iws1N5JmlZfl7T9PwC5+1zSXJKm06m3Vy7Qnxwy+d4ruLtfu/t1OByFc7fST5v1dXUM5CaHTF5ri755Rd44Borz4f6zrm8+6vLieRKNz51swBFSy+Xciw4cIbVcToMDR0jtAyqM6EDGaHDgBLHfDEODAyeIfelGBgdOEPvSjQYHThD70o0RHWhRbJmcBgdaFFsmZ0QHWhRbJqfBgRbFlskZ0YGODZnLaXCgY0PmckZ0oGND5nIaHOjYkLmcER3IGA0O9KzPpRsNDvSsz6UbGRzoWZ9LNxoc6FmfSzdGdGBgXWZyGhwYWJeZnBEdGFiXmZwGBwbWZSZnRAci1FYup8GBCLWVyxnRgQi1lctpcCBCbeVyRnQgYzQ4kICmSzcaHEhA06UbGRxIQNOlGw0OJKDp0o0RHUhQ3UxOgwMJqpvJGdGBBNXN5ObunRZiZitJ30t6JumHTp+sXSnVm1KtUlr1plLrb939fPtk5w3+0xOZLd192suTtSClelOqVUqr3pRq3YUMDmSMBgcy1meDz3t8rjakVG9KtUpp1ZtSrY/0lsEB9I8RHcgYDQ5kjAYHMkaDAxmjwYGM/R8XwMN1enK8wQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.spy(T, markersize=1, precision=1e-6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.459217787080334e-11\n"
     ]
    }
   ],
   "source": [
    "cur_term = np.eye(n)\n",
    "expA_approx = np.eye(n)\n",
    "K = 200\n",
    "denom = 1.\n",
    "conv = [expA_approx]\n",
    "for i in range(1, K):\n",
    "    denom = denom * i\n",
    "    cur_term = U @ T**i @ U.T / denom\n",
    "    expA_approx = expA_approx + cur_term\n",
    "    conv.append(expA_approx)\n",
    "# print(expA_approx.round(4))\n",
    "print(np.linalg.norm(expA_approx - expA) / np.linalg.norm(expA))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoYAAAG5CAYAAAATTvWcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABAVklEQVR4nO3deXRc53nn+d8DcBMJAkWQFEhwB6mNFCUukGXHiU1KoGIncTpjk9KctHs6MyMRcpJJembSouiT02dmzumWIedMd2bmxFwcJz0z7cgi5SS2J45EUqEXKbZFUvtiSQS1kpIogiAJcQXwzB/3LapYKACFwq26tXw/59QB6t6L+z4vq1B88K7m7gIAAADqkg4AAAAA5YHEEAAAAJJIDAEAABCQGAIAAEASiSEAAAACEkMAAABIkiYkHUA1mDVrli9evLioZXz00UeaNm1aUcsoZ9S/dutfy3WXqD/1r93613LdpeLW/+DBgx+6++xc50gMY7B48WIdOHCgqGXs379f69atK2oZ5Yz61279a7nuEvWn/rVb/1quu1Tc+pvZm8OdoysZAAAAkkgMAQAAEJAYAgAAQBKJIQAAAAISQwAAAEgiMQQAAEDAcjUZzKxN0kZJhyStkbTD3XsTDQoAAKBESAyvtMvd10qSmR2QtFPSpmRDAgAAKI2qSwzNbKOkXnffm+PcGkntkroltUnqTl8XzvWkr3X3XjPrKE3UAAAAyauqxDAkcjlb+UI3cZe7b8g4tsvMut09nSj2Zv1Yj5mtcfdDRQwbAACgLFTF5BMzazOz7YqSu55hLuuUtD3r2HZJXeH75iKFBwAAUBGqIjF0925373T3HSNclp5UkulAOC5FCWUq6zzJIgAAqBlVkRiOxsxSytGamJ5xHLqZu5UjEaQbGQAA1IqaSAwVEr4Rlp5py04AQ7I4ZAILAABAtaqqyScjSOV53T1mdp8+XsfwnqJFNAb/9MoH+usXLujRnudlJpkks+icyTKORQej55ZxTThmJgsH6sx01cR6TZ1Ur4bJEzRt8gQ1TJ6g2dMna07TFDVPnaS6Oit9ZQEAQGLM3ZOOIVZmdlhSZ+ZyNWEpmoPuPiTTMTOXtCHX8jajlLNZ0mZJamlpWfvQQw+NL/ARPPbGJf3g8AXJ6iS5Lr9iLnn0JXrqV36fcVl0PON6d2lghJd+Up3U2lCneQ11Wpqq0w3N9ZozzS4nn6XW19enhoaGRMouB7Vc/1quu0T9qX/t1r+W6y4Vt/7r168/6O7tuc7VSouhpGisYVw7mYSJLjskqb293detWxfHbXNaJ+mO/fsVdxmXBgZ19sKA+i726+yFfp250K/jZy7oWO85vdVzTq99cEYvHzujJ45ekCTNbZqi31w5V7+zep5WtDaWNEncX4T6V5Jarn8t112i/tS/dutfy3WXkqt/rSSG3eFrszLWKgyTUjLP15SJ9XVqmlqnpqkTh73G3fXmibN68vAJ/dMvP9B//uc39M2fHtGK1kb9/rpl+tyNc1RPlzMAAFWhJhLDsItJt3IvR9MbFrhGDmamxbOmafGsafrdWxeq9+xF/eC5Y/rWT4/oD759SG2zpumrv3GDOpa3JB0qAAAYp1qZlSxFM4yz+9PXiJnHY5KaOklf/uQi7fmfPqu/+JdrVFdnuvv/PqC7//NTervnbNLhAQCAcajGxLBZuWchb9HQrfI6w3GMUX2d6TdWztU//NGv6f7PX68nXj+h3/jzn+j/e+5Y0qEBAIACVUViaGYpM+sys12KksIuM9tuZuldTdJrGG4J120My9J00Y08PpMm1Onezy7VY//jZ7SspUF/8O1D+tO/e14X+geSDg0AAIxRVYwxTCd9eVx3SEO3xUMMFjRP1cOdn9LXH/2ldvy4W69/0Kcd/027GqcMP7EFAACUl6poMUR5mFhfp6/+xg36T3et0oE3TurObf+s90+fTzosAACQJxJDxO53Vs/Tt37vFr3Vc1Ybtz2p906RHAIAUAlIDFEUn7l2tr59zyd18qNL+pff/JlO9F1IOiQAADAKEkMUzaoFKf3lv27Xu73n9K/+8hc6de5S0iEBAIARkBiiqG5tm6ltX16r1z44oz/4L4fUPzCYdEgAAGAYJIYounXXXa1//zsr9dPXP9S//4eXkw4HAAAMoyqWq0H5u/OWBXrlvTP61hNHdP2c6brrloVJhwQAALLQYoiS+epvXK/PXDtbf/p3L+iFd08lHQ4AAMhCYoiSmVBfpz+/a5VmTpusP/qbp/XRhf6kQwIAABlIDFFSM6ZN0n+8a5WOnPhI/8v3Xkw6HAAAkIHEECX3qaUz9Yfrl2nXwXf0/WePJh0OAAAISAyRiD++/RrdvCClf/f3L7D4NQAAZYLEEImYUF+nr2+8SX0X+vW/fv+lpMMBAAAiMUSCrm2Zrj9cf42+9+xR7X3p/aTDAQCg5pEYIlFfWbdU17VM15/+3QvqY5YyAACJIjFEoiZNqNPXvrRS750+r//z8deSDgcAgJpGYojErV44Q5vWzte3fnpEh4/3JR0OAAA1i8QQZeG+z12vKRPq9b99/yW5e9LhAABQk0gMURZmT5+sP+64Rj969bj2vfxB0uEAAFCTSAxRNv71ryzW0tnT9B9++LL6BwaTDgcAgJpDYoiyMbG+Tvd97np1H/9Iuw6+k3Q4AADUHBJDlJU7lrdozcKU/tPeV3Xu4kDS4QAAUFNIDFFWzEz3f/4GvX/6gv7qySNJhwMAQE0hMUTZ+cSSZt12/dX6xv7D6j17MelwAACoGSSGKEv/9tev05nz/frWT2k1BACgVEgMUZZumNuoz62Yo7964g2dOncp6XAAAKgJJIYoW//D7ct05kK//uoJWg0BACgFEkOUrRWtTdqwvEXf+ukRnb3EbigAABQbiSHK2h/ffo1On+/X3rfoTgYAoNhIDFHWbpzXpNuuv1qPvXGJdQ0BACgyEkOUvXs/u1R9l6TdB99OOhQAAKoaiSHK3i2LZ6itqU47f3JEA4OMNQQAoFhIDFH2zEyfXzJRb/Wc1T++8F7S4QAAULVIDFER1rbUa/HMqdrx48Nyp9UQAIBiIDFERagz03//a2169p1T+vmRnqTDAQCgKpEYomJsXDNfqakT9ddPvJF0KAAAVCUSQ1SMqybV665bFuixl97Tu73nkg4HAICqQ2KIivKvPrlIkvT//uzNhCMBAKD6kBiiosyfMVUblrfob37xls5fYsFrAADiRGKIivN7v7JEvWcv6XvPHE06FAAAqgqJISrOJ9uadV3LdP31k2+wdA0AADEiMUTFMTN9+VOL9NKx03r2nVNJhwMAQNUgMURF+p1VrZo6qV7f/jmTUAAAiAuJISrS9CkT9ds3t+r7zx7T6fOXkg4HAICqQGKIivW7ty7UuUsD+vun3006FAAAqgKJISrWynlNWtHaqP/y87eYhAIAQAxIDFGxzEy/e+tCvfLeGT39dm/S4QAAUPFIDDOY2Roz22xm95nZLjNrSzomjOxfrJqnqZPq9fBTbycdCgAAFY/EMDCzlKR2d9/h7g9K2i5pT7JRYTQNkyfo8zfO1Q+eO6ZzF9kJBQCA8Rh3Ymhmq8xscQyxjFTGRjPrGOZcupWvI/21wGLaJG3JeH5AUltIGFHGNrXPV9+Ffv3ji8eSDgUAgIoWR4vhLklrzOy2GO41REj0dg5zrk1SV2jl2+vuOyR1FtIF7O6HJG3IONQuqdfdewsIGyV065JmLWyeql0H3kk6FAAAKlocieEj7v5dd388hntdZmZtZrZdUUtezzCXdSrq8s20XVJXIWW6e3fWve8p5D4oLTPTxrXz9eThE3q752zS4QAAULEmxHCPw+lvzOweSU+5+zPjvWlI0jrDfbcMc9lGDU0MDyhjbKCZbZa0dISi9rj73swD4We+4+67xxo3kvGltfP1H/e+qkcOvaN/03Ft0uEAAFCR4kgMLy8g5+47zex1M/uau38zhnsPK4z9G9Ka6O69ZiYza3P37tC9PJb7dkjqzk4WUd7mpa7Sp5fO0u6D7+iPbrtGdXWWdEgAAFScOLqS12Y97yp2Uhg0S1EiOMz5MY8zNLM1knrSSaGZbSw4OpTcl9bO0zsnz+ngWyeTDgUAgIoUR2LYaWYDZvaomf2JpKbMk2a2KoYycknFebMwYWWfpINm5mbmKnCsIpJxx/I5umpivf6OLfIAACiIjXcrsTCu8GFJHZLuCl9d0l5FY/3WuvtXxlnGYUmdmd27oXXvoLsP6TMMSd2GYnYHh3GImyWppaVl7UMPPVSsoiRJfX19amhoKGoZ5Szf+m979rye/3BAf75+qiZUUXdyLb/+tVx3ifpT/9qtfy3XXSpu/devX3/Q3dtznRv3GEN3Ty8l80h4yMyWKEoQ75R0u6RxJYYjMbNUEkvKhLGLOySpvb3d161bV9Ty9u/fr2KXUc7yrb/P+UD/7V8/pcGWG7RuxZziB1Yitfz613LdJepP/Wu3/rVcdym5+hdl5xN3P+LuO939Dkn3F6MMSemlZZozD2YsSN0t1JxfvWaWZk6bpL9/5mjSoQAAUHFKsSXemGYF5yu0EnZr6FjDZkULU5MY1qCJ9XX6ws2t2vPy+zp9/lLS4QAAUFGKnhi6+6ki3n6voh1KMq0Jx1Gj/sWqVl3sH9Q/Pv9e0qEAAFBRStFiGIdm5Z6FvEXSpqxjnbpyz2PUmFULUlo0c6q+/xzdyQAAjMWIk0/MbJuiFrnxTV2WTNJhd78r7x+IxgpuVbQeYUpSl5ltULRTyW7p8mLWW8ysS9JT4douupFrm5npN1fO1fYfd+tE3wXNbJicdEgAAFSEERNDd7+3VIHkKLtXebT8ufshSYeKHhAqyhdubtVf7D+sH77wnr78yUVJhwMAQEWolK5kYEyunzNdS2dP0w/oTgYAIG8khqhKZqbfuqlVPz/Sow9On086HAAAKsJoYwy/JmlJTGV1u/vWmO4FjOoLN8/Vn+97Tf/w/DH93qfjehsDAFC9Rtv55AFlLSA9Dj0x3QfIy7Krp+v6OdP1g+dIDAEAyMdok09OSSrmOoRAUf3WTXP1Z4+9qmOnzmlu01VJhwMAQFljjCGq2udXzpUkPfoCi10DADAaEkNUtaWzG3RtS4N+SGIIAMCoxpUYmlmjma0ys8Vm1hhXUECcPnfjXD31Ro8+7LuQdCgAAJS1MSWGZnaPmT1mZj1mdkLSPkkPStoh6ZCZnTCzp8zsT8xscRHiBcbs8zfO0aBLj734ftKhAABQ1kablazQEvhVSasl7ZK0KUxKGelnbpd0v5nNkLTd3R+PI1igENfPma7FM6fqhy8c0+/eujDpcAAAKFujrWO4WtJmSQ+6+5F8b+ru+xS1JqZbGde4+5+NK1KgQGamX79xjv7yJ0fUe/aiUlMnJR0SAABladiuZDNbIqnN3b8ylqQwm7vvlLTTzL5Y6D2A8fr8jXPVP+ja8xLdyQAADGfYxNDdj7j7I3EU4u6n3P27cdwLKMTN85vU2jRFjzLOEACAYbFcDWqCmen2G1r009eP6/ylgaTDAQCgLBUlMTSzbxTjvsB4dCxv0flLg3ry8IdJhwIAQFmKJTHMWM/wNjO7W9KdcdwXiNMn25o1bVK99rz0QdKhAABQlvJZrmaVpGZJbZKWSkqF79PHUtk/IsljjBGIxeQJ9frMtbO17+X3NTh4o+rqLOmQAAAoK6MmhooWsO5QlOyl/yftDo99ipLD7ZJ6JPVKOiLp4bgDBeLQcUOLfvjCe3rh6CndND+VdDgAAJSVURNDd7/DzDYqSg63ZC9ubWbbwpI0mcf2xhsmEI/111+tOpP2vvQ+iSEAAFnyGmPo7rvd/V5Jm0LX8hWnc1x/fwyxAbFrnjZJ7YuatedlxhkCAJBtTJNP3P2bkk6Z2W1FigcouttvuFovHzutd3vPJR0KAABlZcyzksPC14+HWchLihEUUEwdy1skSfteZrFrAAAyFbxcjbs/E7bKO2RmjTHGBBTV0tkNaps1TXvpTgYA4ArjXscwTDyZGdYwbJQk9kVGubv9hqv1z4c/1Jnzl5IOBQCAshHLAtfp7mVJG8xsm6Sdo/0MkKSOG1p0acD1k9fYBQUAgLRYt8Rz90ck7dDQRa+BsrJ20Qylpk7UXsYZAgBwWex7Jbv7IUULXwNla0J9ndZfd7X+6ZUP1D8wmHQ4AACUhdgTw6CzSPcFYtNxQ4tOnr2kQ2/1Jh0KAABloSiJYZitDJS1z1w7SxPrje5kAACCYRNDM1sS1+xiM2s0s7vjuBcQl+lTJurWJTO1/5csWwMAgDRCYhha/Y6Y2TfMbHGhBZjZPZK2hl1TgLLy2Wtn69X3+3SUXVAAABi5K9ndn5Z0v6R7zexRM7s7n8Wsw64o28zsUUmH3X1rTPECsfrMtbMlST9+9XjCkQAAkLwJo13g7qcUJYcysy9J+qaZrZHkknol9UgySW2SZkjqlnRA0vaQWAJl69qWBs1pnKIfv3Zc//UnFiYdDgAAiRo1McwU1il8JP3czJoUJYS9knpCEglUDDPTZ6+drX944Zj6BwY1ob5YE/UBACh/4/1f0BW1Et4uaQ17JqMSfeba2Tpzvl/PvN2bdCgAACSq4MQwTCrplbRHUVfzPkWTVf6reEIDSuNXl81SnTHOEACAghLDMNZwk6S17l7v7svcvU7SBklfMbPb4gwSKKamqRO1euEM/YjEEABQ4wptMbzF3e/Inlzi7ofc/Q5FCSJQMT5zzWw99+4p9Xx0MelQAABITKGJ4VOjnD9Q4H2BRHz2utlyl37yGq2GAIDaVZIpmONZIBsohZXzmjRj6kS6kwEANa3QxHDPcNvlheMHsw53FVgOUBL1daZfvWa2fvzqhxoc9KTDAQAgEWNaxzDDG5JSZpbrf1CTpIxzpmhZG6Csffba2fr+s0f18nuntaK1KelwAAAouUJbDLslLQ0zkrMfdeFx+bkyFsUGytVnrpklSXQnAwBqVqGJ4T3ufmQM1z9QYDlAyVzdOEU3zG3Uj35JYggAqE2FJoYnRzqZPf6QPZNRKT577WwdfPOk+i70Jx0KAAAlV2hiuH2U850F3hdI1Gevna3+QdeTr3+YdCgAAJRcoYnhBjNbn+uEmW2T1FF4SEBy1i6aoWmT6hlnCACoSYUmhjskLcvsMjazVWb2Wnha8ZNNzGy0VlFUoUkT6vSppTP1k9doMQQA1J6Clqtx93slycyWmNndklKStkq60933xRdeMsysQ9Jm0SVekz69bJb2vvyB3u45qwXNU5MOBwCAkil0HcO0E5LulNQmaZO7Pz7+kIYys42Set19b45zayS1K1pCp01Sd67rxlBWSlKPpN5C74HK9itLo2Vr/vnwCRJDAEBNKagr2cy+aGa3KVro+jF3XxYdzr0byniE1rudw5xrk9Tl7jvcfa+775DUGY4XqsPdD43j51Hhrm1p0KyGyXriMN3JAIDaUugYw92KZibf7u5/JkmhC/lpM9tmZq+ONzAzawvj/NoUteDl0qmhM6S3q8At+EISWnBrI6qDmelXls7Uk4dPyJ1NewAAtaPQruRuSe3ufirzoLsfMbP7w/lxcfduhTF+ZrZlmMs2amhieEDSnvQTM9ssaekIRe1x973pLmR37y00ZlSPTy+bqe89e1Svf9Cna1qmJx0OAAAlUWhiuDs7KUxz914ze3gcMeUlJHJDWhND+TKzNnfvDt3L+eiQ1Gxm7eF5KiSVe0OSihqSHmf4xOsfkhgCAGpGQV3J7n6/mS02swfCrGRJkpndbmaN6VnLRdYcYukd5vyYxhm6++4wVnFHOpkM35MU1qAFzVO1oPkqPXH4RNKhAABQMoVOPlmtaCyeKVqqRtLlcYa3mNniOIIbRWrUKwpgZikzuy98f984J7Kggn166Sz9rPuE+gcGkw4FAICSsEIG15vZA+6+NXz/JXd/JOv8F939uzHFKDM7LKkzcxmasEzNQXe3HNe7pA3jWbYmj5g2K1rrUC0tLWsfeuihYhUlSerr61NDQ0NRyyhnSdT/Z8f6te3ZC/p3n5qitqb6kpadrZZf/1quu0T9qX/t1r+W6y4Vt/7r168/6O7tuc4VOsbwQMb3uTLL5gLvO2Zmlkpiwkjobt4hSe3t7b5u3bqilrd//34Vu4xylkT9V5y5oG3P7tX5pkVat25ZScvOVsuvfy3XXaL+1L9261/LdZeSq3+hy9VkJoNDWuwkleJ/0fTYvyuS0DApJfM8UJDZ0yfrupbpevJ1xhkCAGpDoYmhZUw68YyDq8zsUUmPjTuyUYRWwm4NHWvYrGiXFBJDjNuvLJupp97o0YX+gaRDAQCg6AqdlfyIpBlm1iNpp5m9ZmYnJO2TtKtYW+PlsFfRdniZ1ohFqhGTTy+dpQv9gzr0Zm/SoQAAUHSFthjK3b8uaYmivZIfDF/b3P2bMcWWqVm5ZyFvkbQp61hnOA6M261tzaqvMz3J9ngAgBpQ6OQTSVJY5HpfeMQqjBXcqmg9wpSkLjPboGinkt2h/F4z22JmXZKeCtd20Y2MuEyfMlEr5zXpidc/1P98x3VJhwMAQFGNKzEspjCGcNSWP3c/JOlQ0QNCzfr0spna9qNu9V3oV8Pksv2VAQBg3AruSgZqxaeXztLAoOsXR5idDACobiSGwCjWLJqhSRPq9ATL1gAAqhyJITCKKRPr1b5ohp54nQkoAIDqNuKAKTPbpmg5mLHvm5d1K0mH3f2ucd4HSMSn2mbqf9/7qnrPXlRq6qSkwwEAoChGTAzd/d5SBQKUs1vbZspdeuqNk9qwvCXpcAAAKAq6koE83DS/SZMm1Onn3YwzBABULxJDIA9TJtZr9YKUfvFGT9KhAABQNKONMfyaot1N4tDt7ltjuhdQcrcuadb/9U+v68z5S5o+ZWLS4QAAELvRVut9QNF2dHGgqQUV7da2mfo/Hn9dB948qfXXXZ10OAAAxG60ySenJJ0qUSxAWVu9MKUJdaZfHOkhMQQAVCXGGAJ5mjppgm6a38QEFABA1SIxBMbg1raZeu6dUzp3cSDpUAAAiN2wiaGZ3W5mf2JmjaUMCChnn1jSrP5B16G3TiYdCgAAsRs2MXT3fZKOSPqmmT1qZneTJKLWtS+aoTqTfn6EuVQAgOoz2uSTRyQ9Iklm9iVJu83MJe2S9LC7ny5+iED5mD5lola0Ms4QAFCd8h5j6O6PuPsdku5UtPfxbjP7jpl9sWjRAWXo1iXNevrtXl3oZ5whAKC6jHnyibufcvedIUncLGmpmT0WksTb4g8RKC+fWNKsi/2DevZtVnICAFSXcc1KDkni10OSeL+ktWZ2wMy+QZKIavWJJdGa73QnAwCqTWzL1bj7kZAktkt6UNIdGUniqrjKAZKWmjpJ18+Zzr7JAICqM9qWeAVx9yOKWhBlZqsl3WtmM9z9rmKUB5TarUuatevgO7o0MKiJ9SwHCgCoDkX/H83dn3b3e0kKUU0+sWSmzl4c0AvvMs4QAFA9aOoACnB5nCHrGQIAqsi4E0MzW2Vmi2OIBagYs6dPVtvsafoFiSEAoIrE0WK4S9IaZiGj1ty6ZKaeOtKjgUFPOhQAAGIRR2L4iLt/190fj+FeQMW4ZfEMnbnQr1ffP5N0KAAAxCKOxPBw+hszu4elaVAr2hdF4wwPvHky4UgAAIhHHInh5X40d9+paKu8u2O4L1DWFjRfpdnTJ+sg6xkCAKpEHInh2qznXe7+zRjuC5Q1M1P7ohm0GAIAqkYciWGnmQ2Y2aNm9ieSmjJP0rWMarZ20Qy9c/Kc3j99PulQAAAYt1gSQ0nNknZI+oSkr5rZCTP7TuhS7oyhDKAstS8O4wzfoNUQAFD5xp0YuvtOdz/l7o+4+53u3iypXdJeSXdK2jzeMoBytaK1UVMm1unAm4wzBABUvqLsfOLuR0LCeIfCnslANZpYX6eb56d0kHGGAIAqUIot8XaUoAwgMe2LZ+jFo6d19mJ/0qEAADAuRU8M3f1UscsAktS+qFkDg65n3u5NOhQAAMalFC2GQFVbs3CGJOkgE1AAABVuwkgnzWybookk490M1iQddve7xnkfoOw0TZ2oa1saWM8QAFDxRkwM3f3eUgUCVLK1i5r1g+eOanDQVVdnSYcDAEBB6EoGYtC+aIbOnO/Xqx+cSToUAAAKRmIIxKB9cTTOkIWuAQCVbLQxhl+TtCSmsrrdfWtM9wLKysLmqZrVMFkH3zypL39yUdLhAABQkBETQ0kPKNruLg5sDYGqZWZqXzSDHVAAABVttMknpySxDiGQh/bFM/SPL76nD06f19WNU5IOBwCAMWOMIRCTtYvCOEOWrQEAVCgSQyAmK1qbNHlCHRNQAAAVK9bE0MweNrMvZjy/3cwa4ywDKFeTJtTp5gUpHWScIQCgQsXdYrjd3b+bfuLu+xTtnALUhPZFM/Ti0dM6d3Eg6VAAABizuBPDbjNbLElm9piZPSVpQ8xlAGVr7aIZ6h90Pf8uc7YAAJVntOVqRmVmr0naK2mXuz9uZqtCd/J2SXvDzOaKYWYbM5+7++6kYkHlWbUgJUl6+q2T+sSSuFZ6AgCgNMadGEp6RNJTku40s+2SjihaFHtTBSaF9ylaiHu3maUk7ZNEYoi8zWyYrEUzp+rpt3qTDgUAgDEbd2Lo7veHbx+RJDNbIqlD0lfNbI2kPe7+lfGUEVrxet19b45zaxSNY+yW1KYosRtyXR5lpCRtdfcZkuTuvZLWjiNs1KjVC1J68vAJubvMLOlwAADI26iJoZn9W0mHFXULnx7tenc/ImlneKQTxYKZWUe416Yc59okdbn7hoxju8ys2927x1hUu6Ixkhsl9UpaI2l3AfdBjVu9cIb+7pmjOnbqvFpTVyUdDgAAectn8skOSQ9KOmlm38lcjiYfIVEcMzNrC13TbRp+O71ORWMZM22X1FVAkW2KksG9ocVxh6Q9BdwHNW71wpQk0Z0MAKg4+XQlr1W0DM3Xs0+Y2SpJW8PTpyTtyKdVMR+hpa4zlLNlmMs2amhieEAZCZ2ZbZa0dISi9oREsFtRN3RvKL83JKdttBpiLK6f06hJE+r09Fsn9Zs3zU06HAAA8pZPYrgp1xhBM1utKAl7WtJ3JC2TtM/M7nb3Z+MNc6gwJnBIa2JI6JRO6Nx9R563zJX89Y4vStSiSRPqtHJek55+uzfpUAAAGJN8EsPhZhZvVTQh5PIC1iFZe9jMNsbVcjiCZunyJJFc2pQ72cvJ3bvNrNfMUiG5TClqQaS1EGO2ekFK/8/P3tTF/kFNmsDOkwCAypDP/1g+zPEOSQ9fcWGUpN2pwsb4jVWqCPfcJGlr6H7eqhwTXoB8rF44Qxf6B/XKe8X++wgAgPjk02JoZrbY3d/IOHC7pCblmJwRWtsqco2O0Do43HjGK4TkcbMktbS0aP/+/UWMTOrr6yt6GeWs0up//tygJOnhfU+pZ9HEcd+v0uofp1quu0T9qX/t1r+W6y4lV/9RE0N3vz9sb3ePu78ZDqeTp+HWCxyulTF26a7fUpWXFsYu7pCk9vZ2X7duXVHL279/v4pdRjmrtPq7ux58ep/6Js/UunWrx32/Sqt/nGq57hL1p/61W/9arruUXP3zXeD6XkmPmJkr6sJdKmnzCOMIS9FimB7716yMSSJhbGDmeaDkzEyrF8xgAgoAoKLklRiGLtb2MBO5XdIBd38617Vmdpukg/GFOGxMvWbWraFjDZsVTYohMUSiVi9M6R9ffE8n+i5oZsPkpMMBAGBUY9oSLySDQxJCM/uSolnAt0ha4u63xBPeqPYqSlQPZRxbo+G7uIGSWb1whiTpmbd7dfsNLQlHAwDA6GJZR8PdH5G0W9FklCNm9o047puhWblnIW/R0JnDncpzAglQTCvnNam+ztgBBQBQMcbUYjiS7D2SxyuMFdyqqCUyJanLzDYo2qlkdyiz18y2mFmXop1X0nsn042MxF01qV43zJ2uZxhnCACoELElhnELM41Hbflz90O6sisZKBurF8zQ3z79rgYGXfV1FbmKEwCghrAlA1BEqxak1HehX4eP9yUdCgAAoyIxBIpo9cKUJOnpt04mGwgAAHkgMQSKaMmsaWq6aiITUAAAFYHEECgiM9PqhSkSQwBARSAxBIps9YIZevWDMzpz/lLSoQAAMKJxJ4Zm1mRmfxJHMEA1WrUwJXfp+XdOJR0KAAAjiqPFsFnRjicAcrh5fpMk6VkSQwBAmaMrGSiy1NRJWjRzqp57pzfpUAAAGBGJIVACN81P6Vl2QAEAlDkSQ6AEbp7fpKOnzuv4mQtJhwIAwLBIDIESuGl+SpLoTgYAlDUSQ6AEbpzXqDpjAgoAoLyRGAIlMHXSBF3bMp1xhgCAskZiCJTITfOb9Nw7vXL3pEMBACAnEkOgRG6an9LJs5f0zslzSYcCAEBOJIZAidwcJqA8ywQUAECZiisxtJjuA1St6+ZM16QJdYwzBACUrXEnhu5+RNKWGGIBqtqkCXVaPreRmckAgLIVS4thSA4BjOLm+U164d1TGhhkAgoAoPwwxhAooZvmp3T24oAOH+9LOhQAAIYgMQRK6OYFTZKkZxhnCAAoQySGQAm1zWpQw+QJbI0HAChLJIZACdXVmVbOa9JzTEABAJQhEkOgxG5a0KSXj53Whf6BpEMBAOAK40oMzazRzFaZ2WIza4wrKKCa3Tw/pUsDrpePnUk6FAAArjBhLBeb2T2SNklql+SSuiWdDKfbzGxGOPYdSbvd/Y34QgWqw03zowkoz73Tq1ULUskGAwBAhlETw9AS+FVJqyXtkrTJ3UccIGVmt0u6PySK29398TiCBarBvNRVmtUwSc++fUr6VNLRAADwsRETQzNbLWmzpAfHsoi1u++TtC/c4x4zW+PufzauSIEqYWa6aX6KmckAgLIz7BhDM1siqc3dvzKenU3cfaeknWb2xULvAVSbm+Y36fXjfeq70J90KAAAXDZsYujuR9z9kTgKcfdT7v7dOO4FVIOb56fkLj3PsjUAgDLCcjVAAlaGCSgvvEtiCAAoH0VJDM3sG8W4L1AtZjVM1tymKXqexBAAUEbGtFzNcMLM5TZJzeHrnZK+Ese9gWp147wmWgwBAGUln+VqVunjhG+ppJSuTAJT2T+iaI1DACNYOa9Je156X2fOX9L0KROTDgcAgLxaDB+U1KEo2bNwrDs89ilKDrdL6pHUK+mIpIfjDhSoNjfOizYLeunoad3aNjPhaAAAyCMxdPc7zGyjouRwS/bi1ma2LSxJk3lsb7xhAtXnxnnRBJTn3z1FYggAKAt5TT5x993ufq+kTaFr+YrTOa6/P4bYgKp29fQpammczDhDAEDZGNOsZHf/pqRTZnZbkeIBasrKeU3MTAYAlI0xL1cTFr5+3MxWhd1RABToxnlN6v7wI3ZAAQCUhYLXMXT3Z8JWeYfCcjUAxmjlvCa5RxNQAABI2rgXuA4TT2aa2W3pBJF9kYH8rMyYgAIAQNJi2fkk3b0saYOZbZO0c7SfASBd3ThFV09nAgoAoDzEuiWeuz8iaYeGLnoNYBgr2QEFAFAmYt8r2d0PKVr4GkAebpzXpMPH+3T2IhNQAADJij0xDDqLdF+g6tw4r0mDTEABAJSBoiSGYbYygDwwAQUAUC6GTQzNbElcs4vNrNHM7o7jXkC1aWmcrFkNk0kMAQCJG3avZHc/YmYpM/uGpC53f6OQAszsHklt7r61wBhLxszWSGoPT1OSdrt7d3IRoRaYmVbOa2QCCgAgcSN2Jbv705Lul3SvmT1qZnfns5h12BVlm5k9KulwJSSFQYe77wiPByVtSTog1IaV85r0+gdMQAEAJGvYFsM0dz+lKDmUmX1J0jdDy5pL6pXUEy5dKmmGpG5JByRtD4nluJnZRkm97r43x7l0K1+3pDZJ3bmuy1Onme1w996CgwUKkJ6A8vKxM1q7aEbS4QAAatSoiWGmsE7hI+nnZtakKBnrldQTkshYmVmHogWzN+U416aom3tDxrFdZtZdYBdwl6QjZpZuKaTFECWxcn40AeWFd0+RGAIAEjOmxDBbSARjaRXMFpK+LZIO6uNWyWydkrZnHduuKMEbkkiOxt13mFlKHy+3s1dR0gsU1ZzGKZrVMIkJKACARI0rMQzjDdsVunAlHXD3WBZjCy1+naGc4VruNmpoYnhA0p6MGDcr6uYezp5017OZ3RfGFj4Yfm7PKD8LxMLMdCM7oAAAElZwYhhmG6eTsvT4vpNmdre7/20cwY1SfiqUeUVrorv3mpnMrM3du919R57365B0KOM+O8xsqZmtCbu5AEV1Y2uTfvLahzp/aUBTJtYnHQ4AoAYVtMB1mISySdJad69z92XuXidpg6SvmNltcQY5jGYpSgSHOd82xvv1SFqTfZCkEKVy47wmDQy6XjrGDigAgGQU2mJ4i7vfkX0wJFF3mNkDkh4fV2SjS8V5M3c/ZGZtoQs5ff/vxFkGMJL0BJQX3z2lNQuZgAIAKD1z97H/kNmXwgzlgs4XUN5hSZ2Zy9CEZWoOurvluN4lbRjHsjX5xLRZ0mZJamlpWfvQQw8VqyhJUl9fnxoaGopaRjmrhfq7u/7w8bNa2zJB/92Nk684Vwv1H04t112i/tS/dutfy3WXilv/9evXH3T39lznxjX5JF9mtrjQnVPyuHcqiXUHw9jFHZLU3t7u69atK2p5+/fvV7HLKGe1Uv9Vh3+mk+f6tW7dr15xvFbqn0st112i/tS/dutfy3WXkqt/QWMMJe0Zbh/lcPxg1uGuAssZSXqdwuas8lNZ54GKsaK1Sb98/4wuDQwmHQoAoAYV2mL4hqRU6LLNZtLl7tz087H3V48izD7u1tCxhs2KdkkhMUTFWdHaqIv9gzp8vE/Xzxl190kAAGJVaItht6Sl7l6f41EXHpefK2O3lJjtVbSOYqY14ThQcZbPjZLBF99lZjIAoPQKTQzvcfcjY7j+gQLLSWtW7lnIWzR0h5NOsZUdKlTb7AZNmVinF4+SGAIASq/QruSTI500sy+6+3fTz919zNvmhbGCWxWtR5iS1GVmGxTtVLI73LfXzLaYWZekp8K1XXQjo1LV15mun9OoF4+yAwoAoPQKTQy3S/r1Ec53SvruCOdHFWYaj9ryF9ZOZBFqVI0VrY363rNH5e4yG7IaEwAARVNoV/IGM1uf64SZbZPUUXhIQG1b0dqkM+f79c7Jc0mHAgCoMYUmhjskLctcssbMVpnZa+FpsSabAFVveWuYgEJ3MgCgxArqSnb3eyXJzJaY2d2KxgBulXSnu++LLzyg9lw/Z7rq60wvHj2tz904N+lwAAA1ZLw7n5yQdKeiSR+b3L3Y+yMDVW/KxHotnT2NmckAgJIrqCvZzL5oZrcpWuj6MXdfFh3OvRsKgLFZ0dqkl0gMAQAlVugYw92KZibf7u5/JkmhC/lpM9tmZq/GFSBQi1a0Nuq90+d1ou9C0qEAAGrIeHY+ac9enzAsen2/pFnjDQyoZZd3QKHVEABQQgW3GLp7zimTYf3BhwuOCEDGzGQSQwBA6RSUGLr7/Wa22MweCLOSJUlmdruZNaZnLQMoTGrqJM1LXcWSNQCAkip08slqSXslmTL2MA7jDG8xs8VxBAfUshWtjXrpGC2GAIDSKbQr+U53X+bu90s6knkiJIdrxh0ZUONWtDbpyIcf6aML/UmHAgCoEYUmhgcyvvcc55sLvC+AYHlro9ylV96j1RAAUBqFJoaZyaDlOL+swPsCCFYwAQUAUGKFJoaWMenEMw6uMrNHJT027siAGje3aYpmTJ3IQtcAgJIpdFbyI5JmmFmPpJ1m9pqZnZC0T9IutsYDxs/MtKK1iRZDAEDJFLxXsrt/3cx2SGpXtFdyt6QDw61vCGDsVrQ26q+eeEP9K6YkHQoAoAYUnBhKUkgC94UHgJgtb23UxYFBHfso1xwvAADiVegYQwAlkJ6A8ubpgYQjAQDUAhJDoIwtmdWgKRPr9NbpwaRDAQDUABJDoIzV15mun9Oot8+QGAIAio/EEChzy1sb9daZQbkzzhAAUFwjTj4xs22KZh2P938kk3TY3e8a532AmrN8bqO+fUk6euq85qWuSjocAEAVGzExdPd7SxUIgNyWhwkoLx09TWIIACgqupKBMnf9nOkyiR1QAABFR2IIlLmpkyaoZarppWOsHQ8AKK7Rxhh+TdKSmMrqdvetMd0LqCkLG+v08rEzSYcBAKhyo+188oCk5pjK6onpPkDNWTi9Tr947axOn7+kxikTkw4HAFClRpt8ckoS/VdAwhY2RqM+Xjl2Rp9YEtffagAAXIkxhkAFWDg9+lV96Sh/pwEAiofEEKgATZNNM6dN0kvHmJkMACgeEkOgApiZlrc2khgCAIqKxBCoEMvnNurV9/p0aYB9kwEAxUFiCFSI5a2NujgwqO7jHyUdCgCgSpEYAhVi+dywNR4LXQMAioTEEKgQS2ZN06QJdWyNBwAoGhJDoEJMqK/T9XOmMwEFAFA0JIZABVk+t1EvHT0td086FABAFSIxBCrI8tZGnTx7Se+dPp90KACAKkRiCFSQ9ASUl+lOBgAUAYkhUEGuT89MZgIKAKAISAyBCtIweYIWzZzKBBQAQFGQGAIVJj0BBQCAuJEYAhVm+dxGvXHirPou9CcdCgCgypAYAhVmeWs0zvAVupMBADEjMQQqTDoxZGYyACBuJIZAhZnTOEUzpk5kAgoAIHYkhkCFMTPdwAQUAEAR1GxiaGZ7chxrM7P7zKwjfE0lEBowquVzG/XKe2fUPzCYdCgAgCoyIekASs3MOiS1SerIcXqXu68N1x2QtFPSphKGB+RleWujLvQP6siHH+malulJhwMAqBJlkxia2UZJve6+N8e5NZLaJXUrSuq6c12Xj/TPmdn2HGX0ZFzXG5JIoOykJ6C8dOw0iSEAIDZl0ZUcErCdw5xrk9Tl7jvcfa+775DUGY7HqU1Sb9axnpAwAmVl6ewGTaqvYwIKACBWiSaGYUzfdkVJWc8wl3VK2p51bLukrpjDaY75fkDRTKyv07VzGpiAAgCIVaJdye7erSjxk5ltGeayjRqaGB6QdHnyiJltlrR0hKL25NH13CMplXWMZBFl64Y5jXr8lQ/k7jKzpMMBAFSBshljmEuYFTykNTGM/5OZtbl7d+heHq9u5UgE3f1QDPcGYre8tVG7Dr6j42cu6OrGKUmHAwCoAmUxxnAEzVKUCA5zPrZxhtkJYBjDWNAEF6AUls+NJqC8yDhDAEBMzN2TjkGSZGaHJXVmdvmGiR8H3X1IP5mZuaQNY52dHO7ZoWiM4oPK6GbOOHdI0hpJO4ZLSkP39WZJamlpWfvQQw+NJYwx6+vrU0NDQ1HLKGfUf2j9z15y/f6+s9p4zUT91tJJCUVWfLz21J/612b9a7nuUnHrv379+oPu3p7rXFl3JRdDaBk8pCgpHO6cNEprYei+3iFJ7e3tvm7dungDzbJ//34Vu4xyRv1z13/B04/r3FUprVtXvZPnee2pP/Vfl3QYiajlukvJ1b/cu5IlXR5rCCDL8rmNepmuZABATMo9MewOX6+YFJKRKHYLqGE3zG3UkQ8/0tmL/UmHAgCoAmWdGIbxfd3KvYxMb1juBqhZy+c2yl165b0zSYcCAKgCZZ0YBnsVbYeXaY2YMQx8vDUeC10DAGJQTolhs4a2DErSFkmbso51huNATZuXukqNUyawNR4AIBaJzkoOYwW3KlqPMCWpy8w2KFpCZrd0eTHrLWbWJempcG0X3ciAZGZa3tpIiyEAIBZJb4nXqzxa/rKWkQGQYfncJv3NL97SwKCrvo6t8QAAhSunrmQABbhh7nSduzSgN058lHQoAIAKR2IIVDgmoAAA4kJiCFS4a66eron1xgQUAMC4kRgCFW7ShDotu3o6LYYAgHEjMQSqwPK5jbQYAgDGjcQQqALLWxt1/MwFHT9zIelQAAAVjMQQqAI3zJ0uSXqZVkMAwDiQGAJVYPncMDOZxBAAMA4khkAVSE2dpHmpq/Ts271JhwIAqGAkhkCVWHfdbO3/5XGdvdifdCgAgApFYghUiS/c3Kpzlwa09+UPkg4FAFChSAyBKvGJxc1qaZys7z1zNOlQAAAVisQQqBJ1dabfuqlVP3r1A506eynpcAAAFYjEEKgiv31zqy4NuB598b2kQwEAVCASQ6CK3DS/SYtmTtX3nqU7GQAwdiSGQBUxM33hplY9efhDvfDuqaTDAQBUGBJDoMrcdcsCNU+brC/+xZPa/qPDGhj0pEMCAFSICUkHACBeC5qn6tF/82v66t8+rwd++Ip2/qRbLY1TNLNhsibWmSTJLH31x8/Th9LnTNnXJuP48fP6zjsHkw0iQdSf+tdq/Wu17r+/bplWzm9KrHwSQ6AKzWyYrG1fXqvvP3dMP3n1uD7su6Cejy5qwF0eGhAvf5XkfmWr4sfnkm9t/OijQZ1WX9JhJIb6U/9arX+t1v2jhDcpIDEEqpSZ6bdvbtVv39yadCjjsn//fq1b99mkw0gM9af+tVr/Wq57khhjCAAAAEkkhgAAAAhIDAEAACCJxBAAAAABiSEAAAAkkRgCAAAgIDEEAACAJBJDAAAABCSGAAAAkERiCAAAgIDEEAAAAJJIDAEAABCQGAIAAEASiSEAAAACEkMAAABIkszdk46h4pnZcUlvFrmYWZI+LHIZ5Yz61279a7nuEvWn/rVb/1quu1Tc+i9y99m5TpAYVggzO+Du7UnHkRTqX7v1r+W6S9Sf+tdu/Wu57lJy9acrGQAAAJJIDAEAABCQGFaOHUkHkDDqX7tque4S9af+tauW6y4lVH/GGAIAAEASLYYAAAAISAwBAAAgSZqQdAAYmZmtkdQuqVtSm6Rud9+bbFTFYWYbFdVxafi63d13Z5xPSdosabekHknNkjol7an0f5N861aN7wcz65K0R9IBd+8d5pqUqui1D+/13lyx5/MaV/r7YJT6V/3nwHD1r5XPgRHqX9WfBaO9t8M1yf/+uzuPMn2EF3xP1rFdktqSjq0Idd0oaU3G85Skw5I2Z/17eMbjpKSNScce42s9Yt2q9f2g6D8CH+ZxuNpee0kdIf6OYd4HI77Glf4+GKX+Vf85kMfrX9WfA6PUv2o/C8bw3k7895+u5PLWKWl71rHtkroSiKXY2tz9UPqJR38tdmlo/TdImiFpqbvP8Ky/tircaHWr1vdDt6K6r1X0l3T60SlpU8Z1Ff3am1mbmW1X9MHeM8xl+bzGFfk+yLP+Vfs5kGf9pSr9HMiz/tX8WZDPe7s8fv+TzqJ5jPgXxmFl/RWg6K8MTzq2mOuZknRQUirrePovw7aM50P+yqyGRz51q9b3gzL+Yh7ueLW99uG1zNViMuprXA3vg1z1r6XPgRFe/5r4HBih/lX5WTCG93ZZ/P7TYlimwjiKIX9ZeRh3YWZtpY+qOEKd2sIDOVTz+8Hdh6zVZWabcx2vZvm8xlX+PugVnwMjqubXX6rez4J83tvl9PvP5JPy1Sx9/ILn0Kao2b0quPuMHIc7FA1QzqxnWxjAK0X/Rj1eGd0I+RipbjXzfjCzDkkHcpyq5tdeyu81Vh7XVOz7gM8BSXwOXFZNnwWjvbfTSV05/P6TGJavVNIBlIFOSQ9kPO+RJL9yhuIuM1O5fyjkYbS6pZIKLAFr3P3BrGPV/NqnpWK6ptrwOVCbnwNS9X8WZL63U3lcn88140ZXMsqSmW1W9Ffg5Q8Fd+/N0aVQ9oOu81HNdRuL0Aow5C9e/n1qE58DkqqkbmNV7Z8Fud7b5YLEsMyFMQU1JTSpd7r7hjwu71bUrZAqblSJGFK3Kq1npq2S8l2Pqypf+3zqU211zoXPgctq8XNAquLPgpHe2+Xw+09iWL7Sfyk1Zx7MeENUzTiSHLok3Z590Mzuy3FtehBuRQ+6zqNuVf9+CHVZk2v8TDW/9hnyeY2r/n2Qgc+BSE19Dkg18VmQ671dNr//JIZlKvxCdGvomIJmDR2IXTXCOldbsj8Qwl9YXTlmXaV/QSr23yOfutXI+6FDUm/2wWp+7TPl8xrXyPuAz4Er1drngFTFnwXDvbfL6fefxLC87VW07U2mNcq/eb2ihDEXXZlvbjPrMLO2cKwzxxu/Q9KhEWZplb0x1K3a3w+3KPeYoqp97XPI5zWu6vcBnwM1/zkgVelnwUjv7fC0LH7/SQzL2xZdudq7FM1i2pJALEWVsfRAyszWhEeHpE0Zv0Q9mX8phubzTkn3lDbaosinbtX+fhhpR4Rqe+2blXuGYT6vcTW8D3LWv4Y+B4Z7/Wvlc2C4+qdV3WdBnu/tsvj9t7BqNsqURZtl3yXpKUW/LIe8zDcKH6vwi31ymNPd7r4049r0JuQzFX2wdMXVfJ60fOpWze+H0MUid+8c5nxFv/bhfb5VUR3SMy73Ktr3NHPpjVFf40p8H4xW/2r/HBjD61+VnwP51j9cW1WfBWN8byf++09iCAAAAEl0JQMAACAgMQQAAIAkEkMAAAAEJIYAAACQRGIIAACAgMQQAAAAkkgMAVSR9PpnGc83hkVk8/35+8xsT/yRxacSYsxkZm1mtsfMPPv1AVB+SAwBVJPsPVQ3aGz7p96iaHutclYJMV4W9njdEJ5WTEIL1CoSQwBVIewukJ0Eto1lRwR33+TulnXfrswtuEpluHJzxVjuwk4NUnXt5wtUJRJDANWiQ9LB9JOQKPaO54bhHvcp2qS+ZJIqt4g6FG391Zt0IABGRmIIoFps0JUtUu2K9hJF8rJfGwBlisQQQLXI7jYeczISum89fL9LH298vytMnjiZdf1GMzsYzh3MnOgS7nUyfL/dzA6bWSpMxthlZifDz+3J7DIerdzMGDOOrcmY4HHQzDbnqNfJHNetybpuY4jTc8U2Dh0K4wvDv8F94dEVw70BxIjEEEDFG6bbeI27HxrHbe9RlFxKUqekpZLWZpS5UdIuRcnnJkXjG4ckUmZ2UFHr5fbQldqZdf9mZXSBj1ZuthDHwVD+BknfkbQ9xwzgVIh3e7hvW3ievk/6+e5QXrpOqeHKzkfm+MLwOt3p7g+Gem0cz70BxG9C0gEAwFiFpCczAWuWlMpaxqU963mXu+fdgujuvWaWboHsyTGJZaekB919S3i+OySBnZLSx1LhZ9OJnjKuT9flHkkHzWyju+/Oo9xsOyXtcPd0wrnXzA4pSlJ3ZdV5i7vvDuUulXSfmaVCwppO4B4Izw8pShLHq0MfTwq60913pGOJ4d4AYkZiCKDiZCRBki4nilvSkxtCl25bRhISq9AKllKUWN2Xdbon63mnRuDuh8xMGrrUTj5xdIQ4rmgddPe9ZtarqNUvMzHM/P5w1u3S5/aZ2Xck7R1ni2taOinOTArFRBSgPNGVDKAatGUlGtkJUezlha9rFXWJZj42ZV6Yq8UvjOXbFcb5ncw+X0AcuVoVu5WVbI6UjIVz6SSuS1ErZhzrDnYodE/nGtcIoLzQYgig4mR1JefsRpbUFlri0sbUlTyKy4nYWNZJlKQQZ5ui5OtAOHxw+J/IK442RV2/mdoy7p+X8O+zNowF3Cypy8zuC2MCxywjCXwgdJHfp2gc49JC7geg+EgMAVSczK7kkGwcSid9YRLFluzu5gJldwunyz8Uumq79HErWzqe1HAtcyG2DkmbMsb6pfItN4d04tepjC7rMCElpYzJJWMR4n/QzO5StNNKobLXL+xWlMin692R/ncAUB7oSgZQ6W7Jagm8vDTKeGUkNJ2h+zfzvvdI6ghdwh0Z54ddgiW0LvYqaolL7+M8JHkbpdzs6zolbQ5xbMxolds9lhZSM9sclqq5L+M+a/TxMjPpZXbG0hWca8mgzIkoJIVAmSExBFDpUlnP1yre8YUPKko2u5TRXRuSmrWh/D2KZgd3a/TZtvcoajXbFe45XKteznKzhQkdG/Tx8jOdilpMNw33M8N4WNFYwM6M+zyYMWEkpWh5mfYx3LNZGRNjwr9Zd0g6820VBVBC5u6jXwUAZShXt7GZ7clcHgbxMbPDirrB45itDKAM0WIIoJJd0W08zELXiBFJIVDdSAwBVLLsMWyxjS/ElULSzaLUQJUjMQRQyXqyZgC3qbjrF9Ysd+9lsghQ/RhjCAAAAEm0GAIAACAgMQQAAIAkEkMAAAAEJIYAAACQRGIIAACAgMQQAAAAkqT/HwPAxdrAaw38AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x504 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(10, 7))\n",
    "plt.plot([np.linalg.norm(expA - approx, \"fro\") / np.linalg.norm(expA) for approx in conv])\n",
    "plt.yscale(\"log\")\n",
    "plt.xticks(fontsize=20)\n",
    "plt.yticks(fontsize=20)\n",
    "plt.xlabel(\"\\# iterations, $k$\", fontsize=20)\n",
    "plt.ylabel(\"$\\|\\exp(A) - S_k \\|_F / \\|\\exp(A)\\|_F$\", fontsize=20)\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.45921858916028e-11\n"
     ]
    }
   ],
   "source": [
    "import scipy.special as spspec\n",
    "expA_approx = np.zeros((n, n))\n",
    "K = 200\n",
    "denom = 1.\n",
    "conv = [expA_approx]\n",
    "for i in range(K - 1, 0, -1):\n",
    "    cur_term = U @ T**i @ U.T / spspec.factorial(i)\n",
    "    expA_approx = expA_approx + cur_term\n",
    "    conv.append(expA_approx)\n",
    "expA_approx = expA_approx + np.eye(n)\n",
    "# print(expA_approx.round(4))\n",
    "print(np.linalg.norm(expA_approx - expA) / np.linalg.norm(expA))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAn8AAAG5CAYAAAD/F1LPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABBMklEQVR4nO3de4yd933n98937pzr4ZAiJVkXcjhObGtjS8ORvckGCSUN7WaxTBGblNti0aKtyZGLbhEgWVEysNgiQKGMkk2DYrHWDIVt0GYLSaTcbFl0I5N0GOe6K5K+W7ZlDiXralEiZ6gR7+S3fzy/Qx2eOTPnzLnMc3u/gAPNeZ7nPL/vj+fM0Xd+V3N3AQAAIB/a4g4AAAAAq4fkDwAAIEdI/gAAAHKE5A8AACBHSP4AAAByhOQPAAAgRzriDiAt1q9f75s2bWppGR988IH6+vpaWkaSUX/qn9f657nuEvWn/vmtfyvrfvz48Xfd/ZZK50j+arRp0yYdO3aspWUcPXpU27Zta2kZSUb9qX9e65/nukvUn/rnt/6trLuZvbrUObp9AQAAcoTkrwoz22FmM/Pz83GHAgAA0DCSvyrc/aC77xkaGoo7FAAAgIaR/AEAAOQIyR8AAECOkPwBAADkSO6WejGzEUk7JZ2QNCZpxt3nYg0KAABgleQu+ZO03923SpKZHZO0T9KueEMCAABYHalM/sxsp6Q5dz9c4dyYpHFJs5JGJM0WrwvnzhSvdfc5M5tYnagBAADil7rkLyRrFVvrQpfulLtvLzm238xm3b2YDM6VveyMmY25+4kWhg0AAJAIqZnwYWYjZjatKIE7s8Rlk5Kmy45NS5oKPw+3KDwAAIBUSE3y5+6z7j7p7jPLXFacyFHqWDguRUljoew8CSEAAMiN1CR/1ZhZQRVaBYszeUOX8KwqJHt0+QIAgLzITPKnkNQts2zLSHmSFxLCRZNGAAAAsip1Ez6WUajxut1m9qg+XOdvd8siWoEfvDmvP/n+Jb1w5nt138OssRgafHnD5b/55iUdmft+bOU3Xv/G7vDGG5d09NwPGii/oeJlDf4LNFr+669f0l8v/DC28ht9/xp59c9eu6y/u/BSg+XH+/41VP+fXdaLl37UYPnprf8rr17WiSs/ia18qbH6t5tpTVe7ejrb1dvVrjWd7err7tAtA93aMNCttb1damtrNEI0U5aSv5qE1r9iC+CyrX5mtkfSHknauHGjjh492rK4vnv6qk78/Irs9Gt1vd690Qgau0Ezind32VuvxlJ+w/E3Wn6x/m+8Ut/rmxvOystvQgAul712qs7XNlp2vDdwufTqbFzFNyH+Rl/vslMn4ys/5l8gl6STL8cbRAu1mzTUbVrbbbq9v013DLRp02CbRgpt6mwzLSwstPT/r0kWV90zl/yZWaFZO3aEySUzkjQ+Pu7btm1rxm0r2ibpk7ccVSvLSLqjR6k/9d8WdxixyHPdJeofd/29wez36nXXhSvXdPHyNV24ck3nL1/TwqWrOv3+Jf383EW98/4lvXPukt6cu6CX3nlff/XGZUlST2ebtt69ViNdV/Tb235Z6/q7m1GdVInrvc9S8lf8s3lYJWv5hYkgpedXxMx2SNoxOjraSGwAACRSo0MeOttNne1tGuzprOn60+9f0omfndXfz76nb/7ktP7m9GX9Xz86ogc/tkGP/PqItt7NIhytlpnkL+zWMavKS7nMhUWe67nvQUkHx8fHEzE2EACANLtloFufu+dWfe6eW+Xu+tP/9y/0RsftevbFn+kLP/y5Pr15WF/5xx/XvXcW4g41s7I021eKxvCNlx0bEzN6AQBIHDPTnQNteuw3Pqa/eexB/Yt/8gm98u4H+q1/8zf6l//++zp38UrcIWZSWpO/YVWe3btXi7d9mwzH62JmO8xsZn5+vt5bAACAKnq7OvTf/+pmHfmdX9d/88ub9H/8/av6z/7Xb+q7r8/FHVrmpCb5M7OCmU2Z2X5Fid+UmU2bWXH3juIaf3vDdTvDki5T9Xb5hnsedPc9Q0NDjVYBAABUMdDTqf/5N+/R1778KzIz7Xrq7/Rn33oj7rAyJTVj/oqJXQ3XlS7lAgAAUui+u9bq//kf/5H+h393Qr/97Lf12pnz+mcPfTTusDIhNS1/caHbFwCAeKzr79affukz+vzYR/SvDv1E/+boT+MOKRNI/qqg2xcAgPh0trfpD3Z+Sv/5vbfryT//sfZ9s/4F0RFJTbcvAADIp/Y207/a9Sldveb6X/6/l3TXul597p5b4w4rtWj5q4JuXwAA4tfR3qY/+uKn9Kk7hvQ7z31HJ08vxB1SapH8VUG3LwAAydDd0a6v/tOt6upo0+T/eVwLl67GHVIqkfwBAIDUuL2wRv/6v7xPs6cX9C///Q/iDieVSP4AAECq/Mroen152xY9f+J1/dXLp+MOJ3VI/qpgzB8AAMnzzx78qEbW9+kr//f3dP4y3b8rQfJXBWP+AABInp7Odj3x+V/Sa2cu6I++/pO4w0kVkj8AAJBKnxlZp//qM3fp3/7NKf3o7XNxh5MaJH8AACC1Hv3cL6qvu0N/+MKP4w4lNUj+qmDMHwAAyVXo7dIjv75Fh196R8deORN3OKlA8lcFY/4AAEi2//YfbdItA92a+vMfyd3jDifxSP4AAECq9XZ16H966KN68ZWzOvpjln6phuQPAACk3n9x/526a7hXf3z4J7T+VUHyBwAAUq+zvU27f21E33l9XsdePRt3OIlG8gcAADJh59gdKvR26um/mo07lEQj+auC2b4AAKTDmq52/dPP3K2v//DneuXdD+IOJ7FI/qpgti8AAOnxX//y3epoM/3vf3Mq7lASi+QPAABkxobBHv3mpz6i5469rvnzV+IOJ5FI/gAAQKb8d7+6SReuXNOfffuNuENJJJI/AACQKffcPqR7bh/U/uOvxR1KIpH8AQCAzNm19Q59/41zeumtc3GHkjgkfwAAIHN+896PqLPdtP/Y63GHkjgkfwAAIHOG+7o08fGN+rNvv6HLV6/HHU6ikPwBAIBM2jV+h858cFl/8eN34g4lUUj+qmCRZwAA0unXPnqLNgx068Bxun5LkfxVwSLPAACkU0d7m/7JJ2/XX/74tBYuXY07nMQg+QMAAJn1G790qy5fu65v/Iiu3yKSPwAAkFlb71qrWwa69efffyvuUBKD5A8AAGRWW5vpc/ds1F/86LQuXL4WdziJkNvkz8wOxR0DAABovd/4B7fpwpVr+sufnI47lETIXfJnZhNmtkfSRNyxAACA1vvM5mGt7e2k6zfoiDuAIjPbKWnO3Q9XODcmaVzSrKQRSbOVrqtF8XVmNt1AuAAAICU62tu0/RMb9R++97YuXb2m7o72uEOKVSJa/sxsQtK+Jc6NSJpy9xl3P+zuM5Imw3EAAICqfuMf3Kb3L13V3558L+5QYhdr8mdmI6EFbkTSmSUum5RU3ko3LWmqlbEBAIDs+OUt69TT2aa//DHj/mLt9nX3WUXJncxs7xKX7dTi5O+YpBsTNsIYvi3LFHWo3m5iAACQfj2d7fqHI+v0TSZ9JGfMXyVmVlCFVkF3nzMzmdmIu8+GrmAAAIAl/dpHb9Hv/fiHeu3Med053Bt3OLFJxJi/ZQxLUbK3xHnG/QEAgJr8+i/eIkm5X/Il6clfodk3NLMxM3s0/DwVJpsAAICMG1nfp48U1uS+69fcPe4YJElmdlLSZOnYvLDEy3F3twrXu6TtrRzLF8YS7pGkjRs3bn3mmWdaVZQkaWFhQf39/S0tI8moP/XPa/3zXHeJ+lP/1a3/n/zgkv7+zav61w/1qqNtUXqxqlpZ9wceeOC4u49XOpfoMX9FZlZYpuu3ZcJYwhlJGh8f923btrW0vKNHj6rVZSQZ9af+ea1/nusuUX/qv7r1v7j+bR390+Ma2PRJfWZk3aqVW0lc733Su31nw3+HSw+GiSCl51vGzHaY2cz8/HyriwIAAC32K6Pr1N5muR73l+jkL7T2zWrx2L9hRbuBtDz5c/eD7r5naGio1UUBAIAWG+zp1Na71uqbL5P8JdlhRVu7lRoLxwEAAFbkV0bX6QdvntP8hStxhxKLJCV/w6o8u3evpF1lxybD8Zaj2xcAgGz59OZhuUvHX11qc7Fsi3t7t0JYbmW/osRvysymzWxn8ZrQ9bs3XLczLNMytRpdvqF8un0BAMiQ++5cq8520388lc/kL+7t3eZUQwueu5+QdKLlAVVgZjsk7RgdHY2jeAAA0GRrutr1yTsK+k85Tf6S1O2bSLT8AQCQPZ/ePKzvvT6v85evxh3Kqms4+TOze81sUxNiAQAAWBWf3jysq9dd3/rZXNyhrLpmtPztlzRmZg824V4AAAAtN373WrWZcjnurxlj/p5396814T6JxJg/AACyZ6CnU/fcPqT/dOq9uENZdc1o+TtZ/MHMdpvZvU24Z2Iw5g8AgGz69OZhfetnc7p09VrcoayqZiR/fuMH932SDpjZl5pwXwAAgJb59OZhXbp6Xd97PV9r+TYj+dta9nzK3Z9uwn0TgUWeAQDIpvs3DUuSXnzlbMyRrK5mJH+TZnbNzF4ws9+VdFP/aNq7gen2BQAgm4b7unT3ul5957W5uENZVU1J/hRtzTYj6dOSvmJm75nZs6H7d7IJZQAAADTdp+4o6Duvz8UdxqpqOPlz933uPu/uz7v7w+4+LGlc0mFJD0va02gZAAAArfCpOwt6a/6ifn7uYtyhrJqW7PDh7qdCUvhZSY+1ogwAAIBG3XtnQZL07Rx1/a7G9m4zq1BGyzDhAwCA7Lrn9kF1tFmuxv21PPlz91RnTUz4AAAgu3o62/Wx2wZyNe5vNVr+AAAAEutTdxT03dfmdf26V784A5bd3s3MnlI0eaPRfw2TdNLdv9jgfQAAAJrqU3cW9O/+4880++4HGt3QH3c4Lbds8ufuj6xWIAAAAHEoTvr4zmtzuUj+6PYFAAC5tuWWfvV1tedm3B/JXxXM9gUAINva20yfvKOQm+Veqo35+31Jm5tU1qy7P96ke60adz8o6eD4+PjuuGMBAACt8ck7h/Rv//qULl+9rq6ObLeNLZv8SXpC0dZtzXCmSfcBAABoqntuH9KVa66fvrOgT9w+GHc4LVVtwse8JPo7AQBApn3itgFJ0ktvnct88pftdk0AAIAabFrXp+6ONv3wrXNxh9JyJH8AACD3Otrb9LFbB/QSyd/yzGzQzO41s01mlu02UgAAkGkfv21QP3zrnNyzvdNHtQkfNzGz3ZJ26cNdP2YlnQ2nR8xsbTj2rKQD7v5K80IFAABonU/cPqhnXnxNb5+7qNuG1sQdTstUTf5Ci95XJN0nab+kXWEiyHKveUjSYyEZnHb3bzQjWAAAgFb5+G1RJ+ZLb53Lb/JnZvdJ2iPpSXc/VetN3f2IpCPhHrvNbMzd/7ChSGNiZjsk7RgdHY07FAAA0EIfuzWa8fvDN8/pwY9tjDma1llyzJ+ZbZY04u5fXkniV87d90naZ2afr/cecXL3g+6+Z2hoKO5QAABACw30dOqu4V699Nb7cYfSUku2/IWEr+6kr+xe85K+1ox7AQAAtMrHb8v+jF+WegEAAAg+cduQTr33gc5fvhp3KC3TkuTPzL7aivsCAAC00sdvG5C79KO3s9v1u6KlXpYSZgSPKNoHeETSw5K+3Ix7N5uZjSlaqqYg6X5Je919NtagAABAIhS3dvvhm+c0dtfamKNpjVqWerlXHyZ1WxQlTaWJXqH8JYrWAEwcMytIGnf3mfB8QtIhRfUCAAA595HCGg10d+jHOW/5e1LShKKEzsKx2fA4oigBnJZ0RtKcokkiz600EDPbKWnO3Q9XOFdsrZsN5c1Wuq4GI5L2SpoJz48pWpy64O5zddwPAABkiJlpdGO/fvrOQtyhtEzV5M/dPxsSswlFXaQ3LfBsZk+F5VxKj60oMQstcPsU7R5Sfm5E0pS7by85tt/MZlfaXevuJ8xse8mhcUUJ59xK7gMAALLroxv69Y0fnY47jJapacKHux9w90ck7QrdwDedrnD9Y7Xc18xGzGxaUYvcmSUum1TUslhqWtJULWVUiK00YZyUtLue+wAAgGwa3dCvdxcuae785bhDaYkVTfhw96fNbLOZPdiMLdtCIjYpSWa2d4nLdmpx8ndM0Vg9hdfu0fLj9g6VdxOH1zzr7gdWGjcAAMiuj26Idvr46TsLGt80HHM0zbfi2b7FxZ9DC+Cye/w2KkzQWNQq6O5zZiYzG3H32eIEjhXcd0L1jxsEAAAZNrqhX5L0ckaTv7rX+XP3b4dE8ERY6qUVhkNZc0ucH1npDcPkkTPFxC+MZwQAAJAUzfhd09me2UkfDS/yHCZ7rDOzB4tJYBP38S006T6SbkweOSLpuJm5mbnqHDsIAACyqa3NtGVDn17OaPJn7s1bks/MviBpu6Rd7r5uha89KWmytCs2tNIdd3ercL1L2t7KrtswLnCPJG3cuHHrM88806qiJEkLCwvq7+9vaRlJRv2pf17rn+e6S9Sf+iez/tPfvagfn7muP9rW27IyWln3Bx544Li7j1c615QdPorc/XkzO6Umz6CNax2+MJZwRpLGx8d927ZtLS3v6NGjanUZSUb9qX9e65/nukvUn/ons/4/8J/q7174scZ/+VfV393UdOmGuOre9L193f2Eoq7VZiguy3LTaMswEaT0fMuY2Q4zm5mfb+ncFgAAkCDFSR8nM9j12/TkL5hsxk1Ca9+sFo/9G1a0OHPLkz93P+jue4aGhlpdFAAASIiPlsz4zZqWJH9hFnCzHFa0E0epsXAcAACg6e4a7lVXe1smZ/wumfyFxZybMmvXzAbN7EtVLhtW5dm9e7V427fJcLzl6PYFACB/OtrbtHl9n376zvtxh9J0SyZ/JYs5f9XMNtVbgJntlvS4uz9d4VzBzKbMbL+ixG/KzKZL194LXb97w3U7zexRRXv9trzLN5RPty8AADk0urE/ky1/y05fcfdvmdljkh43s/sk7Zf0nLufW+51YfePRyRtVpSo7Vvi/nOqoQUvTCI5Ue26VjCzHZJ2jI6OxlE8AACIyegt/foP33tLF69cU09ne9zhNE3VucvuPi/pMenGOn5Ph/X3XNKcoq3XTNFuG2sVTdA4Jmna3b/VmrBXj7sflHRwfHy8qcvXAACAZBu5pU/XXfrZmfP6hY0DcYfTNCtauMbdn5f0fPG5mQ0pSvrmFG2ZxsA4AACQCZvX90mSTr37QaaSv0Zn+7qi1r6HJI21cI9fAACAVXX3uij5e/W9D2KOpLnqTv7CRI45SYcUdQsfUTRB5LeaE1oyMNsXAIB8GlrTqeG+Lp1693zcoTRVXclfGPu3S9JWd29391F3b1O0r++XzezBZgYZJ2b7AgCQX5vW9eqVd2n5k6T73f2z5RM63P2Eu39WURIIAACQapvW9dHtG7xY5fyxOu+bOHT7AgCQX5vW9+nN+Yu6eOVa3KE0Tav29r1JI4tEx41uXwAA8mvT+uKkj+yM+6s3+Tu01NZv4fjxssNTdZYDAAAQm03reiVJr2So63dF6/yVeEVSwcy8wjmTpJJzpmhJGAAAgFQptvxladJHvS1/s5K2hJm+5Y+28LjxXCULQwMAAKTFYE+n1vV1Zarlr97kb7e7n1rB9U/UWU7smPABAEC+3b2uV69kaK2/epO/s8udLB8PmOY9fpnwAQBAvm1a30fLn6TpKucn67wvAABAomxe16e35i/qwuVsLPdSb/K33cweqHTCzJ6SNFF/SAAAAMlxd5j08bMz2ej6rTf5m5E0Wtq9a2b3mtnL4SkTPAAAQCZsXhclf6cyMuO3rqVe3P0RSTKzzWb2JUkFSY9LetjdjzQvPAAAgHjdvT5ba/3Vu85f0XuSHpY0ImmXu3+j8ZCSxcx2SNoxOjoadygAACAGxeVesrLHb13dvmb2eTN7UNFiz19399HocOVdP9KM2b4AAODO4V69duZC3GE0Rb0tfwcknZT0UHEZF3c/ErqBn5L0oLv/QrOCBAAAiNMda9fo+29kY83fRnb4GC9fvy8s/PyYpPWNBgYAAJAUd6zt1RtzF3T9evp3rK03+Tvg7hXTX3efk/Rc3REBAAAkzJ3Da3Tlmuud9y/FHUrD6kr+3P0xM9tkZk+E2b6SJDN7yMwGi7OBAQAAsuCOtdGM39fOpn+tv3onfNwn6bAkU7TMi6Ro3J+k+81sUzOCAwAASII71q6RJL2egeSv3gkfD4cZvjKzL5SeCBM/Pq9oJjAAAEDqfaQQkr8MzPitd8zfsZKfK418HK7zvoljZjvMbGZ+PhszfAAAwMr1dLZrw0C3Xj+b3+SvNOGzCuczsyIy6/wBAAAp6vrN7Zg/RQs6Fyd6eMnBe83sBUlfbzgyAACABLljbW9+W/7c/XlJa83sjKR9Zvaymb0n6Yik/Vnc5g0AAOTbHWvX6M25C7qW8rX+6t7b193/wMxmJI0r2tt3VtKxpdb/AwAASLM7h3t19brr5+cu6vYwASSN6k7+JCkkekfCAwAAILOKy728duZ8qpO/esf8pZaZjZnZhJntNLNpMxuJOyYAAJB8xYWe0z7ur6GWv5Q6Immzu8+Z2bCk/ZK2xhwTAABIuNsLPTIj+WsaM9spac7dD1c4N6ZobOGswvjCStfVaHPYf1iSztR5DwAAkDPdHe3aONCT+uVeEpH8mdmEpH2SdlU4NyJpyt23lxzbb2az7j670rJKEj9JmpS0d+URAwCAPLpj7ZrUb/G2bPJnZk8panFrdE6zSTrp7l8su/+IouTruJZuhZuUNF12bFrSlCokizUFE5W7U9KhBloQAQBAztyxdo2OvXo27jAasmzy5+6PtLLw0HI3KUlmtlQL3E4tTv6OSTpUfGJmeyRtWaaom5K8UO6TZrbHzA6VtioCAAAs5Y61vTr43bd09dp1dbSnc95sIrp9l2JmBUVj/G5qFQyTNWRmI+4+6+4zNd5vRNJOd38yHHpO0nTxPs2MHQAAZM+dw2t07brrrfmLunO4N+5w6pL0lHVYWjROr9RKl2kZkbSu7PkciR8AAKhFcX2/t+YvxhxJ/aqN+ft9SZubVNasuz++wtcUmlS2JMndD5tZIXQTS9J2SQ81swwAAJBdtw31SJLemk/vci/mvvRcDjMbUmh9a4Izy239ZmYnJU2Wjs0LS7wcd3ercL1L2t7KCRshSdwjSRs3btz6zDPPtKooSdLCwoL6+/tbWkaSUX/qn9f657nuEvWn/umq/4Wrri8fPq+Hf6FT/3ikq6F7tbLuDzzwwHF3H690rtqEj3lJse/Va2aFZbp+WyaMJZyRpPHxcd+2bVtLyzt69KhaXUaSUX/qn9f657nuEvWn/umr/8Bfv6De9R/Rtm33NHSfuOqe9DF/xbF4N7U+hokgpedbxsx2mNnM/HzsOTAAAEiA24Z69OZcert9E538hda+WS0e+zesVZqo4e4H3X3P0NBQq4sCAAApcNvQGr19Lr0TPpZM/szsITP7XTMbXM2AKjisaKHpUmPheMvR8gcAAErdNtST6tm+SyZ/7n5E0ilJT5vZC2b2pRYngsOqPLt3rxbv5LFq27LR8gcAAErdOtSjdxcu6fLV63GHUpdqEz6el/S8JJnZFyQdCLNs90t6zt3PNVJ4GLv3uKL19gqSpsxsu6IdOQ6EGObMbK+ZTUl6MVw7xdp8AAAgDrcPrZG79PNz6VzoueYdPoqJYFj+5WFFieBZSc+6+9fqKTyM6avagufuJySdqKeMRpnZDkk7RkdH4ygeAAAkzK1hrb+3U5r8rXjCh7vPu/s+d/+sojXwtpjZ183sWTN7sPkhxotuXwAAUKq40HNaZ/w2NNs3JIJ/EBLBxyRtNbNjZvbVLCaCAAAAt4Ut3t5O6aSPpi314u6nQiI4LulJSZ8tSQTvbVY5q43ZvgAAoFR/d4cGujtSO+O3Jev8hUTwsZAIzkh6xMyebUVZrUa3LwAAKHfrUE9q9/etecJHvdz9W5IeaXU5AAAAq+W2whq6fQEAAPLitsEevZnX5M/M7jWzTU2IJZEY8wcAAMqleaHnZrT87Zc0ltXZvYz5AwAA5W4v9Mhdeuf99LX+NWPM3/P1LvIMAACQRrcORcu9vDV/UXesTddCz81o+TtZ/MHMdqd5WRcAAIBaFBd6TuNyL81I/vzGD+77FG379qUm3BcAACCRisnf2ylc7qUZyd/WsudT7v50E+6bCEz4AAAA5QZ6OtXf3aE35/LZ8jdpZtfM7AUz+11JN82MSHs3MBM+AABAJWld6LkpyZ+kYUU7eXxa0lfM7D0zezZ0/042oQwAAIBEuXWwR++8fynuMFas4eTP3fe5+7y7P+/uD7v7sKRxSYclPSxpT6NlAAAAJM2GgW69cy6HyV8lYW/ffe7+WUmPtaIMAACAON0y0K3T71+Su1e/OEFWY3u3mVUoAwAAYFXdMtCty9eua/7ClbhDWZGWJ3/unuppssz2BQAAlWwYjJZ7Sdu4v9Vo+Us1ZvsCAIBKNgx0S1Lqxv0tu72bmT2laPJGo53ZJumku3+xwfsAAAAkwo3kL2X7+y6b/Ln7I6sVCAAAQJoUu31P0+0LAACQfX1d7VrT2c6YPwAAgDwwM20Y7E5d8ldtzN/vS9rcpLJm3f3xJt0LAAAgdtFCzxka8yfpCUVbtzXDmSbdBwAAIBE2DPTopbfOxR3GilSb8DEvKdcL3JnZDkk7RkdH4w4FAAAkzC0D3frLn6Sr25cxf1Wwzh8AAFjKhsFuLVy6qvOXr8YdSs1I/gAAAOp0S3+01l+alntpavJnZs+Z2edLnj9kZoPNLAMAACAp0rjFW7Nb/qbd/WvFJ+5+RNEOIQAAAJmTxi3emp38zZrZJkkys6+b2YuStje5DAAAgERI4xZvDSd/ZvaymX3VzB5091OSCqHrd1rSRJLX9jOz6bhjAAAA6bW2t0sdbZa7bt/nJR2W9LCZvSzpSUlTkk6GpWISycwmJO2JOw4AAJBebW2mWwa6U9XtW22R56rc/bHw4/OSZGabJU1I+oqZjUk65O5frnYfM9spac7dD1c4N6Zo7OCspBFFu4Usuq5WZlZQtOj0XL33AAAAkKK1/k4vZCj5M7N/LumkpMPuXnUJ69D1uy88islgtTImwvW7KpwbkTTl7ttLju03s1l3n6127yVMuPsBM6vz5QAAAJENA916/eyFuMOoWS3dvjOKunLPmtmzpUu51CIkgxWZ2UgYdzeipbd/m1Q0frDUtKKu5RULiWbdrYYAAAClbhnoSdU6f7V0+25VtITLH5SfMLN7JRUndLwoaaaW1sGi0HI3Ge61d4nLdmpx8ndM0qGSOPZI2rJMUYfc/XCxu9fd52qNEQAAYDkbBrr13geXdeXadXW2J3//jFqSv12VxuyZ2X2KkrBvSXpW0qikI2b2JXf/TjOCC8naolZBd58zM5nZiLvPuvtMjbeckDRsZsW1BwshcTzcQBcyAADIsQ2D0XIv7y5c0m1Da2KOprpakr+lZuw+rmiCxo1FnEOy9pyZ7VxJC+AyhqUo2Vvi/IiiSSA1cfcDpc/NbHoFiSMAAMAixS3e3jmXjuSvlrZJX+L4hKTnbrowStIeVp3j8SooNOk+NzGzgpk9Gn5+NEwqAQAAWLHiFm9pGfdXS8ufmdkmd3+l5MBDkoZUMu6uKHTJJnoabUhSnwyPJYUu4T2StHHjRh09erSlcS0sLLS8jCSj/tQ/r/XPc90l6k/901//0+evS5L+9sR31fFOZ82vi6vuVZM/d38sbNW2291fDYeLkzOWmjW7VGthXcysEMckjdAlPCNJ4+Pjvm3btpaWd/ToUbW6jCSj/tQ/r/XPc90l6k/901//85ev6p9/8wXdcseItm1bbv7pzeKqe62LPD8i6Xkzc0VdsVsk7VlmXF+zWv6K4/mGVbIgcxhbWHq+Zcxsh6Qdo6OjrS4KAACkUG9Xh9Z0tuu9lCz0XNN85DCjdlxRF+iTkra6+9OVrjWzByUdb0ZwobVvVovH/g0rmmzS8uTP3Q+6+56hoaFWFwUAAFJquK9L731wOe4warKi7d3c/VuKlna5iZl9QdHM2/slbXb3+5sTnqSoa3lc0omSY2NapYWaafkDAADVrO9PT/LXlJUI3f15SQcUTQA5ZWZfreM2w6o8u3evFm/7NqkPxx22FC1/AACgmnX93anp9l1Ry99yyvf0rUUYu/e4olbDgqQpM9uuaEeOA+G+c2a218ymFO0iUtzrl0WZAQBAIgz3demlt5qxxHHrNS35q0cY01e1Bc/dT+jmbt9VQ7cvAACoZl1/l95buCx3V8JXvGtOt2+W0e0LAACqWd/XrcvXruv9S1fjDqUqkj8AAIAGDfd1SZLOLCR/0gfJXxVmtsPMZubnl9riGAAA5N26/ij5e++D5E/6IPmrgm5fAABQzfr+bknSu7T8AQAAZN+Nbt8UrPXXcPJnZkNm9rvNCAYAACCNislfGtb6a0bL37CinT0yiTF/AACgmp7Odg10d9DtmwWM+QMAALUY7u/KR7cvAAAApHV9Xcz2BQAAyIvhvm69R7cvAABAPqzv79J7dPumHxM+AABALdaFMX/Xr3vcoSyL5K8KJnwAAIBaDPd169p117mLV+IOZVkkfwAAAE2wPmzxlvTlXkj+AAAAmmBdX7TFW9IXeib5AwAAaIK0bPHWrOTPmnQfAACAVLrR7Zv15M/dT0na24RYEonZvgAAoBZrU7K/b1Na/kICmEnM9gUAALXobG/T0JrO3HT7AgAA5N66/q7E7/JB8gcAANAk6/u69W4eun0BAAAQzfil2xcAACAnCr2dmrvADh8AAAC5UOjt0tz5y3JP7v6+DSV/ZjZoZvea2SYzG2xWUAAAAGm0trdTV665Prh8Le5QltSxkovNbLekXZLGJbmkWUlnw+kRM1sbjj0r6YC7v9K8UONhZjsk7RgdHY07FAAAkHCF3k5J0tz5y+rvXlGatWqqRhVa9L4i6T5J+yXtcvdlVzw2s4ckPRaSwWl3/0Yzgo2Dux+UdHB8fHx33LEAAIBkK/RGCz3Pnb+iO9bGHMwSlk3+zOw+SXskPbmShZzd/YikI+Eeu81szN3/sKFIAQAAEm5tSfKXVEsmf2a2WdKIu3+5kQLcfZ+ZDZnZ5939a43cCwAAIMmK3b5nzyd3uZclk7/Q0teUbdtCNzGJHwAAyLTSMX9JxVIvAAAATVJYk/xu35Ykf2b21VbctxnMbMrM3MzOmtkhMxuJOyYAAJANXR1t6utq19kEJ39NmYMcZgSPSBoO/31YUkNjBVvopLtb3EEAAIBsKi70nFS1LPVyrz5M6rZIKujmRK9Q/hJFawCuiJntlDTn7ocrnBtTtLbgbChzttJ1AAAAcVvbl+wt3mpp+XtS0oSihK7YYjYbHkcUJWPTks5ImlM0SeS5lQRhZhOS9ilaQLr83IikKXffXnJsv5nNuvvsSsoJCsVEU9J2SU+4+1wd9wEAAFiksKYrnbN9i9z9syFZmpC0t3yBZzN7yt33lR2rqVUuJHZ7JR1XlDxWMqkouSw1LWlKFZLFGswUkz0zO6Mogd1ax30AAAAWKfR26o25C3GHsaSaxvy5+wFJB8zsS2Z2zN2/XXq6wvWP1XjfWUXJncxs7xKX7dTi5O+YpEPFJ2a2R1GX9FIOFbuJS1v53P2EmY2ZWYHWPwAA0AyF3s50j/kr5e5Pm9lmM3twNbZsM7OCom7lm1oF3X3OzGRmI+4+6+4zNd5vTNI+d7+ppY/EDwAANMva3i7NX7ii69ddbW3Jm2O64qVe3P2Uu3/DzO4Nu4C00nAoc26J8ytdpmVWJa2IYazhgboiAwAAqKDQ26XrLr1/8WrcoVRU91Ivxa5fMzthZoPufq5pUX2o0MybhRbD2dBNLEVdxbubWQYAAMi3wpoPt3gbCjt+JIm5r3hVlsU3iVoAN0s65u7n6tnH18xOSposXcIldNMer7Qun5m5pO2tXPIlJIl7JGnjxo1bn3nmmVYVJUlaWFhQf39/S8tIMupP/fNa/zzXXaL+1D979f/2O1f1xycu6V/8wx5tKbQveV0r6/7AAw8cd/fxSueasshzcR9gM/uCmW1XNAu3aXv5xjUhI4wlnJGk8fFx37ZtW0vLO3r0qFpdRpJRf+qf1/rnue4S9af+2av/4M/O6o9P/K1GPv5L2vaLG5a8Lq66N3V7N3d/XlGyVGjSLYvr+A2XHgwTQUrPt4yZ7TCzmfn5+eoXAwCA3Ct2+yZ1xm/T9/Z19xOK1s5rxr3mFCV4hbJTw4p2A2l58ufuB919z9DQUKuLAgAAGbC2t0uSdPaDZO7y0fTkL5hs4r0OK9rardRYON5ytPwBAICVGFzTKTMldou3liR/YQzgSg2rcnfxXi3eyWMyHG85Wv4AAMBKtLeZBnuSu9DzkhM+wgze+1Y6a3eJew1Ketjdny47XpD0uKL1+gqSpsKEkUNhV5Hi8ix7zWxK0ovh2qnV6PIFAACox9reTs2dT2bL35LJn7ufMrOCmX1VUbL1Sj0FmNluSSPu/niFMuZUQwteGEd4op7yG2VmOyTtGB0djaN4AACQQkO9XTqb0Ja/Zbt93f1bkh6T9IiZvRD29h2sdtOw+8dTZvaCpJOVEr+0oNsXAACsVCpb/orcfV5RAigz+4Kkp8Piyy5pTh/uu7tF0lpFs3OPSZoOySMAAECurO3t0snTC3GHUdGKFnkO6/g9X3xuZkOKxuDNSToTEsVModsXAACs1NCaTs1lcakXd59392+5+6ksJn4S3b4AAGDl1vZ26f1LV3Xl2vW4Q1mkoe3dwvi/cUWtf7MKe/s2IzAAAIC0WtsX7fIxf+GK1vd3xxzNzepu+QuzeOcULbb8WPjvKTP7reaEBgAAkE5DCd7ira7kL0z82CVpq7u3ufuou7dJ2i7py2b2YDODjBM7fAAAgJUqbvGWxBm/9bb83e/uny2fzevuJ9z9s4qSwExgzB8AAFipG/v7Zij5e7HK+WN13hcAACD1Cr1Rt28SF3puyd6+5cxs02qUAwAAkASDPVHy9/7FqzFHsli9yd8hM/t8pRPh+PGyw1N1lgMAAJA6/T3RgirnLiSv27fepV5ekVQwM69wziSp5Jwp2g0klVjkGQAArFR7m2mguyNTLX+zkra4e3uFR1t43Hiukl1B0oYJHwAAoB6Dazp17mLyWv7qTf52u/upFVz/RJ3lAAAApNJAT0ciu33rTf7OLneyfDxg+ZIwAAAAWTfYk62Wv+kq5yfrvC8AAEAmDK7p0LkL2Rnzt93MHqh0wsyekjRRf0gAAADpl7WWvxlJo6Xdu2Z2r5m9HJ6mdoJHObZ3AwAA9Rhc05nIMX91LfXi7o9IkpltNrMvSSpIelzSw+5+pHnhxc/dD0o6OD4+vjvuWAAAQHoM9nRo4dJVXb/uamuzuMO5od51/orek/SwpBFJu9z9G42HBAAAkH6Dazp13aUPLl/VQNjxIwnq6vY1s8+b2YOKFnv+uruPRocr7/oBAACQNwPFXT4SttBzvS1/BySdlPRQcRkXdz8SuoGfkvSgu/9Cs4IEAABIm+L+vucuXNFHCmtijuZD9SZ/s5LG3f2mWRDufsrMHgvnAQAAcmtwzYfJX5LUO9v3QHniV+Tuc5KeqzsiAACADLjR8pewbt+6kj93f8zMNpnZE2G2ryTJzB4ys8HibGAAAIC8GlwTdbC+n7C1/uqd8HGfpMOSTNEyL5KicX+S7jezTc0ILglY5w8AANSjdMxfktTb7fuwu4+6+2OSTpWeCAngWMORJYS7H3T3PUNDQ3GHAgAAUiSps33rTf6OlfzsFc4P13lfAACATOhob1NvV3tmWv5KE75KS1aP1nlfAACAzEji/r71Jn9WMtHDSw7ea2YvSPp6w5EBAACk3OCaDp27kIFuX3d/XtJaMzsjaZ+ZvWxm70k6Imk/27wBAAAks+Wv7r193f0PzGxG0riivX1nJR1bav2/JDGznaXP3f1AXLEAAIDsGlzTqdPvX4o7jJvUnfxJUkj0joRHKpjZo5Jm3f2AmRUUxU7yBwAAmm6wp0MnTy/EHcZNGkr+mim0xs25++EK58YUtTDOKrQyVrquhjIKkh5397XSjd1ItjYQNgAAwJIG13QmbrZvIpI/M5uQtE/SrgrnRiRNufv2kmP7zWzW3Ve6h/C4pNlioqloPcIDddwHAACgqoGeDp27eFXuLrNKC6Ssvnpn+zaFmY2Y2bSi1rwzS1w2KWm67Ni0pKk6ihxRlPAdDi2HM5IO1XEfAACAqgZ7OnXtuuv85Wtxh3JDrC1/ocVtUpLMbO8Sl+3U4uTvmEqSNjPbI2nLMkUdCsnerKIu47lQ/lxIQEdo/QMAAM02uCZs8Xbxivq6E9HhunzyZ2ZPKeoqrbSLx0qYpJPu/sUVvSgao7eoVTAkbSombe4+U+MtKyV4cyuJCQAAoFYf7u97VbclZKfYZZM/d39ktQJZwnCIY26J88UlZmri7rNmNmdmhZBAFhS1BNLqBwAAmm5wTZRqvZ+gtf7MvdFGveYws5OSJktn8YZZvsfdfdEISTNzSdtXOus3TCCZlHRSUVfx9FLJX+hO3iNJGzdu3PrMM8+spKgVW1hYUH9/f0vLSDLqT/3zWv88112i/tQ/2/Wfnb+m3/u7i/rtsW7du+HmNrdW1v2BBx447u7jlc4lo/N5FYVEb6nxheXXziiaFKLx8XHftm1bCyOTjh49qlaXkWTUn/rntf55rrtE/al/tut/1+kF/d7f/aU2ffTj2nbfR246F1fdq435+31Jm5tU1qy7P17PC4vdtE2KAwAAYFWUTvhIimotf08ojLtrgqWWcllOsTt2WCUTM8JYvdLzLWNmOyTtGB0dbXVRAAAgYwZ6olQrSQs9V5vwMS8ptr16w6SMWUmFslPDinYDaXny5+4HJR0cHx/f3eqyAABAtnR3tKuns03nLl6NO5QbYl3kuUaHFS03U2osHG85M9thZjPz87HlwAAAIMUGezoTNds3ScnfsBa38EnR5Izybd8mVeOkjUa5+0F33zM0lJDFeQAAQKpE+/smp+Uv1tm+Yeze44rW6ytImjKz7Yp25Dgg3ej63WtmU5JeDNdOsTYfAABIg8GejlRN+GipMIO3agueu5+QdKLlAVXAhA8AANCIgZ5OzZ2/HHcYNySp2zeR6PYFAACNGFzTyYQPAACAvBjs6UjUUi8kf1Uw2xcAADSiv6dDC5do+UsNun0BAEAj+rs6dOnqdV25dj3uUCSR/AEAALRUX3c0v/aDhLT+kfwBAAC0UH/Y4u39hEz6IPmrgjF/AACgEf3Flr/LJH+pwJg/AADQiGLyt0DLHwAAQPYVx/wlZcYvyR8AAEALDfQUJ3xcizmSCMkfAABAC33Y8peMhZ5J/qpgwgcAAGhEf1cx+aPlLxWY8AEAABrR190uiQkfAAAAudDR3qaezjaWegEAAMiL/u5OFnkGAADIi/7udrZ3SwsmfAAAgEb193Swzl9aMOEDAAA0qq+L5A8AACA3+rs7mO0LAACQF/09Hcz2BQAAyIu+7g4mfAAAAOTFQHcHS70AAADkRV93hy5dva4r167HHQrJHwAAQKv1d0f7+yah65fkDwAAoMWKyV8Slnsh+auCRZ4BAECj+ntI/lKDRZ4BAECj+uj2BQAAyI/+7nZJSsSMX5I/AACAFuvv7pQkfXDpWsyRkPwBAAC0XF9o+aPbFwAAIAcGQsvf+yR/q8/MzpqZlz0ejTsuAACQXbT8xcTMCpJ2ubsVH5L2uvuTMYcGAAAyrKO9TT2dbYlY6qUj7gCKzGynpDl3P1zh3JikcUmzkkYkzVa6rhalrwtlHqgvYgAAgNr1d3eQ/BWZ2YSkfZJ2VTg3ImnK3beXHNtvZrPuPruSctx9ruQeBUnDK70HAABAPfq6O7SQ96VezGzEzKYVteadWeKySUnTZcemJU01WPzjkp5r8B4AAAA16e/uSMSYv1hb/kKr26QkmdneJS7bqcXJ3zFJh4pPzGyPpC3LFHWoQjfxhLsvVSYAAEBT9XV3JGK2byK6fZcSumYXtQq6+5yZycxG3H3W3WdWeN+J8nsCAAC00kB3h94+dzHuMBI/23dYunmsXpmROu87JmmpewIAADRdX0K6fc3d445BkmRmJyVNls3GHZN0PCzJUn69S9pez6zfsK7fumrdvqE7eY8kbdy4ceszzzyz0qJWZGFhQf39/S0tI8moP/XPa/3zXHeJ+lP//NT/T35wSSd+flX/24N9klpb9wceeOC4u49XOpfobt9WqXVdv9CdPCNJ4+Pjvm3btlaGpaNHj6rVZSQZ9af+ea1/nusuUX/qn5/6/+35l/T3b79yo75x1T3p3b6Sboz9AwAASK2+rg5dvHJdV69djzWOpCd/xTX4hksPliSDLV+jz8x2mNnM/Px8q4sCAAAZ1t8Tdbh+cOlarHEkOvkLEz1mJRXKTg0r2g2k5cmfux909z1DQ0OtLgoAAGRYf9jf9/1LV2KNI9HJX3BY0dZupcbC8Zaj5Q8AADRDf3enJFr+Sg1rcQufJO3V4m3fJsPxlqPlDwAANENfaPlbiLnlL9bZvmHs3uOK1usrSJoys+2KduQ4IN1Y0HmvmU1JejFcO8WevAAAIE0Gwpi/hZhb/uLe3m1ONbTgufsJSSdaHlAFZrZD0o7R0dE4igcAABnR1x2Sv4vxLvScpG7fRKLbFwAANENfV3G2L8kfAABA5n3Y7Uvyl2jM9gUAAM3Q192hLbf03Zj4EReSvyro9gUAAM3Q2d6mI7+zTV+8/65Y4yD5AwAAyBGSPwAAgBwh+auCMX8AACBLSP6qYMwfAADIEpI/AACAHCH5AwAAyBGSvyoY8wcAALKE5K8KxvwBAIAsIfkDAADIEZI/AACAHCH5AwAAyBGSPwAAgBwh+auC2b4AACBLSP6qYLYvAADIEpI/AACAHDF3jzuGVDCz05JebXEx6yW92+Iykoz6U/+81j/PdZeoP/XPb/1bWfe73f2WSidI/hLEzI65+3jcccSF+lP/vNY/z3WXqD/1z2/946o73b4AAAA5QvIHAACQIyR/yTITdwAxo/75luf657nuEvWn/vkVS90Z8wcAAJAjtPwBAADkCMkfAABAjnTEHQAkMxuTNC5pVtKIpFl3PxxvVK1hZjsV1XFL+O+0ux8oOV+QtEfSAUlnJA1LmpR0KO3/JrXWLYufBzObknRI0jF3n1vimoIy9t6Hz/tcpfhreZ/T/FmoUvfMfw8sVf+8fA8sU//MfxdU+3yHa+L9/Xd3HjE+wht6qOzYfkkjccfWgrrulDRW8rwg6aSkPWX/Hl7yOCtpZ9yxN/G9XrZuWf08KPqy9yUeJ7P43kuaCHWYWOKzsOz7nObPQpW6Z/57oIb3PtPfA1Xqn+nvghV8vmP9/afbN36TkqbLjk1LmoohllYbcfcTxSce/dU3pcX13y5praQt7r7Wy/5iSrlqdcvq52FWUd23KvpruPiYlLSr5LrUv/dmNmJm04q+vM8scVkt73PqPgs11j2z3wM11l/K6PdAjfXP+ndBLZ/v+H//486S8/5Q9BfBSNmxQvTWxB9fE+tZkHRcUqHsePEvvJGS54v+WszCo5a6ZfXzoJK/epc6nsX3PryflVo/qr7Paf8sVKp7nr4Hlnnvc/E9sEz9M/tdsILPd+y//7T8xSiMa1j0F5KHcRBmNrL6UbVGqNNIeKCCLH8e3H3RWlZmtqfS8ayr5X3O6meB74HqsvreF2X5u6CWz3dSfv+Z8BGvYenDN7SCEUVN5Jng7msrHJ5QNCi4tJ4jYcCsFP0bnfH0NPlXs1zdcvN5MLMJSccqnMrye19Uy/usGq5J5WeB7wFJfA/ckLXvgmqf72LiFvfvP8lfvApxB5AAk5KeKHl+RpL85pl/+81MafjFr6Ja3QpxBRaDMXd/suxYlt/7UoUmXZMlfA/k83tAysd3Qennu1DD9bVc0xC6fREbM9uj6K+5G7/47j5Xofk/8QOda5Hluq1E+Gt+0V+t/PvkE98DkjJSt5XKw3dBpc93EpD8JUDo38+V0PQ96e7ba7h8VlEXQKG1UcViUd0yWs9Sj0uqda2qzL73tdQpi/UuxffADXn8HpAy/l2w3Oc77t9/kr94Ff/iGS49WPKGZ2ZcRwVTkh4qP2hmj1a4tjjoNdUDnWuoW+Y/D6EuY5XGsmT5vS9Ty/uc+c9CwPdAJFffA1Juvgsqfb4T8ftP8hej8KGf1eL+/WEtHvycGWEdqL3lv/Thr6SpCjOZir8Aqf33qKVuOfk8TEiaKz+Y5fe+XC3vcx4+C3wP3CRv3wNSxr8Llvp8J+X3n+QvfocVbd9Saky1N4WnShj/MFX64TWzCTMbCccmK3ywJySdWGbmU+KtoG5Z/zzcr8pjfDL73i+hlvc5s58Fvgdy/z0gZfi7YLnPd3ga++8/yV/89urmVc2laGbQ3hhiaamSafsFMxsLjwlJu0p+Sc6U/sUXmrknJe1e3Whbopa6Zf3zsNzK/1l874dVeeZeLe9z2j8LFeueo++Bpd77vHwPLFX/okx+F9T4+Y7999/CqtGIkUWbN39R0ouKfiFOeAo2r16J8Mt7donTs+6+peTa4qbY6xR9eUw1o5k7CWqpW5Y/D6ErRO4+ucT51L/34bP+uKJ6FGczHla0T2fp0hVV3+e0fRaq1T3r3wMreO8z+T1Qa/3DtZn7Lljh5zvW33+SPwAAgByh2xcAACBHSP4AAAByhOQPAAAgR0j+AAAAcoTkDwAAIEdI/gAAAHKE5A9AqhTXByt5vjMsolrr6x81s0PNj6x50hBjKTMbMbNDZubl7w+A5CH5A5A25Xt+btfK9vu8X9FWUUmWhhhvCPuRbg9PU5O0AnlF8gcgNcIK+uWJ3shKVv53913ubmX3rbSRfMstVW6lGJMu7EYgZWv/WSCTSP4ApMmEpOPFJyEZnGvkhuEejyraNH3VxFVuC00o2sJqLu5AACyP5A9AmmzXzS1L44r2vUT8yt8bAAlF8gcgTcq7eFeccISuVg8/79eHG7HvDxMWzpZdv9PMjodzx0snl4R7nQ0/T5vZSTMrhAkQ+83sbHjdodLu3WrllsZYcmysZFLFcTPbU6FeZytcN1Z23c4Qp1eKrQETCuP9wr/Bo+Ex1YR7A2gikj8AqbBEF++Yu59o4La7FSWQkjQpaYukrSVl7pS0X1GCuUvReMNFyZKZHVfUCjkduj0ny+4/rJLu6mrllgtxHA/lb5f0rKTpCjNrCyHe6XDfkfC8eJ/i8wOhvGKdCkuVXYvS8X7hfXrY3Z8M9drZyL0BNF9H3AEAQCUhsSlNsoYlFcqWQBkvez7l7jW3BLr7nJkVWxLPVJg4sk/Sk+6+Nzw/EBK9SUnFY4Xw2mIyp5Lri3XZLem4me109wM1lFtun6QZdy8mlYfN7ISiRHR/WZ33uvuBUO4WSY+aWSEkpcUk7Ynw/ISiRLBRE/pwIs7D7j5TjKUJ9wbQZCR/ABKpJNGRdCMZ3FucUBC6X0dKEo2mCq1ZBUXJ06Nlp8+UPZ/UMtz9hJlJi5epqSWOiRDHTa187n7YzOYUtd6VJn+lP58su13x3BEze1bS4QZbTouKiW9p4icmfwDJRLcvgLQYKUsmypOeppcX/rtVUfdl6WNX6YWVWu7C2Lr9Ydzd2fLzdcRRqXVwVmUJ5XIJVzhXTNSmFLVGNmNdvgmFruRK4wwBJAstfwASqazbt2KXr6SR0KJWtKJu3ypuJFsrWUdQkkKcI4oSrGPh8PGlX1FTHCOKumlLjZTcvybh32drGJu3R9KUmT0axuitWEmi90Tozn5U0bjCLfXcD0DrkfwBSKTSbt+QUJwoJnZh4sLe8q7hOpV34RbLPxG6Vaf0YWtZMZ7CUi1sIbYJSbtKxt4Vai23gmJyN6mS7uUwCaSgkgkdKxHif9LMvqhoR5F6la/vN6soWS/We6L47wAgGej2BZAG95e16N1YVqRRJUnLZOiqLb3vbkkToft2ouT8ksuXhFbCOUUtasV9hxclaFXKLb9uUtKeEMfOkta1Aytp6TSzPWGZl0dL7jOmD5doKS5Rs5Ju20rL7ZRO/iDxAxKG5A9AGhTKnm9Vc8f7PakooZxSSddqSFy2hvIPKZp1O6vqs1h3K2r92h/uuVTrXMVyy4VJFNv14dItk4paPnct9ZolPKdobN5kyX2eLJmkUVC0NMv4Cu45rJLJKOHfbDYklrW2bgJYRebu1a8CgJhU6uI1s0OlS6ugeczspKIu62bMAgaQQLT8AUi6m7p4l1jsGU1E4gdkG8kfgKQrH1PWtPF+uFlIrFmYGcg4kj8ASXembGbtiFq7vl9uufscEzSA7GPMHwAAQI7Q8gcAAJAjJH8AAAA5QvIHAACQIyR/AAAAOULyBwAAkCMkfwAAADny/wOqLXXUgSEKogAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x504 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(10, 7))\n",
    "plt.plot([np.linalg.norm(expA - approx, \"fro\") / np.linalg.norm(expA) for approx in conv])\n",
    "plt.yscale(\"log\")\n",
    "plt.xticks(fontsize=20)\n",
    "plt.yticks(fontsize=20)\n",
    "plt.xlabel(\"\\# iterations, $k$\", fontsize=20)\n",
    "plt.ylabel(\"$\\|\\exp(A) - S_k \\|_F / \\|\\exp(A)\\|_F$\", fontsize=20)\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Главное по семинару\n",
    "\n",
    "- Сдвиги и ускорение методов поиска собственных значений\n",
    "- Метод Релея\n",
    "- Матричные функции для плотных матриц"
   ]
  }
 ],
 "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}