Taylor series as the building block for flow models.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Taylor series as the building block for flow models\n",
    "*Chiel van Heerwaarden (chielvanheerwaarden@gmail.com), 2018*\n",
    "\n",
    "In this assignment, you will learn the basics of Taylor series and how they form the building blocks of numerical approximations for temporal and spatial derivatives in flow models.\n",
    "\n",
    "At the end of this tutorial, you will be able to:\n",
    "1. Understand the meaning and construction of Taylor series.\n",
    "2. Derive numerical approximations for derivatives using Taylor series.\n",
    "3. Understand the concept of errors and convergence.\n",
    "4. Test the error and convergence of numerical approximations of derivatives.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "Many flow models in meteorology and hydrology can be put in the form of an advection-diffusion equation\n",
    "\n",
    "$$\n",
    "\\dfrac{\\partial c}{\\partial t} = - u \\dfrac{\\partial c}{\\partial x} + \\nu \\dfrac{\\partial^2 c}{\\partial x^2},\n",
    "$$\n",
    "\n",
    "where $c$ is an arbitrary variable, $t$ is time, $u$ is the velocity, $x$ is the location, and $\\nu$ is the kinematic viscosity. The equation states that the change in variable $c$ over time (left-hand side) at a certain location, is the result of the net inflow / outflow due to advection (first term right-hand side) and the net change due to diffusion (second term right-hand side).\n",
    "\n",
    "[EXPLAIN THE DERIVATIVES]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculating derivatives using Taylor series\n",
    "Intuitively, the derivative $\\partial c/ \\partial x$ can be approximated using a forward approximation\n",
    "\n",
    "$$ \\left. \\dfrac{\\partial c}{\\partial x} \\right|_x \\approx \\dfrac{c\\left(x + \\Delta x\\right) - c \\left( x \\right)}\n",
    "{\\Delta x}, $$\n",
    "\n",
    "a backward approximation\n",
    "\n",
    "$$ \\left. \\dfrac{\\partial c}{\\partial x} \\right|_x \\approx \\dfrac{c\\left(x \\right) - c \\left( x - \\Delta x \\right)}\n",
    "{\\Delta x}, $$\n",
    "\n",
    "and a centered approximation\n",
    "\n",
    "$$ \\left. \\dfrac{\\partial c}{\\partial x} \\right|_x \\approx \\dfrac{c\\left(x + \\Delta x \\right) - c \\left( x - \\Delta x \\right)}\n",
    "{2 \\Delta x}. $$\n",
    "\n",
    "Formally, each of these approximations can be derived from a Taylor series around $x$. The advantage of doing this, is that it allows for a quantification of the error of the approximation. \n",
    "\n",
    "[INTRODUCE FORMAL DEFINITION TAYLOR SERIES]\n",
    "\n",
    "If we fill in $\\Delta x$ for $a$, we obtain the Taylor series below. The forward approximation, including its error, can be derived from this series\n",
    "\n",
    "$$ c \\left(x + \\Delta x \\right) \n",
    "= c \\left( x \\right) \n",
    "+ \\dfrac{\\left( \\Delta x \\right)  }{1!} \\left. \\dfrac{\\partial   c}{\\partial x  } \\right|_x\n",
    "+ \\dfrac{\\left( \\Delta x \\right)^2}{2!} \\left. \\dfrac{\\partial^2 c}{\\partial x^2} \\right|_x\n",
    "+ \\dfrac{\\left( \\Delta x \\right)^3}{3!} \\left. \\dfrac{\\partial^3 c}{\\partial x^3} \\right|_x\n",
    "+ \\dots $$\n",
    "___\n",
    "<font color=red>\n",
    "**Exercise 1**\n",
    "\n",
    "1. Rewrite the Taylor series above such that $\\partial c / \\partial x$ is on the left-hand side, and all other terms are on the right-hand side.\n",
    "2. Compare the outcome to the intuitive approximation. The difference between the two is the error of the approximation. What is the error (in symbolic math) in this case?\n",
    "3. Derive the backward approximation by filling in $- \\Delta x$ for $a$ in the Taylor series. What is the error?\n",
    "4. Construct the centered approximation by subtracting the calculated Taylor series for $a = -\\Delta x$ from that with $a = \\Delta x$. What is the error?\n",
    "</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Answer**\n",
    "\n",
    "Work out your answers here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "## Numerical implementation of the three approximations\n",
    "Now, you will implement the three approximations \n",
    "\n",
    "[WORK OUT THE TEST]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Global settings.\n",
    "h0 = 1.\n",
    "alpha = 1.\n",
    "total_time = 5."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "class h_euler_forward:\n",
    "    def __init__(self, dt):\n",
    "        self.dt = dt\n",
    "        self.time = np.arange(0, total_time, dt)\n",
    "        self.h_ref = h0*np.exp(-alpha*self.time)\n",
    "        self.h = np.zeros(self.h_ref.size)\n",
    "        self.h[0] = h0\n",
    "        \n",
    "        time = 0\n",
    "        h = h0\n",
    "        for i in range(1, self.time.size):\n",
    "            time += dt\n",
    "            dhdt = -alpha*h\n",
    "            h += dt*dhdt\n",
    "            self.h[i] = h\n",
    "        \n",
    "        self.L1_error = dt*np.sum(abs(self.h - self.h_ref))\n",
    "        self.L2_error = dt*np.sum((self.h - self.h_ref)**2)**.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x10e0e7c18>]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd4VFX+x/H3mZYyCSGBhABpdAgdQu9BkKKgP6XpCroqKAIquorr6rrq6mLbXQVWWHUFFQEBERtIlS6ETqghJBAChBDS68yc3x8TISBClJBLku/reeaZmXvP3Pme+PCZ67nlKK01QgghKheT0QUIIYQoexLuQghRCUm4CyFEJSThLoQQlZCEuxBCVEIS7kIIUQlJuAshRCUk4S6EEJWQhLsQQlRCFqO+uGbNmjoiIsKorxdCiApp+/btqVrrwGu1MyzcIyIiiImJMerrhRCiQlJKJZamnQzLCCFEJSThLoQQlZCEuxBCVEIS7kIIUQlJuAshRCV0zXBXSn2klEpRSu37lfVKKfWuUipOKbVHKdWu7MsUQgjxW5Rmz/1jYMBV1g8EGhU/xgL/uf6yhBBCXI9rhrvWeh2QdpUmQ4E52m0LUF0pVbusCrzc0nUf8qcPB9+ozQshRKVQFmPudYETJd4nFS/7BaXUWKVUjFIq5uzZs7/ry7Yd+oFlluMs+HbW7/q8EEJUBeV6QFVrPUtrHaW1jgoMvObVs1d0a5sHAIg58XVZliaEEJVKWYT7SSC0xPuQ4mU3RPcOA2hYoIgzleoKXCGEqJLKItyXAqOLz5rpDGRorU+VwXZ/VQtbE454aLbt/fFGfo0QQlRYpTkV8nNgM9BEKZWklHpQKfWIUuqR4ibfAfFAHPBfYPwNq7ZYz0bDAPj0h6k3+quEEKJCuuZdIbXWo66xXgOPlVlFpdC3813U3f83Uuw39H8QhBCiwqqQV6iazGZamMI46FlEUkqC0eUIIcRNp0KGO0CPhnfhUIrZX79qdClCCHHTqbDhPrjbGAIcLmJzZMIPIYS4XIUNd4vFSrPCAOJ8i8jIvtoFtEIIUfVU2HAH6BY+mDyTiS9/nG50KUIIcVOp0OF+Z69H8Xa5WBf/jdGlCCHETaVCh7uP3Y8muV4csWdz6nSy0eUIIcRNo0KHO0BUYE/SzSa2HFhqdClCCHHTqPDhPjz6SSxas+Hol0aXIoQQN40KH+7BNUNplm9jD0kkJiQYXY4QQtwUKny4AzT1aMFpq4lV2xYbXYoQQtwUKkW4j4yejNKao9kbjC5FCCFuCpUi3BuHt6FRoZl9jiMUFBQYXY4QQhiuUoQ7QKS1CfE2mPrui0aXIoQQhqs04X57h4cAyLYfMbgSIYQwXqUJ944t+hNWCInmY0aXIoQQhqs04Q7QxlqPQx5FfLFMznkXQlRtlSrcB7UahVMpFq5+G6fTaXQ5QghhmEoV7l3bjKCmw4UlIh+llNHlCCGEYSpVuCuTiShzHQ555bMr4bjR5QghhGEqVbgDDGx6BwUmxaz5fyY7O9vocoQQwhCVLtx7tH8AX6eLLO8jbNi40ehyhBDCEJUu3K1Wb6JUDRJ8c6lWv6XR5QghhCEqXbgDDKjfn0yziR82f2h0KUIIYYhKGe69o8Zh05r4tB/493vTjC5HCCHKncXoAm4Eb3sg7V12jvqk4Z2sjS5HCCHKXaXccwfoF9qDFKuJgDCX0aUIIUS5q7Th3rf9I5i05lT6d5w6e87ocoQQolxV2nAPCGhIS5eNFO+TRA25H4fDYXRJQghRbiptuAP0qxXFcQ9F14GdKCoqMrocIYQoN6UKd6XUAKXUIaVUnFJqyhXWhyml1iildiql9iilBpV9qb9ddOs/AuBt30Zh5Tx2LIQQV3TNcFdKmYHpwEAgEhillIq8rNlfgAVa67bASGBGWRf6e4SGdKah00SWPYG35q3g7NmzRpckhBDlojR77h2BOK11vNa6EJgHDL2sjQaqFb/2A5LLrsTr09c/kjhPJ19u+JHZs2cbXY4QQpSL0oR7XeBEifdJxctKegn4g1IqCfgOmFgm1ZWBWyJHoZWideMEbrvnQaPLEUKIclFWB1RHAR9rrUOAQcAnSqlfbFspNVYpFaOUiimvIZImjW6jtgtMvgf5cvfpcvlOIYQwWmnC/SQQWuJ9SPGykh4EFgBorTcDnkDNyzektZ6ltY7SWkcFBgb+vop/I2Uy0de/Bfu9nGza9ANTnnu+XL5XCCGMVJpw3wY0UkrVU0rZcB8wXXpZm+NAXwClVDPc4X7THL2Mbv0ghSZFpH0Z3+5MkNMihRCV3jXDXWvtACYAy4EDuM+KiVVKvayUGlLc7CngYaXUbuBz4H6t9U1zU5e24b2pjpki36N0HTkeq9VqdElCCHFDlerkb631d7gPlJZc9mKJ1/uBbmVbWtmxmCz0CmrPaucW2sV9x8m01tSu7o3JVKmv4RJCVGFVJt2iI+8hy2yilfdKOo2YxPLly40uSQghbpgqE+5d63bDS5lJ8T1D09ZtCAgIMLokIYS4YapMuHtaPOka3Ik13l4M8YyBwIZGlySEEDdMlQl3gOgGg0mxWGjtvYHPNhzmp59+MrokIYS4IapUuPcK6YUZE7vsDnIPLOPuUX/A6XQaXZYQQpS5KnWrRD8PP6KC27PK8RMTLetp9ObHmM1mo8sSQogyV6X23AH6hPXlmMVEfc+9HEhKNbocIYS4IapcuEeHRgOw1tuDyJRveODJv3DgwAGDqxJCiLJV5cK9tk9tmgU0Y5V/LUZa1rLsUDqbN282uiwhhChTVS7cAfqG9WWvKsTLksrADhEMv3e00SUJIUSZqpLhHh3mHppZ6RvAHa6VLN6RRF5ensFVCSFE2amS4d6wekNCfUP5MSiMgeZtfPj1Oho1akx+fr7RpQkhRJmokuGulKJvWF9+cpwnXzmJVluJvmc8BQUFRpcmhBBlokqGO7iHZhzaybq6kfzBuhZz02j8/PyMLksIIcpElQ33VjVbEeAZwJoatWnACdKPbObrtT+xe/duo0sTQojrVqWuUC3JbDLTJ7QP3x/7ngKrN6Nca3n0nUQ6WhJZvHix0eUJIcR1qbJ77uAemsl15LK1SW+GmrdQq1U33nv/v0aXJYQQ161Kh3un2p3wtniz2q8GHjqPfnozPxzJNLosIYS4blU63D3MHnSv2501aftwBjblj17rmLXmMAMGDiQlJcXo8oQQ4ner0uEO7qGZc/nn2NvsVho7DuOXE8/+83DkyBGjSxNCiN+tyod7z5CeWEwWVnt7oM02HvBaR6+HX6Jbt5t2vm8hhLimKh/uvjZfOgZ3ZFXyJmh6G0NN69l6JJm4M1kcO3bM6PKEEOJ3qfLhDu7bAB/POs7Rpv3xdGRxmzWGP/5jNl27dpV7zgghKiQJd6B3aG8AVrsywT+CR303csa7Hi/+/Q2sVquxxQkhxO8g4Q7UsteiVc1WrD6xBtreR8PcndTRp8gL6YjFUmWv8xJCVGAS7sX6hPUh9lwspxv3A2Xm2aBtfLwpgVn/+4QvvvjC6PKEEOI3kXAv9vM93lefj4XGt3JL/goK83N47/udfPrppwZXJ4QQv42Ee7H6fvWp51eP1SdWQ9eJWPLP8ULwVqwtbmXeF4uMLk8IIX4TCfcSokOjiTkdQ0ZwcwjvzrCCReTk5rIgJon8/HwcDofRJQohRKlIuJcQHRaNUztZl7QOej2DLS+FpwO3MX31YRo0asxnn31mdIlCCFEqEu4ltKjZgkCvQFYfXw31ekJoJ+5zLuZ8dh5t7xhLw4YNjS5RCCFKpVThrpQaoJQ6pJSKU0pN+ZU2w5VS+5VSsUqpuWVbZvkwKRPRYdFsTN5IvrMAej2DZ+4pJgRsJTe8O527dDW6RCGEKJVrhrtSygxMBwYCkcAopVTkZW0aAc8B3bTWzYEnbkCt5SI6NJo8Rx6bkzdDg75Qpx0Ps4TjqRl8uS2eadOm4XQ6jS5TCCGuqjR77h2BOK11vNa6EJgHDL2szcPAdK31eQCtdYW9X26H4A74WH3cZ80oBb2exTs3iYf9Ynhn2X4mTpzIjz/+aHSZQghxVaUJ97rAiRLvk4qXldQYaKyU2qiU2qKUGnClDSmlxiqlYpRSMWfPnv19Fd9gVrOVXqG9WJW4iszCTGh8KwS3YrzlK87kaf777Saio6ONLlMIIa6qrA6oWoBGQG9gFPBfpVT1yxtprWdpraO01lGBgYFl9NVl7/7m95NVlMXH+z527733/BO+OYmM9t3OskSN1pqCggKjyxRCiF9VmnA/CYSWeB9SvKykJGCp1rpIa30MOIw77CukpgFNGRAxgE8PfEpqXio0vQ2CInnctpTtied4eeYCGjRoQFpamtGlCiHEFZUm3LcBjZRS9ZRSNmAksPSyNktw77WjlKqJe5gmvgzrLHfj24ynwFnAh3s/BJMJej5N9Zx47vXdxfrMAHr07ElhYaHRZQohxBVdM9y11g5gArAcOAAs0FrHKqVeVkoNKW62HDinlNoPrAH+pLU+d6OKLg/1/OoxtMFQ5h+az6nsUxB5B9RszNOeSzmals99z71NcHCw0WUKIcQVlWrMXWv9nda6sda6gdb678XLXtRaLy1+rbXWk7XWkVrrllrreTey6PLySOtHAJi5ZyaYzNDjaapnHWa0/37eWn6I+IREPvjgA4OrFEKIX5IrVK+ijk8dhjcZzpK4JSRkJECLuyCgPk95fEV8ajZT3v+SSZMmkZJSYc/8FEJUUhLu1/BQy4ewmW3M2DUDzBbo8RTV0mN5qNYREnybs2tvLEFBQUaXKYQQl5Bwv4aaXjW5t9m9fJ/wPYfSDkGrEVA9jMetS0jJKmDtSRcA+fn5BlcqhBAXSbiXwv3N78fX6su0ndPAbIXuk/FN3cX40BNMX3OUJ555jh49euByuYwuVQghAAn3UvHz8OOBFg+wNmktu1J2QZt7oFpdHjMtIiOvkMzaHenXrx9FRUVGlyqEEICEe6nd2+xeAjwDeHfnu2izDbo/if3MNp5slMLmNC+eeO5FPDw8jC5TCCEACfdS87Z683DLh9l2ehtbTm2BtveBTzBjXQspcrqYtjqOmJgY3nvvPaNLFUIICfffYliTYQTbg3lv53toiwd0m4TXyY38KfI8n289zrv/+5ypU6eSk5NjdKlCiCpOwv038DB78GjrR9mbupc1J9ZA+wfAHsj9RV9gNZswt7ub2NhY7Ha70aUKIao4CfffaEiDIYRXC+e9ne/htHhAlwl4JK7h5ah81hxJY9eZQrTWJCQkGF2qEKIKk3D/jSwmCxPaTCAuPY7vE76HDg+BVwB3Zs2lXk07f1say6QnnqRz585kZGQYXa4QooqScP8d+kf0p4l/E2bsmkGR1QO6jMd8ZDlvdnEQn5qDf6f/45VXXsHX19foUoUQVZSE++9gUiYmtp3IiawTLIlbAh3HgT2IqNi/c0uTGiw6mMuQEfdhMsmfVwhhDEmf36lnSE9aB7bm/d3vk2+xwa2vQfIOptbbSZFT84/vD/LNN99wzz33oLU2ulwhRBUj4f47KaV4vN3jpOSmMP/QfGh5N9TrSY0tr/N452p8ufMkMQlpxMbGcrPOFyuEqLwk3K9Dh+AOdK7dmQ/3fkiOIxcGvQ2FuYwr+Jjafp5sdUawdVuM3DVSCFHuJNyv06S2kzhfcJ5P9n8CgY2h2+NY9i3gnY5Z7D+VxRc7kikoKOCTTz6R4RkhRLmRcL9OLQNbEh0azezY2aTnp0PPp6F6OJ0PvEb3etV464dDvP+/Txg9ejRbtmwxulwhRBUh4V4GJrSdQE5RDh/FfgRWLxj0Fir1EP8M3UBWvoOTNTrw448/0qVLF6NLFUJUERLuZaCRfyMG1x/M5wc+52zuWWjcH5reRuCOf/N0R08WbE9C1W4GwJkzZ2R4Rghxw0m4l5HxrcfjcDmYtWeWe8HAqaBMjM2ZSf2adqYs2sumrdtp0KABCxYsMLZYIUSlJ+FeRkKrhXJnoztZeGQhSVlJ4BcCvadgPrKMWR1Pk5yRx7JkG2PHjpXhGSHEDSfhXobGtRqHCRP/2f0f94LOj0JQJA23v8pDHWsxZ8txRk78C2FhYcYWKoSo9CTcy1Atey1GNR3FN/HfEJ8e755vdfA7kHGCZ7yWEBrgxbOL9pB6PoPRo0fL8IwQ4oaRcC9jD7Z8EE+zJ9N2TXMvCO8Cbf6Adet/eDfag2OpOczcmERcXBzJycnGFiuEqLQk3MuYv6c/o5uPZkXiCmLPxboX9nsZPHxpu+dVRnUI5cONCbw3dylPPPGEscUKISotCfcbYHTkaPw8/HhvZ/F8qvYacMtLcHwTL4btJsjXk+e+jKXA4WTbtm1MmzbNyHKFEJWQhPsN4Gvz5cEWD7Lx5EZiTse4F7YdDSEd8VrzV94YHMLhM9lMWx3HzJkzefvtt2XeVSFEmZJwv0FGNh1JoFegezJtrcFkgtvegbx0eh6fwV3tQpi+Jo77Jr/E9u3bZd5VIUSZknC/QbwsXoxtNZYdKTvYmLzRvTC4JXR6BLZ/zCtRuYQFePPcV4cwefrgcrlYuHChXL0qhCgTpQp3pdQApdQhpVScUmrKVdrdpZTSSqmosiux4rqr0V3U9anLuzvexaVd7oV9ngPf2ngv/xPvjWjF2ewCnl20hyVLvmLYsGF8/fXXxhYthKgUrhnuSikzMB0YCEQCo5RSkVdo5ws8DvxU1kVWVFazlfFtxnMg7QArE1e6F3r4woDX4cxeWp6cxzO3NmV57Blyarfl66+/5vbbbze2aCFEpVCaPfeOQJzWOl5rXQjMA4Zeod0rwFQgvwzrq/AG1xtMfb/6TNs1DafL6V4YeQc06g8r/8aDDTLp1TiQV749QMOoniilOHXqFOfPnze2cCFEhVaacK8LnCjxPql42QVKqXZAqNb62zKsrVIwm8xMbDuRYxnH+Cb+G/dCpeCO/4B3DUwLx/D27RFU87Qyce5Ozmfm0KVLFx566CFjCxdCVGjXfUBVKWUC3gGeKkXbsUqpGKVUTFWaV7RvWF8ia0QyY9cMCp2F7oX2mjDsY8hIoubKJ/jn8FYcScnmjZXxvP7667z00ktGliyEqOBKE+4ngdAS70OKl/3MF2gBrFVKJQCdgaVXOqiqtZ6ltY7SWkcFBgb+/qorGKUUk9pOIjknmUVHFl1cEdYJ+r0Ch76lx9l5jOtVn8+3HsevRW9atmwJwKlTpwyqWghRkZUm3LcBjZRS9ZRSNmAksPTnlVrrDK11Ta11hNY6AtgCDNFax9yQiiuornW60r5We2btmUWeI+/iis6PQrMhsPIl/tT0HK1Dq/Psoj3En81m9uzZNGzYkL179xpXuBCiQrpmuGutHcAEYDlwAFigtY5VSr2slBpyowusLH7ee0/NS+Xzg5+XXAFDp4N/BJZFD/KfoSFYzSYenhNDj+j+jB8/nkaNGhlXuBCiQlJGXTQTFRWlY2Kq3s79oysfZc/ZPSy7axm+Nt+LK07vgw/6QkgHNnf7iD/8L4Y+TQKZdV8UJpMiPz8fk8mEzWYzrnghhOGUUtu11te8lkiuUC1nE9tOJLMwk9mxsy9dEdzCfe/3hPV0OT6TFwY3Y+WBFP658jAFBQX06tVL7iIphCg1CfdyFlkjkv7h/Zmzfw7n8s5durLtvdD2Plj/FmNqHmJ4VAjvrY5j9eE0Bg0aRL9+/YwpWghR4Ui4G+Cxto9R4Cy4OB1fSYPehFotUV+O45XefrQNq85TX+xm+NgnufPOOwHIyMgo54qFEBWNhLsB6vvVZ1TTUcw/NJ9FhxddutLqBcNng3bhsfgB3h/ZAh8PC2M/ieF8TiHr168nIiKCH3/80ZjihRAVgoS7QZ6KeopudbrxypZX2Hhy46UrazRwn0GTvINam19h5n3tOZNRwGNzd9A0sjm33347zZo1M6ZwIUSFIOFuEKvJytu936Zh9YZMXjuZQ2mHLm0QOQS6TIBt/6VtxipevbMFm46eY+rqE8yePZugoCC01qSmphrTASHETU3C3UB2q53pfafja/Nl/KrxnM45fWmDW16C0E6wdBLDw3OZ1LcRC7cn8dYP7h+CJ554gm7dupGVlVXutQshbm4S7garZa/F9L7TySnK4bFVj5FdmH1xpdnqvv+M1QsWjObJnnUY1TGU6WuOMntTAsOGDWP06NH4+PgYVr8Q4uYk4X4TaBLQhHd6v0N8ejyT106myFV0cWW1OnDXB3D2EOqbJ3llSHNuaVaLl76OJdOvAc8//zxKKZKSksjNzTWuE0KIm4qE+02ia52uvNjlRTaf2swrm1+5dLq9Bn2gz59h7wIs617nvVFtaRfmzxPzdrH56Dny8/Pp2bMn999/v2H1CyFuLhajCxAX3dnoTk5mn2TmnpnU9anLuNbjLq7s8TSkH4d1b+Jl9ebDMRO4+/3NjJ0Tw4JHuvDKK6/IGTRCiAtkz/0m81ibx7i9/u1M2zWNr4+WmE/VZILb/w0th8Gqv1F9z0fM/mNH7B4Wxny0lZ4D76Rdu3YALFu2jPx8mRBLiKpMwv0mo5Tib13/Rsfgjry46UW2ntp6caXJ7J7BqeltsOxZ6sZ/wew/dqTA4WLkrC2cSMvlyJEjDB48mDfeeMO4TgghDCfhfhOymq38s88/CfcN54k1T3A0/ejFlWYr3P0RNLwFlk6iScoyPnuoE9kFDkbO2oJnjbosXbqUZ555xrgOCCEMJ+F+k6pmq8aMW2bgYfFg/MrxpOaVuFjJ4gHDP4HwbvDlOFpkrr8k4Ft16YOnpyf5+fn85S9/IScnx7iOCCEMIeF+E6vjU4dpfadxvuA8j616jNyiEqc62rzhnnlQtx0sfIAWedsuCfgTabn8+OOPTJ06Ve5DI0QVJOF+k2teozlv9nyTg2kHeWbdMzhcjosrPXzh3oUQ2ATm3UuLor2XBHxkhx4cOnSIQYMGAeByuQzqhRCivEm4VwC9QnvxXMfn+DHpR/6x9R+XngPvVR3uWwL+ETB3BC1chy8JeGv1YAB27NhB69atOXz4sDGdEEKUKwn3CmJk05Hc3/x+5h+a/8tZnOw1YfRXYA+Ez+6ihSnxQsDf/f4mDp7ORCmFp6cn3t7exnRACFGuJNwrkCfbP0n/8P68vf1tlicsv3SlbzCMWQo2X/jkDlpYTzF/XGcAhr2/mUK/MLZu3UpISAhaaw4cOGBAD4QQ5UXCvQIxKROv9XiNNoFt+PP6P7MzZeelDaqHuQPeZIE5Q2lqPcuiR7sS5OvBfR9tZXms+66T77//Pi1btqQqTlAuRFUh4V7BeJg9eDf6XWr71GbS6kkkZiZe2qBGA/cQjbMQ5gwlRKWy8JGuNK9TjUc/28EnWxK59957eeONN2jfvr0xnRBC3HAS7hWQv6c/M/rOQKF4dOWjpOWnXdogqBmMXgL5mfBBP/zT9zH3oc5ENwnihSX7+GDLKZ588kmUUqSkpDBmzBjOnTt35S8TQlRIEu4VVFi1MN6NfpeU3BQmrp5IvuOye8nUbg1/XAYWG/xvIF6Hv2Lmfe0ZHhXCu6vjeGbhHgocTmJiYli6dCmJiYlX/iIhRIUk4V6BtQlqw+s9Xmfv2b38ecOfcenLzmOvFQkPrYbabWDhA1jWv8HU/2vJpL6N+GJ7Evf89yeiekSTkJBw4aZjEvJCVA4S7hVcv/B+PBX1FCsSV/B2zNu/bOAT6D7I2noUrH0dtehBJvcOZfo97difnMmQ9zaSkOk+b3716tU0bNiQJUuWlHMvhBBlTcK9EhgdOZpRTUcxZ/8c5h6Y+8sGFg/33SRv+RvEfgn/G8TgeopFj3bFbFIMm7mZxTuS6NSpE1OmTKF///7l3wkhRJmScK8ElFI82+FZeof2Zuq2qaw5vuZKjaD7EzDyMzh7CGb1IZJ4vp7YnXZh1Zm8YDf/XJPIX1/6G97e3hQVFTF06FBWrFhR/h0SQlw3CfdKwmwyM7XHVJoFNOOZdc+wL3XflRs2HQwPLgdlgo8GEJD4PZ882In7u0bwwYZj3P+/bZzNKiAlJYX4+HjOnz9fvh0RQpQJdcl9SspRVFSUlotoyl5qXip/+O4P5Dny+GzQZ4T4hly5YdYZmH8vJG2D6L9Aj6dZEJPEC1/tw9fTyj9HtKZzRHWsVivgHo9v0KAB4eHh5dgbIcTllFLbtdZR12one+6VTE2vmszoO4MiVxHjV40noyDjyg19a8GYb6DlcFj9Kiwey/A2gSyd0B1/byv3fbiVt1bGUehwUVBQwOjRo3n00UfLtzNCiN+tVOGulBqglDqklIpTSk25wvrJSqn9Sqk9SqlVSinZvTNQ/er1+Xeff5OUlcQ9397DrpRdV25o9YT/m+Xec9+7AGbfRhOfPJZO6M49ncKY+WM8w2Zu5ky2kw0bNjB9+nQAsrOzSU5OLsceCSF+q2uGu1LKDEwHBgKRwCilVORlzXYCUVrrVsBCQCbwNFiH4A7M6jcLh8vBmGVjeHfHuxQ5i37ZUCno+ScYPgdO74OZvfA6vpbX7mzJjHvbcexsNoPeXc/udCv16tUD4LnnnqNVq1akp6eXc6+EEKVVmj33jkCc1jpea10IzAOGlmygtV6jtf55mqAtwK8M9IryFBUcxaIhi7i9/u38d+9/uee7ezhy/siVG0cOdR9o9fCBT/8Plk5iUCM73z3egybBvjw+bxePzd1BanYBEydO5NVXX6V69eoA5OXllWOvhBClUZpwrwucKPE+qXjZr3kQ+P56ihJlx8fmw6vdX+Vfff5FSm4KI74Zwcf7Psbpcv6yce3WMG4ddJ0EOz+B/3QlJO0n5o/tzNP9G7Mi9gz93vmRA7l2xo0bB0BsbCxhYWGsXr26nHsmhLiaMj2gqpT6AxAFvPkr68cqpWKUUjFnz54ty68W19A3rC+Lhyyme93uvL39bR784UGSspJ+2dDqBf1fgT8ud1/89MkdWL6bzIRuwXw7qTvhNew8Pm8XD8/ZzpnMfDw9PenVqxetWrUCwOm8wo+GEKLcXfNUSKVUF+AlrfWtxe+fA9Bav35Zu1uA94BeWuuUa32xnAppDK01Xx396sJ0fVM6TuGOhneglPpl46I895k0m6eDXygMfQ9nRC/+t/GyiFLUAAAXq0lEQVQYby4/hM1i4oXbIhnWPgSlFFprhgwZQosWLXj99dd/uT0hxHUry1MhtwGNlFL1lFI2YCSw9LIvawvMBIaUJtiFcZRS3NHwDhYPWUzzms15cdOLTFo9idS81F82tnrBrX8v3ou3wZyhmL97ioc6BrHsiZ40q12NZxbuYeSsLRw8nYnD4SAiIoLatWtf2IRMyi2EMUp1EZNSahDwL8AMfKS1/rtS6mUgRmu9VCm1EmgJnCr+yHGt9ZCrbVP23I3n0i4+O/AZ/9r+L+xWO3/t8lf6hve9cuOSe/HVQ2HodFzhPfh823HeXH6IrHwH93UO58l+jfHzcl/4tGzZMp5//nkWL14sFz8JUUZKu+cuV6gKjqYf5bn1z3Eg7QBDGgxhSscp+Np8r9z4+BZY8iikxUOHh+GWlzjvsPH2ikPM/ek4/t42nh3QlLvbh7Bs2fe89dZbLFu2DJvNRmFhITabrVz7JkRlI+EufpMiVxEzd8/kg70fEOgdyKvdXqVT7U5XblyYC6tfgS3/cc/bOnAqNB7AvuRMXloaS0zieVqHVuevt0fSLswfAIfDQfv27Rk+fDjPP/98OfZMiMpFbj8gfhOrycqEthOYM3AOnmZPHvrhIaZunfrLGZ4AbN4w4HV44DswW+HzkfDRAFo4D/DFI13454jWJKfn8X8zNjF2TgxHzmRRUFBAr169aN68OQCFhYXIGVNC3Diy5y5+Ic+Rxz+3/5PPD35OPb96vN79dZrXbH7lxs4i9znxa6dC9mloPBD6vkhO9cZ8tOEYs9bFk1Po4K52ITzRrzF1q3sBMGvWLCZPnsyuXbto2LBhOfZOiIpNhmXEdduUvIkXNr7AubxzjGs1jodaPYTVZL1y48Ic+Ol92PBvKMh0z/zU5znSrMH8Z20cszcngob7uoQzvncDUk8mMHfuXF566SWUUqxbt47IyEhq1qxZvp0UooKRcBdlIqMgg39s/QffxH9Dixot+HuPv1Pfr/6vfyA3DTa8Az/NArT7oGuPp0gu8ubfK4/wxfYTeFrN3NspjId61KdWNU+KiooIDQ2la9euLF68uNz6JkRFJOEuytQPCT/w8paXyXfk82T7JxnVdBQmdZVDNhlJsPZ12DUXbD7uWxp0fpS4DJixJo6vdidjVophUSE80qsBWacTcLlctGjRgnPnzjFp0iT++te/0rhx4/LrpBAVgIS7KHNnc8/y101/Zf3J9XSq3YlXu71KsD346h9KOeg+s+bgN2APgl7PQPv7OZ5exPvrjrIwJgmn1gxpXYfxvRvQqJYvq1at4q677mLjxo00b96cnJwcvLy8MJnk+L8QEu7ihtBas+jIIt7Y9gYWZeG5Ts9xW/3brnz7gpJObIWVL0HiRvetDDo8BO1Gc7rImw/Wx/PZT8fJK3LSs3EgD3SNoGOYD3ZvbwAmTpzIqlWr2LNnDxaL5cZ3UoibmIS7uKFOZJ7g+Y3PszNlJ80CmjGy6UgG1huIl8Xr1z+kNcSthI3/hoT1YPGCVsOg4zjSfBvz2ZZEPtmSSEpWAfVq2hnTJZy72oewatm3HDx4kClT3PPE/Otf/6Jr16507NixnHorxM1Dwl3ccE6Xk8Vxi5l7YC5x6XH42nwZ2mAow5sMp55fvat/+Ews/DQT9iwARx6Ed4dO4yhsOIBlB1L538Zj7Dyejo+Hhbvbh3BvpzAa1fIlOzub8PBwxo0bx2uvvYbWmrS0NGrUqFE+nRbCYBLuotxordmRsoP5B+ez4vgKHC4HnWt3ZkSTEfQO7Y3FdJWhlNw093nyWz+AjOPFQzYPQrsx7DpnYvamBL7Zk0yRU9M2rDojO4TSu341bCaNv78/W7dupVu3bnz77bf079+//DothEEk3IUhUvNS+fLIlyw4vIDTOacJ8g7i7sZ3c3ejuwn0Dvz1D7qccOh797nyCevB4gkth0GncZzzacyXO08yb9sJ4lKy8baZua1VbUZ0CKWGzmTGjBk8//zz+Pn5sXTpUlasWMFrr72Gr++v3B9HiApMwl0YyuFysD5pPfMPzWdj8kYsykJ0WDQjmoygQ3CHqx+APRMLW2fB7vnFQzbdoPVIdNPb2XEWFmw7wdd7ksktdBJRw5shreswpE0dGgb5MnXqVD766CMOHDiAyWRixYoVBAQE0L59+/LrvBA3kIS7uGkczzzOF4e/4Mu4L8koyKC+X32GNxnOkAZDfv3uk1A8ZPMpxHwE54+ByQoN+0KLu8iJ6Me3h7NZuiuZTUdTcWloVrsaQ1rXYVDzIMID3dtt3bo1/v7+rF27FoDExETCwsKufXaPEDcpCXdx08l35LM8YTnzD81nb+pevCxeDK4/mBFNRtA0oOmvf1BrSN4J+xZB7BLITHIP2zTqDy3+j5Tavfj2QAZf7Upm14l0ANqEVufW5sF0DvHEy5FF06ZNycvLIygoiEceeYQ333TPBCm3IRYVjYS7uKnFpsYy/9B8vj/2PfnOfNoEtmF4k+H0j+iPh9nj1z/ockHS1otBn5MCVjs0HQQt7uK4f2e+jj3Hsn2n2XsyA4CGQT7c2rwWvRr4s3/DMlq0aE5UVBTHjh2jVatWfPbZZwwZctW5ZYS4aUi4iwohoyCDpUeXMv/QfBIzE/H38OfORncyrPEwQnxDrv5hlxMSNkDsYtj/FeSdBw8/aHY7NLuNZP8O/HAkix/2n+GnY2k4XZpa1Tzo1TiQXo2DCLPlMOPfb/P000/ToEEDli1bxrPPPsvChQtp1KhR+fwBhPiNJNxFheLSLn469RMLDi1gzYk1uLSLHiE9GNFkBN3qdMNsMl99A84iiF8L+xa7b3VQkOkeow/tBA36kFm3Jz+cD2bNoVTWHTlLVr4Ds0nRNrS6O+ybBHLmYAxvvfkmS5YswcvLi/fff58lS5awZMkSPD09y+XvIMS1SLiLCut0zmkWHVnEwsMLSc1LpY69Dt3rdqdDcAeigqOo6XWN2wI7CtzTAR5d7X6c3uNe7hUA9XvjrB/NPs92rDxpYe2hsxeGb6p5WuhUvwZd6tegS4MarFs6jyVLvuT7778H4OmnnyYpKYl58+YB7vP75cCsKG8S7qLCK3IVsfr4apYeXcr2M9vJKcoBIKJaBFHBUXSo5Q77IO+gq28o+6x7r/7oKnfYZ59xLw9sCg2iyajTg/WFjdiQmMfm+HMknssFIMBuo1O9ADpEBNA+3J+ls6dz5vQppk2bBsCgQYOIiIhgxowZAOTk5GC322/I30KIn0m4i0rF4XJwMO0gMadj2HZmGzvO7CC7KBuA8GrhRNWKIio4iqhaUVe/U6XWkLL/4l594iZw5LuHcIJbQkgH0gJa81NhfVae8mLLsTROpucB4GU10zrUj/bh/kSFB/DN7GmEBdfg8ccfR2tNnTp1uPfee3nrrbcA2L59O82aNcO7+AZoQpQFCXdRqTldTg6ed4d9zOkYtqdsJ6swC4AQn5ALQzgdanWgtk/tX99QUR4c3wzH1kFSDJzcDkXuPXfsgRDSgcyabYg1NWZNVghbkgqITc7E6XL/uwkN8KJVSHWa1/bh4MYfiG7TgDtuG0h6ejr+/v68+uqrPP/88xQWFvLhhx8ycOBAIiIibvBfR1RmEu6iSnG6nBw+f5iYM+6wjzkTQ2ZhJgB1fepe2LPvENyBuj51r7IhB5w9AEnb4MQ29/O5I+51ygRBzSmq055Er0h2FoWxIc2f7cm5JJ13790rBfVq2mkcZMeSdYrOTULp16EZJw7vo0OHDixYsIBhw4Zx9OhRnn32WV544QVat26Nw+HAZDLJPevFNUm4iyrNpV0cOX/kkrBPL3Bf4FTbXtu9Z18c+CE+IVc/MJqbBid3uIM+aRucjIF890FYTBao0ZCCgKac9KhHrCOELdlBrD/rzfHzBRc2EWC3Uc/fStPafkSGBFBwJoHXpkxk4dw5tGnThq+//poRI0awefNmWrduTXx8PPv27eOWW26RYR1xCQl3IUpwaRdH04+y7fS2C4F/vuA8ALW8a9EuqB31q9cnoloE4dXCCa8Wjrf1V0LV5YJzcXBmL5zZ7x7DPxML6YkX29h8cNRsQpq9IcdM4ewqqMvGzJrEpJrJLXRdaBbk60GjWj5UU/mcPLCTR+65kxbhQSz65L9MfvIJzpw5Q1BQEPPmzWPu3LnMnTsXHx8fTp06hcPhICTkGj9MotKRcBfiKrTWxGfEXwj73Wd3czrn9CVtgryDLgn7n1/X9a2L1WT95UYLstzTCqbEXhr6eWkXv9fmS5FfOOmeoZxUwRx2BLErJ4CfMvyIz/cF3EFtNStqeZtoElKD8AA7yYd3s2nF1yyeM5PQAG+en/Is06ZNIzc3F5PJxOzZs4mNjeWNN94A4PTp03h6elK9evUb9jcUxpBwF+I3ynPkcTzzOImZiSRmJpKQmXDh9c9DOgBmZSbEN+SSwP/5Ocg76NI9aa0hO8Ud+GcPQdoxSIt3P9ITweW42NTiRZ5vGGm2uiSZanOkKJCDuX7szfLhWJE/WXgBCqUgwNOEjyqkTeMw6lb3YvOqb4nbvZXFn3xAsJ8nD435Azt27ODIEffxgrfeeou8vDxeeOEFAA4cOICPjw+hoaHl8rcVZUfCXYgylJ6ffknY//z6eOZx8p35F9p5Wbwu7On/HPoR1SII9wunmq3apRt1OiDjhDvozx+7NPjTjoGz4NLmFjs5nrVIswRxhgBOOPyJK6zOwZxqJLn8OaVrkIN7mkNPM/hZXTQODaK2nydb1izHlZvOC09PIsjXg0fuvwdvVcSaVSsBGDNmDHXq1OH1118HYNGiRQQHB9OtWzcACgoK8PC4yj1/RLmRcBeiHLi0i5TcFHfYZ1wM/YTMBE5mn8SlL46vV7NVo4ZXDWp41rjic4BnwIX3niYbZCVDxkn3XTAzTkLmSchIgsxk9+ufL8YqwWHxJscaQIYpgFT8OO2sRlKRD4n5dlK1H2e1H6n4kar9yMWTALuNmj42ziYdo7qHid5d2lPTx4Opf/sLLRvX4+XnnyHAbqNnx7YM7NeHWTNnAnD//fcTHR3N6NGjAfjiiy9o3rw5kZGRAGRnZ2O32+V4wA0g4S6EwYqcRZzIPkFihntvPyk7ibT8NM7lnbvwnFWUdcXP2q32S8K/ZPDX8Cp+WKtRw1GAd855VGay+0cg67R7GCjn7MXnEmP+l9Rn8iTL4k+mqka6tpPqtJPi8Oas05t07eN+cPE5AztWewB+dk/iD+4jJMifjq2bU93LyhuvvsSt0T24b8Rd+NhM3NqnJxPG/pG/PDMZm3IycMCtTJgwgeHDh5Ofn8+bb77J4MGDadeuHYWFhWzdupUmTZoQGBjIz5kkPwxXJuEuRAVQ4CwgLS+Nc/nnLoZ+8etzeecuvs4/d8m4f0meZs8L4V/dozp2qx271Y6P1cf92uKJj0tjdzqwFxXgU5SHPT8be0EW9rwM7LkZWPLT3XfVzEtD56WjuHIuuFAUmLzJwZtsvMjQXqQ7vcjUnmRpbzLxJlt7kYU3WdqLbNyvc51mlM2OV7UaOE0ebNq6gxYtmtO6WROcBTl88P40ht1xO326dyXr/FkmTxzPn595mtsH9iMt5RRTnnqSF5+fQv8+PTmVfJK3336bsWPHEhkZSUpKCsuXL6dfv34EBweTlZXFsWPHaNiwId7e3jidTpRSleYagjINd6XUAODfgBn4QGv9j8vWewBzgPbAOWCE1jrhatuUcBfit3G4HJzPP39J4P/8I/Dzj0J6QTq5RblkF2WTU5RDniOvVNv2NHu6fxBsPtgtduxmG3aTDR9lxq4Vdg0+LifeTgfejiJsjgI8ivLxcOTjUZiHtTAHW0E2ngU5eDoL8NAam9bFzz+fA3RZfzCTiyc52pMc7UEOnuThQY72JA8bBdjI0zbysZGHB/naRj5WCvAgt0hj9vLF6lWN3CLNkcQkGjRsSvWaQZxPz2LNug30HzCIOmERHEtI4tNPP2Xy4xNpVD+Cwwf28dmcj/nLn6dQPzyU/Xv3sHD+57zw5ymE1q3N/r27+f7br/nT5CepFViDQwdi2bhhA2PGjMFutxMfH8/Bgwe55ZZbsNlsnDlzhpSUFJo3b47JZCI3N5fCwkL8/PxuyP99lFm4K6XMwGGgH5AEbANGaa33l2gzHmiltX5EKTUSuFNrPeJq25VwF+LGc7gc5DpyySnMIaco50LoZxdlX/gRKPk6p/Bim5LtsoqycJQ4s+e38lAWbCYLHsqMhzLhgQkPFDYNHsU/AlaXC5vTic3pwOpyYnE5sLgcWH9+aBcWDRa0+1nrS15bKV52oY3GUrzM5AIXVrQ249IWXNqCU5txYsGFGae24MSCw+V+dmLBoc04sODQ7teYbWiTlQKnJievgGr+NcBkIys3j9TzmYTVa4DJbOVMahpJyadp16ETR6t1ZGvsMWL37eXhPz6AxaxITzvHfX1a0i7M/3f9LUsb7pZSbKsjEKe1ji/e8DxgKLC/RJuhwEvFrxcC05RSShs15iOEAMBislDNVu2XZ+r8DoXOQrKLssl35FPgLKDQWUiBs+CSx4VljuL3rhJtHJe2KXQWku/Mp9BZSE6Jz+c783G4HBS5inC4HBce+leGiq6fCygsflyZSWtMgFmDCY25eJn7GcwRmi3sda/30ZjD4ST7sJ/SnLc2xz+yG9/uPUWR00WRw0mnZlm/O9xLqzThXhc4UeJ9EtDp19porR1KqQygBpBaspFSaiwwFiAsLOx3liyEMILNbCPAHGDY9ztdTpza+cvg145LfgQuWX+VdS7twqmdlzz/2mun6/JlTvcyVxFOlwOXy3Hps3bidDmY3PEuQvwbGPL3Kk24lxmt9SxgFriHZcrzu4UQFZvZZMaMGZtZJjQvjdIcPj4JlLyMLaR42RXbKKUsgB/uA6tCCCEMUJpw3wY0UkrVU0rZgJHA0svaLAXGFL++G1gt4+1CCGGcaw7LFI+hTwCW4z4V8iOtdaxS6mUgRmu9FPgQ+EQpFQek4f4BEEIIYZBSjblrrb8Dvrts2YslXucDw8q2NCGEEL9X5bhkSwghxCUk3IUQohKScBdCiEpIwl0IISohw+4KqZQ6CyRes+GV1eSyq1+riKrab6i6fZd+Vy2l6Xe41jrwWhsyLNyvh1IqpjQ3zqlsqmq/oer2XfpdtZRlv2VYRgghKiEJdyGEqIQqarjPMroAg1TVfkPV7bv0u2ops35XyDF3IYQQV1dR99yFEEJcRYULd6XUAKXUIaVUnFJqitH1lAel1EdKqRSl1D6jaylPSqlQpdQapdR+pVSsUupxo2sqD0opT6XUVqXU7uJ+/83omsqTUsqslNqplPrG6FrKi1IqQSm1Vym1SylVJvOPVqhhmdLM51oZKaV6AtnAHK11C6PrKS9KqdpAba31DqWUL7AduKMK/PdWgF1rna2UsgIbgMe11lsMLq1cKKUmA1FANa31bUbXUx6UUglAlNa6zM7tr2h77hfmc9VaFwI/z+daqWmt1+G+lXKVorU+pbXeUfw6CziAe0rHSk27ZRe/tRY/Ks5e2HVQSoUAg4EPjK6loqto4X6l+Vwr/T92AUqpCKAt8JOxlZSP4qGJXUAKsEJrXSX6DfwLeAb3rNVViQZ+UEptL55r+rpVtHAXVZBSygdYBDyhtc40up7yoLV2aq3b4J7WsqNSqtIPxymlbgNStNbbja7FAN211u2AgcBjxUOx16WihXtp5nMVlUjxmPMi4DOt9WKj6ylvWut0YA0wwOhaykE3YEjx+PM8IFop9amxJZUPrfXJ4ucU4EvcQ9DXpaKFe2nmcxWVRPGBxQ+BA1rrd4yup7wopQKVUtWLX3vhPoHgoLFV3Xha6+e01iFa6wjc/7ZXa63/YHBZN5xSyl58wgBKKTvQH7juM+MqVLhrrR3Az/O5HgAWaK1jja3qxlNKfQ5sBpoopZKUUg8aXVM56Qbch3sPblfxY5DRRZWD2sAapdQe3Ds0K7TWVea0wCqoFrBBKbUb2Ap8q7Vedr0brVCnQgohhCidCrXnLoQQonQk3IUQohKScBdCiEpIwl0IISohCXchhKiEJNyFEKISknAXQohKSMJdCCEqof8HsYT42l1Hv1YAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "ef1 = h_euler_forward(0.1)\n",
    "plt.plot(ef1.time, ef1.h_ref, 'k:')\n",
    "plt.plot(ef1.time, ef1.h, 'C0-')\n",
    "\n",
    "ef2 = h_euler_forward(0.2)\n",
    "plt.plot(ef2.time, ef2.h, 'C1-')\n",
    "\n",
    "ef4 = h_euler_forward(0.4)\n",
    "plt.plot(ef4.time, ef4.h, 'C2-')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "dt = np.array([ef1.dt, ef2.dt, ef4.dt])\n",
    "L1_error = np.array([ef1.L1_error, ef2.L1_error, ef4.L1_error])\n",
    "L2_error = np.array([ef1.L2_error, ef2.L2_error, ef4.L2_error])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEACAYAAABRQBpkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xl4lOXV+PHvnUw2yB5AAiFkAcMSwk4CsiggWEBAFhUC+hNFtO5WrKK+data7WuLb6WtWrWaICCLgoi4YMUq+45sQsKSECRkhewzc//+eCZDQNAkhDwzk/O5Li7NZGZyZgjPmec5576P0lojhBCiafMyOwAhhBDmk2QghBBCkoEQQghJBkIIIZBkIIQQAkkGQgghkGQghBACSQZCCCGQZCCEEAJJBkIIIQCL2QH8khYtWuiYmBizwxBCCLexZcuWU1rrlnV9nEsng5iYGDZv3mx2GEII4TaUUkfq8zi5TCSEEEKSgRBCCEkGQgghkGQghBACFy8gCyFEU/LRtmxeWb2f44VltAkNYPbIBMb3bNsoP1uSgRBCuICPtmXz+NJdlFXZAMguLOPxpbsAGiUhyGUiIYRwAa+s3u9MBNXKqmy8snp/o/x8OTMQQggTHcsvZd2hPLILyy74/eMXub2hSTIQQohGlF1YxrpDeaw7lMf6jLNJwEuBXf/8/m1CAxolLpdMBkqp64HrO3ToYHYoQghxSXKKylifYRz812XkcSzfOPiHNvMhJTaCOwfHkRIXwZ7jRcxZtvucS0UBPt7MHpnQKHG6ZDLQWq8AVvTp02em2bEIIURdnCwuZ11GnjMBHM4rBSDY30JyXAS3DYilf3wECVcE4eWlnI9LaB2EUkq6iYQQwh3lnq4wDvyOBJCRWwJAkL+F5NhwpqW0JyUugs6RwXjXOPhfyPiebRvt4H8+SQZCCFEHeWcq2JCZ77zsc/DkGQAC/Sz0jQnj5r7tSImLoGubkF89+LsSSQZCCPELCkoq2ZCZx/oMIwHs/+k0AM18vekTE87EXlH0j48gsU0wFm/37daXZCCEEDUUlVaxIbP6sk8++04UozX4+3jRNyacsT3akBIXQVJUCD5ufPA/nyQDIUSTVlxexaYal3325BgHfz+LF73bh/Hw8CvpHx9BUlQovhbPOfifT5KBEKJJOVNhZVNmvrPouzu7CLsGX4sXvaJDeWBYR/rHRdAjOhQ/i7fZ4TYaSQZCCI9WUmFl85EC5yKvXdlF2OwaH29Fz3Zh3Du0Iylx4fSKDsPfp+kc/M8nyUAI4VHKKm1sOVLAuoxTrM/IZ8exQqx2jcVL0b1dKHcPiSclLoLe7cMI8G26B//zSTIQQri18iobW48UOC/7bD9WSJVN4+2l6NY2hJmD4+jvOPg395ND3sXIOyOEcCsVVhvbjhY6V/huO1ZIpdWOl4JubUOYcVUsKfER9I0JJ1AO/rUm75QQwqVVWu3syCp0XvPfcqSACqsdpaBrm2Bu7W+s8O0bG06wv4/Z4botSQZCCJdSZbOzM6vI+cl/85F8yqvsAHSODCY1uT394yPoFxtOSIAc/BuKJAMhhKmsNju7souMFb4ZeWw+nE9ppbFzZ6fWQdzcN5qUuAiSY8MJa+5rcrSeS5KBEKJR2eyaH44XOS/7bDpcwJkKKwAdWwUyqXcU/eOMT/4RgX4mR9t0SDIQQlxWNrtmb06x87LPxsx8TjsO/vEtmzOuRxv6x0eQHBtByyA5+JtFkoEQokHZ7Zp9J047Wz03ZuZTVFYFQGyL5ozp3oaUuHD6x0XQKtjf5GhFNUkGQohLorXmwE9nWHfIWOS1ITOPglLj4B8d3ozrurY2PvnHhRMZ0jgjHEXdSTIQQtSJ1ppDuWecG7ttyMgnr6QSgLahAQzrfAX94yJIiY+gbSPN7xWXTpKBEOIXaa3JOFXivOa/PiOfU2cqAGgT4s+QhJakxEXQPy6CduHNTI5W1JckAyHEObTWHMkrdY5xXJ+Rx0/FxsH/imA/BnaIoH98BClxEUSHN0Mp95nmJS5OkoEQgmP5pc7LPusz8sgpKgegRaAf/eONT/0pceHEtmguB38PJclAiCYou7DMOPg7ev2zC8sAiGjuS4rjU3//uAjiW8rBv6mQZCBEE5BTVOa85r8uI49j+cbBP6yZDylxEdw5OI7+8RF0bBUoB/8mSpKBEB7oZHG585LPukN5HM4rBSAkwIfk2HBjZ8+4CBKuCMLLSw7+QpKBEB4h93SFc5HX+ow8MnJLAAjyt5AcG860FGNzt86tg+Xg76LsdjtHjx4lJibGlJ8vyUAIN5R3poINNYa4Hzx5BoBAPwv9YsOZ4tjcrUubYLzl4O+SKisr2b59O3379kUpxZw5c5g7dy5FRUX4+jb+hnySDIRwAwUllWzINHr81x3KY/9PpwFo5utN35hw5+ZuXdsEY/H2MjlacSGlpaV8//33JCcnExQUxDvvvMNdd93FwYMHiY+PZ/LkyXTq1Am73W5KfJIMhHBBRaVVbMisvuyTz74TxWgNAT7e9IkJY6xjc7dubUPwkYO/SyotLeXzzz+nZ8+etG/fng0bNnDttdeycuVKRo0axejRo/nwww9p1aoVAL1796Z3796mxau01qb98F/Tp08fvXnzZrPDEOKyKy6vYlONyz57coyDv5/Fiz4xYY4+/wiSokLxtcjB3xVVVlayaNEiOnbsSHJyMtnZ2URFRfHqq6/y0EMPUVZWxjfffMPAgQMJDAy8bHEopbZorfvU+XGSDIRofGcqrGzKzHcWfXdnF2HX4Gvxold0KP3jWtA/PoLu7ULws3ibHa64iLfeeouwsDAmTpyIzWYjPDycW2+9lddeew2ADRs20L17d/z9G293VkkGQriwkgorm48UOBd57couwmbX+Hp70SM61LnIq2d0KP4+cvB3Vf/6178oKCjgkUceAYxLO+3bt2fp0qUAZGZmEh0djbe3eX+HLp8MlFJxwBNAiNZ6Um0eI8lAuKuyShtbjhSwLsPY1nnHsUKsdo3FS9GjXahzb59e0WEE+MrB31WlpaXx7bff8s9//hOAqVOnkpWVxdq1awEoKCggNDTUpRbq1TcZ1KqArJR6GxgDnNRaJ9a4/TpgLuANvKW1fuliz6G1zgBuV0otrmuQQriKj7Zl88rq/RwvLKNNaACzRyYwvmdbyqtsbD1S4Lzss/1YIVU2jbeXIikqxLnCt3f7MJr5St+Gq1q2bBlz587lyy+/xGKxcPjwYbZs2UJVVRU+Pj78+9//xsfHx3n/sLAwE6NtWLU6M1BKDQbOAO9VJwOllDdwALgWyAI2AVMwEsOL5z3FDK31ScfjFsuZgXBHH23L5vGluyirsjlvs3gp2oc341hhGZVWO14KurUNIcWxuVufmHAC/eTg72q01iilWLNmDffddx+fffYZ7dq1Y8mSJcydO5eFCxcSGRnpvJ87uaxnBlrrtUqpmPNu7gccdHziRym1ABintX4R4yxCCI9RXmXj+ZV7zkkEAFa75mhBKbddFUtKXDh9Y8IJ8ve5yLMIM9jtdqxWK76+vuzcuZMxY8bw9ttvM3z4cFq0aEH79u05fdpYtzFx4kQmTpzofKy7JYJLcSkfWdoCx2p8nQUkX+zOSqkI4I9AT6XU446kcaH73QncCRAdHX0J4QlRfyUVVrYcKWBDpjHJa0eWcdnnQqw2zZxRnRs5QnExVVVVlJSUEBoayqlTp0hISODZZ5/lnnvuISYmhpSUFIKCggBISkri008/NTli19Bo569a6zzgrlrc7w3gDTAuE13uuIQAKCqrYvPhfDZkGn92O7p9vL0U3dqGMGNgLIs3ZznHO9bURkY7mqqkpIRTp07Rvn177HY7UVFR3HTTTbz22mtERERwyy230LVrVwCCg4NZtGiRyRG7pktJBtlAuxpfRzluE8Ll5ZdUsjEz3/nJf69jha+vtxc92oXy26vj6RcbTq/oMJo7rvl3bh38s5pBgI83s0cmmPUymqS8vDyOHDlCr169ABg0aBAtW7Zk9erVeHl58dRTT9GpUyfAuMzzl7/8xcxw3calJINNQEelVCxGErgZmNogUQnRwE4Wlzs+9eexMTOfAz8ZG7v5+3jRKzqMB4ddSXJcOD3aXbzPf3zPtgAX7CYSl8/x48fZvn07o0aNAuCuu+5i06ZNHD58GICnn36a5s2bO+9/7733mhGm26ttN9EHwNVAC+An4A9a638ppUYBf8XoIHpba/3HhgxOuolEfWUXlrHR8al/Q2Y+maeMLZ2b+3rTJyacfrHhpMSF062tbO/gao4dO8Znn33G7bffjpeXF0888QQvv/wyhYWFNG/enE2bNlFeXs7AgQObVIG3tlx+0VldKKWuB67v0KHDzB9//NHscISL01pzNL+UDRn5rHd88s8qMCZ5BfsbWzonx0aQHBdOl0jZ1dPVZGdns2DBAqZPn06rVq147733uPXWW9m5cyfdunXj8OHDFBQUkJSUZOrKXnfhUcmgmpwZiAvRWnMo94xx2Scjn42Z+ZwoNga4RzT3pV9suDMBdGotk7xczalTp/jb3/7G+PHj6dGjBxs2bCAlJYXly5dz/fXXU1hYyKlTp4iPj5dP/vVwWdcZCGEmu12z/6fTbMjIY0OmcfCv7uppFeRHclwEyY7LPvEtZYavqykvL+fJJ59k8ODBjB07Fi8vL55//nkiIyPp0aMHvXv35vjx40RGRgIQGhpKaGioyVE3PZIMhMux2uzsySlmY2Y+6zPy2XQ4n6KyKgDahgYwJKElyY5P/u0jmsnB3wU9+OCDtG/fnoceegg/Pz8+/PBDQkJCGDt2LOHh4RQWFjq3cbZYLM5EIMwjyUCYrtJqZ1d2kbPTZ/PhAs5UWAGIbdGc67q2JjnOuPQTFdbM5GjFhTz++OPk5+c7N3Q7dOgQXl5GbUYpRUZGxjnX+y/nfv6ifiQZiEZXXmVj+7FCZ5//1iOFzt79jq0CGd+zDf1ijUs/VwQ33j7w4pfZ7XbnAf7FF19kxYoVfP/994BRx6k5rnHFihXnPFYKv67PJZNBjW4is0MRDaC00srWI4XGAq/MfLYfK6TSakcpYyHXTX3bOff1iQj0Mztc4VBaWkpAQABKKd544w0ef/xxsrOz8ff3JzIykk6dOjl383zppYtuWCzchHQTiQZ3uryKzYcLnIu8dmUVYXVs7ZDYJpjkuAj6xRgH/5Bmsqmbq8jPzycgIICAgAA+/vhjJk2axK5du+jUqRNr167lww8/5OmnnyYiIsLsUMUvkG4iYZqCkko2Ofb12ZiZzw/HjRGOPt6KpKhQ7hwcR3KcsZe/bOfsOo4dO4a/vz8tW7Zk8+bN9O3bl48++ohx48bRvXt3HnnkEZo1M2o0gwcPZvDgwSZHLC4nOTMQdZZ7uoKNmfnGCt/MfPadMLb/9bN40TM61FjgFRtOT5ni5TK01uzduxdvb28SEhLIzc2lVatWvPzyy8yePZvy8nJeffVVbrzxRuTyrIl2LoKvnoWiLAiJgmH/A0k31ukpZNGZuGxyisqcbZ4bM/M4lGts7dDM15ve7cOMNs+4CJKiZHi7K9myZQulpaUMGjQIu91OixYtmDBhAm+99RYA77zzDgMHDqRjx44mRyoAIxGsuB+qys7e5hMA179Wp4QgyUA0CK01WQVlrM/Ic3T75HM0vxSAID8LfWPDSXas8E1sG4KPbO3gMtavX8+RI0e46aabAEhJScFisfDf//4XgC+++IL4+Hji4uLMDFNciNbw5yuh5OTPvxfSDh7aXeunkmQg6kVrTcapEuPA70gAx4uMrR1Cm/nQLybcucK3c2Qw3rK1g8tYv349a9eu5dFHHwXg1ltvZfXq1eTk5KCUYseOHbRs2ZI2bdqYHKm4qNMnYMcHsC0N8g5e5E4Kni6s9VN6VDKQjeouH7td8+PJM842z42Z+eSergCgRaAfyXHhztW9HVsFyr4+LmTLli28/fbb/OUvf8HX15cXXniBZ599lpMnTxIcHExWVhaBgYGylYOrs1bCj6uNBPDjF6BtED0AcvdBWf7P7y9nBnJm0BBsds3enGLHpm55bDqcT0GpsbVDZIi/83p/cmw4sS2ay9YOLmTPnj28+OKLPPfcc8TExLB48WJmzJjBhg0b6Ny5M8XFxfj6+uLvLwvz3MLJvUYC2LEASk9BYGvoMRV6pEKLDqbXDKTPz8NU2ezszi5yfurfdDif0+XG1g7R4c0Y3vkKx17+EUSFBcjB34VkZWVx3333cf/993PNNddgs9n4/PPPmTFjBjExMYwbN478/HwsFuOfbXBwsMkRi19VXgS7lxhJIHsLePlAwm+g5zSIHwbeNQ7B1Qf8S+wmqi9JBm6uwmpjx7EiZ5vnliMFlFYaWzvEt2zOmKQ2pDj29YkMkVm9rkBrjVKKsrIyxo4dy6RJk5g1axZhYWHs3buXvLw8ABITEzlx4oQzYfv4yAI9t2C3w5H/Gglgz3KwlkHLzjDyBUi6CZq3uPhjk25stIP/+SQZuJmyShvbjhaw3tHnv+1oIRVWY0+YTq2DmNw7in6xEfSLDadlkGzt4Aqqt2wAuO6660hISGDu3LkEBATg5+fn/KTfvHlz9u3b53ycnLW5mcJjZ4vBhUfALxh6TDHOAtr0Ahf/+5Rk4OLOVFjZfDjf2ea5M6uQKpvGS0HXNiFMS2nvbPUMbeZrdrgCKCwsdBZxb7rpJvLy8vjyyy8B6NatG9HR0c77fvLJJ6bEKBqItQL2rYRt78OhrwENsYNh6JPQaQz4us8uu5IMXExRaZVjawejzXP38WJsdo3FS9EtKoTbB8aRHBtO75gwgv3lsoHZtNYcPnyY2NhYwBjGvmzZMrKyslBKMXToUEpKSpz3f+WVV8wKVTSknJ3GGcCuRVBWAMFRMORRoyAcFmN2dPUiycBkeWcqnJ/6ja0ditEafL296BEdym+vjic5NoJe7UNp5it/XWazWq1s27aNnj17YrFYeOmll3jyyScpKCggODiYsWPH0rFjR2w2GxaLhVmzZpkdsmgopfmwa7FxFnBiJ3j7QecxxmWg2CHg5d6r713y6OLJW1j/VFzubPPcmJnPjyfPAODv40Xv9mE8NPxK+sWG06NdKP4+7v3L5QnKy8tZt24dSUlJREREsHjxYqZMmcLmzZvp3bs3Y8eOpVWrVs59/keMGMGIESNMjlo0GLsNMv5jnAXs+wRslRDZHUb9GRInQrNwsyNsMLLO4DLLKih1Dm3fkJnH4Txja4dAPwt9YsKcg9u7tQ3B1yJbO5itvLycL774goSEBK688kq2bt1K7969SU9PZ+rUqeTm5rJmzRpGjhwpi7s8WX4mbJ9v/CnOgoAwoxOoRypEJpkd3S+SRWcuQGvNkbxSY3VvhnHZJ7vQWEASEuBD35hwZ5tnl8hgLLKvj+lsNhuLFi2iTZs2DBkyhOLiYsLCwnjyySd55plnsNlsfPbZZwwcOJCQkBCzwxWXU2Up7F1hXAY6/C2gIH4o9JoOCaPA4h7deZIMTKC15uDJM442T+PSz0nH1g4RzX2Nub2OvX0SrgiSrR1cxHvvvYeXlxfTpk1Da01UVBTDhg3jvffeA2Dz5s107dqVgABZl+HxtIbsrUYC2L0EKoqNAnCPaUZbaEiU2RHWmaxAbgR2u2bfidPOT/4bD+eTX1IJwBXBfqTERTj39olvGSh94i4iPT2djIwMnnrqKQDeffddlFJMmzYNpRTfffcdUVFn/9H36VPnf0fC3ZScMraF2JYGuXvBEgBdxhnF4PZXgVfTO2uXM4NfYLXZ+eF4sfN6/8bMfIodWztEhQU4h7gkx4UTHd5MDv4uYsmSJSxbtoy0tDQA7r77btavX8/WrVtRSlFYWEhISIj8fTU1Nisc/BK2p8H+VWC3Qts+RgJInAD+nnEZUM4MGkCl1c6u7ELWO673bzmcT4lja4fYFs0Z1S3SuPQTG0HbULmE4CpWr17Nc889x+rVq2nevDnHjx9n586dnD59mqCgIF577bVztnKQwm8Tc+qgkQC2fwBnTkCzFpB8l5EEWnU2OzqX0aSTQXmVjW1HC52f/LceLaC8ytja4corArmhV1vnp/9WwbIzpNmq9/TZsGEDM2fOZP78+SQmJuLt7Y3Wmp9++om4uDjuvfde7rvvPufjZE+fJqjiDOz5yLgMdHQdKG/oOMJIAFeOBG/5nTifxyWDj7Zl88rq/RwvLKNNaACzRyYwvmdbAEoqrGw9WuDo9Mljx7EiKm12lILOrYOZ0i+a5NgI+saEERHoHp0DnkprTWVlJX5+fmRkZDBs2DBeeeUVJk2axBVXXEFkZCQVFUaxfvjw4QwfPtz5WLn800RpDcc2OIrBy6CqBCI6wvBnoPvNENTa7Ahdmkclg4+2ZfP40l2UVRmXdrILy3h08U6W78gmv6SK3dlFWO0aby9FYtsQ/t9VMSTHhtMnJpyQAPmkYCar1UpxcTHh4eGUlZURExPDAw88wJw5c4iKiqJfv360aGHs9hgTE8Pq1atNjli4jPOnhfkGGjWAntOhXT+X3yDOVbhkMqjvCuRXVu93JoJqlTY7a/bl0qd9GLOGxDm2dggj0M8lX3qTUVZWxokTJ5x7+nTq1IkBAwbw3nvvERAQwIwZM+jduzcAvr6+LFy40Mxwhau54LSw/jDwIegyHvwCzY7Q7XhUN1HsYyu50KtRQOZLoxssLlF3+fn5HDx4kH79+gFwzTXXUFpayoYNGwB46623aNOmDaNGjTIzTOHqLjgtbIqxLqCF521fUx/STQS0CQ1wrvg9/3bRuE6cOMGGDRsYN24cAI8++ihLly7l1KlTeHl58dhjj1Hzg8gdd9xhVqjC1f1sWpjFMS1s+s+nhYl686h3cfbIhHNqBgABPt7MHplgYlRNQ05ODitXrmTatGn4+/vz3nvv8fvf/56cnBxat27Nfffdx/Tp050JYOTIkSZHLFzapUwLE/XiUcmgumvoYt1EouGcPHmS+fPnM2HCBKKjo9m4cSMzZ84kISGBQYMGkZqaytChQ51F3+7du5scsXAL1dPCtqdDwWG3mxbmzjyqZiAun+LiYubOncvw4cPp378/+/bto3PnzqSlpZGamsqZM2fIzs7myiuvlNZOUTcXmxbWc7rbTQtzBVIzEA3KbrczZ84cunfvzpQpU/D19eWFF17A39+f/v37k5CQQFZWFm3bGmddgYGBJCTI5ThRBx44LcydSTIQTk888QT+/v489dRTeHl5sXLlSux2O1OmTMHf35/c3FwCA42WPaWUMxEIUWs/mxbmC52v95hpYe5MkkET9sc//pE9e/aQnp4OwOHDh2nevLnz+9u3b8fb++w/zupEIESdNKFpYe5MkoGHs9lszgP6a6+9xptvvsnOnTtRSjmv7Vfv+VOdFKrVTARC1Nn508L8Q6H3bcZZgItPC2uKJBl4mLKyMnx9ffH29mb+/PncddddZGZmEhERQZs2bejVqxclJSUEBgYyZ84cs8MVnuZi08JGPGdMC/ORDR9dlSQDN1dYWIjFYiEwMJA1a9Zw3XXXsXbtWlJSUkhISGD69OnODd0mTZrEpEmTTI5YeJyLTQu75km3nRbWFLlkMqjv3kRNQXZ2Nl5eXkRGRnLgwAE6derEO++8w6233kpiYiIPPfSQs7e/d+/ezv19hGhwMi3Mo8g6Axe3f/9+rFYrXbt2pbS0lJCQEGbPns0LL7yA3W7nT3/6E+PHj6dzZxnSIRpBE5kWZqYff/yR6Oho/Pzqt41+fdcZSDJwMTt27CAvL4+hQ4cCxnbNvXv3ZsmSJYAxz7dXr15y8BeN60LTwrrfLNPCGsjJkydZuHAhaWlpbNy4kWXLljF+/Ph6PZcsOnNTW7ZsYc+ePUyfPh0wNnQ7fvw4u3btAuCdd945Z1h7amqqKXGKJuiXpoV1HAEWX7Mj9Ai//e1veeONN7DZbHTv3p2XX36Z5OTkRo9Dzgwa2datW1m5ciVPPvkkSinuv/9+3nnnHQoKCrBYLPzwww+EhISckwCEaDQXnBbWwdgaQqaFXTKr1cqXX37J4sWLef311/Hz8+P1118nKyuL1NRUEhMTL/lnyGUiF7Vr1y7+8Y9/8OKLLxIcHMy8efN48MEHOXLkCJGRkeTk5ODv709YWJjZoYqm7ELTwrreINPCGoDWmk2bNpGens6CBQs4efIkoaGhfP311/To0aPBf54kAxeRkZHBM888w+zZs0lMTOSrr75i3LhxfP311/Tt25czZ87g5eVFs2ay+ZYw2cWmhfWcJtPCGoDVasVisbBp0yb69euHn58fY8aMYdq0afzmN7+pd4H410jNwCR5eXnMmjWLW265hbFjx+Lj48Onn37KpEmTSExMZMiQIc61ACBbOggXcKFpYVfdL9PCGkBubq6zENy7d29ef/11+vTpw/vvv8+YMWMIDQ01O8SLkmRQS9VbNthsNsaOHcs111zDI488QkhICAcOHKCgoACAdu3acfLkSedWD9VJQAhTybSwy+qjjz7izTffZPXq1c5CcFKSseWGUopp06aZHOGvk9+Ai6iqqsLHxwcwVu4GBgby7rvv4u3tTWBgoPMUz2KxsHPnznMeK/v5C5cg08IuG6vVytq1a7nmmmtQSrFq1Sp2797N7NmzG6wQ3NgkGTgUFBQ4i7gzZ85k27ZtVNcrEhMTCQg4O0d54cKFpsQoRK0UZRnrAbanybSwBnShQvC6detISUnhz3/+M3//+9/xcuNV100yGWitycjIID4+HjD28f+///s/CgoK8Pb2ZvDgwcTExDgvDT399NPmBizEr3FOC0uDQ2twTgu75gmZFtYA9u3bx9ixY/nxxx/PKQT37NkTgKCgIJMjvHRNIhlYrVZ27NhBly5dCAgIYN68edx7770cO3aMqKgorrvuOlq0aEFlZSUBAQHOBWBCuDyZFnZZVBeCg4KCuPXWW4mJiSEhIYHf//73TJw40aULwfXlka2lFRUVrF+/niuvvJLIyEhWrVrFqFGj+Oqrrxg6dCiHDh1izZo1TJ482SP/UoWHu9C0sE5jjMtAcVfLtLB6Kikp4eOPPyY9Pd1ZCJ4wYYJzKxh3IesMajh06BAdOnRg3rx53H2Fn4zFAAAYNklEQVT33RQXF7Ny5UpGjBhBRETEZYhUiMvsQtPCWicZ3UDdJsm0sHqy2+3O6/yTJ09m8eLFtGvXjqlTp5Kamkq3bt1MjrDuPCoZ1NjCeuaPP/5Y58drrVm1ahX9+/eXlb3CvV1oWljSTdAz1RgdKepMa83mzZtJS0tj0aJFrF+/nvbt27Np0yZKS0sZNGiQWxeCPWrRmdZ6BbCiT58+M+vzeKUUo0aNauCohGgkMi3sssjLy+P1118nLS3tnEJwWVkZAH379jU5QnO5ZDIQosnRGo5vha0yLawh5ebmkpubS5cuXbDZbDz//PMMHDiQxx57jAkTJkjNsAZJBkKYqeQU7Fxo1AJO7pFpYQ2gtLT0nELwwIED+frrr2nVqhU5OTlSN7wISQZCNDabFQ59ZVwGqjktbMxfZVrYJXr++ef505/+xJkzZ2jXrh2/+93vzpkBIong4iQZCNFYLjQtLPkumRZWT9WF4Pnz5/PMM88QHBxMZGQkU6ZMITU11e0LwY3NJbuJqrnjFtZCnONn08K8HNPCpsu0sHo6dOgQ6enpzkKwr68vq1atco6Kbeo8qptICLd2sWlhw5+GpJshONLsCN1O9dYwGRkZdOhgbLN99dVXe/SK4MYmyUCIhnL+tDCf5pBYPS0sWTaIq6OaheCIiAj+/e9/ExcXx1tvvcWIESNo166d2SF6FEkGQtTFzkXw1bPGzqAhUXDNHPAL+vm0sIEPybSwevr222958803Wbp0KSUlJURFRXHHHXc4v3/77bebGJ3nkmQgRG3tXAQr7ocqY5ESRcfgo7uN/5dpYfWmtWbLli306NEDi8XCypUrWb58ubMQPHjwYCkENwIpIAtRW692geLsn9/erAX8br9MC6ujjIwMZyH4wIEDfPbZZ4wcOZLCwkL8/f3x95eV1vUhBWQhLofqofE7Flw4EQCU5kkiqIOcnBwmTpzIunXrAKMQPHv2bJKTkwGkGGwS+Q0W4nzVW0Ns/wB2LzbmBDRvBb5BUHn65/eXrSJ+UUlJCcuXL6esrIwZM2bQqlUrAgMD+dOf/sSUKVOkEOwiJBkIUa0oy9gaYscCOHUAvP2g02joPsXYKO6HpefWDAB8AmDY/5gXs4uyWq189dVXpKenOwvB/fr1Y8aMGXh7e/P555+bHaI4jyQD0bRVnDF2CN3xAWSuBbTRDXT9XKMbKKDGJYukG43/1uwmGvY/Z29v4qrXAgA88MADzJs3j5CQkHMKwcJ1SQFZND12Oxxea5wB7FluLAoLbW+cAXS/CcLjzI7QrVQXgtPT01m0aBFJSUls27aNzMxMRo0aJYXgRiYFZCF+Te4B4wxg5yJjUIxfMHSbCN2nQnSKLAqrg5KSEt59913S09OdheAhQ4Y4ZwP07NnTOSxeuAdJBsKzleYb8wF2fADZW4y9geKHwYhnHYNiAsyO0G2UlpZy7NgxEhIS0Frz6KOPEhcXx0svvcSUKVOIjo42O0RxCSQZCM9jrYQfPzcSwIHVYK+CKxJhxB+h22QIusLsCN2G1WplzZo1zkJwhw4d2LZtG4GBgezfv5+oKOmk8hSSDIRnqG4H3bEAdi2GsnyjHTR5FnS/GVq732Bzs73xxhv84Q9/4MSJE4SEhHDzzTeTmprqLBRLIvAsLpkMlFLXA9dX704oxEUVZddoB93vaAcdZdQB4ofKYrA6qC4E33HHHURGRhIYGMiAAQNITU2VQnATIN1Ewv1UlpxtB834BtDQLsWYFXx+O6j4Rbm5uSxatOicQvAHH3zAzTffbHJkor6km0h4NrsdDn/raAf9+Gw76JDfSztoPRUUFBAVFUVlZSWJiYlSCG7iJBkI13bqR+MMYMfC89pBpxhnA7KbZa3ULASXl5ezcOFCwsLCmDt3LgMGDCApKcnsEIXJJBkI13OxdtBrnzG2h5B20FrbtWsXb7/9NgsWLDinEFxdBL7rrrvMDlG4CEkGwjVYK+HgF0YC2P+Z0Q7aqiuMeN7RDtra7AjdRkZGBpGRkQQEBPDJJ58wb948Ro8eTWpqKqNHj5ZCsLggKSAL82gNx7cZdYDdi42toJu3hG43nm0HlVXBtXLq1ClnIfj77793FoELCwvRWhMWFmZ2iKKRSAFZuI+ibNi1yEgCuftqtIM6dgf19jE7Qrdx+vRppk6dymeffYbVanUWggcNGgTIbABRe5IMROOoLIG9n8CO+ee2g475K3S9QdpBa6m6EHzs2DFuv/12AgMDqaio4OGHHyY1NVUKwaLe5DKRuHzsdjjy37PtoJVnIDTaOANIugki4s2O0C1ordm6dStpaWnOQnB0dDSZmZkyG1j8jFwmEq6juh105yJjaLxvkPHpv/sUY1aAHMDq5LnnnuMPf/gDvr6+5xSCJRGIhiTJQDSM0nxjEtj2DyB7s6MddCgMf9rYHdS3mdkRuoWaheAXX3yRwYMHM2HCBCIjI5k0aZIUgsVlI8lA1J+1Eg5+adQBDqwGWyW06gLXPmdM/5J20Fqpqqpi6dKlpKWlOQvBXbt2paSkBIDExEQSExNNjlJ4OkkGom60hpztZ4fFl+ZBsxbQ9w7jMpC0g9aKzWbjyJEjxMUZ22jcc889+Pv789BDDzkLwUreR9GIJBmI2ik+fnZ30Nx94O1rXP7pPgU6DJN20FrQWrNt2zZnIdjPz4+MjAx8fHxYt24dcXFxeHt7mx2maKIkGYiLqyyBfSth+3zI+A9GO2iyox10PATI9evaWrZsGXPmzGHfvn3nFILtdjve3t507NjR7BBFEyfJQJzLbocj3xndQDXbQQfPNlYFSztoreTl5bFo0SJGjhxJXFwcFouFVq1a8fDDD0shWLgkWWcgDKcOOtpBF9ZoBx3nGBYv7aC1UVZWxvLly0lPT2fVqlVYrVb++te/8sADD5gdmmhCZJ2BqLvqdtAdCyBrk9EOGneNtIPWQ0VFBdHR0Zw6dYq2bdueUwgWwh1IMmhqbFXwo2N30AOfndsO2m0yBEeaHaHLqy4Ep6enk5mZydKlS/Hz8+Ppp5+mS5cuDB48WArBwu1IMmgKqttBdyyAXR+e1w56M7ROknbQWjh69Cjvv/8+aWlp7Nu3Dx8fH0aPHk1lZSW+vr7cc889ZocoRL1JMvBkxTk12kH3OtpBf2PUAaQdtFby8vLw8/MjMDCQTz/9lCeffJJBgwbx0EMPMWnSJMLDw80OUYgGIQVkT1PdDrrjA6MdVNuNdtDuNzt2B5Uull9TVlbGihUrSEtLY9WqVfztb39j1qxZFBcXU1BQQPv27c0OUYiLkgJyU+ZsB10Aez4y2kFDomHQI9IOWgdWq5WZM2eyZMkSTp8+TZs2bXjwwQcZMmQIAMHBwQQHB5scpRCXhyQDd3bqIOxcYAyLLzrqaAcd79gddIC0g/4KrTXbt29n+/bt3HbbbVgsFnJycpg0aRLTpk1jyJAhUggWTYZcJnI3ZQWwu7oddOPZdtDuU4xh8dIO+qsOHz7M/PnzSUtLY+/evQQFBfHTTz8REBDgHBQvhLuSy0SezFbl2B30A9i/ymgHbdkZrn3WmBcs7aC19uabb3LnnXcCMGjQIP75z38yadIkAgICACQRiCZLkoGr0hpydtRoBz1ltIP2ud2oA0R2l3bQX1GzEDxr1ixGjx7N0KFDeeGFF5g6daoUgoWoQZKBqynOOTss/uSeGu2gU6DDcGkH/RV2u501a9aQnp5+TiG4oKAAgPj4eB5//HGToxTC9UgycAWVpTXaQb822kGj+sHoVyFxgrSD/gqtNTk5ObRp0walFDNnziQ/P18KwULUgSQDs9jtcPR7IwH88DFUnna0g/7OOAuQdtBflZmZyfz580lPTyc3N5fjx4/j4+PDJ598QlxcnLMOIIT4dZIMGlveIeMS0M4FUHgUfAOhy3ijDtD+KmkHrYVvvvmGJ554gu+++w6AgQMHcv/992Oz2fDx8aFr164mRyiE+2m0ZKCUGg+MBoKBf2mtP2+sn226sgL4YZkxKtLZDno1DH0KOo2RdtBfUVZWxieffEJiYiKdO3dGKUV+fj5//OMfmTp1KjExMWaHKITbq9U6A6XU28AY4KTWOrHG7dcBcwFv4C2t9Uu1eK4w4M9a69t/7b5uvc7AVgUHv6rRDloBLTsZl4CSboTgNmZH6NJsNhv/+c9/SEtLcxaCH3vsMV588UWqf2elDVSIn7vc6wzeBf4GvFfjB3oDrwPXAlnAJqXUcozE8OJ5j5+htT7p+P8nHY/zPFrDiZ3GGYCzHTQC+txmJAFpB60VrTXdunVj7969BAcHn1MIBkkCQlwOtUoGWuu1SqmY827uBxzUWmcAKKUWAOO01i9inEWcQxn/gl8CVmmtt15K0C7nQu2gV14HPaZKO2gtVBeCN23axLJly1BKcc8999CqVSvGjBkjhWAhGsGl1AzaAsdqfJ0FJP/C/e8DhgMhSqkOWut/XOhOSqk7gTsBoqOjLyG8y6yyFPZ/6hgWX90O2tdoB+16AzSTrY1/SX5+PosWLSItLe2cQnBRURGhoaEyG0CIRtZoBWSt9WvAa7W43xvAG2DUDC53XHVit8PRdbBjfo120HZGO2jSzdCig9kRurSysjKsVitBQUF88cUX3H333XTu3FkKwUK4gEtJBtlAuxpfRzlu8zzSDlpv1YXg9PR0Fi9ezJw5c3jssccYO3YsW7dupUePHlIDEMIFXEoy2AR0VErFYiSBm4GpDRKVKygrPDss/tgGQNVoBx0Nvs1NDtC1aa15/PHHef/99zl+/DhBQUFMnDiRq6++GoCAgAB69uxpbpBCCKdaJQOl1AfA1UALpVQW8Aet9b+UUvcCqzE6iN7WWv9w2SJtDLYqOLTGqAPUbAcd/oy0g9bCkSNH+Pbbb5k2bRpKKfbu3UufPn1ITU3l+uuvl0KwEC5M5hlUt4NW7w5akmu0g3ab7NgdtIe0g/6C/Px8PvzwQ9LS0vjvf/+Ll5cXJ06coGXLltjtdrzkEpoQjcqj5hkopa4Hru/Q4TIWZE+fgJ3V7aA/nG0H7T4FOl4r7aC18PHHHzN58mSqqqro0qWLsxDcsmVLAEkEQriRpnVmUFV2dnfQQ2vOtoN2vxm6TpB20F9QsxA8cuRIbrrpJnJycvjf//1fpk2bRvfu3aUQLIQL8Kgzg0uycxF89SwUZUFIlFHwDYkyEsCej6Gi2GgHHfiwkQRadDQ7Ype2fft20tLS+OCDD5yF4MREY0eSyMhI/vznP5scoRCiIXhWMti5CFbcb5wBABQdg2WzAO1oBx3naAcdKO2gv6CgoICwMGOGwm233cbu3bsZNWqUFIKF8GCelQy+evZsInDSEBAOD+2WdtBfULMQvH37dnJycggMDOTdd98lKiqKiIgIs0MUQlxGnvXxuCjrwreXFUgiuIgdO3Ywfvx4WrduzV133UVeXh6PPfYYNpsNgO7du0siEKIJcMkzg3p3E4VEGZeGLnS7AIxC8DfffEPLli3p1q0bABs3buS+++5j2rRpsiJYiCbKs7qJzq8ZAPgEwPWvGYvGmiitNTt27CA9PZ358+dz/Phxbr/9dt566y201tjtdpkRLISHkG4iOHvAr9lNNOx/mnQiABg5ciRffPEFFovlnEIwGLMBJBEIITzrzEA4C8GrV6/mww8/xNvbm7///e8opZg8ebJc/xfCw8mZQRNWPSM4PT2dTz/9lKqqKjp37kx2djbR0dHcfffdZocohHBxkgzclM1mo7y8nObNm7N27VpuvPFGIiMjpRAshKgXSQZu5PxC8PTp03nppZcYNmwYX375JVdffbVc/xdC1ItLJoNG2ajOzcydO5c333yTH374AYvFwm9+8xsGDx4MgMViYdiwYSZHKIRwZy656ExrvUJrfWdISIjZoZgmPz+fhQsXOr/euHEjISEhzJs3j5ycHJYvX86oUaNMjFAI4Ulc8sygqSovL3cWgleuXElVVRU9evQgISGBd999Fx8f2VZbCHF5SDJwEd999x2jR4+mqKjonELwlVdeCSCJQAhxWUkyMIHWmp07d5Kenk5iYiK33HIL3bp1Y8KECUydOpVrrrlGCsFCiEYlyaARHT16lPnz55Oens7u3buxWCw8/PDDAAQHB/P222+bHKEQoqmSZHCZlZaW0qxZMwBuvfVW/vOf/3DVVVcxb948Jk+eTIsWLUyOUAghJBlcFjULwZ9//jmZmZm0atWKV199ldDQUGJjY80OUQghzuGSycBd1xlkZmby/PPPs2TJEoqKimjdujWzZs1yzgbo2bOnyREKIcSFyUZ1l6C6EKy1pkePHhw5coSkpCRuuOEGUlNTGTp0qBSChRCNSjaqa0TnF4LHjRvHRx99RPv27Tl58iR+fn5mhyiEEHUiyaCOZsyYwTvvvAPAgAEDeP3117nxxrPzEiQRCCHckSSDX1BeXs7KlStZsmQJb7/9Nv7+/vTv35+4uDimTp1KXFyc2SEKIUSDkGRwHrvdzjfffEN6ejqLFy92FoIPHDhAUlISM2fONDtEIYRocJIMMArBFRUV+Pv7s3nzZoYOHUpgYCATJ06UQrAQoklo0smgZiF4wIAB/POf/6Rv374sW7aMESNGOBeLCSGEp2uSyeCDDz7gH//4B2vXrgWMQvBVV10FGAPix48fb2Z4QgjR6FxynoFS6nql1BtFRUUN8nzl5eWsWLGC6jUV3377LT/99BPPPfcchw4d4rvvvuOWW25pkJ8lhBDuyGMXnV2oEPzdd98xYMAAysrK8Pf3lxnBQgiPI4vOatizZw8jR44kKyuLwMBAbrjhBqZNm0a/fv0ACAgIMDlCIYRwLR6ZDOLj4+nfvz8TJkxg7NixUggWQohf4ZHJwM/Pj0WLFpkdhhBCuA2XLCALIYRoXJIMhBBCSDIQQgghyUAIIQSSDIQQQiDJQAghBJIMhBBC4KLJoKH3JhJCCPHLXHpvIqVULnCkng9vAZxqwHCEqCkEkE8rl1dTfo8v5bW311q3rOuDXDoZXAql1Ob6bNYkRG0opd7QWt9pdhyerCm/x2a8dpe8TCSEG1hhdgBNQFN+jxv9tcuZgRBCCI8+M3jD7ACEEMJdeOyZgRBCiNrz5DMDIYQQtSTJQAiTKKXilFL/UkotNjsWT9WU3+O6vnZJBqJJUEq1U0p9rZTao5T6QSn1wCU819tKqZNKqd0X+N51Sqn9SqmDSqnHful5tNYZWuvb6xuHq1FK+SulNiqldjje42cu4bnc8j1WSnkrpbYppT65hOcw5bU3yWTQlD8tNGFW4Hda6y5ACnCPUqpLzTsopVoppYLOu63DBZ7rXeC6829USnkDrwO/AboAU5RSXZRS3ZRSn5z3p1XDvCyXUgEM1Vp3B3oA1ymlUmreoQm8xw8Aey/0DVd/7W6XDC6WNV3504Iwn9Y6R2u91fH/pzH+wbY9725DgI+UUn4ASqmZwP9d4LnWAvkX+DH9gIOO369KYAEwTmu9S2s95rw/Jxvu1bkGbTjj+NLH8ef8DhWPfY+VUlHAaOCti9zFpV+72yUDLpA13eDTgnAhSqkYoCewoebtWusPgdXAQqVUKjADmFyHp24LHKvxdRY/Tzg144hQSv0D6KmUerwOP8dlOS6TbAdOAl9orZvSe/xX4FHAfqFvuvprt9QhEJegtV7r+MdckzNjAiilqjPmi8CYxo1QuDKlVCCwBHhQa118/ve11i87fn/+DsTX+KTb4LTWecBdl+v5zaC1tgE9lFKhwDKlVKLWevd59/G491gpNQY4qbXeopS6+hficdnX7o5nBhfiyp8WhItQSvlgJIJ0rfXSi9xnEJAILAP+UMcfkQ20q/F1lOO2JkdrXQh8zYWvfXvie3wVMFYpdRjj8s1QpVTa+Xdy5dfuKcmgTrTWeVrru7TW8Y6zB+HhlFIK+BewV2v96kXu0xNj5fo44DYgQin1fB1+zCago1IqVinlC9wMLL+0yN2HUqql44wApVQAcC2w77z7eOR7rLV+XGsdpbWOwYhpjdZ6Ws37uPpr95Rk4KqfFoTruAqYjvGJbbvjz6jz7tMMuFFrfUhrbQdu4QJbqCulPgDWAQlKqSyl1O0AWmsrcC/GdeG9wCKt9Q+X7yW5nEjga6XUTowD1xda6/NbLJvye+zSr90tt6Nw1Aw+0VonOr62AAeAYRhJYBMw1Y1+SYQQwlRud2Zwoazp5p8WhBDCdG55ZiCEEKJhud2ZgRBCiIYnyUAIIYQkAyGEEJIMhBBCIMlACCEEkgyEEEIgyUAIIQSSDIQQQiDJQAghBPD/AfBujAv4+KJ1AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1071fb6a0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.loglog(dt, L1_error, 'C0o-')\n",
    "plt.loglog(dt, L2_error, 'C1o-')\n",
    "plt.loglog(dt, 0.3*dt, 'k:')\n",
    "plt.loglog(dt, 0.15*dt**1.5, 'k--');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.6.5"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "navigate_num": "#000000",
    "navigate_text": "#333333",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700",
    "sidebar_border": "#EEEEEE",
    "wrapper_background": "#FFFFFF"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "30px",
    "width": "252px"
   },
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 4,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false,
   "widenNotebook": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}