tutorials/3_convergence.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial 4: Convergence rates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The problem\n",
    "\n",
    "It is reasonably easy to check that $p=\\mathrm{e}^{\\mathrm{i}kx_0}$ is the solution of\n",
    "\n",
    "$$\n",
    "\\Delta p + k^2p=0 \\quad \\text{in }\\Omega,\n",
    "$$\n",
    "$$\n",
    "\\frac{\\partial p}{\\partial \\mathbf{n}} = \\mathrm{i}kn_0\\mathrm{e}^{\\mathrm{i}kx_0}\\quad \\text{on }\\Gamma,\n",
    "$$\n",
    "where $\\mathbf{x}=(x_0,x_1,x_2)$ is a point in $\\Omega$ and $\\textbf{n}=(n_0,n_1,n_2)$ is the normal to the surface $\\Gamma$.\n",
    "\n",
    "In this tutorial, we solve this problem using BEM and compare the approximate solution to this knows actual solution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## BEM formulation\n",
    "\n",
    "### Representation formula\n",
    "\n",
    "$$\n",
    "p = \\mathcal{D}p-\\mathcal{S}\\frac{\\partial p}{\\partial \\mathbf{n}}\n",
    "$$\n",
    "\n",
    "where $\\mathcal{S}$ is the single layer potential operator and $\\mathcal{D}$ is the double layer potential operator.\n",
    "\n",
    "### Boundary integral equation\n",
    "\n",
    "$$\n",
    "(\\mathsf{D}+\\tfrac{1}{2}\\mathsf{I})p=\\mathsf{S}\\frac{\\partial p}{\\partial \\mathbf{n}},\n",
    "$$\n",
    "\n",
    "where $\\mathsf{S}$ is the single layer boundary operator; $\\mathsf{D}$ is the double layer boundary operator; and $\\mathsf{I}$ is the identity operator."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solving with Bempp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can solve this problem using the code below. In this example, we use a cuboid 1.5 units long, 1 unit wide, and 1 unit high: we import this mesh from the file `cuboid.msh`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import bempp.api\n",
    "from bempp.api.operators.boundary import helmholtz, sparse\n",
    "from bempp.api.operators.potential import helmholtz as helmholtz_potential\n",
    "from bempp.api.linalg import gmres\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "k = 4.\n",
    "\n",
    "grid = bempp.api.import_grid(\"cuboid.msh\")\n",
    "\n",
    "space = bempp.api.function_space(grid, \"DP\", 0)\n",
    "\n",
    "identity = sparse.identity(space, space, space)\n",
    "single_layer = helmholtz.single_layer(space, space, space, k)\n",
    "double_layer = helmholtz.double_layer(space, space, space, k)\n",
    "\n",
    "@bempp.api.complex_callable\n",
    "def lambda_callable(x, n, domain_index, result):\n",
    "    result[0] = 1j * k * np.exp(1j * k * x[0]) * n[0]\n",
    "\n",
    "lambda_fun = bempp.api.GridFunction(space, fun=lambda_callable)\n",
    "\n",
    "p_total, info = gmres(double_layer + 0.5 * identity, single_layer * lambda_fun, tol=1E-5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now compare our solution to the solution we are expecting. To do this, we create a `GridFunction` representing the actual solution, and take the $L^2$ norm of the difference between our approximation and the actual solution.\n",
    "\n",
    "0.4 isn't too big, so it looks like we've got a pretty good approximation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.41326500769117885\n"
     ]
    }
   ],
   "source": [
    "@bempp.api.complex_callable\n",
    "def actual_solution_f(x, n, domain_index, result):\n",
    "    result[0] = np.exp(1j * k * x[0])\n",
    "\n",
    "actual_solution = bempp.api.GridFunction(space, fun=actual_solution_f)\n",
    "\n",
    "print((p_total - actual_solution).l2_norm())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get a better understanding of the accuracy of our BEM approximation, we can look at how the $L^2$ norm changes as we change the number of elements in our mesh.\n",
    "\n",
    "The following code uses a for loop to calculate the error for a range of values of $h$. The $h$ parameter controls the size of each triangle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.5 0.7974278842777981\n",
      "0.3535533905932738 0.4639766101213504\n",
      "0.25 0.3615649671805228\n",
      "0.1767766952966369 0.19117350563272814\n",
      "0.125 0.11249838362456548\n",
      "0.08838834764831845 0.05736033728592835\n",
      "0.0625 0.03268628009237414\n",
      "0.04419417382415922 0.01696026668525607\n"
     ]
    }
   ],
   "source": [
    "import bempp.api\n",
    "from bempp.api.operators.boundary import helmholtz, sparse\n",
    "from bempp.api.linalg import gmres\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "k = 4.\n",
    "\n",
    "@bempp.api.complex_callable\n",
    "def lambda_callable(x, n, domain_index, result):\n",
    "    result[0] = 1j * k * np.exp(1j * k * x[0]) * n[0]\n",
    "\n",
    "@bempp.api.complex_callable\n",
    "def actual_solution_f(x, n, domain_index, result):\n",
    "    result[0] = np.exp(1j * k * x[0])\n",
    "\n",
    "h_values = []\n",
    "errors = []\n",
    "ndofs = []\n",
    "for i in range(2, 10):\n",
    "    h = 2 ** -(i/2)\n",
    "\n",
    "    grid = bempp.api.shapes.cuboid(length=(1.5, 1, 1), h=h)\n",
    "    space = bempp.api.function_space(grid, \"DP\", 0)\n",
    "\n",
    "    identity = sparse.identity(space, space, space)\n",
    "    single_layer = helmholtz.single_layer(space, space, space, k)\n",
    "    double_layer = helmholtz.double_layer(space, space, space, k)\n",
    "\n",
    "    lambda_fun = bempp.api.GridFunction(space, fun=lambda_callable)\n",
    "\n",
    "    p_total, info = gmres(double_layer + 0.5 * identity, single_layer * lambda_fun, tol=1E-5)\n",
    "\n",
    "    actual_solution = bempp.api.GridFunction(space, fun=actual_solution_f)\n",
    "    \n",
    "    h_values.append(h)\n",
    "    ndofs.append(space.global_dof_count)\n",
    "    errors.append((p_total - actual_solution).l2_norm())\n",
    "    print(h, errors[-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then make a convergence plot showing how the error decreases as we decrease $h$.\n",
    "\n",
    "We have used some features of matplotlib to make this plot as informative as possible:\n",
    "\n",
    "- We use `plt.xscale` and `plt.yscale` to plot $\\log(h)$ against $\\log(\\text{error})$. The error will be approximately polynomial: plotting it on log-log axes will turn it into a straight line where the gradient is the polynomial order.\n",
    "- We use `plt.axis` to make the axes equal aspect. The gradient is an important feature of the plot, and it is much easier to just the gradient on a equal aspect plot.\n",
    "- We use `plt.xlim` to reverse the $h$-axis. This one is due to my own personal preference (and you may prefer otherwise), but I think it makes more sense for the number of elements to increase to the right (and so for $h$ to decrease to the right)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEcCAYAAADQqlM0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAfV0lEQVR4nO3deZgU1b3/8fd3GEBGEVBARWAmrpEr4IIajUZjFNEEt5hcSftLNAjiHqOoCKKCCBpicEEJLtHEEbfcuEXvlegl0bjigmjU6zYDisoi4DKgAuf3x6mJzdgzdM109+mu/ryep5/pOtVd9e2umfpMbafMOYeIiEi2KkIXICIipUXBISIisSg4REQkFgWHiIjEouAQEZFYFBwiIhKLgkOklcz7g5ktN7Nn8zyvlJk9ks95NJnffmb2RqHmF5eZ1ZnZQc2MO8DM3it0TeVEwZFQZvYzM5trZp+Z2Qdm9rCZ7Ru6roTZFzgY6O2c2zOfM3LO1TrnBmfzWjM73syeaOP8HnfO7diWaUhyKTgSyMx+DUwDLgO2APoC1wFHhKwrnZlVhq4hB6qBOufc56ELyaW2Lhszaxdy/lIAzjk9EvQAugCfAT9p4TUd8cGyKHpMAzpG4w4A3gPOBhYDHwAnROO+A3wItEub1lHAy9HzCuB84G1gGXAXsFk0rgZwwHBgAfCPqP3nQH30+guBOuCgGNP7RTS9pcDYtLraARdE7/0UeB7oE437NjAb+Bh4A/hpC99VL+D+6LVvASOi9uHAamBt9H1fkuG92wKPRbUvBWqBrmnj64BzgJeBlcCdwEbN1HE88ETasANGAW8Cy4HpgAE7NalrRdoynxp9Vx8BM4BOTZb5edHy/VNjW9r8dgLmACuAV4HD08bdAlwPPAR83rj8mtQ/B5gMPBt91vuy+N04PJrXiuj9OzX57sYA/4o+/x8av7sMtfcC/gwsAd4FzkgbdzFwN3Bb9HsyH9ghmvZiYCEwOPTfdbE9ghegR44XKAwB1gCVLbxmAvA00BPoATwJTIzGHRC9fwLQHjgMaAC6RePfBg5Om9bdwPnR819F0+0drah+D8yKxjWuHP4IbAx0AvpFK7d9gQ7Riu0rvg6ObKZ3QzStgcAXjSsXYHS0EtgRv0IdCGwezXshcAJQCeyGX6n/RzPf1d/xW2sbAbtEK58fROOOJ21lnuG92+F3ZXWMvud/ANPSxtfhV6S9gM2A14BRzUxrvXlFn/1BoCt+i3IJMKS5uvD/HNwfzacz8AAwuckyvzyqtRNpK9/o9+AtfBB3AA7Er2R3jMbfgg+D7+LD/hvhh1/xvw/sHC2DPwO3tfC7sQM+hA6O5n9uVEOHtO/uFaBP9Jn+CVya9nkaa6/A/9MwPqp9G+Ad4JBo/MX4oD0k+n34Iz5cxkbzHQG8G/rvutgewQvQI8cLFFLAhxt4zdvAYWnDh+B3uTT+0a0iLXjw/3l9J3p+KXBz9Lxz9MddHQ2/RrRSjYa3wgdBZdrKYZu08eOJgiAargK+5OvgyGZ6vdPGPwscGz1/Azgiw2f/T+DxJm2/By7K8No++P/cO6e1TQZuiZ4fTwvBkWF6RwIvpg3XAcelDV8BzGjmvevNK/rs+6YN38XXAd70tRYtp23T2vZuXCFGy/xL0lb4rL/y3Q+/JVKRNn4WcHH0/Bbgjxv47HOAKWnD/aJ5tmvmd+NC4K604Qp88ByQ9t2NSht/GPB2htr3AhY0qWUM8Ifo+cXA7LRxQ/H/zLRL+x13pG0p6uHQvsTkWQZ0N7NK59yaZl7TC797qFF91PbvaTR5bwOwSfT8duBJMzsZOBp4wTnXOK1q4C9mti7tvWvxx1kaLWxSx7+HnXMNZrYsbXw20/uwmTr74AOyqWpgLzNbkdZWid8901Qv4GPn3KdpbfXAoAyv/QYz6wlcjV/xdsav/JY3eVnT+nuRveY+e1M98KH8vJn9uzz8SrvREufc6mbe3wtY6JxLXw71wNZpwwvZsPTX1OP/o+/ezPj1fkedc+vMbGEL82z6O9yoGujVZHm3Ax5PG/4o7fkqYKlzbm3aMPjvNn0aZU0Hx5PnKfym95EtvGYR/g+qUd+obYOcc//C/5EeCvwMHySNFgKHOue6pj02cs69nz6JtOcf4HdDAWBmnfC7k+JMrzkL8ccYMrX/vck0N3HOnZzhtYuAzcysc1pbX/x/vtmYjP+8A5xzmwLH4VfY+eaaDC/FrwD/I+0zd3HObdLCe9ItAvqYWfr6oun30NL7G/Vp8v6votoyTWO931HzidenyTybTi/T7/BC/JZV+vLu7Jw7LIt6pRkKjoRxzq3E7wKabmZHmlmVmbU3s0PN7IroZbOAcWbWw8y6R6+/LcZsbgfOAL6HP8bRaAYwycyqAaLpt3Qm1z3AUDPbx8w6AJew/oo17vTS3QhMNLPto+stBpjZ5vjjAjuY2f+Lvpf2ZraHme3UdALOuYX44z+TzWwjMxuAP4Bbm2UNnYkOUJvZ1vjjLoXwEdA7+k6JthRuAH4XbQVhZlub2SFZTu8Z/K6uc6Pv6wD8Lp07YtZ1nJn1M7Mq/DG0e9L+s2/qLuCHZvYDM2uPP1njC/zyaHSqmfU2s83wx1/uzDCdZ4FPzOw8M+tkZu3MbGcz2yNm7ZJGwZFAzrkrgV8D4/AHTRcCpwH3Ri+5FJiLP5tnPvBC1JatWfj9yI8559L/Y7wKfwD2ETP7FH9ge68W6nwVOB2/AvoAf8B1MX4FEXt6TVyJX/k8AnwC3IQ/i+hTYDBwLP4/1A/5+qBwJsPw++AXAX/BHwuZnWUNl+APvq8E/gr8V5bva6vH8GcjfWhmjcvnPPzB5afN7BPgb/gTBzbIOfcl/gynQ/FbCNcBP3fOvR6zrj/hj4d8iD/Z4IwW5vkGfgvtmmieQ4GhUS2Nbscv33eixzd+h6NgGoo/seHdaFo34s8+lFay6ACQSHBm1rgfeXvn3Luh65HcMbM5+LOobgxdi7SdtjgkKDMbGu1O2xh/Ou58/BkzIlKkFBwS2hF8fSHi9vjTabUZLFLEtKtKRERi0RaHiIjEouAQEZFYEn/lePfu3V1NTU3oMkRESsrzzz+/1DnXI9O4xAdHTU0Nc+fODV2GiEhJMbP65sZpV5WIiMSS2OCIrg+YuXLlytCliIgkSmKDwzn3gHNuZJcu6llARCSXEhscIiKSHwoOCaq2FmpqoKLC/6zNtt9ZEQkm8WdVSfGqrYWRI6GhwQ/X1/thgFQqXF0i0jJtcUgwY8d+HRqNGhp8u4gUr8QGh86qKn4LFsRrF5HikNjg0FlVxa9v33jtIlIcEhscUvwmTYKqqvXbqqp8u4gULwWHBJNKwcyZUF0NFt1pfOhQHRgXKXYKDgkqlYK6Oli3Do4+Gu67D958M3RVItISBYcUjWuugY4d/Sm5ur+YSPFScEjR6NULrrgC5syBm28OXY2INCexwaHTcUvTiSfC974H55wDH34YuhoRySSxwaHTcUtTRYU/YL5qFZxxRuhqRCSTxAaHlK4dd4Tx4+Huu/3BchEpLgoOKUqjR0P//nDKKaC9jSLFRcEhRal9e7jxRn+cY8yY0NWISDoFhxStPff0xzmuvx6eeCJ0NSLSSMEhRW3iRH9l+YgR8MUXoasREVBwSJHbZBOYMQNefx0uuyx0NSICCQ4OXceRHEOGwHHHweTJ8MoroasRkcQGh67jSJbf/Q66dPEXCK5dG7oakfKW2OCQZOneHaZNg2eegeuuC12NSHlTcEjJ+NnP/G6rCy7QXQJFQlJwSMkw86fmrlsHJ5+sHnRFQlFwSEmpqfF3CHzoIbjzztDViJQnBYeUnNNPhz328BcHLlsWuhqR8qPgkJLTrp3vjmT5cjj77NDViJQfBYeUpAED4Lzz4NZbYfbs0NWIlJfEBocuAEy+ceNghx3gpJOgoSF0NSLlI7HBoQsAk2+jjeCGG+Ddd2GrrfxNoGpqoLY2dGUiyVYZugCRtli4ECor4ZNP/HB9PYwc6Z+nUuHqEkmyxG5xSHkYOxbWrFm/raHBt4tIfig4pKQ1dwW5riwXyR8Fh5S0vn0zt3fp4q8wF5HcU3BISZs0Caqq1m9r1w5WrIDBg2HRojB1iSSZgkNKWioFM2f6uwSa+Z+33urbnnzSX+/xwAOhqxRJFnMJ7ylu0KBBbu7cuaHLkABefx2GDYOXXoJTT4Xf/AY6dQpdlUhpMLPnnXODMo3TFock1re/DU8/Db/+NUyf7vu3mj8/dFUipU/BIYnWsSP89rfw8MOwZIkPj2uvVZfsIm2h4JCyMGQIvPwyHHig71338MN9kIhIfIkNDvVVJU1tsQX89a9w1VXwyCP+wLk6SBSJL7HBob6qJBMzfx+PZ5+Fbt38KbujR8OXX4auTKR0JDY4RFoycCDMnQujRsHUqbDPPv5YSE2NOksU2RCdjitl79574bjj4PPP12+vqvLXg6izRClHOh1XpAVHHgldu36zXZ0limSm4BCh+a5J1FmiyDcpOERovrPE5tpFypmCQ4TMnSV26ODbRWR9Cg4RvtlZYocOPkiOOip0ZSLFR8EhEkmloK7O38fj0Ud91+zTpoWuSqT4KDhEMth3X98tyZQp6ppEpCkFh0gzJk/213ZcemnoSkSKi4JDpBn9+sHw4XD99fDOO6GrESkeCg6RFlx8MVRW6kJAkXQKDpEW9OrlbwR1xx3w3HOhqxEpDokNDnWrLrly7rnQvTucd55uACUCCQ4OdasuubLppjB+PPzv/8J//3foakTCS2xwiOTSSSfBttv6rY61a0NXIxKWgkMkC43dj8yfD7fdFroakbAUHCJZ+slPYNAgGDcOVq0KXY1IOAoOkSxVVMAVV8B778E114SuRiQcBYdIDN//Phx2mL+q/OOPQ1cjEoaCQySmKVNg5Uq47LLQlYiEoeAQial/f/jFL/zuqrq60NWIFJ6CQ6QVJkzwxzwuvDB0JSKFp+AQaYU+feDMM6G2Fl56KXQ1IoWl4BBppfPPh27d/EWBIuVEwSHSSl27+l5zH3kEZs8OXY1I4Sg4RNrg1FOhpsZvdaxbF7oakcJQcIi0QceO/g6BL74Is2aFrkakMBQcIm00bBjsuqvviuSLL0JXI5J/Cg6RNqqogMsv99d0XHdd6GpE8k/BIZIDBx8Mgwf73VYrVoSuRiS/FBwiOXL55bB8ue+SRCTJFBwiObLLLpBKwVVXwcKFoasRyZ/EBofuOS4hTJwIX30FO+3kj33U1Piry0WSJLHBoXuOSwj//CeYweefg3NQXw8jRyo8JFkSGxwiIYwdC2vWrN/W0ODbRZJCwSGSQwsWxGsXKUUKDpEc6ts3XrtIKVJwiOTQpElQVbV+W7t2vl0kKRQcIjmUSsHMmVBd7Q+Sd+0Ka9dChw6hKxPJHQWHSI6lUr77kXXrYMkS2H1334vu0qWhKxPJDQWHSB5VVsLNN/sryn/1q9DViOSGgkMkzwYM8Kfj1tbCgw+Grkak7bIODvP65LMYkaS64ALYeWcYNQrUmYGUuqyDwznngHvzWItIYnXo4HdZffABjB4duhqRtom7q+ppM9sjL5WIJNwee8DZZ8MNN8Bjj4WuRqT1zG9IZPlis38BOwD1wOeA4TdGBuSnvLYbNGiQmzt3bugyRABYtQoGDvTdksyfDxtvHLoikczM7Hnn3KBM4ypjTuvQHNQjUrY6dYIbb4T99/cHzKdNC12RSHyxdlU55+qBrsDQ6NE1ahORLH3ve/66jquvhiefDF2NSHyxgsPMzgRqgZ7R4zYzOz0fhYkk2eTJ0KcPDB8Oq1eHrkYknrgHx4cDeznnxjvnxgPfAUbkviyRZOvc2R8kf/11mDAhdDUi8cQNDgPWpg2vjdpEJKbBg+H44+GKK+CFF0JXI5K9uMHxB+AZM7vYzC4GngZuynlVImXiyiuhRw/45S/9LWdFSkGsK8eBu4ETgI+B5cAJzjmdFyLSSt26wfXXw7x5cPnloasRyU7Wp+M655yZ3euc2x3QhrVIjhx5JPz0pzBxIhx9NPTrF7oikZbpynGRInDNNf6A+S9/6e/fIVLM4gbH94GnzOxtM3vZzOab2cv5KEyknPTs6a/reOYZuOqq0NWItCzrXVXRMY5R+O5GRCTHhg2DWbNg3Dg4/HDYbrvQFYlkFrd33N855+qbPvJYn0jZMIMZM6B9exgxwt9BUKQY6RiHSBHZemv47W9hzhx/73KRYtSaYxxP6xiHSP4MHw4/+IG/1Wzv3lBRATU1/g6CIsVAveOKFBkz+NGP4NFH4f33fVt9PYwc6Z+nUuFqE4H4WxwLgP2AX0THNhywRc6rEilzmbpbb2jwXbGLhBY3OK4D9gaGRcOfAtNzWpGIsGBBvHaRQoobHHs5504FVgM455YDHXJelUiZ69s3XrtIIcUNjq/MrB1+FxVm1gPQSYMiOTZpElRVrd/WqZNvFwktbnBcDfwF6Glmk4AngMtyXlULzGwbM7vJzO4p5HxFCimV8qfjVlf7g+Xg+7HSgXEpBnFvHVsLnAtMBj4AjnTO3Z3t+83sZjNbbGavNGkfYmZvmNlbZnb+Bmp4xzk3PE7dIqUolYK6On8h4IEHwuzZ8NlnoasSib/FgXPudefcdOfctc6512K+/RZgSHpDtOtrOv5U337AMDPrZ2b9zezBJo+ecesVSYKJE2HxYrj22tCViLQiONrCOfcP/L080u0JvBVtSXwJ3AEc4Zyb75z7UZPH4mzmY2YjzWyumc1dsmRJjj+FSOHtsw8cdpi/W+DKlaGrkXJX0OBoxtbAwrTh96K2jMxsczObAexqZmMyvcY5N9M5N8g5N6hHjx65rVYkkAkTYPnyzNd4iBRSMQRHpnuWu+Ze7Jxb5pwb5Zzb1jk3OY91iRSV3XeHo47yt5v9uOl2u0gBbTA4zOxgM7vBzHaJhkfmuIb3gD5pw72BRTmeh0giXHIJfPopTJ0auhIpZ9lscZwCjAaOM7MDgV1yXMNzwPZm9i0z6wAcC9yf43mIJEL//nDssf5mT4uzOuInknvZBMcS59wK59w5wGCg1d2qm9ks4ClgRzN7z8yGO+fWAKcB/wO8BtzlnHu1tfNIm9dQM5u5UkcSJWEuughWr4YpU0JXIuXK/P2ZWniB2RHOufvShk93zl2T98pyZNCgQW7u3LmhyxDJqRNO8HcLfPttfw8PkVwzs+edc4MyjdvgFkd6aETDJRMaIkk1fjysXQuXFbTfBhGvzWdVmdk/c1GIiGTvW9+CE0+EG27wV5eLFFIuTsftlYNpiEhMY8f6uwNOnBi6Eik3WQWHmV0TXY29t5l1bjK65YMkIpIXvXvDqFFw663w5puhq5Fyku0Wx3xgADAFqDOzd83s/qiH3KZBUhR0VpWUg/PPh44d/fUdIoWSVXBEXXic5pzb3zm3Of72sTOAT/Cn0RYd59wDzrmRXbp0CV2KSN5suSWcfjrcfju82uaT2EWys8HTcUudTseVpFu2zB8sHzwY7tFdaiRH2nQ6rogUt803h7POgj//GV58MXQ1Ug4UHCIJcNZZ0K2bv75DJN8UHCIJ0LUrnHMOPPggPP106Gok6RQcIglxxhnQvTtceGHoSiTpEhscOh1Xys0mm8CYMfC3v8Hf/x66GkkynVUlkiCrVsG228J22/nwsEy3SRPJgs6qEikTnTr5rkgefxxmzw5djSSVgkMkYU48Efr2hXHjIOE7FCQQBYdIwnTs6E/Lfe45f5aVSK4pOEQS6Oc/98c5LrwQ1q0LXY0kjYJDJIHat/e3mJ03z19RLpJLCg6RhBo2DHbayQfI2rWhq5EkSWxw6DoOKXft2sGECfDaa7DFFv6mTzU1UFsbujIpdYkNDnWrLgKrV/trOZYt82dY1dfDyJEKD2mbxAaHiGQ+JbehwV/rIdJaCg6RBFuwIF67SDYUHCIJ1rdvvHaRbCg4RBJs0iSoqlq/rWNH3y7SWgoOkQRLpWDmTKiu9gfJKyuhSxf48Y9DVyalTMEhknCpFNTV+SvIH34YFi+GKVNCVyWlTMEhUkYOOshfGDh5Mvzf/4WuRkqVgkOkzFx5pe9+/ZRT1HuutE5ig0NXjotktuWWcNll8OijMGtW6GqkFOkOgCJlaO1a2HtvfyX5G29A166hK5JiozsAish62rWDGTNg6VJdRS7xKThEytRuu8Fpp8H118Ozz4auRkqJgkOkjE2c6I95jBoFa9aErkZKhYJDpIxtuilcdRW8+CJMnx66GikVCg6RMnfMMXDIIf42s++/H7oaKQUKDpEyZ+a3Nr78Es46K3Q1UgoUHCLCttv6e3fcfbfvlkSkJQoOEQFg9GjYcUd/ptWqVaGrkWKm4BARwHe3ft118M47/spykeYkNjjU5YhIfAceCMcdB5dfDq+/HroaKVaJDQ7n3APOuZFdunQJXYpISZk6FTbeGE4+WZ0gSmaJDQ4RaZ0ttvD365gzB267LXQ1UowUHCLyDSNGwF57wdlnw/LloauRYqPgEJFvqKjwnSAuWwZjxoSuRoqNgkNEMtplFzjzTPj97+Gpp0JXI8VEwSEizbrkEth6a3+gXJ0gSiMFh4g0q3Nn3wnivHnQs6ffhVVTA7W1oSuTkCpDFyAixW31ah8YjQfJ6+th5Ej/PJUKV5eEoy0OEWnR2LGwbt36bQ0NunNgOVNwiEiLFiyI1y7Jp+AQkRb17RuvXZJPwSEiLZo0Caqq1m/r2NG3S3lScIhIi1IpmDkTqqv9TZ8qK2HzzeHYY0NXJqEoOERkg1IpqKvzB8lra2HRIp2SW84SGxzqVl0kP445BnbbDS66yN9uVspPYoND3aqL5EdFhb/RU12d34Ul5SexwSEi+TN4MOy/P0ycCJ99FroaKTQFh4jEZgaTJ8Pixb5LEikvCg4RaZW994bDD4ff/AY+/jh0NVJICg4RabVJk+CTT/w9yqV8KDhEpNV23tmfqnv11fD++6GrkUJRcIhIm1xyib9Xx8SJoSuRQlFwiEibbLMNnHQS3HQTvPVW6GqkEBQcItJm48ZBhw4wfnzoSqQQFBwi0mZbbunvTz5rFrz0UuhqJN8UHCKSE6NHQ9euusFTOVBwiEhOdOsG558PDz0ETzwRuhrJJwWHiOTM6afDVlvBmDHgXOhqJF8UHCKSM1VVcOGFfovj4YdDVyP5ouAQkZwaPtyfonvBBf7+HZI8Cg4RyakOHfzFgPPmwZ13hq5G8kHBISI5d+yxMGCA32311Vehq5FcU3CISM5VVPgOEN9+G26+OXQ1kmsKDhHJix/+EL77XZgwAVatCl2N5JKCQ0TyovFmT4sWwbXXhq5GcimxwWFmQ81s5sqVK0OXIlK29tsPDj3UB8iKFaGrkVxJbHA45x5wzo3s0qVL6FJEytqkSbB8OUydGroSyZXEBoeIFIddd/VnWU2bBh99FLoayQUFh4jk3YQJsHq13/qQ0qfgEJG82357f0X5jBnw7ruhq5G2UnCISEGMH+87Puzf31/nUVMDtbWhq5LWqAxdgIiUhzlz/M/PP/c/6+th5Ej/PJUKUpK0krY4RKQgxo6FNWvWb2to0I2fSpGCQ0QKYsGCeO1SvBQcIlIQffvGa5fipeAQkYKYNMnf6CldVZVO0S1FCg4RKYhUCmbOhOpq349VdbUf1oHx0qOzqkSkYFIpBUUSaItDRERiUXCIiEgsCg4REYlFwSEiIrEoOEREJBZzzoWuIa/MbAlQ34q3dgFyefvAXE+vrboDS5u0NVdjS7Vv6HNl+7njfD9t+S7zsRyKadlmWq6Q+2Wr5VpYhVqu6eOrnXM9Mr7COadHhgcws5inl4N65mZbY0u1b+hzZfu543w/bfku87EcimnZZlqu+Vi2Wq7JXK7Zfm7tqmreA0U+vXxorsaWat/Q58r2c8f5ftryXeZjOZTjstVyLQ5B/mYTv6tKMjOzuc65QaHrkNzSck2mYluu2uIoXzNDFyB5oeWaTEW1XLXFISIisWiLQ0REYlFwiIhILAoOERGJRcEhAJjZxmZ2q5ndYGbq+DqBzGwbM7vJzO4JXYvkjpkdGf3d3mdmgwsxTwVHgpnZzWa22MxeadI+xMzeMLO3zOz8qPlo4B7n3Ajg8IIXKy2KuSwzcs6945wbnt9KJY4cLdd7o7/b44H/zGO5/6bgSLZbgCHpDWbWDpgOHAr0A4aZWT+gN7AwetnaAtYo2bmFLJelmfU3swebPHoWvmTJwi3kbrmOi96Xd7oDYII55/5hZjVNmvcE3nLOvQNgZncARwDv4cPjJfQPRdGJsyydc5OBHxW2QmmNXCxXMzNgCvCwc+6F/FbsaQVRfrbm6y0L8IGxNfBfwI/N7HpKo6sFaX5ZZmRmm5vZDGBXMxuT7+Kk1WItV+B04CDgGDMblc/CGmmLo/xYhjbnnPscOKHQxUibZFyWzb3YObcMKMiKRdok7nK9Grg6f+V8k7Y4ys97QJ+04d7AokC1SNtoWSZT0S9XBUf5eQ7Y3sy+ZWYdgGOB+wPXJK2jZZlMRb9cFRwJZmazgKeAHc3sPTMb7pxbA5wG/A/wGnCXc+7VkHXKhmlZJlOpLld1cigiIrFoi0NERGJRcIiISCwKDhERiUXBISIisSg4REQkFgWHiIjEouAQEZFYFBwiIhKLgkMkADM7yMz+FLoOkdZQcIiEMRB4MXQRIq2h4BAJYyCwpZk9bmYfmtlBoQsSyZaCQySMgcBS59x+wClAKnA9IllTcIgUmJm1BzYDpkZNlcCKcBWJxKPgECm8fsA859y6aHgA8ErAekRiUXCIFN5AYF7a8ADg5UC1iMSm4BApvIGsHxQ7oy0OKSG6kZOIiMSiLQ4REYlFwSEiIrEoOEREJBYFh4iIxKLgEBGRWBQcIiISi4JDRERiUXCIiEgs/x9mnElSjiqVMwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "plt.plot(h_values, errors, \"bo-\")\n",
    "plt.xlabel(\"$h$\")\n",
    "plt.ylabel(\"$L^2$ error\")\n",
    "plt.xscale(\"log\")\n",
    "plt.yscale(\"log\")\n",
    "plt.axis(\"equal\")\n",
    "plt.xlim(plt.xlim()[::-1])\n",
    "\n",
    "plt.title(\"Convergence of an interior problem\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, we could plot the error against the number of degrees of freedom (DOFs). In the code above, we used `space.global_dof_count` to get these values and stored them in the list `ndofs`. This line has a different gradient to the plot above as the number of DOFs scales like $h^{-2}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEaCAYAAAAG87ApAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd5xU1fnH8c+zLCAriBowEWF3xYoaYlmjqBQVpShiLLGsCkhCNIklxdhQQENAk1iwxKACIS4YiYoFLCSIiGJZ0B/YsC5FQcEoFhQQzu+Pc1cuw2wZdmbvzJ3v+/Wa184tc+8zc2fvM+eee84x5xwiIiL1VRB1ACIikluUOEREJCVKHCIikhIlDhERSYkSh4iIpESJQ0REUqLEIbKVzBtvZp+a2YsZ3le5mT2ZyX0k7K+rmS1qrP2lysyqzKxnDct6mNmyxo4pnyhxxJSZnWlmlWb2pZktN7PHzOyIqOOKmSOAY4D2zrkfZ3JHzrkK59yx9VnXzAaa2ZwG7u8Z59xeDdmGxJcSRwyZ2W+Bm4A/Ad8HioHbgf5RxhVmZoVRx5AGJUCVc+6rqANJp4YeGzNrEuX+pRE45/SI0QNoDXwJnFrLOs3xieXD4HET0DxY1gNYBvwO+BhYDgwKlh0KrACahLb1E2BB8LwAuAx4F/gEuA/YMVhWCjhgMLAEmB3MPwdYHKx/FVAF9ExhewOC7a0CrgzF1QS4InjtF8A8oEOwbG9gBvA/YBHw01o+q3bAw8G67wA/D+YPBr4BNgSf94gkr90NmBnEvgqoALYPLa8Cfg8sAFYD/wK2qSGOgcCc0LQDzgPeBj4FbgMM6JQQ12ehY/6X4LP6CLgDaJFwzC8Nju8/q+eF9tcJmAV8BrwGnBBaNgH4GzAd+Kr6+CXEPwsYBbwYvNeH6vHdOCHY12fB6zslfHaXA68H73989WeXJPZ2wP3ASuB94MLQsuHAFOCe4HuyENgz2PbHwFLg2Kj/r7PtEXkAeqT5gEJv4FugsJZ1rgGeB3YC2gLPAdcGy3oEr78GaAr0BdYAOwTL3wWOCW1rCnBZ8PziYLvtgxPV34HJwbLqk8NEYFugBbBPcHI7AmgWnNjWsylx1Gd7dwbb+hGwtvrkAlwSnAT2wp9QfwR8L9j3UmAQUAgciD+p71vDZ/U0vrS2DbB/cPI5Olg2kNDJPMlrd8dfymoefM6zgZtCy6vwJ9J2wI7AG8B5NWxrs30F7/1RYHt8iXIl0LumuPA/Dh4O9tMKeAQYlXDMrwtibUHo5Bt8D97BJ+JmwFH4k+xewfIJ+GRwOD7Zb5H88Cf+D4D9gmNwP3BPLd+NPfFJ6Jhg/38IYmgW+uxeBToE7+lZ4I+h91MdewH+R8PVQewdgfeAXsHy4fhE2yv4PkzEJ5crg/3+HHg/6v/rbHtEHoAeaT6gUA6sqGOdd4G+oele+Esu1f90XxNKPPhfXocGz/8IjAuetwr+uUuC6TcITqrB9M74RFAYOjl0DC2/miARBNNFwDo2JY76bK99aPmLwOnB80VA/yTv/TTgmYR5fweGJVm3A/6Xe6vQvFHAhOD5QGpJHEm2dyLwcmi6CjgrNH09cEcNr91sX8F7PyI0fR+bEnjiuhYcp91C87pUnxCDY76O0AmfzU++XfElkYLQ8snA8OD5BGBiHe99FjA6NL1PsM8mNXw3rgLuC00X4BNPj9Bnd15oeV/g3SSxHwIsSYjlcmB88Hw4MCO0rB/+x0yT0HfcESop6uHQtcT4+QRoY2aFzrlva1inHf7yULXFwbzvtpHw2jVAy+D5JOA5MzsfOAmY75yr3lYJ8KCZbQy9dgO+nqXa0oQ4vpt2zq0xs09Cy+uzvRU1xNkBnyATlQCHmNlnoXmF+MszidoB/3POfRGatxgoS7LuFsxsJ2AM/sTbCn/y+zRhtcT421F/Nb33RG3xSXmemX0XHv6kXW2lc+6bGl7fDljqnAsfh8XALqHppdQtvM5i/C/6NjUs3+w76pzbaGZLa9ln4ne4WgnQLuF4NwGeCU1/FHr+NbDKObchNA3+sw1vI6+pcjx+5uKL3ifWss6H+H+oasXBvDo5517H/5P2Ac7EJ5JqS4E+zrntQ49tnHMfhDcRer4cfxkKADNrgb+clMr2arIUX8eQbP7TCdts6Zw7P8m6HwI7mlmr0Lxi/C/f+hiFf7+dnXPbAWfhT9iZ5hKmV+FPgPuG3nNr51zLWl4T9iHQwczC54vEz6G211frkPD69UFsybax2XfUfMbrkLDPxO0l+w4vxZeswse7lXOubz3ilRooccSMc241/hLQbWZ2opkVmVlTM+tjZtcHq00GhppZWzNrE6x/Twq7mQRcCHTD13FUuwMYaWYlAMH2a7uT699APzM7zMyaASPY/MSa6vbC7gKuNbM9gvYWnc3se/h6gT3N7Ozgc2lqZgebWafEDTjnluLrf0aZ2TZm1hlfgVtRzxhaEVRQm9ku+HqXxvAR0D74TAlKCncCNwalIMxsFzPrVc/tvYC/1PWH4PPqgb+kc2+KcZ1lZvuYWRG+Du3foV/2ie4DjjOzo82sKf5mjbX441HtV2bW3sx2xNe//CvJdl4EPjezS82shZk1MbP9zOzgFGOXECWOGHLO3QD8FhiKrzRdCvwamBqs8kegEn83z0JgfjCvvibjryPPdM6FfzHejK+AfdLMvsBXbB9SS5yvARfgT0DL8RWuH+NPEClvL8EN+JPPk8DnwN34u4i+AI4FTsf/Ql3BpkrhZM7AX4P/EHgQXxcyo54xjMBXvq8GpgEP1PN1DTUTfzfSCjOrPj6X4iuXnzezz4H/4G8cqJNzbh3+Dqc++BLC7cA5zrk3U4zrn/j6kBX4mw0urGWfi/AltFuCffYD+gWxVJuEP77vBY8tvsNBYuqHv7Hh/WBbd+HvPpStZEEFkEjkzKz6OvIezrn3o45H0sfMZuHvoror6lik4VTikEiZWb/gctq2+NtxF+LvmBGRLKXEIVHrz6aGiHvgb6dVMVgki+lSlYiIpEQlDhERSYkSh4iIpCT2LcfbtGnjSktLow5DRCSnzJs3b5Vzrm2yZbFPHKWlpVRWVkYdhohITjGzxTUt06UqERFJSWwTR9A+YOzq1aujDkVEJFZimzicc48454a0bq2eBURE0im2iUNERDJDiSPLVVRAaSkUFPi/FfXtl1VEJENif1dVLquogCFDYM0aP714sZ8GKC+PLi4RyW8qcWSxK6/clDSqrVnj54uIRCW2iSMOd1UtWZLafBGRxhDbxBGHu6qKi5PPb9Mm+XwRkcYQ28QRByNHQlHR5vPMYOVKGDgQcrgwJSI5TIkji5WXw9ixUFLiE0ZJCYwfD0OHwj33wA9/CDNnRh2liOQbJY4sV14OVVWwcaP/O2AAXHstPPsstGgBRx8NF18MX38ddaQiki+UOHLUIYfAyy/DBRfAzTfDgQfCSy9FHZWI5AMljhxWVARjxsCMGfDll9ClCwwbBuvXRx2ZiMRZbBNHHG7Hra+ePWHhQjjzTLjmGjj0UHj99aijEpG4im3iiMPtuKnYfnuYOBHuv9+38zjwQLjxRl83IiKSTrFNHPnqpJPg1VehVy/47W/hqKN8pbqISLooccTQ978PU6fCuHEwfz507uyfOxd1ZCISB0ocMWUGgwbBggX+stXgwdC/P3z0UdSRiUiuU+KIudJS30jwxhvhySdhv/18PYiIyNZS4sgDBQW+keD8+b71+SmnwNlnw2efRR2ZiOQiJY48ss8+MHeub+sxebLvsmTGjKijEpFcE9vEkU/tOFLRtCkMH+4TSMuWcOyxvvV54rgfIiI1iW3iyLd2HKk6+GB/6erii+HWW+GAA+CFF6KOSkRyQWwTh9StRQtfaT5zJnzzDRx2mO95d926qCMTkWymxCEceaS/bfecc/wYIIcc4hsRiogko8QhALRu7cf6mDoVPvgADjoI/vIX2LAh6shEJNsocchm+vf3pY3jjoNLLvGlkffeizoqEckmShyyhZ128o0E//EP+L//812W3HmnuiwREU+JQ5Iy83UeCxf6Oo8hQ+D442H58qgjE5GoKXFIrYqLfSPBMWP83Vf77Qf33Rd1VCISpdgmDjUATJ+CAt9I8OWXYbfd4LTT/KBR//tf1JGJSBRimzjUADD99t4bnnvOjzI4ZYrvsuSJJ6CiwnemWFDg/1ZURB2piGRSYdQBSG4pLISrroK+fX0dSO/eft633/rlixf7+hCA8vLo4hSRzIltiUMy66CDYN48aNVqU9KotmYNXHllNHGJSOYpcchW22Yb+PLL5MuWLGncWESk8ShxSIMUFyef37IlrFrVuLGISONQ4pAGGTkSioo2n1dY6Esiu+8Of/6z70BRROJDiUMapLwcxo71Iwua+b8TJvhuS444Av7wB+jUCe69Vy3PReJCiUMarLwcqqpg40b/t7zcjzb46KO+8WDr1nDGGdCli7+dV0RymxKHZFTPnv7uq3HjfIX54YfDqafCu+9GHZmIbC0lDsm4Jk1g0CB4+20/bO306f7y1e9+B59+GnV0IpIqJQ5pNNtuC8OG+QRy9tl+9MHddoObbtKogyK5JLaJQ31VZa927eDuu+GVV6CsDH7zG9h3X3jgAVWgi+SC2CYO9VWV/Tp39n1dTZ8OzZrBySdD9+7w0ktRRyYitYlt4pDcYAZ9+vgBo+64AxYtgh//2N+ZtXixOlAUyUbmYn5toKyszFVWVkYdhtTT55/DddfBDTfA+vU+sYT7wioq8u1G1IGiSGaZ2TznXFmyZSpxSFbZbjvfGv2tt3xfWOpAUST7KHFIVurQwSeJZNSBoki0lDgka9XUgWKHDo0bh4hsTolDslayDhTB3867dm3jxyMinhKHZK3EDhSLi/14588/D716qdW5SFSUOCSrhTtQXLzY97JbUeE7S+zaFZYujTpCkfyjxCE558wz4fHHfdLo0gUWLow6IpH8osQhOemoo+CZZ/zzI46AmTOjjUcknyhxSM7q3BnmzvV3WfXuDZMmRR2RSH5Q4pCc1qEDzJkDhx3m60Ouv14dJYpkmhKH5Lztt/edJf70p3DppXDhhbBhQ9RRicRXYdQBZIqZ9QP67b777lGHIo2geXOYPBnat/f9XH3wgb/7qkWLqCMTiZ/YljjUrXr+KSiAv/7VDxA1dSoccwx88knUUYnET2wTh+Sviy+Gf/0LKiv9GOdVVVFHJBIvShwSS6eeCjNmwEcf+bYeL78cdUQi8aHEIbHVtSs8+yw0bQrduvkKdBFpOCUOibV99vF9W3XsCMcfDxMmRB2RSO5T4pDYa9fOtzLv3h0GDfK97qqth8jWU+KQvLDddjB9Opx1FgwdCuefv+XogiJSP7FtxyGSqFkzmDjRt/UYPRo+/NC3/dh226gjE8ktKnFIXjGDUaPgtttg2jTfWeLKlVFHJZJblDgkL/3yl/DAA7Bgge/n6t13o45IJHcocUje6t/fd8f+6ae+rceLL0YdkUhuUOKQvNali2/r0bIlHHmkv3wlIrVT4pC8t9defijaTp18KeSuu6KOSCS7KXGIAD/4Acya5TtG/PnPYdgwtfUQqYkSh0igZUt4+GHfSPCaa2DwYFi/PuqoRLKP2nGIhDRtCnffDcXFMGIELF8OU6b4pCIinkocIgnMYPhwuPNO38Nu9+6wYkXUUYlkDyUOkRr87Gfw0EPw5pv+7qtFi6KOSCQ7KHGI1OK443yl+Vdf+UGh5s6NOiKR6MU2cZhZPzMbu3r16qhDkRx38ME+Yeywg++i5De/gdJSP1Rtaakf21wkn5iL+T2HZWVlrrKyMuowJAZWrvSXrBK7JykqgrFjobw8mrhEMsHM5jnnypIti22JQyTd2raFdeu2nL9mDVx5ZePHIxIVJQ6RFCxblnz+kiWNG4dIlJQ4RFJQXJx8/i67NG4cIlFS4hBJwciRvk4j0caNsHhx48cjEgUlDpEUlJf7ivCSEt9QsKTE12989RUceijMnx91hCKZp8QhkqLycqiq8qWMqir44x991+zNmkG3bvDYY1FHKJJZShwiabDvvr6txx57QL9+6ppd4k2JQyRN2rWD2bOhZ0/fNfvVV6trdomneicO8zpkMhiRXNeqFTzyiO+S/dprYcCA5G0/RHJZvROH803Mp2YwFpFYaNrU96x77bXwz39C376gnm8kTlK9VPW8mR2ckUhEYsQMhg6Ff/wDnn4aunatufGgSK5JNXEcCcw1s3fNbIGZLTSzBZkITCQOzjnH32VVVeVv112g/xaJgVRHAOyTkShEYqxnT5gzx1+yOuIIuP9+P7a5SK5KqcThnFsMbA/0Cx7bB/NEpBadO8Pzz/tu2Pv2hQkToo5IZOullDjM7CKgAtgpeNxjZhdkIjCRuGnfHp55Bnr0gEGD4JprdLuu5KZUL1UNBg5xzn0FYGbXAXOBW9IdmEgctW4N06b5dh7Dhvledf/2N38nlkiuSDVxGLAhNL0hmCci9dSsmb9UVVrqSx3LlsGUKb4NiEguSDVxjAdeMLMHg+kTgbvTG5JI/JnBiBG+m/Zf/ML3cTVtmm99LpLtUmo5DkwBBgH/Az4FBjnnbspQbCKxN3iwTxjvvOOHpX3ttagjEqlbyi3HnXPznXNjnHM3O+dezmBsInmhVy/fx9X69XD44fDUU1FHJFI7tRwXyQIHHOB7191lF59IJk2KOiKRmqnluEiWKCnxDQUPP9yP+TFqlG7XlexU78rxoI7jPEAN/kQyZIcd4PHH4dxz4Yor/O26t9wChanexiKSQfX+OjrnnJnd6Jw7KJMBieS75s19r7olJb7UsXQp3HsvtGwZdWQinuo4RLJQQQH86U++ceBjj/nW5itWRB2ViLc1dRzPq45DpHGcdx489BC88Ya/XffPf/YNBwsK/N+KiqgjlHxkLoXaNzMrSTY/mzs6LCsrc5WVlVGHIdIglZVw1FHwxRebzy8qgrFjfWW6SDqZ2TznXFmyZamWOJYAXYEBQbJwwPcbGJ+I1KGsLHmXJGvWwJVXNn48kt9STRy3A12AM4LpL4Db0hqRiCS1fHny+UuWNG4cIqkmjkOcc78CvgFwzn0KNEt7VCKyheLi1OaLZEqqiWO9mTXBX6LCzNoCG9MelYhsYeRIX6eR6KijGj8WyW+pJo4xwIPATmY2EpgD/CntUdXCzDqa2d1m9u/G3K9I1MrLfUV4SYnvXbdDB9h/fxg/Hm6+OeroJJ+kOnRsBfAHYBSwHDjROTelvq83s3Fm9rGZvZowv7eZLTKzd8zssjpieM85NziVuEXiorwcqqpg40Zft/HCC3DyyXDxxTB6dNTRSb5IuSMD59ybwJtbub8JwK3AxOoZwaWv24BjgGXAS2b2MNAEn6DCznXOfbyV+xaJnWbNfKvyAQPg8svhm2/8yIKm4dUkgxq1Bxzn3GwzK02Y/WPgHefcewBmdi/Q3zk3Cjh+a/ZjZkOAIQDFqjmUmCsshIkTfVclI0b45DFqlJKHZE6qdRyZsAuwNDS9LJiXlJl9z8zuAA4ws8uTreOcG+ucK3POlbVt2za90YpkoSZN4K674Pzz4brr/KUr9awrmZINfW4m+11U41feOfcJvpdeEQkpKIDbbvMlj5tugrVr4fbb/XyRdKrzK2Vmx5jZnWa2fzA9JM0xLAM6hKbbAx+meR8iecEMbrjB13f8/e9+aNoNG6KOSuKmPiWOX+LHGR9qZjsC+6c5hpeAPcxsV+AD4HTgzDTvQyRvmPk2Hy1awNVX+zqPiROhadOoI5O4qE8hdqVz7jPn3O+BY4Gt7lbdzCYDc4G9zGyZmQ12zn0L/Bp4AngDuM8599rW7iO0r35mNnb16tUN3ZRIzjGDq67y9R333gunnQbr1kUdlcRFnb3jmll/59xDoekLnHO3ZDyyNFHvuJLvxoyBiy6C446Df/8bttkm6ogkFzSod9xw0gimcyZpiAhceKGv75g+Hfr18z3qijREg++3MLNn0xGIiGTOkCG+a5KZM6FPny3H9RBJRTpu1GuXhm2ISIYNGACTJsGzz8Kxx8Jnn0UdkeSqeiUOM7vFzIaYWRczSxxORs2MRHLEaafBlCkwbx4cfTR88knUEUkuqm+JYyHQGRgNVJnZ+2b2cNBDbpJxyaKnu6pEkvvJT2DqVHjtNTjySPhYvb9JilIac/y7F5m1xyeSHwI/dM6dle7A0kV3VYkk95//wAkn+G7a//tfaKeLzhKSzjHHAXDOLXPOTXfOXZfNSUNEatazJzz+OCxbBt26aQhaqT/1YiOSx7p1gxkzYNUq//y996KOSHKBEodInjv0UH+p6osvfPJYtCjqiCTbKXGICAcdBE895bsl6d7dV5yL1ESJQ0QA6NwZnn7ad8Peowe88krUEUm2im3i0O24Iqnr1Almz/Y96x55JLz4YtQRSTaKbeJwzj3inBvSunXrqEMRySm77+6Tx447+juvnlWnQpIgtolDRLZeaalPHjvvDL16+foPkWpKHCKS1C67+DqP0lLo29e3+RABJQ4RqcUPfgCzZsHee0P//vDww1FHJNlAiUNEatWmje+O/Uc/gpNP9p0kSn5T4hCROu2wg+/b6pBD4PTT4Z57oo5IoqTEISL1st12vp6je3c45xy4++6oI5KoxDZxqB2HSPq1bAnTpvk7rX72M3/LbkGBr0CvqIg6OmkssU0caschkhktWvjLVU2awKefgnOweLEfnlbJIz/ENnGISOYMGwYbNmw+b80auPLKaOKRxqXEISIpq2nsDo3pkR+UOEQkZcXFyec3b+67Z5d4U+IQkZSNHAlFRZvPa9YM1q71d10tXx5NXNI4lDhEJGXl5TB2rB+v3Mz/HTfO33H11lvQpQu8+WbUUUqmmHMu6hgyqqyszFVWVkYdhkjeqKyE446Db7+FRx6Bww6LOiLZGmY2zzlXlmyZShwiklZlZTB3Lnzve3D00TB1atQRSbopcYhI2nXs6MfxqO7f6vbbo45I0im2iUMtx0Wi1bat7xzxuOPgV7+Cyy/3jQUl98U2cajluEj0iorggQd8q/LRo2HAAFi3LuqopKEKow5AROKtsBDuuAM6dICrrvK36t5/v+80UXJTbEscIpI9zGDoUBg/3g9Dq7YeuU2JQ0QazcCB8Oij8PbbauuRy5Q4RKRR9e7txzL/+ms4/HB/95XkFiUOEWl0Bx20qa1Hz57w4INRRySpUOIQkUh07AjPPbeprcett0YdkdSXEoeIRKZNG9/Wo18/uOACuOwy2Lgx6qikLkocIhKpoiJ/e+5558F116mtRy5QOw4RiVxhoe+WpEMHP4rgihVq65HNYlviUJcjIrnFDK64wrf1mDULunWDDz+MOipJJraJQ12OiOSm6rYe77zj23q88UbUEUmi2CYOEcldvXr5th5r1/q2HnPmRB2RhClxiEhWqm7r0batb+vxwANRRyTVlDhEJGvtuqtvWX7ggXDKKWrrkS2UOEQkq7VpA//5D5xwgm/rcemlausRNSUOEcl61W09zj8frr8ezj5bbT2ipMQhIjmhSRO47TYYORImTYK+feGuu6C0FAoK/N+KiqijzA9qACgiOaO6rUf79v623ZkzNw1Hu3ixH2kQoLw8shDzgkocIpJzzjnH322VOIb5mjW+5blklhKHiOSklSuTz1+ypHHjyEdKHCKSk4qLU5sv6aPEISI5aeRIf7dVotNPb/xY8o0Sh4jkpPJyGDsWSkp8pXn79r533RtvhKlTo44u3pQ4RCRnlZdDVZVvELh0KbzyChxwgG9lPmlS1NHFV2wTh7pVF8k/O+4IM2ZA165w1lm+RCLpF9vEoW7VRfJTq1YwfTr07g2/+AXccEPUEcVPbBOHiOSvFi18Pccpp8DvfgfXXLNlmw/Zemo5LiKx1KwZTJ4M224Lw4bBF1/4fq7Moo4s9ylxiEhsFRbCuHHQsiX85S/w5Ze+v6sCXWtpECUOEYm1ggK45RafPK67zieP8eN9UpGto49ORGLPDEaPhu22831ZffWVv4zVvHnUkeUmFdhEJG9ccQXcdBM8+CD07+87RZTUKXGISF656CK4+2548kl/y+7nn0cdUe5R4hCRvHPuuf5S1dy50LMnfPJJ1BHlFiUOEclLp50GDzwACxZAjx6wYkXUEeUOJQ4RyVv9+sG0afDee9Ctm8byqC8lDhHJa0cf7fu3+vhj38fV229HHVH2U+IQkbx32GF+/PI1a3zJ49VXo44ouylxiIgABx4ITz/t23x07w6VlVFHlL2UOEREAvvsA3Pm+IaCRx0FzzwTdUTZSYlDRCSkY0efMNq1g169fHsP2ZwSh4hIgvbtYfZs2HNPf+eVhqLdnBKHiEgSO+0ETz21aSjaioqoI8oeShwiIjXYYYdNQ9GefbaGoq2mxCEiUovqoWj79NFQtNVimzjMrJ+ZjV29enXUoYhIjmvRwveoe+qpfijaESPyeyja2CYO59wjzrkhrVu3jjoUEYmBZs1g0iQYOBCGD4dLLsnf5KGBnERE6qmw0HfJ3rIl/PWvfjTB22/Pv6FolThERFJQUABjxvjkMXq0H00w34aizaO3KiKSHmYwapSvOM/HoWjzrIAlIpI+V1wBN9/sK87LyqC42JdISkvj3e5DJQ4RkQa48EJYuBDuumvTvMWLYcgQ/7y8PJq4MkklDhGRBpoxY8t5a9b4y1hxpMQhItJANY0cGNcRBZU4REQaqLg4tfm5TolDRKSBRo6EoqLN5xUV+flxpMQhItJA5eW+A8SSEn+rbkmJn45jxTjorioRkbQoL49vokikEoeIiKREiUNERFKixCEiIilR4hARkZQocYiISErMxXwkEjNbCSwOJlsDqQwJWJ/161qntuU1LUs2vw2wqo5YGlOqn2Wmt5nKa7PpuEJ2Hdu4H9e61tNx3aTEOdc26RLnXN48gLHpXr+udWpbXtOyZPOByqg/v4Z8lpneZiqvzabjmm3HNu7HdWuPnY7r5o98u1T1SAbWr2ud2pbXtCzVOKOQiRgbss1UXqvjWrO4H9e61tNxrYfYX6qKCzOrdM6VRR2HpJ+ObTzF+bjmW4kjl42NOgDJGB3beIrtcVWJQ0REUqISh4iIpESJQ0REUqLEISDCOQUAAAeoSURBVCIiKVHiyFFmdqKZ3WlmD5nZsVHHI+lhZp3M7A4z+7eZnR91PJI+Zratmc0zs+OjjqWhlDiyiJmNM7OPzezVhPm9zWyRmb1jZpcBOOemOud+DgwETosgXKmnFI/rG86584CfArG8lTMuUjmugUuB+xo3ysxQ4sguE4De4Rlm1gS4DegD7AOcYWb7hFYZGiyX7DWBFI6rmZ0AzAH+27hhSoomUM/jamY9gdeBjxo7yExQ4sgizrnZwP8SZv8YeMc5955zbh1wL9DfvOuAx5xz8xs7Vqm/VI5rsP7DzrnDgDwZTy43pXhcjwQOBc4Efm5mOX3u1dCx2W8XYGloehlwCHAB0BNobWa7O+fuiCI42WpJj6uZ9QBOApoD0yOISxom6XF1zv0awMwGAquccxsjiC1tlDiynyWZ55xzY4AxjR2MpE1Nx3UWMKtxQ5E0Snpcv3vi3ITGCyVzcrq4lCeWAR1C0+2BDyOKRdJHxzWe8uK4KnFkv5eAPcxsVzNrBpwOPBxxTNJwOq7xlBfHVYkji5jZZGAusJeZLTOzwc65b4FfA08AbwD3OedeizJOSY2Oazzl83FVJ4ciIpISlThERCQlShwiIpISJQ4REUmJEoeIiKREiUNERFKixCEiIilR4pC0MjNnZn8NTf/ezIanadsTzOyUdGyrjv2camZvmNlT2RBPuplZWzN7wcxeNrOuCcu6mtlrZvaKmbVI8357mNmj6dymREOJQ9JtLXCSmbWJOpCwoLvr+hoM/NI5d2Sm4gkzs8buM+5o4E3n3AHOuWcSlpUDf3HO7e+c+7p6Zoqfn8ScEoek27fAWOA3iQsSf6Gb2ZfB3x5m9rSZ3Wdmb5nZaDMrN7MXzWyhme0W2kxPM3smWO/44PVNzOzPZvaSmS0ws1+EtvuUmU0CFiaJ54xg+68GXdRjZlcDRwB3mNmfE9Y3M7vVzF43s2nATqFlBwXvYZ6ZPWFmOwfzDw5imhvE+Gowf6CZTTGzR4Ang3mXhN7DiNC2zwo+i1fM7O/B+20SfJ6vBu8h2eddYmb/Dbb3XzMrNrP9geuBvomlCjP7GX4AqavNrCLZ55cslmD+scF7nB+8r5bB/N5m9qaZzcH3+lu9rx3NbGoQ2/Nm1jmYP9zM/mFmT5pZlZmdZGbXB+/xcTNrmvg+JQLOOT30SNsD+BLYDqgCWgO/B4YHyyYAp4TXDf72AD4DdsZ3J/4BMCJYdhFwU+j1j+N/8OyB71BuG2AIMDRYpzlQCewabPcrYNckcbYDlgBt8b1EzwRODJbNAsqSvOYkYAbQJHj9Z8ApQFPgOaBtsN5pwLjg+avAYcHz0cCrwfOBQfw7BtPH4hOuBe/vUaAb0Al4BGgarHc7cA5wEDAjFNv2SeJ9BBgQPD8XmBra9601HL/vjlHi51dLLG2A2cC2wfxLgauDY7M0OFaGH/3u0WCdW4BhwfOjgFeC58Pxg1g1BX4ErAH6BMserD5GekT7ULfqknbOuc/NbCJwIfB1XesHXnLOLQcws3cJfoXjf+mGLxnd5/xYBm+b2XvA3viTbudQaaY1/mS1DnjROfd+kv0dDMxyzq0M9lmBP1FPrSXGbsBk59wG4EMzmxnM3wvYD5hhZuATy3Iz2x5o5Zx7LlhvEhAeb3qGc656IKBjg8fLwXTL4D10xieJl4JttwA+xp/AO5rZLcC00OcV1oVNv/L/iS9ppCr8+R1dQyyH4ke7ezaY3wzfh9PewPvOubcBzOwefJIHX6o7GcA5N9PMvmdmrYNljznn1pvZQvxn+XgwfyFQuhXvQdJMiUMy5SZgPjA+NO9bgsuj5s8wzULL1oaebwxNb2Tz72li52oO/2v2AufcE+EF5gdF+qqG+JKNm1AfyTp3M+A151yXhP3vUMe2wrEZMMo59/eEbVwA/MM5d/kWOzX7EdAL+BX+EtO5WxF7XRJj3CIWM+uHT4JnJMzfv5Z91jZuxVoA59xGM1vvguIGW34XJCKq45CMCH5J34evaK5Whf/FCn44za25Xn2qmRUE9R4dgUX4nkjPr77+bWZ7mtm2dWznBaC7mbUJrtOfATxdx2tmA6cH9Qs7s6kktAhoa2Zdgv03NbN9nXOfAl+Y2aHBeqfXsu0ngHNDdQO7mNlO+HHHTwmeV9cNlJi/+aDAOXc/cBVwYJJtPhfaZzn+ElBDJI0FeB443Mx2D+YXmdmewJvArrapjiqcWGYHMVUn+FXOuc8bGJ80EmVvyaS/4ruYrnYn8JCZvYg/CdVUGqjNIvwJ/vvAec65b8zsLvwljPlBSWYlcGJtG3HOLTezy4Gn8L9+pzvnHqpj3w/ir8cvBN4K4sA5ty64TDYmuNxSiC9xvYZPnHea2Vf4upPVNcTzpJl1AuYGl3u+BM5yzr1uZkOBJ82PU70eX8L4Ghhvm8au3qJEgr9UOM7MLgk+k0F1vL9a1RSLc+5580OiTjaz5sHqQ51zb5nZEGCama3CJ679guXDg/gX4OsxBjQkNmlc6lZdJIPMrKVzrvruscuAnZ1zF0UclkiDqMQhklnHBSWbQmAx/o4mkZymEoeIiKREleMiIpISJQ4REUmJEoeIiKREiUNERFKixCEiIilR4hARkZT8P2ufK0QabhiBAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "plt.plot(ndofs, errors, \"bo-\")\n",
    "plt.xlabel(\"Number of degrees of freedom\")\n",
    "plt.ylabel(\"$L^2$ error\")\n",
    "plt.xscale(\"log\")\n",
    "plt.yscale(\"log\")\n",
    "plt.axis(\"equal\")\n",
    "\n",
    "plt.title(\"Convergence of an interior problem\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this final plot, we add a trend line to indicate order 1.5 convergence (ie $\\text{error} = ch^{1.5}$).\n",
    "\n",
    "We do this by plotting $ch^{1.5}$ for appropriate values of $h$ and a value of $c$ chosen to put the trend near our error line. As the trend line will be a straight line on a log-log plot, we only need to tell matplotlib the first and last coordinates on the line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEcCAYAAADQqlM0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deZgU5bn38e89MAij0KKiccBhFDcQRY2KJm6JygyeF4159UQdj2ERXqJGc4waEUEEEQJhXKKGTUEjeGL0xLjExBiTuERlUQEDGpRdxYjIuCD78/7x1EgzzMA001XVXf37XFdfM13dXffdXT11T9WzlDnnEBERaayiuBMQEZH8osIhIiIZUeEQEZGMqHCIiEhGVDhERCQjKhwiIpIRFQ6RXWTeFDP71MxmhByrysyeDTNGnXinmNk7UcXLlJktMbMzG3jsdDNbEXVOhUSFI6HM7GIzm2VmX5jZh2b2jJmdHHdeCXMycBbQwTl3QpiBnHPTnHM9GvNcM+ttZi81Md6LzrnDmrIOSS4VjgQys2uAO4DbgP2AMuBe4Nw480pnZs3jziELOgJLnHNfxp1INjV125hZszjjSwScc7ol6AakgC+AC3bwnN3wheWD4HYHsFvw2OnACuCnwL+BD4E+wWMnAiuBZmnrOg+YG/xeBNwAvAd8AjwC7BU8Vg44oB+wDHghWH4psDR4/hBgCXBmBuv7YbC+VcDgtLyaATcGr/0cmA0cEDx2OPBnYDXwDvCfO/isSoEngue+C/QPlvcD1gGbg8/7lnpe2wl4Psh9FTAN2DPt8SXAtcBcoAb4DdCygTx6Ay+l3XfAQGAh8ClwD2BA5zp5rUnb5r8IPquPgPFAqzrb/GfB9v117bK0eJ2BvwFrgH8C56Q9NhX4FfAH4Mva7Vcn/78Bo4AZwXv9fSO+G+cEsdYEr+9c57MbBMwP3v+U2s+untxLgceAj4HFwFVpjw0Dfgs8FHxP5gGHBuv+N7Ac6BH333Wu3WJPQLcsb1CoBDYBzXfwnOHAq8C+QDvgH8CI4LHTg9cPB4qBs4G1QNvg8feAs9LW9VvghuD3nwTr7RDsqCYADweP1e4cHgR2B1oBXYKd28lAi2DHtpGthaMx65sUrKsbsL525wJcF+wEDsPvULsBewexlwN9gObAsfid+hENfFZ/xx+ttQSODnY+ZwSP9SZtZ17Paw/Gn8raLficXwDuSHt8CX5HWgrsBSwABjawrm1iBe/9KWBP/BHlx0BlQ3nh/zl4IojTGngSGFVnm/88yLUVaTvf4HvwLr4QtwC+i9/JHhY8PhVfDL6NL/bbFT/8jv99oGuwDR4DHtrBd+NQfBE6K4h/fZBDi7TP7i3ggOA9vQzcmvZ+anMvwv/TMDTI/SBgEVARPD4MX2grgu/Dg/jiMjiI2x9YHPffda7dYk9AtyxvUKgCVu7kOe8BZ6fdr8Cfcqn9o/uKtMKD/8/rxOD3W4H7g99bB3/cHYP7Cwh2qsH9/fGFoHnazuGgtMeHEhSC4H4JsIGthaMx6+uQ9vgM4MLg93eAc+t57z8AXqyzbAJwcz3PPQD/n3vrtGWjgKnB773ZQeGoZ33fA95Iu78EuCTt/hhgfAOv3SZW8N5PTrv/CFsLeN3nWrCdOqUtO6l2hxhs8w2k7fDZdud7Cv5IpCjt8YeBYcHvU4EHd/Le/waMTrvfJYjZrIHvxhDgkbT7RfjCc3raZzcw7fGzgffqyb07sKxOLoOAKcHvw4A/pz3WC//PTLO077gj7UhRN4fOJSbPJ8A+ZtbcObepgeeU4k8P1VoaLPt6HXVeuxbYI/h9OvAPM/sR8H3gdedc7bo6Ar8zsy1pr92Mb2eptbxOHl/fd86tNbNP0h5vzPpWNpDnAfgCWVdHoLuZrUlb1hx/eqauUmC1c+7ztGVLgePqee52zGxf4C78jrc1fuf3aZ2n1c2/lMZr6L3X1Q5flGeb2dfp4XfatT52zq1r4PWlwHLnXPp2WAq0T7u/nJ1Lf85S/H/0+zTw+DbfUefcFjNbvoOYdb/DtToCpXW2dzPgxbT7H6X9/hWwyjm3Oe0++M82fR0FTY3jyfMK/tD7ezt4zgf4P6haZcGynXLOzcf/kfYELsYXklrLgZ7OuT3Tbi2dc++nryLt9w/xp6EAMLNW+NNJmayvIcvxbQz1Lf97nXXu4Zz7UT3P/QDYy8xapy0rw//n2xij8O/3KOdcG+AS/A47bK7O/VX4HeARae855ZzbYwevSfcBcICZpe8v6n4OO3p9rQPqvH5jkFt969jmO2q+4h1QJ2bd9dX3HV6OP7JK396tnXNnNyJfaYAKR8I452rwp4DuMbPvmVmJmRWbWU8zGxM87WHgJjNrZ2b7BM9/KIMw04GrgFPxbRy1xgMjzawjQLD+HfXkehToZWbfMrMWwC1su2PNdH3pJgMjzOyQYLzFUWa2N75d4FAz+6/gcyk2s+PNrHPdFTjnluPbf0aZWUszOwrfgDutkTm0JmigNrP2+HaXKHwEdAg+U4IjhUnA7cFREGbW3swqGrm+1/Cnuq4PPq/T8ad0/ifDvC4xsy5mVoJvQ3s07T/7uh4B/sPMzjCzYnxnjfX47VHrCjPrYGZ74dtfflPPemYAn5nZz8yslZk1M7OuZnZ8hrlLGhWOBHLOVQPXADfhG02XA1cCjwdPuRWYhe/NMw94PVjWWA/jzyM/75xL/4/xTnwD7LNm9jm+Ybv7DvL8J/Bj/A7oQ3yD67/xO4iM11dHNX7n8yzwGXAfvhfR50AP4EL8f6gr2dooXJ+L8OfgPwB+h28L+XMjc7gF3/heAzwN/G8jX9dUz+N7I600s9rt8zN84/KrZvYZ8By+48BOOec24Hs49cQfIdwLXOqcezvDvH6Nbw9Zie9scNUOYr6DP0L7ZRCzF9AryKXWdPz2XRTctvsOB4WpF75jw+JgXZPxvQ9lF1nQACQSOzOrPY98iHNucdz5SPaY2d/wvagmx52LNJ2OOCRWZtYrOJ22O7477jx8jxkRyVEqHBK3c9k6EPEQfHdaHQaL5DCdqhIRkYzoiENERDKiwiEiIhlJ/MjxffbZx5WXl8edhohIXpk9e/Yq51y7+h5LfOEoLy9n1qxZcachIpJXzGxpQ4/pVJWIiGQksYUjGB8wsaamJu5UREQSJbGFwzn3pHNuQCqlmQVERLIpsYVDRETCocIhGZs2DcrLoajI/5zW2LliRSQREt+rSrJr2jQYMADWrvX3ly719wGqquLLS0SioyMOycjgwVuLRq21a/1yESkMiS0c6lUVjmXLMlsuIsmT2MKhXlXhKCvLbLmIJE9iC4eEY+RIKCnZdlmrVn65iBQGFQ7JSFUVTJwIHTuCBVcH/9a31DAuUkhUOCRjVVWwZAls2QLXXAN/+Qu8+GLcWYlIVFQ4pEmGD/djOfr3h3Xr4s5GRKKgwiFNsvvuMGECvPOO2jlECkViC4e640anRw+49FIYPRrmzYs7GxEJW2ILh7rjRqu6Gtq29aesNm+OOxsRCVNiC4dEa++94c474bXX4J574s5GRMKkwiFZc+GF0LMn3Hijn8NKRJJJhUOyxgx+9Sv/+49+BM7Fm4+IhEOFQ7KqY0e47TZ45hl4+OG4sxGRMKhwSNZdcQV07w5XXw2rVsWdjYhkmwqHZF2zZjB5MqxZAz/9adzZiEi2JbZwaBxHvLp2hUGD4MEH4dln485GRLLJXMJbMI877jg3a9asuNMoSOvWwdFHw/r18NZbfpS5iOQHM5vtnDuuvscSe8Qh8WvZEiZN8hMiDh0adzYiki0qHBKqU06BgQPhjjtg5sy4sxGRbFDhkNCNHg3f+IafjmTjxrizEZGmUuGQ0KVScO+9MGcOjBsXdzYi0lQqHBKJc8+F88+HYcNg4cK4sxGRplDhkMjcdZdvMO/f3189UETykwqHRGb//eEXv4C//x3uvz+eHJxzzFQrvUiTJLZwaABgburXD04/Ha69Fj78MPr4t99+OyeeeCIzZsyIPrhIQiS2cOhCTrnJDCZOhC++gIMPhqIif83yadOiid+vXz9KS0vp3bs363SRdJFdktjCIblrxgxfMNau9VOvL10KAwZEUzxSqRSTJk1iwYIF3HLLLeEHFEkgFQ6J3ODB24/nWLvWL49CZWUlffv2ZcyYMTplJbILmsedgBSeZcsyWx6G6upqFixYwFdffRVdUJGEUOGQyJWV1X9p2eJi+Ne/4NBDw88hlUrx8ssvY2bhBxNJGJ2qksiNHAklJdsu2203XziOPRamTInmsrNmxoYNGxg6dKi66IpkQIVDIldV5XtWdezoe1l17Aj33Qdvvw3HHw99+8JFF/kLQYXtq6++YsqUKeplJZIBFQ6JRVWVn259yxb/s6oKOnSA556DUaPgscegWzd46aVw80ilUkyePJn58+erl5VII6lwSE5p1gxuuAFefhmaN4fTTvPzW23aFF7MiooK+vXrp15WIo2kwiE56YQT4I034JJL4JZb/Gjz+hrUs2XcuHGUlpYyYMAAkn5VTJGmSmzh0JQj+a9NG3jgAT8wcN48f+rqN78JJ1YqlWL69OlMnTpVPa1EdkLXHJe8sHgxXHwxvPoq9O7tZ9pt3Tq8eF988QV77LFHeAFEcpyuOS5578AD4YUXYMgQePBB32131ix/NFJent05rwYNGkT37t3Vy0qkASockjeKi2H4cPjrX2HdOt8O0qePb/vI5pxXp512mnpZieyACofknVNPhblzoVWrcOa80lxWIjumNg7JW0VF9Y8wN2v6FQZramro2rUrbdq0Yfbs2bRs2bJpKxTJM2rjkEQqK8tseSZqp19fsWIFc+fObfoKRRJEhUPyVn1zXpWU+OXZUFlZyZIlSzjhhBOys0KRhFDhkLyVPudVrR//2C/PlrZt2+KcY/r06eplJRJQ4ZC8Vjvn1bp1vsvuH//Y9PaNul577TWqqqrUy0okoMIhibDbbv4U1Zw52b8E7Yknnvj1XFaafl1EvaokQbZs8WM7Pv4Y3nkHstkRSr2spNCoV5UUhKIiGDPGX4L2nnuyu25Nvy6ylQqHJMp3vwuVlf601aefZnfdFRUVDB48mNNPPz27KxbJMzpVJYkzdy4cfTRce60/AhGRzBXkqSpNq164jjoK/uu//Ay6y5Zlf/3OOYYOHcqQIUOyv3KRPJDYwuGce9I5NyCVSsWdisRgxAj/c+jQ7K/bzPjggw+47bbb1MtKClJiC4cUtrIyuOoqPwV7GDOG1F4xsHfv3hoYKAVHhUMSa9Ag2HNP+NnPsr/u2rms5s+fz/Dhw7MfQCSHqXBIYrVtCzfe6EeTP/989tdfO/367bffzkcffZT9ACI5Sr2qJNHWrYPDDoN27WDGDD/WI5tqampYvHgxRx99dHZXLBKzguxVJQJ+9Pitt8Ls2fDII9lffyqV+rpoLFy4MPsBRHKQCockXlUVdOvmT1utXx9OjKlTp9K5c2f1spKCoMIhiVdUBD//OSxeDOPHhxPjvPPOY//991cvKykIKhxSEHr0gDPO8OM7whgTqrmspJCocEhBMPPTj3zySXjTkFRUVGj6dSkIKhxSMI49Fi6+GG6/Hd5/P5wY48aNo3PnzqxYsSKcACI5QIVDCsqtt8LmzXDzzeGsP5VKMWfOHM4777xwAojkABUOKSgHHgiXXw5TpsD8+eHEaNasGc45JkyYwIwZM8IJIhIjDQCUgrNqFXTqBKedBk88EU6Mzz//nC5duuiKgZK3NABQJM0++/h5rJ58El54IZwYrVu3Vi8rSSwVDilIV18N7dvD9ddDWAfd6b2sdMpKkkSFQwpSq1YwfDi89ho89lh4cWqnX+/Tpw+bNm0KL5BIhNTGIQVr82Y/Fcn69b6hvLg4nDjPP/88mzZtokePHuEEEAmB2jhE6tGsmZ+K5N13Yb/9/NQk5eUwbVp243z3u9/9umhs3LgxuysXiUFiC4euOS6NsWaNLxiffurbOpYuhQEDsl88AO68806OP/54zWUleS+xhUPXHJfGGDwYtmzZdtnatX55th1++OHMmTNHvawk7yW2cIg0xrJlmS1vCs1lJUmhwiEFrawss+VNVdvLStOvSz5T4ZCCNnIklJRsu6xFC788DLXTry9cuJCXX345nCAiIWsedwIicaqq8j8HD/anp3bbzTeSn3pqeDErKipYtGgRHTp0CC+ISIh0xCEFr6oKlizxjeTz5/tuuv/v/4U3ohz4umg899xzOmUleUeFQyTNgQfCqFHwzDPw0EPhxpo3bx5nnXWWellJ3lHhEKnjyivhW9/y81mtXBlenCOPPJK+fftqLivJOyocInUUFcF99/nxHFdeGW6s6urqr+ey0ikryRcqHCL1OPxwGDbMT4AY5iSItb2sNP265JNGFw7zDggzGZFccu21/jrlV1wBn3wSXpyKigr++7//m8MOOyy8ICJZlNHsuMFsid8MMZ+s0+y40hRz5sBxx8FFF8GDD8adjUh0sjk77qtmdnwWchLJC926+asF/vrX8Ic/hB9vwoQJjBgxIvxAIk2QaeH4DvCKmb1nZnPNbJ6ZzQ0jMZFcMXgwdOnix3Z89lm4sWbMmMGwYcPUy0pyWqanqjrWt9w5tzRrGWWZTlVJNrz2mu+i278/jB8fXpyamhq6du1KmzZtmD17Ni1btgwvmMgOZO1UVVAg9gR6Bbc9c7loiGRL9+7wk5/AhAnwt7+FFyeVSjFx4kT1spKcllHhMLOrgWnAvsHtITP7cRiJieSaESOgUyfo1w++/DK8OD179qRv376MHTuWJUuWhBdIZBdl2sbRD+junBvqnBsKnAj0z35aIrmnpMQPDFy0CIYMCTfWuHHjeOaZZygvLw83kMguyLRwGLA57f7mYJlIQTjtNBg4EO64A159Nbw4e+65J2eddRYAK8Oc90RkF2RaOKYAr5nZMDMbBrwK3Jf1rERy2M9/Dh06QN++sH59uLEef/xxysvL1ctKckpGI8eB3wJ9gNXAp0Af59wdIeUmkpPatIGJE2HBAt/uEabvfOc7tGvXTnNZSU5pdOFwvt/u4865151zdznn7nTOvRFibiI5q7ISLr0URo+GN98ML456WUku0shxkV10++2wzz7+lNXGjeHFqe1lpenXJVdo5LjILtprL7j3XnjjDRg7NtxY1dXVlJeXM2fOnHADiTRCo0eOB20cpwDbDfjL5UGAGjkuYbvgAnjiCX/KqnPn8OKsW7dOI8klMlkZOR60cdzunFta95a1TEXy0N13wx57wDnnQMeO/kJQ5eUwbVp249QWjaeffhr9MyRxUhuHSBPttx/853/Cu+/CsmXgHCxdCgMGZL94fPXVVwwcOJBLL71UvawkNrvSxvGq2jhEtlXflOtr1/qZdbOpVatWTJo0iQULFqiXlcRGs+OKZEFRkT/SqMsMtmzJfrx+/foxdepUXnnlFU444YTsB5CCl80LOS3DN5D/MCgWDtivifmJ5L2yssyWN1V1dTWlpaX07t2b9WEPXxepI9PCcS9wEnBRcP9z4J6sZiSSh0aO9JMgpisp8cvDkEqluP/++7nuuuto0aJFOEFEGtA8w+d3d84da2ZvADjnPjWzSL+1ZnYQMBhIOefOjzK2SEOqqvzPwYN9wzjAjTduXR6G2kkQAZxz+B7zIuHL9Ihjo5k1w5+iwszaAY0+g2tm95vZv83srTrLK83sHTN718xu2NE6nHOLnHP9MsxbJHRVVbBkCaxZA23bwiuvRBN3+vTpnHTSSeplJZHJtHDcBfwO2NfMRgIvAbdl8PqpQGX6gqAQ3QP0BLoAF5lZFzM70syeqnPbN8N8RSKXSsH118PTT0dTPNq2bctrr72mXlYSmYx6VQGY2eHAGfjrcPzFObcgw9eXA08557oG908ChjnnKoL7gwCcc6N2sp5HGzpVZWYDgAEAZWVl31y6NGc7fUlCffEFHHQQdOsGf/5z+PHUy0qyLZu9qnDOve2cu8c5d3emRaMB7YHlafdXBMvqZWZ7m9l44JjaIlNPjhOdc8c5545r165dFlIUycwee8CgQfDcc+Feo7xWbS8rTb8uUci4cISgvha9Bg+DnHOfOOcGOuc67eyoRCROAwdCaam/zGyGB/YZS59+/Q/1jUYUyaJcKBwrgAPS7ncAPogpF5GsadUKbroJXnoJnn02/Hg9e/bkrbfe4vvf/374waSg5ULhmAkcYmYHBl17LwSeiDknkazo189PfBjFUQfAEUccAcCbb76pU1YSmp0WDjM7y8wmmdnRwf0BuxrMzB4GXgEOM7MVZtbPObcJuBL4E7AAeMQ5989djZEWq5eZTaypqWnqqkR2WYsWMHQozJwJTz4ZTcz33nuP4447Tr2sJDQ77VVlZr/DX2f8JuAPwPnOucsjyC0rNFeVxG3TJn+djpISf9GnogiO89XLSpqqqb2qPnbOrXHOXQv0ADStukgGmjeHYcNg7lx49NFoYqbPZaVTVpJtjSkcT9f+4py7AXgwvHREkunCC6FLF7j5Zti8Ofx4qVRK069LaHZaOJxzv69z/5fhpSOSTM2awfDh8PbbMH16NDErKysZOHAgrVu3jiagFIyMR45vtwKzl51z385SPlmnNg7JFVu2wDe/CZ995gtIcXH4MTX5oeyqrI4cr0dpFtaRdepVJbmmqAhGjIBFi+CBB6KJWVs0nnrqKcaOHRtNUEm8RhUOM/ulmQ0ws5PMrO5xbwS90zPnnHvSOTcglUrFnYrI1/7jP6B7d3/aKsrrLz3++OPccMMNzJgxI7qgkliNPeKYBxwFjAaWmNliM3simCFXJ1BFGskMbr0Vli+HSZOiiztu3Dj1spKsaVThCCYNvNI5d5pzbm/85WPHA5/hB+6JSCOdcQaceqq/OuDatdHETKVSTJ48mQULFjBs2LBogkpi7VIbh3NuhXPuD865nzvnLsl2UiJJZubbOlauhF/9Krq4FRUV9OvXj7FjxzJ//vzoAkvi5MJcVSIF59RToUcPGD0aPv88urjjxo3jwQcfpHPnztEFlcRR4RCJyYgRsGoV3HVXdDFTqRRVVVWYGZ9HWbEkURJbONQdV3LdCSdAr17wi1/465RH6a9//SsHHHAAM2fOjDawJEJiC4e640o+GD7cF43q6mjjHnvssbRu3Vq9rGSXJLZwiOSDo4+GCy6A22/3p62iUjuX1fz58zWXlWRMhUMkZsOGwZdfwpgx0catrKykb9++jBkzRgMDJSMqHCIx69IFqqrg7rt9F90oVVdX0759e5577rloA0teU+EQyQE33wwbNsCoUdHGTaVSzJ07lxtvvDHawJLXVDhEcsDBB0Pv3jB+vJ+OJEp77rknALNmzeKNN96INrjkJRUOkRwxZIi/zGyXLn4m3fJymDYtmtgbN27kggsu4JJLLlEvK9mpxBYOjeOQfPPSS75gfPEFOAdLl8KAAdEUj+LiYsaPH69eVtIoTb6QU67ThZwkX5SX+2JRV8eOsGRJNDlcdtllTJkyhVdffZXjjz8+mqCSk3Z0IScVDpEcUVTkjzTqMvNXD4xCTU0NXbt2pU2bNsyePZuWLVtGE1hyzo4KR/OokxGR+pWV1X/EUVYWXQ6106/PnDmTZs2aRRdY8kpi2zhE8s3IkVBSsu2y4mK/PEoVFRXcdNNNFEdxUXTJSyocIjmiqgomTvRtGmbQqpX/eeqp8eTz7LPPcuaZZ6qXlWxHhUMkh1RV+YbwLVvgn/+EZs3g6qvjycU5x1/+8hf1spLtqHCI5KgDD4ShQ+F3v4Onnoo+fu0VA8eMGaPp12Ub6lUlksM2bIBjjvGTIM6fv30bSNjUy6pw7ahXVWKPODQAUJKgRQt/XfKlS/0VA6OWPv36ww8/HH0CkpN0xCGSB/r0gYcegjffhCOOiD7+Sy+9xLe//W3MLPrgEouCPOIQSZIxY6BNG/jRj+ofJBi2k08+GTNjyZIl6mUlKhwi+aBdO188XnwRHnggnhw++OADjjzySPWyEhUOkXzRpw9861tw7bXwySfRxy8tLeUHP/iBelmJCodIvigq8tfrWLMGbrghnhzGjRtHaWkpvXv31imrAqbCIZJHjjwSrrkGJk+Gl1+OPn7tXFaafr2wqXCI5Jmbb/YTHw4cCBs3Rh+/oqKCyy67jM8//5yk98qU+ml2XJE8s/vucNdd8L3vwR13wHXXRZ/DhAkTKCrS/52FSlteJA+dey6ccw4MGwbLlkUfv7ZozJw5k7vvvjv6BCRWKhwieequu/zPq66KL4dJkyZx9dVXM2PGjPiSkMgltnBoyhFJuo4d/RHH73/vb3EYO3YspaWl9OnTR72sCkhiC4dz7knn3IBUKhV3KiKh+clPoGtXf9Tx5ZfRx1cvq8KU2MIhUgiKi/0kiMuWwfDh8eSQPv36G2+8EU8SEin1qhLJcyefDP36QXU1XHKJH+sRtXHjxtGlSxe6du0afXCJnGbHFUmATz6Bww6Dww+HF17wo8zjsmHDBlq0aBFfApIVmh1XJOH23hvGjvWjyadOjS+P119/nU6dOmkuq4RT4RBJiB/+EE45xQ8IXLUqnhw6deoEoLmsEk6FQyQhiop8Q/lnn8H118eTg3pZFQYVDpEEOeIIP+36lCnwjW/4YlJeDtOmRZdDei8rnbJKJhUOkYQ59FAwg48+8lcLXLoUBgyItnjUTr8+LcqgEhn1qhJJmPJyXyzq6tgRliyJLo/333+f0tJSXac8T6lXlUgBaWjSw6gnQ2zfvj1mxuLFi5k3b160wSVUGgAokjBlZfUfcZSVRZ/Lli1bOPvssykqKmL27Nm0bNky+iQk63TEIZIwI0dCScm2y1q18sujVlRURHV1tXpZJYwKh0jCVFXBxIm+TaO2eeE73/HL49CzZ0/69u3LmDFjNP16QiS2cdzMegG9Dj744P4LFy6MOx2R2PTpAw8/DO++Cx06xJNDTU0NXbt2pU2bNjpllScKsnFc06qLeMOG+W65cc2eC1sHBp5yyils3rw5vkQkKxJ7xFFL3XFF4Oqr4Z57YP58P85DZGcK8ohDRLYaPBhatoQhQ+LOBGbNmsX555+vuazymAqHSAHYd1+45hp45BF4/fV4c1m1ahWPPfaYelnlMYdFPJoAAAcUSURBVBUOkQLx05/CXnv5o484VVZWft3LSnNZ5ScVDpECkUrBoEHwxz/C3/8eby7V1dWUlpZq+vU8pcIhUkCuuALat/cFJM5+MalUikmTJjF//nwmTJgQXyKyS1Q4RApIq1YwdCi88go89VS8uVRWVvLEE09w+eWXx5uIZEzdcUUKzMaN/rodLVvCm2/Ge33yWqtXr6akpEQDA3OIuuOKyNeKi2HECJg3z48oj9vq1avp2rWrelnlERUOkQJ0wQVw9NF+XMeGDfHmstdee9GzZ0/NZZVHVDhEClBREdx2GyxeDJMnx53N1l5Wffr0US+rPKDCIVKgKivh1FP9aasvv4w3l1QqxcSJEzX9ep5Q4RApUGYwahSsXAm//GXc2Wydfn3BggVs2bIl7nRkB9SrSqTA9eoFL70EixZB27bx5rJ+/XpatGih65TnAPWqEpEGjRwJNTUwZkzcmcBuu+329XXKp0yZEnc60gAVDpECd9RRcPHFcOed8OGHcWfjjRkzhssuu0y9rHKUCoeIcMstfmDgrbfGnYk3evRozWWVw1Q4RIROnaB/f3+t8kWL4s5m61xWCxYsUC+rHJTYwmFmvcxsYk1NTdypiOSFIUP8qPKhQ+POxNP067krsYVD1xwXycz++/tLzE6f7qcjyQXV1dVcd911HH744XGnImnUHVdEvvbpp3DQQXDKKfDEE3Fns60tW7ZQlAszMhYIdccVkUZp2xauvx6efBL+8Y+4s9nqX//6F8ccc4x6WeUIFQ4R2cZVV8F++8V/sad0++23H6tXr9ZcVjlChUNEtrH77r6h/IUX4E9/ijsbT3NZ5RYVDhHZTv/+sM8+cM45fibd8nKYNi3enGrnstL06/FT4RCR7fz2t/DZZ35QoHOwdCkMGBB/8aidfv3OO++MN5ECp15VIrKd8nJfLOrq2BGWLIk6m20tXLiQ8vJyiouL400k4dSrSkQysmxZZsujdMghh1BcXMzq1at5++23406nIDWPOwERyT1lZfUfcZSVRZ9LfZxzVFZW8uWXXzJ79mxatmwZd0oFRUccIrKdkSOhpGTbZSUlfnkuMDNGjBjB/PnzmThxYtzpFBwdcYjIdqqq/M/Bg/3pqbIyXzRql+eCiooKnnrqKSoqKuJOpeCocVxERLajxnEREckaFQ4REcmICoeIiGREhUNERDKiwiEiIhlR4RARkYyocIiISEZUOEREJCOJHwBoZh8D9cy6k5EUUJOFdLJhH2BVnWUN5bejvHf2nhrznjP5XJryGWb789f23PXn7Mpzs/naKNbXFEnanh2dc+3qfbZzTred3ICJceeQlsusxua3o7x39p4a854z+Vya8hlm+/PX9tT21PZs2meoU1WN82TcCexEQ/ntKO+dvafGvOdMPpemfIbZ/vy1PXf9Obvy3Gy+Nor1ZVvitmfiT1UljZnNcg3MHyP5R9szWQple+qII/9oDulk0fZMloLYnjriEBGRjOiIQ0REMqLCISIiGVHhEBGRjKhw5DEz293MHjCzSWaWQxf1lGwws4PM7D4zezTuXKTpzOx7wd/q782sR9z5NIUKR44xs/vN7N9m9lad5ZVm9o6ZvWtmNwSLvw886pzrD5wTebLSoAy3Y72cc4ucc/3CzVQaI0vb8/Hgb7U38IMQ0w2dCkfumQpUpi8ws2bAPUBPoAtwkZl1AToAy4OnbY4wR9m5qTRyO5rZkWb2VJ3bvtGnLDswlextz5uC1+Wt5nEnINtyzr1gZuV1Fp8AvOucWwRgZv8DnAuswBePN9E/ATklk+3onBsF/J9oM5RMZGN7mpkBo4FnnHOvh5txuLSzyQ/t2XpkAb5gtAf+F/i/ZvYrcn/aBWl4O9bLzPY2s/HAMWY2KOzkJGMZbU/gx8CZwPlmNjDMxMKmI478YPUsc865L4E+UScju6ze7djQk51znwB5vYNJuEy3513AXeGlEx0dceSHFcABafc7AB/ElIvsOm3HZCnY7anCkR9mAoeY2YFm1gK4EHgi5pwkc9qOyVKw21OFI8eY2cPAK8BhZrbCzPo55zYBVwJ/AhYAjzjn/hlnnrJj2o7Jou25LU1yKCIiGdERh4iIZESFQ0REMqLCISIiGVHhEBGRjKhwiIhIRlQ4REQkIyocIiKSERUOERHJiAqHSAzM7Ewz+3XceYjsChUOkXh0A96IOwmRXaHCIRKPbsA3zOxFM1tpZmfGnZBIY6lwiMSjG7DKOXcKcDlQFXM+Io2mwiESMTMrBvYCfhEsag6siS8jkcyocIhErwswxzm3Jbh/FPBWjPmIZESFQyR63YA5afePAubGlItIxlQ4RKLXjW0LRVd0xCF5RBdyEhGRjOiIQ0REMqLCISIiGVHhEBGRjKhwiIhIRlQ4REQkIyocIiKSERUOERHJiAqHiIhk5P8Dnz7V3OLknG4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "plt.plot(h_values, errors, \"bo-\")\n",
    "plt.plot([0.2, 0.02], [0.5, 0.5 * (0.02 / 0.2) ** 1.5], \"k--\")\n",
    "plt.xlabel(\"$h$\")\n",
    "plt.ylabel(\"$L^2$ error\")\n",
    "plt.xscale(\"log\")\n",
    "plt.yscale(\"log\")\n",
    "plt.axis(\"equal\")\n",
    "plt.xlim(plt.xlim()[::-1])\n",
    "\n",
    "plt.title(\"Convergence of an interior problem\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What next?\n",
    "\n",
    "After reading this tutorial, you should attempt [exercise 3](../exercises/3_iterations.ipynb)."
   ]
  }
 ],
 "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}