{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<h1>Jacobi vs Conjugate Gradient</h1>\n",
        "\n",
        "Lets consider solving the discretized 2D Poisson equation with iterative methods. "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "collapsed": true
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import numpy.linalg as la\n",
        "from matplotlib import pyplot as pt"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<matplotlib.image.AxesImage at 0x7f11e955ae48>"
            ]
          },
          "execution_count": 2,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQYAAAD8CAYAAACVSwr3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAETBJREFUeJzt3V2oXXV+xvHvU+N4MQpqPQ1pjNWR9CJz0RgPVhgZLMKo\nuYlCkXihYbBkLiIoTC/izMV4OS3VAaEVIoaJxWoFFXNh23GCIL1Q5xzJxLzUmhkjJsTkTC0qHbA1\n/nqx/8fsnNe93t+eDxz2Omuvlf3L4uznrPX/rf0/igjMzMb9QdMFmFn7OBjMbBEHg5kt4mAws0Uc\nDGa2iIPBzBZpPBgk3SHpPUnHJe1uup6FJJ2Q9K6kg5Jm0rorJb0m6f30eEUDde2VdFbS4bF1S9al\nkSfSMT4kaUsLan1U0ql0XA9K2jr23COp1vck3V5jnRskvS7pqKQjkh5K61t1XFeos7xjGhGNfQEX\nAb8BvgV8A/g1sKnJmpao8QRw1YJ1fwvsTsu7gb9poK7vAluAw6vVBWwF/gUQcDPwVgtqfRT46yW2\n3ZR+Di4Brks/HxfVVOc6YEtavgz4z1RPq47rCnWWdkybPmO4CTgeEb+NiP8Fnge2NVzTJLYB+9Ly\nPuCuuguIiDeATxasXq6ubcAzMfImcLmkdfVUumyty9kGPB8RX0TEB8BxRj8nlYuI0xHxTlr+HDgG\nrKdlx3WFOpeT+Zg2HQzrgY/Gvj/Jyv/BJgTwC0mzknamdWsj4nRa/hhY20xpiyxXV1uP84PpFHzv\n2OVYK2qVdC1wA/AWLT6uC+qEko5p08HQBbdExBbgTmCXpO+OPxmjc7XW3Vfe1rrGPAlcD2wGTgOP\nNVvOeZIuBV4EHo6Iz8afa9NxXaLO0o5p08FwCtgw9v3VaV1rRMSp9HgWeJnRKdiZ+VPG9Hi2uQov\nsFxdrTvOEXEmIs5FxFfAU5w/tW20VkkXM3qzPRsRL6XVrTuuS9VZ5jFtOhh+BWyUdJ2kbwDbgf0N\n1/Q1Sd+UdNn8MvA94DCjGnekzXYArzRT4SLL1bUfuD+Not8MfDp2atyIBdfidzM6rjCqdbukSyRd\nB2wE3q6pJgFPA8ci4vGxp1p1XJers9RjWsco6iojrFsZjar+Bvhx0/UsqO1bjEZzfw0cma8P+EPg\nAPA+8EvgygZqe47R6eL/MbpmfGC5uhiNmv99OsbvAtMtqPUfUy2H0g/uurHtf5xqfQ+4s8Y6b2F0\nmXAIOJi+trbtuK5QZ2nHVGknM7OvNX0pYWYt5GAws0UcDGa2iIPBzBZxMJjZIpUFQ9ZPTY7dbtxq\nXakTulOr6yxf0VorCQZJFzHq797J6JNd90ratMpuXTnoXakTulOr6yxf+4KB7n5q0sygmhucJP0l\ncEdE/FX6/j7gzyPiwbFtdnI+1W6cX3/jjTcyqdnZ2Uzbl2Fubo6pqalaXzOvrtTqOsu3VK2zs7O/\ni4iJ/gNrKqlqAhGxB9gDMD09HbOzswDMzMxk+nck4bs3zVYn6cNJt63qUiLzp7m+vkdbYvQZkcnM\nh0KWfcxsZVUFQ+FPTWZ9o8+HipkVV8mlRER8KelB4N8Yzeu4NyKOTLjv12/wrJcJ8/v60sKsmMrG\nGCLiVeDVnPsWCgczK6a1dz6Ov8GzXiL4ksKsmNYGA+T/7e/xBrNiWh0M4G6FWRNaHwwLuVthVr3O\nBEORMQeHg1k2nQkGKB4OZjaZTgUDFHuD+6zBbDKtCYasA4sekDSrTmuCocg4gMcczMrVmmCA7JcJ\nHpA0q0arggGKDSo6HMzK0bpgyPNmdbfCrFytCwYoHg5Z+azB7EKtDAbI10Fwt8KsHK0NhnnuVpjV\nr/XBAO5WmNWtE8EA9XcrzIasM8FQd7fCZw02ZJ0JBqi3W+FLChuyTgUDFO9WVP1aZn3QuWCYl+dM\nIE8rc3w/s6HobDBAsd/kDgez5XU6GHz7tFk1Oh0M4G6FWRU6HwzgboVZ2XoRDODPVpiVqTfBMK/O\nT1n67MH6qnfBAJ7sxayoXgaDuxVmxfQyGMCTvZgV0dtgAA9ImuXV62CY58lezLIZRDCAJ3sxy2JN\nkZ0lnQA+B84BX0bEtKQrgX8GrgVOAPdExH8XK7MckjIFRJE3uAckrcvKOGP4i4jYHBHT6fvdwIGI\n2AgcSN+3Qt4BybyfyvRZg3VVFZcS24B9aXkfcFcFr5Fbnaf5vqSwrioaDAH8QtKspJ1p3dqIOJ2W\nPwbWFnyN0rlbYbayQmMMwC0RcUrSHwGvSfqP8ScjIiQtebGdgmQnwDXXXFOwjHzm3+h5xgPyjld4\n7MG6oNAZQ0ScSo9ngZeBm4AzktYBpMezy+y7JyKmI2J6amqqSBmFuFthtljuYJD0TUmXzS8D3wMO\nA/uBHWmzHcArRYusmm+fNrtQkUuJtcDL6Y2xBviniPhXSb8CXpD0APAhcE/xMquV5zR//Ld/1n19\nSWFtlzsYIuK3wJ8tsf6/gNuKFNWEouFQ9WuZ1Wkwdz5Owt0KsxEHwxI82YsNnYNhGZ7sxYbMwbAM\nT/ZiQ+ZgWIEne7GhcjCswgOSNkQOhgl5shcbEgdDBr592obCwZCRuxU2BA6GjNytsCFwMOTgboX1\nnYMhJ3crrM8cDAW5W2F95GAogbsV1jcOhpJ4QNL6xMFQkrq7FT5rsCo5GEpUZ7fClxRWJQdDydyt\nsD5wMFSkznEAnz1Y2RwMFarzz+E5HKxMDoYKFX2zulthTXEwVMzdCusiB0MN3K2wrnEw1MTdCusS\nB0PNPDW9dYGDoQGe7MXazsHQAE/2Ym3nYGiIJ3uxNnMwNMgDktZWDoYW8GQv1jYOhpbwZC/WJg6G\nFnG3wtrCwdAidQ9Iulthy3EwtEzecMh7BuCzBlvKqsEgaa+ks5IOj627UtJrkt5Pj1ek9ZL0hKTj\nkg5J2lJl8X2Vt4PgboWVZZIzhp8DdyxYtxs4EBEbgQPpe4A7gY3payfwZDllDpO7FdaUVYMhIt4A\nPlmwehuwLy3vA+4aW/9MjLwJXC5pXVnFDpG7FdaEvGMMayPidFr+GFibltcDH41tdzKtW0TSTkkz\nkmbm5uZyljEMvn3a6lZ48DFGP0mZf5oiYk9ETEfE9NTUVNEyes2TvVjd8gbDmflLhPR4Nq0/BWwY\n2+7qtM4K8mQvVqe8wbAf2JGWdwCvjK2/P3UnbgY+HbvksIL82Qqry5rVNpD0HHArcJWkk8BPgJ8C\nL0h6APgQuCdt/iqwFTgO/B74fgU1D17RT1lm2X8+VDz2MCyrBkNE3LvMU7ctsW0Au4oWZavL+wYv\nsq/DYTh852NHebIXq5KDocPcrbCqOBg6zt0Kq4KDoQfcrbCyORh6xFPTW1kcDD3jyV6sDA6GnnG3\nwsrgYOghT01vRTkYesoDklaEg6HnPNmL5eFgGABP9mJZORgGwt0Ky8LBMBDuVlgWDoYBcbfCJuVg\nGBh3K2wSDoaBcrfCVuJgGLAi3Yo8r+Vw6A4Hw8AV+XN4HpDsLwfDwBX9Te7JXvrJwWCe7MUWcTAY\n4G6FXcjBYBfwZC8GDgZbgm+fNgeDLeLbp83BYEvy1PTD5mCwZblbMVwOBluRuxXD5GCwibhbMSwO\nBpuYuxXD4WCwiblbMRwOBsvEk70Mg4PBMvOAZP85GCw3T/bSX6sGg6S9ks5KOjy27lFJpyQdTF9b\nx557RNJxSe9Jur2qwq0dPDV9P01yxvBz4I4l1v8sIjanr1cBJG0CtgPfTvv8g6SLyirW2sndiv5Z\nNRgi4g3gkwn/vW3A8xHxRUR8ABwHbipQn3WAuxX9U2SM4UFJh9KlxhVp3Xrgo7FtTqZ11nPuVvRL\n3mB4Erge2AycBh7L+g9I2ilpRtLM3NxczjKsTdyt6I9cwRARZyLiXER8BTzF+cuFU8CGsU2vTuuW\n+jf2RMR0RExPTU3lKcNayt2K7ssVDJLWjX17NzDfsdgPbJd0iaTrgI3A28VKtC5yt6Lb1qy2gaTn\ngFuBqySdBH4C3CppMxDACeAHABFxRNILwFHgS2BXRJyrpnRrO0mZAmL8DZ5nXyuP2nBAp6enY2Zm\npukyrAJZ3+Dz+8zLsm+e1xoSSbMRMT3Jtr7z0SrlyV66ycFglSvaraj6tWwxB4PVJs+ZQJE/h+dw\nyM/BYLWq88/hORzyczBYrXz7dDc4GKx2npq+/RwM1gh3K9rNwWCN8Wcr2svBYI3z1PTt42CwVvBk\nL+3iYLBWcLeiXRwM1hqe7KU9HAzWKh6QbAcHg7WSJ3tploPBWsuTvTTHwWCtVueYgwckz3MwWKvl\nHZDM+6lMnzWMOBis9eo8zfclxYiDwTrB3Yp6ORisU9ytqIeDwTrH3YrqORisk3z7dLUcDNZJnuyl\nWg4G6yxP9lIdB4N1mrsV1XAwWC94spdyORisNzzZS3kcDNYbnuylPA4G6xVP9lIOB4P1jgcki3Mw\nWG/59un8HAzWa759Oh8Hg/WeuxXZrRoMkjZIel3SUUlHJD2U1l8p6TVJ76fHK9J6SXpC0nFJhyRt\nqfo/YbYSdyuym+SM4UvghxGxCbgZ2CVpE7AbOBARG4ED6XuAO4GN6Wsn8GTpVZtl5G5FNqsGQ0Sc\njoh30vLnwDFgPbAN2Jc22wfclZa3Ac/EyJvA5ZLWlV65WUbuVkwu0xiDpGuBG4C3gLURcTo99TGw\nNi2vBz4a2+1kWmfWCu5WrG7iYJB0KfAi8HBEfDb+XIziMdN5l6SdkmYkzczNzWXZ1awwdytWNlEw\nSLqYUSg8GxEvpdVn5i8R0uPZtP4UsGFs96vTugtExJ6ImI6I6ampqbz1m+XmAcnlTdKVEPA0cCwi\nHh97aj+wIy3vAF4ZW39/6k7cDHw6dslh1hqe7GV5aybY5jvAfcC7kg6mdT8Cfgq8IOkB4EPgnvTc\nq8BW4Djwe+D7pVZsVqL5cMjyGz3vpUGe12rKqsEQEf8OLHcUblti+wB2FazLrDbjHYRJ37QLuw55\n9mtzQPjOR7Okzjdq2wckHQxmY+r8c3htDgcHg9mYom/WvnQrHAxmC7hb4WAwW9LQp6Z3MJgtY8if\nrXAwmK1iiFPTOxjMJjC0yV4cDGYTGNpkLw4GswkNabIXB4NZBkMZkHQwmOXQ98leHAxmOfV5shcH\ng1kBfe1WOBjMCqh7QLKuboWDwaygvOGQ9wygjrMGB4NZCfJ2ENrarXAwmJWoL90KB4NZyfrQrXAw\nmFWg67dPOxjMKtD1yV4cDGYV6fJkLw4Gswp19bMVDgazGnRtshcHg1lNunT7tIPBrCZdmuzFwWBW\no650KxwMZjXrQrfCwWDWgCa6FVk4GMwa1NY5IR0MZg2rc8xhUg4Gs4bVPSA5CQeDWQvUPRPUahwM\nZi1R54DkalYNBkkbJL0u6aikI5IeSusflXRK0sH0tXVsn0ckHZf0nqTbS6vWbADqnOxlOWsm2OZL\n4IcR8Y6ky4BZSa+l534WEX+3oLBNwHbg28AfA7+U9KcRca6Uis0GIM9kL/OhIKnwZcaqZwwRcToi\n3knLnwPHgPUr7LINeD4ivoiID4DjwE2FqjQboCYHJDONMUi6FrgBeCutelDSIUl7JV2R1q0HPhrb\n7SQrB4mZLaHJbsXEwSDpUuBF4OGI+Ax4Erge2AycBh7L8sKSdkqakTQzNzeXZVezwWiqWzFRMEi6\nmFEoPBsRL6UXPxMR5yLiK+Apzl8unAI2jO1+dVp3gYjYExHTETE9NTVV5P9g1mtNdCsm6UoIeBo4\nFhGPj61fN7bZ3cDhtLwf2C7pEknXARuBtzNXZmYXqLNbMUlX4jvAfcC7kg6mdT8C7pW0GQjgBPAD\ngIg4IukF4CijjsYudyTMylGkW5GF6vpbeCsWIc0B/wP8rulaJnAV3agTulOr6yzfUrX+SURMdN3e\nimAAkDQTEdNN17GartQJ3anVdZavaK2+JdrMFnEwmNkibQqGPU0XMKGu1AndqdV1lq9Qra0ZYzCz\n9mjTGYOZtYSDwcwWcTCY2SIOBjNbxMFgZov8P87xQmqjOEe2AAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<matplotlib.figure.Figure at 0x7f11eb8c00b8>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "n=16\n",
        "T = 2*np.eye(n)-np.diag(np.ones(n-1),1)-np.diag(np.ones(n-1),-1)\n",
        "A = np.kron(np.eye(n),T)+np.kron(T,np.eye(n))\n",
        "pt.spy(A)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Define a right-hand side and solve the resulting system directly."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {
        "collapsed": true
      },
      "outputs": [],
      "source": [
        "h = 1/(n-1)\n",
        "b = h*np.arange(0,n*n)\n",
        "x = la.solve(A,b)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Split the matrix into its diagonal and strictly lower/upper triangular parts."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "0.0"
            ]
          },
          "execution_count": 4,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "d = np.diag(A)\n",
        "D = np.diag(d)\n",
        "L = np.tril(A,-1)\n",
        "U = np.triu(A,1)\n",
        "la.norm(A-(D+L+U))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Jacobi iteration proceeds using \n",
        "$$\\boldsymbol x^{(i+1)} = \\boldsymbol D^{-1}(\\boldsymbol b- (\\boldsymbol L+\\boldsymbol U)\\boldsymbol x^{(i)}).$$"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEACAYAAACznAEdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl0lfW97/H3NztzwhAyMCRAEkAQwYmIA6JW6ym2Xmm9\n1oLWoSKorZ2uPT12tbenXee0PV22t6e2VotAUdui1mNbalFasQgoVoITICAxIASEhJkkhEzf+0e2\nGCKQkOzsZ+/k81prL9m//Tz7+bjcyw/P+DN3R0REJCHoACIiEhtUCCIiAqgQREQkTIUgIiKACkFE\nRMJUCCIiAqgQREQkTIUgIiKACkFERMJUCCIiAkBi0AFORU5OjhcWFgYdQ0QkrqxevXq3u+e2t1xc\nFUJhYSGlpaVBxxARiStm9l5HltMhIxERAVQIIiISpkIQERFAhSAiImEqBBERAQIsBDM73cweMrOn\nzOyuoHKIiEiLiF52ambzgKuBSncf12p8CvBzIATMcff/cvf1wJ1mlgA8CjwYySytbdtby4adh0hN\nSiAtKURq+JWWHCItqeWVkphAQoJ1VwQRkZgX6fsQ5gO/pOV/8ACYWQh4ALgSqABWmdlCd3/bzK4B\n7gIei3COY6wo2823nl7T7nIpiQmkJYdITWwpi5biSDhaGqlHP2tTLEcLpmU8JenDovng+1JbrZMU\n0pE6EYk9ES0Ed19mZoVthicCZe5eDmBmjwNTgbfdfSGw0Mz+Cvw+kllam3LGIMYN6UddYxOH65s4\n3NBEXfh1uL6Jusbmln82fPjZ4YaWsSPhdQ4cbgiv08zho+s14X7qeRIT7MPiSP6wcFoXSZ/URLIy\nkumfnkRWejJZ6Un0T08+5s/JiSoWEYmcaNypnA9sa/W+AjjfzC4DrgVSgEUnWtnMZgGzAIYNG9ap\nAFkZyWRlJHdq3ZNxd440Nh+3KOrCxXO41Wd1rUrnw1JqPqak9tfWs6O+iUN1jeyrredIY/MJt5+Z\nkni0MD5aHEkt/97hEukffp+RHMJMh8ZE5KMCe3SFuy8FlnZgudnAbICSkpJO/H28+5jZ0cNG3eVw\nfRN7a+vZV1PP/toG9tXWs7+2nn1H/9zyz321DWzdW8u+mnoO1jWe8PuSQna0MI4WR3oy/dOTGZDx\n0b2QrPQk+qUlkajDXCI9XjQKYTswtNX7gvCYdEBacoj85DTy+6d1eJ3Gpmb2H274sDhqji2OlvGW\nP2/eXcPqmv3sr62nsfnEfdv36CGslpIYkP7hn7MzUyjMSac4J5OBfVO0ByISp6JRCKuAUWZWREsR\nTANuiMJ2e63EUAI5mSnkZKZ0eB13p/pI40eLo6Z1ibR8trv6CJt2VbO/tp6a+qZjvic9OURhdgZF\nuRkU52RQ1OrVPz3yh+1EJHIifdnpAuAyIMfMKoB/d/e5ZnY3sJiWy07nufu6SG5Xus7M6JOaRJ/U\nJIYOSO/wekcam6g6dIT39tRSvruGzVU1bN5dzbrtB3hu7U6aWu11DMhIpjA7naKcTIpzPyyKwuwM\n0pK777CbiHSMeWcukwlISUmJ6/HX8aO+sZlt+2rDJVHTUhi7q9m8u4ZdB48cs+yQfqkUHS2JzKN7\nFwVZaTp/IdJFZrba3UvaWy6u5kOQ+JKcmMCI3ExG5GZ+5LOaI41s3t1SFFt2f1gYC9/YccxJ8cQE\nY1h2eqvDT5kU5WRQnJtBXh+drxCJJBWCBCIjJZFx+f0Yl9/vmHF3D5/srqY8vGfxwWv5pt3HXIZ7\novMVxTmZ9EtPiva/kkjcUyFITDEzBmQkMyBjABOGDzjms+Zm5/2DdUfPU5SHi2Lt9gM8u+Z9Wl8k\nNSAj+ZgT2sU5LcVRmJ3RrZcJi8QzFYLEjYQEI79/yyW4F4/KOeaz+sZmtu6tPebw0+bd1SzfVMVT\nqyuOLmcGZwzpy+VjBnLFmDzG5/fTM6xEwnRSWXq86iONR4uirLKal8p289rWfTQ75GSm8LHRuVxx\neh4Xj8olM0V/R5Kep6MnlVUI0ivtq6nnxXeqWLKhkqUbKzlU10hSyLigOJvLx+Rx+Zg8hmdnBB1T\nJCJUCCId1NDUzOr39vHChkqWrN/Fu1U1AIzMyzxaDhOGZ+kptRK3VAginfTenhpe2FDJCxsqeaV8\nDw1NTt/URC4dnccVY/K49LTcbnlYokh3USGIRED1kUZWbKpiyfpK/rGxkt3V9SQYnDssi8tPz+OK\nMQM5bWCm7oeQmKZCEImw5mbnre0HwnsPu1i7/SAA+f3TuOL0lkNLFxRn67JWiTkqBJFutvNAHf/Y\nWMmS9ZWsKKuirqGZtKQQk0bmHC2IgX1Tg44pokIQiaa6hiZWlu/hhfUt5x627z8MwLh83fMgwVMh\niATE3XlnVzVLNuzihfWVuudBAqdCEIkRe2vqefGdlkNLL75TpXseJOpUCCIxSPc8SBBUCCJxYMvu\nD+95+Ofmj97zcNnoXM00J12mQhCJM8e75yElMYEvXjaSOy4t1uWs0mlxUQhmVgx8G+jn7te1t7wK\nQXqL5mbnzYr9zFmxmb++9T4FWWl89+qxXDl2oG6Ck1PW0UKI+IFKM5tnZpVmtrbN+BQz22hmZWZ2\nL4C7l7v7jEhnEIl3CQnGOcOyeOCGc/n97eeTlhRi1mOrufU3qyivqg46nvRQ3XHmaj4wpfWAmYWA\nB4CrgLHAdDMb2w3bFulxLhqZw6KvTub/Xj2W197bxyf+exn/9ewGao40tr+yyCmIeCG4+zJgb5vh\niUBZeI+gHngcmBrpbYv0VEmhBGZcXMSSb1zKNWfl89CL73LFT19k4Zs7iKfzgBLbonVtWz6wrdX7\nCiDfzLLN7CHgHDP71vFWNLNZZlZqZqVVVVXRyCoSs/L6pPLT68/if+66iJw+yXxlwetMm/0KG3Ye\nDDqa9ACBXuzs7nvc/U53H+HuPzrBMrPdvcTdS3Jzc6MdUSQmTRiexZ+/dDE/+Mw4Nu46xKfuX8H3\nFq7jwOGGoKNJHItWIWwHhrZ6XxAeE5FOCiUYN54/nH/ccxnTJw7lkZVbuPwnS3mydBvNzTqMJKcu\nWoWwChhlZkVmlgxMAxZGadsiPVpWRjL/+enx/OXuiynMyeCbT73FtQ++zFsV+4OOJnGmOy47XQCs\nBEabWYWZzXD3RuBuYDGwHnjS3ddFetsivdm4/H48deeF/L/rz6Ji32GmPvAS33r6LfbW1AcdTeKE\n7lQW6YEO1TXw8+c3Mf/lLWSkJHLPv5zGDROHkahnJPVKgd2YJiLB65OaxHeuHsuzX53MuPy+fPfP\n6/hfv3yJVVvaXhEu8iEVgkgPNmpgH34743x+deO5HKit57MPreRrj7/OroN1QUeTGKRCEOnhzIxP\njh/M8/dcypcvH8miNTu5/CdLmb3sXeobm4OOJzFEhSDSS6QnJ3LPv4zmb1+/hAuKs/nhog1c9fNl\nLN+kGz6lhQpBpJcpzMlg7q3nMe/WEhqbnZvmvsqdj62mYl9t0NEkYCoEkV7q8jEDWfy1S/jXT4zm\nxXequOKnL/Lz5zdR19AUdDQJiApBpBdLTQrxpY+NZMk9l/LxsQP52fPvcOXPXuTvb+/SQ/N6IRWC\niDCkf9rRuRdSE0PMfLRUcy/0QioEETnqeHMv/Pg5zb3QW6gQROQYbedeeHCp5l7oLVQIInJcmnuh\n91EhiMhJHW/uhe//RXMv9EQqBBFpV9u5F+a/vIUrfqq5F3oaFYKIdFjruReGZ2vuhZ5GhSAip+xE\ncy8cqNVhpHimQhCRTjEzrj23gH9841JmTCriD6UV3Dr/Vd3pHMdUCCLSJR/MvfDLG87h9a37+eZT\nb+ny1DgVaCGYWbGZzTWzp4LMISJdN2XcYL45ZTQL39zB/UvKgo4jndDpQjCzeWZWaWZr24xPMbON\nZlZmZvee7DvcvdzdZ3Q2g4jElrsuHcH/PreAnz3/Dgvf3BF0HDlFiV1Ydz7wS+DRDwbMLAQ8AFwJ\nVACrzGwhEAJ+1Gb929y9sgvbF5EYY2b88NpxbNtbyzf+8CYFWWmcOywr6FjSQZ3eQ3D3ZUDbCVon\nAmXhv/nXA48DU919jbtf3ealMhDpgVISQzx00wQG9U1l1qOlmmchjkT6HEI+sK3V+4rw2HGZWbaZ\nPQScY2bfOsEys8ys1MxKq6o0s5NIPBiQkcy8W0s40tjM7Y+UUq2H48WFQE8qu/sed7/T3Ue4e9tD\nSh8sM9vdS9y9JDc3N9oRRaSTRub14Vc3nsumymq+suB1mnRHc8yLdCFsB4a2el8QHhORXmjyqFy+\nf80ZvLChkh8uWh90HGlHV04qH88qYJSZFdFSBNOAGyK8DRGJI5+/YDjvVlUzd8VminMzuPH84UFH\nkhPoymWnC4CVwGgzqzCzGe7eCNwNLAbWA0+6+7rIRBWRePWdT43lY6Nz+e6f17Fi0+6g48gJWDzd\nUVhSUuKlpaVBxxCRTjhU18B1D65kx4HD/PGLkxiZlxl0pF7DzFa7e0l7y+nRFSISFX1Sk5hzSwkp\niQnMeGQV+2rqg44kbagQRCRqhg5I59c3lfD+gTru+O1q6hubg44kragQRCSqJgzP4r7rzuTVzXv5\n9h/X6EF4MSTSVxmJiLRr6tn5lFfV8PMlmxiRl8mdl44IOpKgQhCRgHzt46Mo313Dj5/bQGF2BlPG\nDQo6Uq+nQ0YiEggz477rzuSsgv58/Yk3WLv9QNCRej0VgogEJjUpxOybJzAgI5kZj6xi54G6oCP1\naioEEQlUXp9U5txSQnVdI7c/uoraej0ILygqBBEJ3OmD+/KLG87h7R0H+T9PvEmzHoQXCBWCiMSE\ny8cM5NufGstz63byk79tDDpOr6SrjEQkZtw2qZB3q6r51dJ3Kc7N5LoJBUFH6lW0hyAiMcPM+P41\nZzBpZDbfevotXt3cdlJG6U4qBBGJKUmhBH51wwSGDkjnjsdKeW9PTdCReg0VgojEnH7pScy75Twc\nuG3+Kg4cbgg6Uq+gQhCRmFSYk8FDn5/A1r213P3712ho0oPwupsKQURi1gXF2fzgM+NZvmk331u4\nTg/C62a6ykhEYtr1JUN5t6qaX79Yzsi8TL4wqSjoSD1WoIVgZqcDXwVygCXu/mCQeUQkNv3bJ8aw\nuaqG/3jmbQqzM/jYmLygI/VIXZlTeZ6ZVZrZ2jbjU8xso5mVmdm9J/sOd1/v7ncC1wOTOptFRHq2\nhATjv6edzemD+/LlBa+zYefBoCP1SF05hzAfmNJ6wMxCwAPAVcBYYLqZjTWz8Wb2TJtXXnida4C/\nAou6kEVEerj05ETm3FJCenKIGfNLqTp0JOhIPU6nC8HdlwFt7xqZCJS5e7m71wOPA1PdfY27X93m\nVRn+noXufhVwY2eziEjvMLhfGnNuKWFPzRFmPVZKXUNT0JF6lEhfZZQPbGv1viI8dlxmdpmZ3W9m\nv+YEewhmNsvMSs2stKqqKrJpRSTunFnQn59dfzavb93PN596S1ceRVCgJ5XdfSmwtJ1lZgOzAUpK\nSvRfXkS4avxg/vUTo7lv8UZG5Gby1Y+PCjpSjxDpQtgODG31viA8JiISUV+8bATvVlXzs+ffoSg3\ng2vOGhJ0pLgX6UNGq4BRZlZkZsnANGBhhLchIoKZ8aNrx3NeYRbf+MObvLZ1X9CR4l5XLjtdAKwE\nRptZhZnNcPdG4G5gMbAeeNLd10UmqojIsVISQ/z6phIG9U1l1qOlVOyrDTpSXLN4OiFTUlLipaWl\nQccQkRhTVnmIz/zqZfL7p/HUXReRmaKHMLRmZqvdvaS95fQsIxGJeyPz+vDADeeyqbKaryx4nSZN\nwdkpKgQR6REuOS2X711zBi9sqOSHi9YHHScuab9KRHqMmy4YzruV1cxdsZni3AxuPH940JHiigpB\nRHqU73zqdLbsqeG7f17H8AEZXDwqJ+hIcUOHjESkR0kMJfCL6ecwIjeDu363mrLK6qAjxQ0Vgoj0\nOH1Sk5h7y3kkhxKY8cgq9tXUBx0pLqgQRKRHGjogndk3T+D9A3Xc8dvV1DdqCs72qBBEpMeaMHwA\n9113Jq9u3su3/7hGD8Jrh04qi0iPNvXsfN6tquH+JZsYkZfJnZeOCDpSzFIhiEiP9/WPj6K8qpof\nP7eBwuwMpowbFHSkmKRDRiLS45kZP/nsWZxZ0J+vP/EGa7cfCDpSTFIhiEivkJoU4uGbJ5CVnsSM\nR1ax80Bd0JFijgpBRHqNvD6pzL31PKrrGrnjsVKa9cyjY6gQRKRXOX1wX/7j0+N4s+IAz6/fFXSc\nmKJCEJFe55qzhlCQlcbDy8uDjhJTVAgi0uskhhK4bVIRq7bs00xrragQRKRXuv68ofRNTWSO9hKO\nCrQQzOwyM1tuZg+Z2WVBZhGR3iUzJZEbLxjOc2t3snWPpt6Ers2pPM/MKs1sbZvxKWa20czKzOze\ndr7GgWogFajobBYRkc649aJCQgnGvJc2Bx0lJnRlD2E+MKX1gJmFgAeAq4CxwHQzG2tm483smTav\nPGC5u18F/Bvw/S5kERE5ZQP7pnLNWfk8sWob+2v1RNROF4K7LwP2thmeCJS5e7m71wOPA1PdfY27\nX93mVenuHzx+cB+Q0tksIiKdNfOSIg43NPG7f24NOkrgIn0OIR/Y1up9RXjsuMzsWjP7NfAY8MsT\nLDPLzErNrLSqqiqiYUVExgzqyyWn5TL/5S0caWwKOk6gAj2p7O5Pu/sd7v45d196gmVmu3uJu5fk\n5uZGOaGI9AazJhdTdegIf35jR9BRAhXpQtgODG31viA8JiISsyaNzGbMoD48vKy8V8+ZEOlCWAWM\nMrMiM0sGpgELI7wNEZGIMjNmXVLMpspqlr7Tew9Nd+Wy0wXASmC0mVWY2Qx3bwTuBhYD64En3X1d\nZKKKiHSfq88cwqC+qb36RrVOT5Dj7tNPML4IWNTpRCIiAUhOTOALkwr50bMbWLv9AOPy+wUdKer0\n6AoRkbBpE4eRkRzqtXsJKgQRkbB+aUlMmziMv7z1Pjv2Hw46TtSpEEREWvnCpEIA5r+8JdAcQVAh\niIi0UpCVzqfGD+b3/9zKwbqGoONElQpBRKSNmZOLqT7SyBOvbmt/4R5EhSAi0sb4gn5cUDyAeS9t\npqGpuf0VeggVgojIccy6pJj3D9Tx17feDzpK1KgQRESO47LT8hiZl8nDy3vP4yxUCCIix5GQYMyc\nXMS6HQdZ+e6eoONEhQpBROQEpp6dT05mMrN7yY1qKgQRkRNITQpxy4WFLN1YxTu7DgUdp9upEERE\nTuLzFwwnNSmhVzzOQoUgInISWRnJXF8ylD+9voPKg3VBx+lWKgQRkXbcNqmIhuZmHlm5Jego3UqF\nICLSjsKcDD4xdhC/fWUrtfWNQcfpNioEEZEOmHlJMQcON/CH0oqgo3QbFYKISAdMGJ7FhOFZzFlR\nTlNzz7xRLdBCMLPJZvaQmc0xs5eDzCIi0p6Zk4vYtvcwi9ftDDpKt+jKnMrzzKzSzNa2GZ9iZhvN\nrMzM7j3Zd7j7cne/E3gGeKSzWUREouHKsYMYnp3O7GU983EWXdlDmA9MaT1gZiHgAeAqYCww3czG\nmtl4M3umzSuv1ao3AL/vQhYRkW4XSjBuv7iIN7btZ/V7+4KOE3GdLgR3XwbsbTM8EShz93J3rwce\nB6a6+xp3v7rNqxLAzIYBB9y9598GKCJx77oJQ8lKT2L2sp53o1qkzyHkA61nlKgIj53MDOA3J/rQ\nzGaZWamZlVZVVUUgoohI56Ulh/j8BcP5+/pdlFdVBx0nogK/ysjd/93dT3hC2d1nu3uJu5fk5uZG\nM5qIyHHdfGEhSQkJzF2xOegoERXpQtgODG31viA8JiLSY+T2SeHac/N5anUFe6qPBB0nYiJdCKuA\nUWZWZGbJwDRgYYS3ISISuNsnF3GksZnfvrI16CgR05XLThcAK4HRZlZhZjPcvRG4G1gMrAeedPd1\nkYkqIhI7Rub14fIxeTy6cgt1DU1Bx4mIrlxlNN3dB7t7krsXuPvc8Pgidz/N3Ue4+w8iF1VEJLbM\nnFzMnpp6nn6tZxwZD/yksohIvLqgeADj8/sxZ3k5zT3gcRYqBBGRTjIzZl5STPnuGl7YUBl0nC5T\nIYiIdMEnxw0iv39aj5h3WYUgItIFiaEEvjCpkFc37+WNbfuDjtMlKgQRkS6aNnEYfVITeTjO9xJU\nCCIiXZSZksgN5w/j2TXvs21vbdBxOk2FICISAV+4qIgEM+a9FL+Ps1AhiIhEwKB+qVxz1hCeWLWN\nA7UNQcfpFBWCiEiE3D65mNr6Jn736ntBR+kUFYKISISMHdKXyaNymP/SFuobm4OOc8pUCCIiETRz\ncjGVh46w8M0dQUc5ZSoEEZEImjwqhzGD+vBwHM67rEIQEYkgM+P2ycVs3HWIZZt2Bx3nlKgQREQi\n7JqzhjCwbwpz4uxGNRWCiEiEJScmcOtFRSzftJu3dxwMOk6HqRBERLrBDROHkZ4ciqu9BBWCiEg3\n6JeexOfOG8rCN3fw/oHDQcfpkMAKwczGmtmTZvagmV0XVA4Rke5y26Qimt2Z/9KWoKN0SKcKwczm\nmVmlma1tMz7FzDaaWZmZ3dvO11wF/MLd7wJu7kwOEZFYNnRAOp8cP5jf/3Mrh+pi/3EWnd1DmA9M\naT1gZiHgAVr+Rz8WmB7eCxhvZs+0eeUBjwHTzOw+ILvz/woiIrFr5uRiDh1p5IlV24KO0q7Ezqzk\n7svMrLDN8ESgzN3LAczscWCqu/8IuPoEX/WlcJE83ZkcIiKx7qyh/ZlYNIDfvLSFWy4qJCkUu6du\nI5ksH2hdgRXhseMys0Izmw08Ctx3kuVmmVmpmZVWVVVFLKyISLTMmlzM9v2HWbTm/aCjnFRgVeXu\nW9x9lrvf6O4rTrLcbHcvcfeS3NzcaEYUEYmIy8fkUZybwcPLY/txFpEshO3A0FbvC8JjIiK9WkKC\ncfvFxazdfpBXyvcGHeeEIlkIq4BRZlZkZsnANGBhBL9fRCRuXXtuPtkZyTE973JnLztdAKwERptZ\nhZnNcPdG4G5gMbAeeNLd10UuqohI/EpNCnHzhYW8sKGSsspDQcc5rk4VgrtPd/fB7p7k7gXuPjc8\nvsjdT3P3Ee7+g8hGFRGJbzddOJyUxATmLI/NeZdj9/onEZEeZkBGMp8tKeDp17ZTeagu6DgfoUIQ\nEYmiGRcX09DczGMrY2/eZRWCiEgUFeVkcOXpA3nslfeorW8MOs4xVAgiIlE265Ji9tc28D+rK4KO\ncgwVgohIlE0YnsU5w/ozZ8Vmmppj50Y1FYKISJSZGTMnF/Penlr+/vbOoOMcpUIQEQnAJ84YxNAB\nacxeFjs3qqkQREQCEAo/zuK1rftZ/V5sPM5ChSAiEpDPlhTQLy2Jh5fFxo1qKgQRkYCkJyfy+QuG\nsfjtnWzZXRN0HBWCiEiQbrmwkKSEBOauCH4vQYUgIhKgvL6pfPqcIfxh9Tb21dQHmkWFICISsNsn\nF1PX0MxvXwn2cRYqBBGRgJ02sA+Xjc7lkZVbqGtoCiyHCkFEJAbMmlzM7up6/vR6cBNNqhBERGLA\nhSOyOWNIXx5eXk5zQI+zUCGIiMQAM2PWJcW8W1XD0ncqA8kQtUIws2Izm2tmT51sTESkt/rk+MEM\n7pca2OMsOlQIZjbPzCrNbG2b8SlmttHMyszs3pN9h7uXu/uM9sZERHqrpFACt00q4pXyvbxVsT/q\n2+/oHsJ8YErrATMLAQ8AVwFjgelmNtbMxpvZM21eeRFNLSLSQ02bOJQ+KYk8HMC8yx0qBHdfBrR9\n+tJEoCz8t/x64HFgqruvcfer27yCOSAmIhJn+qQmMf38YSxa8z4V+2qjuu2unEPIB7a1el8RHjsu\nM8s2s4eAc8zsWycaO856s8ys1MxKq6qquhBXRCQ+3HpRIQb85qUtUd1uYrQ25O57gDvbGzvOerOB\n2QAlJSWxM7WQiEg3GdI/javPHMzjr27lK1eMol9aUlS225U9hO3A0FbvC8JjIiLSRbdPLqamvokF\nr26N2ja7UgirgFFmVmRmycA0YGFkYomI9G7j8vsxaWQ281/aQn1jc1S22dHLThcAK4HRZlZhZjPc\nvRG4G1gMrAeedPd13RdVRKR3mTm5mJ0H63jmrR1R2V6HziG4+/QTjC8CFkU0kYiIAHDpabmcNjCT\n2cvK+cw5+ZhZt25Pj64QEYlRZsbtk4vZsPMQK8p2d/v2VAgiIjFs6tlDGDYgnc1RmGIzapediojI\nqUtJDLHknktJCnX/39+1hyAiEuOiUQagQhARkTAVgoiIACoEEREJUyGIiAigQhARkTAVgoiIACoE\nEREJM/f4mWLAzKqA96K82X7AgTjaTme/51TW6+iy7S3X2c9zgO6/jz9yovEbiqffT0eX7+oyJ/ss\nnn5DkfhvO9zdc9tdyt31OskLmB1P2+ns95zKeh1dtr3lOvs5UBr07yKI/7bR2EY0fj8dXb6ry7Tz\nWdz8hqL1/yB31yGjDvhLnG2ns99zKut1dNn2luvq5/EiGv8e8fT76ejyXV1Gv59TFFeHjEQAzKzU\n3UuCziHxS7+h49MegsSj2UEHkLin39BxaA9BREQA7SGIiEiYCkFERAAVgoiIhGnGNIl7ZvZp4FNA\nX2Cuu/8t4EgSR8zsdOCrtNystsTdHww4UmC0hyAxyczmmVmlma1tMz7FzDaaWZmZ3Qvg7n9y95nA\nncDngsgrseUUfz/r3f1O4HpgUhB5Y4UKQWLVfGBK6wEzCwEPAFcBY4HpZja21SLfCX8uMp9T+P2Y\n2TXAX4FF0Y0ZW1QIEpPcfRmwt83wRKDM3cvdvR54HJhqLX4MPOvur0U7q8SeU/n9hJdf6O5XATdG\nN2ls0TkEiSf5wLZW7yuA84EvAx8H+pnZSHd/KIhwEvOO+/sxs8uAa4EUevkeggpB4p673w/cH3QO\niU/uvhRYGnCMmKBDRhJPtgNDW70vCI+JdIR+P+1QIUg8WQWMMrMiM0sGpgELA84k8UO/n3aoECQm\nmdkCYCUw2swqzGyGuzcCdwOLgfXAk+6+LsicEpv0++kcPdxOREQA7SGIiEiYCkFERAAVgoiIhKkQ\nREQEUCFFvksmAAAAH0lEQVSIiEiYCkFERAAVgoiIhKkQREQEUCGIiEjY/weENUt5d+LD5QAAAABJ\nRU5ErkJggg==\n",
            "text/plain": [
              "<matplotlib.figure.Figure at 0x7f11e92ed588>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "def jacobi(niter,x0):\n",
        "    xi = x0\n",
        "    for i in range(niter):\n",
        "        xi = np.diag(1./d)@(b-(L+U)@xi)\n",
        "    return xi\n",
        "\n",
        "niters = np.asarray(2**np.arange(4,12),dtype=int)\n",
        "\n",
        "x0 = np.random.random(n*n)\n",
        "\n",
        "jacobi_results = []\n",
        "err = []\n",
        "\n",
        "for niter in niters:\n",
        "    jacobi_results.append(jacobi(niter,x0.copy()))\n",
        "    err.append(la.norm(jacobi_results[-1]-x))\n",
        "\n",
        "pt.plot(niters,err)\n",
        "pt.yscale('log')\n",
        "pt.xscale('log')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 16,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEACAYAAACznAEdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4lPW9/vH3JzsECFvYEvY9BK01olJAqwWDImi1LXZv\nqR5t1Zb219aeI1q32vb0lCMejv5opXSxolULiApYtcYFlUBVwh42SWQJBAKEJYR8zh8ZIESQkEzm\nmUzu13XNFeebZ7nxSnLn+3yfzJi7IyIiEhd0ABERiQ4qBBERAVQIIiISokIQERFAhSAiIiEqBBER\nAVQIIiISokIQERFAhSAiIiEqBBERASAh6ABno2PHjt6rV6+gY4iINClLly7d6e7pZ9quSRVCr169\nyM/PDzqGiEiTYmab67KdLhmJiAigQhARkRAVgoiIACoEEREJUSGIiAigQhARkZAmddtpfe0/XMkb\n63bSoVUS7Vom0SE1ibQWicTFWdDRRESiRrMohE07y7n5L0tPGoszaNcyiXapSbRPTaJ9yyTat6ou\ni3Ytk46XR/vUE4+UxPiA/gUiIo2vWRRCv06tmH/bCErLK9h9oIJd+0MfyyvYXV79cX3JfpZsqh6v\n8lMfp2VS/EkFUbNI2tcqj/apSbRJ0SxERJqOZlEIKYnxZGek1Wnbqiqn7OARSg9UUFp+ojxKy09+\n7Npfwbrt+yktr+DgkaOnPFZ8nNGuZeKpSyS1enbSITWZdqmJxz8mJ2gWIiLBaBaFcDbi4ox2oR/W\nfc/4yh/VDlYcpfTAidnGsY+l5YcpLT9CaflhdpcfYc22few+cITdByrw08xCWiUn1CiLky9fdUg9\ncYnr2H+3SUnATLMQEWk4FUIYtEiKJyOpBRltW9Rp+6PHZiE1CqP2x13lFezYd4jVW/eyq7yCw5VV\npzxWQqjAjpVHzctXtYukQ2oSbVsmkZSgm8tE5ONUCAGIj7Pjl4/q6kBF5ccuWx17HFsXKS2vYNXW\nvZSWV7DnwJHTHqt1SsLHLl/VXgfp0CqZvumptE5JDMc/WUSaABVCE9EyKYGWSQlktmtZp+0rj1ax\n5+CRj13GOv4xtC6ytewQK0OzkIpTzEL6dExlaGYaQzPSyM5IY0i3NioJkRilQohRCfFxdGyVTMdW\nyfSvw/buzoGKo8dnHSX7DrN6214+KCpjycZS5r73EQBm0LtjKkMz0o4/hmSk0SpZX0oiTZ2+iwUA\nMyM1OYHU5AS6t6+ehXwuq/Pxz+/cf5jlxWUUFJWxvLiMd09REueEZhEqCZGmKdDvWDO7BrgKaAM8\n5u6Lgswjp9exVTKfHdiJzw7sdHzsWEksD5XEOxtLmVOjJPocm0lktq0uiW5tSFVJiEQt89Pd/1jf\nA5rNBMYBO9w9u8Z4LvAQEA/83t1/WeNz7YDfuPukTzp2Tk6O6x3TolvJvsMUFFcXxAdFZRQUl7Ft\n7yGguiT6prc6vh5xTmYaWV1VEiKNzcyWunvOGbdrhEIYBewH/nSsEMwsHlgLjAaKgCXADe6+MvT5\n/wIed/dln3RsFULTtGPfoeqSKNpbPaMo3sP2vYeBEyVx/HKTSkIk7OpaCGH/rnP3PDPrVWt4GFDo\n7htC4WYDE8xsFfBL4MUzlYE0XZ1ap3DZoBQuG3RiTeJYSRybRby5fifP/qsYqH6dqWMziWN3OGV1\na0PLJJWESGOK1HdYBrClxvMi4ELgNuBzQJqZ9XP3R2vvaGY3ATcB9OjRIwJRJRJOWRJ7D4VmENXr\nEm8UnqIkMk/c3aSSEAmvQL+b3H0aMO0M28wAZkD1JaNI5JJgdGqTwuVtUrh88ImS2L730PFF64Li\nMl5ft5Nnl50oiX6dWlWvRxy/3JRGiyS9HpRIfUSqEIqB7jWeZ4bGRD5R5zYpdM5KOekW2GMl8UGo\nJPLWnlwS/Tu1Dt3+2oahmW3J6tpGJSFSB5EqhCVAfzPrTXURTAS+HKFzS4ypXRLuzva9h2tcbtrD\na2tLeGZZEXCiJGr+xbVKQuTjwl4IZvYEcCnQ0cyKgLvd/TEzuxVYSPVtpzPdfUW4zy3Nk5nRJS2F\nLmkpjK5VEh8U7Tl+G+w/1+zg6aXVJREfZ/Q/drkp80RJ6E2QpDkL+22njUm3nUpDuDvbaqxJHFu8\n3lVeAZwoiWN3N6kkJFYEdtupSLQyM7qmtaBrWgvGDOkCVJfE1rJDxxetPygq45XVO/hbrZnE0Boz\nicEqCYlRKgRp1syMbm1b0K1tC66oVRLH/kZieXEZL9cqiQGdW1cvWodemmNQl9YqCWnyVAgitdQs\nidzsEyXxUdmxy017WF68l3+s2sFT+dUlkRBnfPGC7vz7lYP1on7SZOkrV6QOzIyMttXvilezJIr3\nHKy+9XXdTp5490NeW1PCr68/h8/06xhwYpGzp/dSFKknMyOzXUtys7vyi2uH8vTNF5OUEMdXfv8O\nU+YUUH64MuiIImdFhSASJuf3bM8Lt49k0oje/OWdzeQ+lMfbG3YFHUukzlQIImHUIimeKeOyePKm\ni4kzY+KMt/n5vBUcqNBsQaKfCkGkEQzr3Z4Xvz+Sbw7vxay3NjH2odd5d2Np0LFEPpEKQaSRtExK\n4Ofjh/DEjRdR5c6XZizmvvkrOVhxNOhoIqekQhBpZBf37cCC74/iqxf25LE3NnLltNdZulmzBYk+\nKgSRCEhNTuC+a7J5/DsXUlFZxfWPLuYXL6zi0BHNFiR6qBBEIugz/Tqy4AcjmXhBD2bkbeCqaa/z\nrw93Bx1LBFAhiERc65REHvz8UP707WEcrDjKdY+8xS9fXK3ZggROhSASkFED0lkweRRfzOnOo6+t\n5+qH3+D9LXuCjiXNmApBJEBtUhL55XXn8IdvXcC+Q5V8/pG3+M3CNRyu1GxBIk+FIBIFPjuwEwsn\nj+La8zL4n1cLGf/wmxQUlwUdS5oZFYJIlEhrkchvvnAuM7+Zw+4DFUyY/ia/fWktFZVVQUeTZkKF\nIBJlLhvUmZcmX8KEc7sx7eV1TJj+Jis/2ht0LGkGVAgiUSitZSK//dKnmPG18ynZd5jx//MG015e\nx5Gjmi1I4wm0EMws1cz+aGa/M7OvBJlFJBqNGdKFlyaP4qpzuvLbl9Zy7f++yeptmi1I4wh7IZjZ\nTDPbYWYFtcZzzWyNmRWa2R2h4c8DT7v7jcD4cGcRiQXtUpN4aOJ5PPrVT7N1zyGufvgN/ueVdVRq\ntiBh1hgzhFlAbs0BM4sHpgNjgSzgBjPLAjKBLaHNdJ+dyCfIze7KosmjGDOkC79ZtJbPP/IWa7fv\nCzqWxJCwF4K75wG1X7lrGFDo7hvcvQKYDUwAiqguhUbJIhJrOrRKZvqXP830L3+aot0HGTftDR75\n53rNFiQsIvVDOIMTMwGoLoIM4FngOjN7BHjuVDua2U1mlm9m+SUlJY2fVKQJuOqc6tnCZYM68asF\nq7n+0cUU7tgfdCxp4gL9rdzdy939W+5+i7s/fpptZrh7jrvnpKenRzqiSNTq2CqZR776aabdcB6b\ndpVz5bTXmZG3nqNVHnQ0aaIiVQjFQPcazzNDYyLSAGbG+HO7sWjyKC4ZkM4vXljNFx59iw0lmi3I\n2YtUISwB+ptZbzNLAiYC8yJ0bpGY16l1CjO+dj5Tv3Qu60vKGfvQ6/z+9Q2aLchZaYzbTp8AFgMD\nzazIzCa5eyVwK7AQWAU85e4rwn1ukebMzLj2vEwWTR7FiH4duf/5VUycsZhNO8uDjiZNhLk3nd8g\ncnJyPD8/P+gYIlHP3XlmWTH3PLeCI0eruCN3EF+/uBdxcRZ0NAmAmS1195wzbadbPUVikJlx/fmZ\nvDT5Ei7q04GfP7eSG373Nh/uOhB0NIliKgSRGNYlLYU/fPMCfn3dOaz8aC+5D+Xx58WbqNLagpyC\nCkEkxpkZX7ygOwsmj+L8nu2YMncFX33sHbaUarYgJ1MhiDQTGW1b8KdvD+MX1w7l/S17yP3vPB5/\nZzNNaR1RGpcKQaQZMTO+fGEPFk4exad6tOU//l7A12e+S/Geg0FHkyigQhBphjLbteQvky7k/muy\nWbp5N1dMzePJJR9qttDMqRBEmikz46sX9WThD0aRndGGnz6znG/+YQlbyzRbaK5UCCLNXPf2Lfnr\ndy7invFDeHdjKWOm5vG3/C2aLTRDKgQRIS7O+MbwXiz4wUgGd2nDj5/+gEl/zGf73kNBR5MIUiGI\nyHE9O6Qy+6aLuGtcFm+t38no377Gs8uKNFtoJlQIInKSuDjj2yN688LtI+nfuTU/fOp9bvzTUnbs\n02wh1qkQROSU+qS34ql/u5g7rxpM3roSxkzNY+57xZotxDAVgoicVnyc8Z2RfXjh9pH06pDK92e/\nx81/WUrJvsNBR5NGoEIQkTPq16kVz9wynDvGDuLV1SWMmfoa8z/4KOhYEmYqBBGpk/g44+ZL+vL8\n7SPo0b4lt/71X3zv8WXs2q/ZQqxQIYjIWenfuTXP3DKcH18xkEUrtzFmah4vLt8adCwJAxWCiJy1\nhPg4vvfZfjx32wi6tk3hlseXcdsT/2J3eUXQ0aQBVAgiUm+DurTh79/9DD8aPYAFBVsZPTWPRSu2\nBR1L6kmFICINkhgfx22X92fu90bQqXUyN/15KT+Y/S/2HNBsoakJtBDM7Boz+52ZPWlmY4LMIiIN\nk9WtDXO+9xm+f3l/5n9QPVt4edX2oGPJWah3IZjZTDPbYWYFtcZzzWyNmRWa2R2fdAx3n+PuNwI3\nA1+qbxYRiQ5JCXFMHj2AOd/7DB1Sk5j0x3x+9NT7lB08EnQ0qYOGzBBmAbk1B8wsHpgOjAWygBvM\nLMvMhprZ/FqPTjV2vTO0n4jEgOyMNObdOoLbLuvHnPeKGTP1NV5dvSPoWHIG9S4Ed88DSmsNDwMK\n3X2Du1cAs4EJ7r7c3cfVeuywar8CXnT3ZfX/Z4hItElKiONHYwby9+8Op01KIt+atYSfPP0+ew9p\nthCtwr2GkAFsqfG8KDR2OrcBnwOuN7ObT7WBmd1kZvlmll9SUhK+pCISEedktmX+7SO45dK+PL20\niCum5pG3Vt/L0SjQRWV3n+bu57v7ze7+6Gm2meHuOe6ek56eHumIIhIGyQnx/DR3EM/cMpyWSfF8\nfea7/C5vQ9CxpJZwF0Ix0L3G88zQmIgI5/Vox/O3j+Rzgzvzn4vWsHlXedCRpIZwF8ISoL+Z9Taz\nJGAiMC/M5xCRJiwlMZ4Hrs0mKT6Ou+au0MtpR5GG3Hb6BLAYGGhmRWY2yd0rgVuBhcAq4Cl3XxGe\nqCISKzq3SWHy6AG8traEhfrL5qhhTamdc3JyPD8/P+gYIhIGlUerGPfwG+w9eISXfngJqckJQUeK\nWWa21N1zzrSdXrpCRAKREB/H/ddk81HZIaa9si7oOIIKQUQClNOrPV84P5PHXt/Iuu37go7T7KkQ\nRCRQd4wdRGpyAnfOKdACc8BUCCISqA6tkvlJ7kDe2VjK3Pf0tpxBUiGISOAmXtCDc7u35f7nV+mF\n8AKkQhCRwMXHGfdPyGZX+WGmvrQ26DjNlgpBRKLC0Mw0vnZRT/60eBMFxWVBx2mWVAgiEjV+NGYg\n7VOTuHNOAVVVWmCONBWCiESNtBaJ/PuVg3lvyx6ezN9y5h0krFQIIhJVrj0vg2G92/OrBaspLdf7\nMkeSCkFEooqZcd+EbPYfquRXL64OOk6zokIQkagzsEtrvj2iN0/mb2Hp5t1Bx2k2VAgiEpW+f3l/\nurRJ4c45BVQerQo6TrOgQhCRqJSanMBdV2exaute/vz25qDjNAsqBBGJWmOzuzBqQDr/tWgtO/Ye\nCjpOzFMhiEjUMjPuHT+EiqNVPPDCqqDjxDwVgohEtV4dU7n5kr7Mfe8j3ircGXScmKZCEJGo991L\n+9KjfUumzC2golILzI1FhSAiUS8lMZ57xg9hfUk5v39jQ9BxYlbghWBmqWaWb2bjgs4iItHrs4M6\nccWQzjz8ciFFuw8EHScm1bsQzGymme0ws4Ja47lmtsbMCs3sjjoc6qfAU/XNISLNx11XDwHg3udW\nBpwkNjVkhjALyK05YGbxwHRgLJAF3GBmWWY21Mzm13p0MrPRwEpgRwNyiEgzkdG2Bbdf3p9FK7fz\nyurtQceJOQn13dHd88ysV63hYUChu28AMLPZwAR3fxD42CUhM7sUSKW6PA6a2QvurhUjETmtSSN6\n88yyIu6et4LhfTuSkhgfdKSYEe41hAyg5mvWFoXGTsnd/8PdfwD8FfjdqcrAzG4KrTHkl5SUhDmu\niDQ1SQlx3DthCFtKD/K/rxYGHSemBL6oDODus9x9/mk+N8Pdc9w9Jz09PdLRRCQKDe/bkQmf6saj\nr21g487yoOPEjHAXQjHQvcbzzNCYiEhY/ceVg0lOiOOuuQW4693VwiHchbAE6G9mvc0sCZgIzAvz\nOURE6NQmhR+OGcDr63byYsG2oOPEhIbcdvoEsBgYaGZFZjbJ3SuBW4GFwCrgKXdfEZ6oIiIn+9pF\nPcnq2oZ7n1vJ/sOVQcdp8updCO5+g7t3dfdEd89098dC4y+4+wB37+vuD4QvqojIyRLi47jvmmy2\n7T3EtJfXBR2nyYuKRWURkfo6v2c7Jl7Qncfe2MiabfuCjtOkqRBEpMn7Se4gWqckMGWOFpgbQoUg\nIk1e+9Qk7sgdxLubSnl2mW5srC8VgojEhC/mdOe8Hm158MVVlB04EnScJkmFICIxIS7OuG9CNqXl\nFfxm0Zqg4zRJKgQRiRnZGWl8/eJe/OWdzSwvKgs6TpOjQhCRmPLDMQPokJrMnXOWc7RKC8xnQ4Ug\nIjGlTUoid141mPeLypi95MOg4zQpKgQRiTkTPtWNi/q059cL1rBr/+Gg4zQZKgQRiTlm1QvM5Ycr\n+eWLq4OO02SoEEQkJvXv3JrvjOzD35YWkb+pNOg4TYIKQURi1u2X96NbWgp3zimg8qjejPFMVAgi\nErNaJiVw19VDWL1tH7Pe2hR0nKinQhCRmHbFkM5cOjCdqS+tZVvZoaDjRDUVgojENDPjnvFDOFLl\n3P/8yqDjRDUVgojEvJ4dUvnepf2Y/8FW3li3M+g4UUuFICLNwr9d0oeeHVpy19wCDlceDTpOVFIh\niEizkJIYzz3jh7BhZzm/f31j0HGiUqCFYGZxZvaAmT1sZt8IMouIxL5LB3ZibHYXpr28ji2lB4KO\nE3XqXQhmNtPMdphZQa3xXDNbY2aFZnbHGQ4zAcgEjgBF9c0iIlJXU8ZlER9n3PPciqCjRJ2GzBBm\nAbk1B8wsHpgOjAWygBvMLMvMhprZ/FqPTsBA4C13/yFwSwOyiIjUSbe2Lfj+5f35x6odvLRye9Bx\nokq9C8Hd84Dafw8+DCh09w3uXgHMBia4+3J3H1frsYPqWcHu0L5a5RGRiPj2iN7079SKn89bwcEK\n/eg5JtxrCBnAlhrPi0Jjp/MscIWZPQzknWoDM7vJzPLNLL+kpCR8SUWk2UqMj+O+a7Ip3nOQ6a8W\nBh0nagS6qOzuB9x9krvf5u7TT7PNDHfPcfec9PT0SEcUkRh1UZ8OfP68DP5/3nrWl+wPOk5UCHch\nFAPdazzPDI2JiESdn105mJTEeO6euwJ3vbtauAthCdDfzHqbWRIwEZgX5nOIiIRFeutkfnzFQN4o\n3Mn8D7YGHSdwDbnt9AlgMTDQzIrMbJK7VwK3AguBVcBT7q57u0Qkan3lwp5kZ7Thvvkr2XfoSNBx\nAtWQu4xucPeu7p7o7pnu/lho/AV3H+Dufd39gfBFFREJv/i46ndXK9l/mP/+x7qg4wRKL10hIs3e\neT3aMfGCHsx6axOrtu4NOk5gVAgiIsBPrhhIWotEpswpoKqqeS4wqxBERIB2qUnckTuI/M27eWZZ\n83wlHRWCiEjI9edncn7Pdjz44mr2HKgIOk7EqRBERELiQgvMew5U8J8L1wQdJ+JUCCIiNWR1a8M3\nh/fmr+9+yPtb9gQdJ6JUCCIitUwe3Z/0VsncOaeAo81ogVmFICJSS+uURO4cl8Xy4jL++s7moONE\njApBROQUrj6nK8P7duDXC9dQsu9w0HEiQoUgInIKZsa9E7I5dOQoD764Kug4EaFCEBE5jX6dWnHj\nyD48u6yYdzbsCjpOo1MhiIh8gtsu609G2xZMmVvAkaNVQcdpVCoEEZFP0CIpnruvzmLt9v3MenNT\n0HEalQpBROQMRmd15vJBnZj6j7VsLTsYdJxGo0IQETkDM+Pn44dwtMq5f37sLjCrEERE6qB7+5bc\n+tl+PL98K3lrS4KO0yhUCCIidXTTJX3o3TGVu+YWcOjI0aDjhJ0KQUSkjpIT4rl3whA27TrAjLwN\nQccJOxWCiMhZGNk/navO6cr0Vwv5cNeBoOOEVaCFYGY9zGyOmc00szuCzCIiUldTrsoiIc64e14B\n7rHz4nf1LoTQD/EdZlZQazzXzNaYWWEdfsgPBZ52928D59U3i4hIJHVJS2Hy6AG8uqaERSu3Bx0n\nbBoyQ5gF5NYcMLN4YDowFsgCbjCzLDMbambzaz06AW8Dk8zsFWBBA7KIiETUN4b3YmDn1tz73EoO\nVFQGHScs6l0I7p4HlNYaHgYUuvsGd68AZgMT3H25u4+r9dgBfAu4290vA66qbxYRkUhLjI/j/muz\nKd5zkIdfKQw6TliEew0hA9hS43lRaOx0FgC3m9mjwKZTbWBmN5lZvpnll5TE5r2/ItI0XdCrPdd9\nOpPfv76Bwh37go7TYIEuKrt7gbtf7+43u/v/O802M9w9x91z0tPTIx1RROQT/ezKQbRIjGfKnBVN\nfoE53IVQDHSv8TwzNCYiEpM6tkrmx7mDWLxhF/Pe/yjoOA0S7kJYAvQ3s95mlgRMBOaF+RwiIlHl\ny8N6cE5mGvc/v4q9h44EHafeGnLb6RPAYmCgmRWZ2SR3rwRuBRYCq4Cn3H1FeKKKiESn+Djj/muy\n2bn/MFNfWht0nHpLqO+O7n7DacZfAF6odyIRkSbonMy2fOXCHvzxrU1cf34mQ7qlBR3prOmlK0RE\nwuTHYwbRrmUSU+YUUFXV9BaYVQgiImGS1jKRn105mGUf7uFvS7eceYcoo0IQEQmj6z6dwQW92vHL\nF1ezu7wi6DhnRYUgIhJGZsZ912Sz91Alv164Oug4Z0WFICISZoO6tOFbw3sxe8kWln24O+g4daZC\nEBFpBD8YPYBOrZOZMqeAo01kgVmFICLSCFolJzBlXBYrPtrLX97eHHScOlEhiIg0kquGdmVk/478\nZuEaduw7FHScM1IhiIg0EjPjnvFDOFxZxYMvRP8CswpBRKQR9Ulvxb9d0oe//6uYxet3BR3nE6kQ\nREQa2Xcv7UdmuxbcNbeAisqqoOOclgpBRKSRtUiK557xQ1i3Yz8z39wYdJzTUiGIiETA5YM7Mzqr\nMw/9Yx0f7TkYdJxTUiGIiETI3Vdn4Tj3Prcy6CinpEIQEYmQzHYtue2y/ixYsY1X1+wIOs7HqBBE\nRCLoxpF96JOeyt1zV3DoyNGg45xEhSAiEkFJCXHcNyGbD0sP8Mg/1wcd5yQqBBGRCPtMv45cfW43\nHnltPZt2lgcd5zgVgohIAO68ajBJ8XHcPW8F7tHx4ncRKwQz62Nmj5nZ0zXGUs3sj2b2OzP7SqSy\niIgErXObFCaPHsBra0tYuGJb0HGAOhaCmc00sx1mVlBrPNfM1phZoZnd8UnHcPcN7j6p1vDngafd\n/UZg/FklFxFp4r5xcU8GdWnNPc+tpPxwZdBx6jxDmAXk1hwws3hgOjAWyAJuMLMsMxtqZvNrPTqd\n5riZwLE3Ho2u5XYRkUaWEB/HA9dms7XsENNeWRd0nLoVgrvnAaW1hocBhaHf/CuA2cAEd1/u7uNq\nPU53w20R1aVQ5ywiIrHk/J7t+WJOJo+9vpG12/cFmqUhP4QzOPHbPVT/cM843cZm1sHMHgXOM7Of\nhYafBa4zs0eA506z301mlm9m+SUlJQ2IKyISnX6aO4jU5ASmzCkIdIE5Yr+Vu/sud7/Z3fu6+4Oh\nsXJ3/5a73+Luj59mvxnunuPuOenp6ZGKKyISMR1aJfPT3EG8s7GUOe8VB5ajIYVQDHSv8TwzNCYi\nImdp4gXdObd7Wx54fjVlB48EkqEhhbAE6G9mvc0sCZgIzAtPLBGR5iUuzrh/Qjal5Yf57aI1wWSo\ny0Zm9gSwGBhoZkVmNsndK4FbgYXAKuApd1/ReFFFRGLb0Mw0vnpRT/789mYKissifn6Llr+Qq4uc\nnBzPz88POoaISKMpO3iEy//rn2S2a8mztwwnLs4afEwzW+ruOWfaTrd6iohEkbQWifz7lYN5b8se\nnszfcuYdwkiFICISZa49L4NhvdvzqwWrKS2viNh5VQgiIlHGzLj/mmz2H6rkVy+ujth5VQgiIlFo\nQOfWTBrRmyfzt7B0c+0XimgcKgQRkSh1++X96ZqWwp1zVlB5tKrRz6dCEBGJUqnJCdw1LotVW/fy\np8WbG/18KgQRkSiWm92FSwakU1iyv9HPldDoZxARkXozM2Z8/XySE+Ib/VyaIYiIRLlIlAGoEERE\nJESFICIigApBRERCVAgiIgKoEEREJESFICIigApBRERCmtQb5JhZCdCQv99OAxrzbYjCffxwHK8h\nx6jPvmezT0dg51kevzlr7K/fxhBk5kicu6l8z/d09/Qz7unuzeYBzGhKxw/H8RpyjPrsezb7APlB\nf000pUdjf/3GWuZInDvWvueb2yWj55rY8cNxvIYcoz77Nvb/4+asKf6/DTJzJM4dU9/zTeqSkcQW\nM8v3OrzPq4hERnObIUh0mRF0ABE5QTMEEREBNEMQEZEQFYKIiAAqBBERCVEhSNQwsz5m9piZPR10\nFpHmSIUgjcrMZprZDjMrqDWea2ZrzKzQzO4AcPcN7j4pmKQiokKQxjYLyK05YGbxwHRgLJAF3GBm\nWZGPJiI1qRCkUbl7HlBaa3gYUBiaEVQAs4EJEQ8nIidRIUgQMoAtNZ4XARlm1sHMHgXOM7OfBRNN\npPlKCDqAyDHuvgu4OegcIs2VZggShGKge43nmaExEQmQCkGCsATob2a9zSwJmAjMCziTSLOnQpBG\nZWZPAIvxei7DAAAAVElEQVSBgWZWZGaT3L0SuBVYCKwCnnL3FUHmFBG9uJ2IiIRohiAiIoAKQURE\nQlQIIiICqBBERCREhSAiIoAKQUREQlQIIiICqBBERCREhSAiIgD8HxZ1TVY9r4zmAAAAAElFTkSu\nQmCC\n",
            "text/plain": [
              "<matplotlib.figure.Figure at 0x7f11e901f668>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "def cg(A,b,niter,x0):\n",
        "    rk = b - A @ x0\n",
        "    sk = rk\n",
        "    xk = x0\n",
        "    for i in range(niter):\n",
        "        alpha = np.inner(rk,rk)/np.inner(sk, A @ sk)\n",
        "        xk1 = xk + alpha * sk\n",
        "        rk1 = rk - alpha * A @ sk\n",
        "        beta = np.inner(rk1,rk1)/np.inner(rk,rk)\n",
        "        sk1 = rk1 + beta*sk\n",
        "        rk = rk1\n",
        "        xk = xk1\n",
        "        sk = sk1\n",
        "    return xk\n",
        "\n",
        "iters = np.asarray(2**np.arange(2,7),dtype=int)\n",
        "\n",
        "x0 = np.random.random(n*n)\n",
        "\n",
        "cg_results = []\n",
        "err = []\n",
        "\n",
        "for niter in iters:\n",
        "    cg_results.append(cg(A,b,niter,x0.copy()))\n",
        "    err.append(la.norm(cg_results[-1]-x))\n",
        "\n",
        "pt.plot(iters,err)\n",
        "pt.yscale('log')\n",
        "pt.xscale('log')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": true
      },
      "outputs": [],
      "source": []
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.5.2"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 2
}