{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Numerical Methods \n",
    "\n",
    "# Lecture 2: Numerical Differentiation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Learning objectives:\n",
    "\n",
    "* Learn about finite difference approximations to derivatives.\n",
    "* Be able to implement forward and central difference methods.\n",
    "* Calculate higher-order derivatives.\n",
    "* Solve simple ODEs using the finite difference method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# some imports we will make at the start of every notebook\n",
    "# later notebooks may add to this with specific SciPy modules\n",
    "\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Differentiation (or taking a derivative) - the continuous definition\n",
    "\n",
    "\n",
    "The classical definition of the derivative of a function $f$ at a point $x_0$ is of course given by:\n",
    "\n",
    " $$f'(x_0)=\\lim_{h\\rightarrow 0} \\frac{f(x_0+h)-f(x_0)}{h}. $$\n",
    " \n",
    " \n",
    "**Notation** \n",
    " \n",
    "1. Note that the following are all equivalent mathematical ways of writing the derivative of the function $f$ with respect to (w.r.t.) $x$ and evaluated at the location $x_0$:\n",
    " \n",
    " $$  f'(x_0) = \\frac{df}{dx}(x_0) = \\left.\\frac{df}{dx}\\right|_{x_0}. $$\n",
    " \n",
    " \n",
    "2. We're using $h$ here to denote a small (potentially infinitesimally small) increment to the $x$ coordinate, as is common. But note that in the literature $\\Delta x$ is commonly used to mean the same thing. Also, of course for finite $h$ there is significant overlap here with the mesh spacing in a numerical approximation, as we shall see below, and so $\\Delta x$ is also used."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Finite differences \n",
    "\n",
    "Finite differences are a class of approximation method for estimating/computing derivatives of functions.\n",
    "\n",
    "\n",
    "In this lecture we shall largely focus on the case of functions of a single spatial dimension $x$ (e.g. $\\;f\\equiv f(x)\\;$ or $\\;u\\equiv u(x)$), but these ideas extend for time derivatives (initial-value problems - IVPs - at end of this lecture) and multiple spatial dimensions (PDEs).\n",
    "\n",
    "\n",
    "Approximations to the derivatives of a function can be computed by using weighted sums of function evaluations at a number of points.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## The forward difference\n",
    "\n",
    "For example, we can turn the formal definition of a derivative given above into an approximation rule by replacing the limit as $h$ approaches zero (i.e. the $\\text{lim}_{h\\rightarrow\\infty}$) with a small but finite $\\Delta x$ value:\n",
    "\n",
    "$$ f'(x_0)\\approx \\frac{f(x_0+\\Delta x)-f(x_0)}{\\Delta x},\\;\\;\\;\\; \\Delta x>0. $$\n",
    "\n",
    "Since this approximate gradient method uses values of $x$ greater than $x_0$ ($\\Delta x>0$), this algorithm is known as the **forward difference method**. \n",
    "\n",
    "The figure below illustrates this approximation. \n",
    "\n",
    "In the figure the derivative is approximated by the slope of the red line, while the true derivative is the slope of the blue line -- if the second (and/or higher) derivative of the function is large then this approximation might not be very good, unless you make $\\Delta x$ very small.\n",
    "\n",
    "You can see this yourself by varying the value of $\\Delta x$ in the code below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcYAAAGTCAYAAACyIbJtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3gUVffA8e8FEkJC7xAQFASp0pQOUeHVV35YAEEkNAUBQUEp0rs0kWJBOqEEiYDCiyCRFhSpERCRGKQjBAioFCkJyfn9Mcu6CSkbSLLZcD7PMw/MzL0zZzbJnJ07984YEUEppZRSliyuDkAppZTKSDQxKqWUUg40MSqllFIONDEqpZRSDjQxKqWUUg40MSqllFIONDGqDMsY08kYI4lMTVwdX1owxvjZjs/vHuufMMYEOMzf+QxLOyzzMcYsNsZcsK2bZlte1xizyxjzj215tfs7GpXR2H6uI10dR0aXzdUBKOWEV4A/4i075IpA3NBaoC4Q4bCsJ9AWeB047LBuHnADaA5ct61T6oGjiVG5g/0iciS1N2qMyS4it1J7uxlpvyISCUTGW1wBOCsiixxiygKUBz4Qkc2psW9Xfb5K3S9tSlVuzxhT3hjztTHmb2PMDWPMTmPMc/HKjLQ1I1U2xgQbY64BXxpj+tmaDj0dyq6M31xrjOlqjLltjMltm3/CGLPCGPOHbZ/hxphxxpgc8fYbYozZZoxpbozZZ4y5BbxlW1fIGLPUGHPFFvsiIG8Kjru3ren0pjEm1BjTMIEycZpSjTECdAJKOjRLdwJisM4Hw2zLTjhso7ExZpMx5qrtswo2xlROwXFmM8YMMsb8Zoy5ZYw5a4z5yBjj5VC/tG2/3Ywxo40xEbbPZI0xpkQCx9XVGLPX9tn/ZYzZaoyp57De2xgz0Rhz3BgTZft3iO0LQHKfa7J1bU3RfxtjSjksK2aMiTTGLHdY9h9jzDrb8Vw3xhw0xvQ1xmSNt88Txpglxpj2tt+lG8aYH4wxjxqr6XuWMeaSMea87bPL5lD3TvN7S2NMgO3zuGKMCTTGFHDieB83xvzPVu+GMebHhH6XHigiopNOGXLCOoEL1pVMNocpq0OZ4lhXRMcAf6xmwPVYJ/r/OpQbadvWUWAw8DTgB9SwLW9kK2eAi1hNieMc6n8B7HKYbwkMBf4PaIyVBM4By+IdQwhwATiO1XTpB1S1rfsBuAL0Ap4F5gOnbfH4JfPZvGErtwB4zraNP4DLQEACn2Fp23wd2+cTYft/HaAUUN9Wbq5tWXVb+WbAbWA18KJt2g78BZR08jiXAf8Aw4EmwNvA38BKh/qlbfs/ASwF/gt0tP0stsY79skOsTa3xTgGeNW2Ppvts70E9AGeAYYAN4GPkvlcnaoL5AKOAD8CWbG+VGwCTgJ5Hcp1B/rajucpoD9wFZgQb78ngFPADuAloDVwFjgArLIdc1PbcQrwlkNdP9uy0w6/D2/b9rMl3n4EGOkwX8P2s9kGtAKeB/4H3AJquvoc4LJzj6sD0EmnxCb+PanHn7Y5lJlsO3GXdViWFQgH9josG2mr2zvePrIAfwIjbPPVgFhgKrDDoVxE/JOZwzpjO6H62+oWcFgXYltWLV6dprZ4Xo23/FuSSYy2mE8D6+Mtb2OrG5DAZ1jaYdkS4ES8utninzRty48Am+Ity42VsKY5cZwNbdvtEG95O9vyarb50rb5+Emwn215cdt8WawvPVOS+Hza4/Blx2H5ECAKKJwadYEnbctGYX3Zug00TGLbd35PhmB9scjisO6E7fcwj8Oyd2yxzI23nb04JDz+TYzxfx/ufMbPOCyLnxg3AWGAZ7y/nzBgVVr9bWf0SZtSlTt4GXjCYXrDYV0jYKc43IMUkRisK7xqd5o+HXztOCMiscD3WFeQ2P49AHwJ1DLG5DLGVASKAvZ7b8aY3LbmtqNY366jgcVYJ79H4+3zhIjsj7esLtYJfmW85csS/ATiKmGbvoy3fCXWyTlVGGMeBcoAgbbm0Gy2JrzrWFc2jeJVSeg4n8NKHivjbeM72/r421gbb/4X278P2f5tgvXFYHYSoT+HdeW2PYF9emBdEd93XRHZjXUVPAQrOX4gIj84bszWvDrLGHPS9jlEA2OxmswLx9v3DhG57DD/m+3f4HjlfgNKJhB7/N+H5VhfVuomdKDGavZvfKecw7EaYCN3/2weGNr5RrmDg5J455v8wL4Elp/D+gPPh9VceUdEAmU3Ax/aThRPAVuAPVjNZw2Bh7FOaD861FmAdZIeDuzHao56EvgM8CKuhPZZDPhLRKLjLT+fQNmE6t5VVkRuG2MuOVHfWXdO3PNsU3yn4s0ndJyFAU/gWiL7iH8P7M9483c679z5TO+Uj99LOf4+S2H9zJzZ5/3UXcq/zZufOa6w3ZP8H1Zz/0ishHYDq6l0CHf/nvwVbz4qieXx68Ldvw9Rxpi/AN+ED4X8WFeHw2zTXYwxWWxfHh8omhiVu/sT62ouvqJYJ6v4J9qE3rO2Bevk3cg2zbYlmR+wriAfBnaLyD8Atk4jL2I1SU2/sxFjTJVEYkxonxFAPmOMR7zkWCSRbcSve1dZ27f9ZDtbpMCdJDsI6woivqh48wkd5yX+/YKRkLMpjOmi7V9frObyhFzCutfZOpH1J5LYvtN1bYlvIVaSzgXMBFo4lC0D1ALai8gSh3rNk9j//Yj/++CJ9cXwTCLl/8a6ovwMWJRQgQcxKYImRuX+tgJ9jDGlReQEgK3HXxtgn4hcdWIbB7E68PQHfLCaVsG6kmyH1Ww1w6F8dqxv2vGvKjqlIO4dtm20JG7z6atO1P0D6x5ja6wOO3e0JHX/psOxEkElEZlwj9tYD7yPde9sUyrEtBHrZP4mVqeWxPbZErgmIr8lUiYxKak7CGiA9WWqELDaGNNNRGbZ1nvb/rX/nhhjPLB+p9JC/N+HV7CanXckVFhE/rF9+Xsc6378A5kEE6KJUbm7qVgJaYMxZgRWs+lbQDms3orJEhExxoRgnUj2ONzn2QJ86PD/O+UvG2N2An2NMRFYVzGvk3iTVUL73GCM2QbMMsYUBH7HSuaVk65pfYs3xowC5hpjFmAl1rJYJ+orSVZOAdvn0hPrhO+JdQ/rItaVST3glIhMSWYbIcaYL4AVxpgpwG6sxFYaqwfk+yLi9IMEROSoMWYq8J4xJhdWU2UMVjP2byISBAQCnYFNxpiPgJ+xWgTKAC8AL4nI9UR24VRdY0xtrObRUSKyA8AYMwOYYoz5XkTCsDqwnAQ+MMbEYCXId5091ntQyeH3oRzwAVZnpqS+kLyH9UUw2BgzD6s1oiBWb9WsIjIwDePNsLTzjXJrInIW61v7r8DnwAqseyfNRGR9CjZ1J/E5Dm7fh3V/5xZ3f+tuC/yE1QwVgHVPs3cKw28BrAPGA0FYX1R7OVNRROZhDSd4GmsoRWesq83496Pui4isw7oi8sEaHhEMTMJqqk7wSiQB/lhJpJUt1hVYx/k7zt1TjR9TP6wvP3WwOhwFYt0bPmVbH401/GUO1pXlOluZjlhDTeI3ATtuO9m6toS81DY/zqF6X6xhQ18Y6+EGUVj3E89hNVV+hpWE7vXqOzm9se6rB9ni+gbrM0+UiOzF6tB2CfgYq5PRdKAK/7acPHCMrXuuUkopN2Ss5+puAZqKSEL3glUK6RWjUkop5UATo1JKKeVAm1KVUkopB3rFqJRSSjnQxKiUUko50HGMbq5gwYJSunRpV4ehlFJu5aeffrooIoUSWqeJ0c2VLl2a0NBQV4ehlFJuxfZg9wRpU6pSSinlQBOjUkop5UATo1JKKeVAE6NSSinlQBOjUkop5UB7pWZyV65c4cKFC0RHJ/ZCcqVURuTh4UHhwoXJnTu3q0N54GhizMSuXLnC+fPn8fX1JUeOHBhjXB2SUsoJIsKNGzc4c+YMgCbHdKZNqZnYhQsX8PX1xdvbW5OiUm7EGIO3tze+vr5cuHDB1eE8cDQxZmLR0dHkyJHD1WEope5Rjhw59DaIC2hizOT0SlEp96V/v66hiVEppZRyoIlRqXTQq1cv/Pz8UlQnNDQUYwwnTpwAICQkBGMMFy9etJdZvXo1jz76KNmyZaNTp06JLlOJmzBhAmXLlr3v7Rw5cgRjDPv370+FqJQraWJUGU6nTp0wxtw1PegnnHr16hEREUGBAgXsy7p06ULLli05efIk06dPT3SZSnsPP/wwERERVK5c2dWhqPukwzVUhtSkSRMWL14cZ1nBggXveXvR0dF4eHjcb1hJio2NRUTImjVrmmzf09OTokWL2uf//vtvLl68yLPPPouvr2+iyzK6qKgoPD09XR3GfblzDI4/H+W+9IpRZUjZs2enaNGicaZs2azvcbdu3aJPnz4UKVIELy8v6tSpw7Zt2+x17zQ5rlu3jieffBJPT0+Cg4Px8PBg165d9nIlSpSgQoUK9vkNGzbg4+Nj7wU4ZcoUqlatio+PD76+vnTp0oW///7bXj4gIICcOXOybt06KleujKenJ2FhYcTExNCvXz/y5ctHvnz56NOnDzExMcke8/r163nsscfw8vKiYcOGHD58OM56x6bUkJAQ8uXLB8DTTz+NMSbRZc6IjIzk1VdfpUSJEnh7e1O5cmUWLVoUp0yDBg3o1asXvXr1Il++fOTPn5+BAwcSGxsb5zMdM2YMr732Gjlz5qRYsWJMnTrVvv727dsYY5g5cyYvvvgiPj4+DB8+3H58Tz75JF5eXhQtWpR+/foRFRUFwPnz5ylSpAjjxo2zb2vfvn1kz56dVatWOXWMd4wfP54iRYqQK1cuOnXqxPXr1+8qM3fuXCpUqICXlxfly5fn448/RkSSPAbHptSYmBiKFSvG559/Hme7YWFhGGP45ZdfUhSzSmciopMbTzVr1pTEHDp0KNF1GVnHjh2lWbNmia5/5513pGjRovLNN9/IoUOHpEuXLuLj4yNnz54VEZEtW7YIIJUrV5bg4GA5evSoXLhwQZ588kkZP368iIgcPnxYcubMKV5eXvZ6gwcPliZNmtj3M3XqVNm0aZMcP35cQkJCpEqVKuLv729fv2DBAsmaNavUrVtXtm3bJuHh4XLlyhWZOHGi5M6dW4KCgiQsLEx69eoluXLlksaNGyd6TKdOnZLs2bNLr169JCwsTIKCgsTX11cAOX78eJzjioyMlFu3bsmvv/4qgKxcuVIiIiISXeaMkydPyocffij79u2To0ePyueffy7ZsmWTLVu22MvUr19fcubMKb1795awsDD54osvJFeuXDJ9+nR7GV9fX8mVK5eMGzdOwsPDZcaMGZItWzZZvXq1iIhER0cLIIULF5Z58+bJ0aNH5fjx43Lq1Cnx8vKSHj16yKFDh2T16tVSqFAhGTBggH3b33zzjXh4eMiuXbvkn3/+kccee0y6dOni1PHdERgYKJ6enjJ79mwJDw+XUaNGSa5cuaRMmTL2MjNmzJBixYrJihUr5NixY7Jq1SopVKiQfP7550kew++//y6A7Nu3T0RE+vTpI/Xr14+z/8GDB0vVqlVTFLO7/h1ndECoJHJedfmJXad0Toy9e4s0bpy+U+/eicaYkI4dO0rWrFnFx8fHPj333HMiInLt2jXx8PCQhQsX2svfvn1bHnnkERkyZIiI/JtAVqxYEWe7AwYMkP/85z8iIjJ79mx57rnnpFGjRrJ06VIREalXr56MHTs20bi+/fZb8fT0lJiYGBGxEiMgoaGhccoVK1YsznZiYmLk0UcfTTIxDho0SB599FGJjY21LxszZkyiiVFEJDIyUoA4ySuhZfeqZcuW0q1bN/t8/fr1pUKFCnFiHDFihJQqVco+7+vra/9Z3dGxY0f7sd9JKn369IlTZsCAAVK+fPk4254zZ45kz55dbty4YV/Ws2dPKVu2rHTo0EHKlSsn165dS9ExPfHEE9K9e/c4yxo3bhwnMRYvXtz+O3HHhx9+KFWqVEnyGOInxtDQUDHG2H9+sbGxUrp0aZk0aVKKYtbEmDaSSozalKoypEaNGrF//377NHfuXACOHj1KdHQ09evXt5fNmjUrdevW5dChQ3G2UatWrTjzfn5+/Pjjj0RHRxMSEsJTTz2Fn58fISEhXL9+nT179sTpObp582aaNm1KiRIlyJUrFy1atCAqKopz587Zy2TLlo1q1arZ5y9fvkxERAR169a1L8uSJQu1a9dO8njDwsKoU6dOnHFrjttIa7dv32bMmDFUrVqVAgUKkDNnTlavXs2pU6filEsoxpMnT/LPP/8kGrczP5uwsDDq1q0bZ9sNGjTg1q1bHDt2zL5s8uTJGGMIDAwkMDAQHx+fFB3nnf3Ej++OiIgIzp49yxtvvEHOnDnt09ChQzl69GiSxxBfzZo1eeyxx1i6dCkA27dv59SpU7z22mspilmlP+1886CZNs3VETjF29s7wS701he9hAc+x18W/6TZsGFDbt26xZ49e9i6dSt9+vTh2rVrdOvWjR9//BEPDw+efPJJAE6ePEmzZs3o2rUro0ePpkCBAuzdu5e2bdva73uBdS80NTrb3DkuV5k4cSLTp09n2rRpVK5cmZw5c/L+++9z+fLlNNlf/J+NiCQ6mN1x+bFjx/jjjz8AOH78eLLJKaXu3C+dM2fOXV9mkvv9Ski7du0IDAxk8ODBBAYG4ufn5zadoh5kesWo3ErZsmXx9PSM09kmJiaGHTt2ULFixSTr5syZkxo1ajB79myuXr1KjRo1qFu3LqdOnSIwMJB69erZe66GhoYSFRXF1KlTqVu3LuXKlePs2bPJxpcnTx6KFSvGzp077ctEhN27dydZr2LFiuzatStOgnTcRlrbtm0bL774Iv7+/lSrVo0yZcoQHh5+V7n4Me3cuZOHHnooTpJIqIxjJ6eEVKxYke3bt8c5/m3btpE9e3YeeeQRwOr52a5dO1q2bMmECRPo1q2bPUk6q0KFCgnGd0fx4sUpUqQIx44do2zZsnGmMmXKpGhfAP7+/oSFhbFnzx6WL1+Ov79/ireh0p8mRuVWfHx86NGjBwMHDmTdunWEhYXRo0cPzp8/z1tvvZVsfT8/P5YsWULDhg3JmjUrXl5e1K5dmyVLlsRpRn300UeJjY1l2rRpHD9+nC+++IJpTl5t9+7dm0mTJrFixQrCw8Pp06cPERERSdbp3r07J06coE+fPoSHh7NixQpmzpzp1P6Sc+rUKR577DHWrFmTaJly5cqxYcMGtm/fbv9MT58+fVe506dP89577xEeHs6XX37JlClTePfdd+OU+fHHH5k0aRK///47M2fOJDAw8K4y8fXq1YuTJ0/Sq1cvfvvtN9asWcPgwYPp3bs32bNnB2Dw4MH89ddffPbZZ/Tt25caNWrQoUMHezLdsWMHjz32GHv37k10P71792b+/PnMmzePw4cPM3bsWH766Sf7emMMI0eOZPz48UyfPp3w8HAOHjzIwoULmThxYpLHkJBSpUpRv359unbtyrVr12jZsmWc9a+++ipdu3ZN8XZV2tLEqNzOxIkTad26NZ07d6ZatWocOHCA9evXU6xYsWTrPvXUU8TExMRJggktq1q1KtOnT2fKlClUrFiRuXPnMnnyZKfi69u3L507d6ZLly7Url2b2NhY2rVrl2Sdhx56iK+++or169fz+OOPM3XqVCZMmODU/pITFRVFeHh4ks2iI0aMoGbNmjz77LM0btyYvHnz0qZNm7vKdejQgRs3bvDkk0/SrVs3unXrxjvvvBOnTN++ffnpp5+oVq0aI0eOZNy4cbz00ktJxliyZEm+/fZb9uzZw+OPP06XLl1o3749Y8aMAWDTpk1Mnz6dJUuWkDt3bowxLFq0iAMHDth/Lv/88w/h4eEJDr+4o127dgwdOpRBgwZRo0YNfvvtN3r37h2nTPfu3Zk9ezYBAQE8/vjjNGrUiLlz5/Lwww8neQyJad++PT///DMvvPDCXa+POnHiRIJfQJRrGVff21D3p1atWhIaGprgurCwsGSbsJRyVoMGDahVq1aSV84lSpSgX79+9OnTJx0jy9z07zhtGGN+EpEEb1LrFaNSSinlQBOjUkoplxo1apTT9/DTgw7XUEo5xbEncGJS2ktUqdOnTzNu3LgM1WNXrxiVUkq5zAcffICIMGzYMFeHYqeJUSmllEscO3aMefPm0bVrV0qXLu3qcOw0MSqllHKJ0aNHky1bNoYMGeLqUOLQxKiUUirdhYWFsXjxYt566y2KFy/u6nDi0MSolFIq3Q0bNgxvb28GDRrk6lDuoolRKaVUugoNDWXlypX07duXggULujqcu2hiVCqVTZ48OVU6Epw4cQJjDIk92Sg1+fn50atXrzTfj1IAQ4cOJX/+/Lz33nuuDiVBmhhVhtOpUyeMMXdNderUSbcYRo4cSeXKldNtfwkpWbIkERERcd73eL8CAgLImTPnXcu/+uorxo8fn2r7USoxW7duJTg4mEGDBt317NiMQgf4qwypSZMmLF68OM4yT09PF0WT/qKiovD09KRo0aLpsr/8+fOny37Ug01EGDx4MMWLF6dnz56uDidResWoMqTs2bNTtGjRONOdk/fWrVvx8PAgJCTEXn7mzJnkzp3b/rb39evX07BhQ/Lly0f+/Pl59tlnCQsLi7OPs2fP0q5dOwoUKIC3tzfVqlVjy5YtBAQEMGrUKH799Vf71WpAQECisU6aNImiRYuSM2dOOnTowLVr1+4qs2DBAipWrIiXlxflypVj6tSp9pfigvW6o88++4wWLVrg4+PD4MGD4zSlxsbGUqJECT755JM42z18+DDGGPbt2wfAlClTqFq1Kj4+Pvj6+tKlSxf+/vtvAEJCQujcuTP//POP/bhGjhwJxG1KHTRoEDVr1rzrGOrVqxfnTRTJHZNS8a1bt47t27czfPhwcuTI4epwEiciOrnxVLNmTUnMoUOHEl2XkXXs2FGaNWuWZJlBgwZJiRIl5NKlSxIWFibe3t4SEBBgX79ixQpZsWKFHD58WH7++Wd55ZVXpEyZMnLr1i0REbl27ZqULVtW6tWrJ1u3bpUjR47IypUrZfPmzXL9+nXp27evlC9fXiIiIiQiIkKuX7+eYBxBQUHi4eEhM2fOlPDwcBk7dqzkypVLSpUqZS8ze/ZsKVq0qCxfvlyOHTsm//vf/6RIkSLyySef2MsAUqhQIZkzZ44cPXpUjh07JsePHxdA9uzZIyIi/fr1k9q1a8fZ//Dhw6VixYr2+alTp8qmTZvk+PHjEhISIlWqVBF/f38REbl165ZMmzZNvL297cd19epVERFp3Lix9OzZU0REfv31VwEkLCzMvt1jx44JILt373b6mFTqcNe/4/hiYmLk8ccflzJlykhUVJSrwxEgVBI5r7r8xK5T+ibG3r1FGjdO36l370RDTFDHjh0la9as4uPjE2caMGCAvUxUVJTUqlVLXn75Zalevbq0bt06yW1eu3ZNsmTJIj/88IOIWCf2nDlzSmRkZILlR4wYIZUqVUo21rp160qXLl3iLHvmmWfiJMaSJUvKokWL4pSZOnWqVKhQwT4PSK9eveKUiZ8Yf/75ZwHk999/t5cpW7asjBs3LtH4vv32W/H09JSYmBgREVmwYIH4+PjcVc4xMYqIVKtWTYYOHWqfHzNmjJQrVy5Fx6RSR2ZJjMuWLRNAAgMDXR2KiCSdGPUeo8qQGjVqxOzZs+Msy5s3r/3/Hh4eLF26lEqVKlG4cGE2b94cp+zRo0cZNmwYu3btIjIyktjYWGJjYzl16hQA+/bto2rVqvfdVTwsLIwuXbrEWVa3bl2OHDkCQGRkJKdPn6Zbt2706NHDXub27dvWN1MHtWol+Go4u6pVq1KlShWWLl3K8OHD2bVrF0ePHuW1116zl9m8eTPjx48nLCyMy5cvExMTQ1RUFOfOnUvRIGp/f39mzJhhf1FwYGCg/SHPKTkmpQCio6MZNmwYVapU4dVXX3V1OMnSxPiAyUBvdkmSt7c3ZcuWTbLMzp07iY2N5e+//yYyMjJO4mzevDm+vr7MmjULX19fsmXLRsWKFYmKigJItxP4nXtuM2fOpF69ekmW9fHxSXZ77dq1Y/78+QwfPpzAwEAaNmxIqVKlADh58iTNmjWja9eujB49mgIFCrB3717atm1rP25nvfbaawwYMIAdO3aQPXt2fvvtN9q1a5fiY1IKYOHChfz++++sXr2aLFkyfteWjB+hUgk4ceIEvXr14rPPPqNp06a0a9eO27dvA3Dp0iXCwsIYPHgwTZo0oUKFCly9etW+HqBGjRocOHCAixcvJrh9T09PYmJiko2jQoUK7Ny5M84yx/kiRYrg6+vL0aNHKVu27F1TSrVr144jR46wc+dOgoKC4ryqJzQ0lKioKKZOnUrdunUpV64cZ8+evafjKlasGE8//TSBgYEEBgZSr149HnnkkTQ5JpW53bx5k1GjRlGnTh2aN2/u6nCcoleMKkO6desW586di7Msa9asFCpUiJiYGPz9/WncuDHdunWjVatWVKlShVGjRjFmzBjy5ctHwYIFmTNnDiVLluTMmTP079+fbNn+/XV/7bXXmDBhAi+99BLjx4+nRIkS/PLLL+TKlYunnnqK0qVLc/LkSfbu3ctDDz1Erly5yJ49+11x9u7dmw4dOvDEE0/g5+fHihUr2LVrV5zhDyNHjuTtt98mb968PP/880RHR7N3717OnDmT4sdhlShRgkaNGtG9e3cuX77MK6+8Yl/36KOPEhsby7Rp02jRogU7d+686+WvpUuX5ubNm2zYsIHq1avj7e2Nt7d3gvvy9/enX79+eHp6MnTo0DjrUvOYVOY2c+ZM/vjjDxYtWoQxxtXhOCexm486uceUWXulAndNvr6+IiIyevRoKVKkiFy4cMFe57vvvpNs2bLZO9ds2rRJKlWqJNmzZ5dKlSrJ+vXrxcfHRxYsWGCvc/r0aWndurXkyZNHcuTIIdWqVZMtW7aIiMjNmzelZcuWkjdvXgHi1Itv3LhxUqhQIfHx8ZG2bdvKiBEj4nS+ERFZunSpVK9eXbJnzy558+aV+vXryxdffGFfD8jy5cvj1Inf+eaOefPmCSAtWrS4K5bp06dL8eLFxcvLS55++mkJCgoSQHM8ldwAACAASURBVI4fP24v0717dylQoIAAMmLECBG5u/ONiMjVq1fF29tbPDw85OLFi3ftK7ljUqnDXf+ORUSuXLkiBQsWlCZNmrg6lLuQROcbI3qz3K3VqlVLEntkWFhYGBUqVEjniJRSqcmd/47HjBlj7yj25JNPujqcOIwxP4lIgj3e9B6jUkqpVHfx4kUmT57MSy+9lOGSYnI0MSqllEp1Y8eO5dq1a4wbN87VoaSYJkallFKp6vjx48yYMYPXX3/dLZuBNTEqpZRKVUOHDiVbtmz2Z/G6G02MmZx2rlLKfbnj3+++fftYunQpffr0wdfX19Xh3BNNjJmYh4cHN27ccHUYSql7dOPGDTw8PFwdRoq8//775M+fn/fff9/VodwzTYyZWOHChTlz5gzXr193y2+eSj2oRITr169z5swZChcu7OpwnLZhwwY2bNjA0KFDyZMnj6vDuWc6jtHNJTWOEeDKlStcuHCB6OjodIxKKXW/PDw8KFy4cIZ9y318sbGx1KpViz///JPw8PAEnxSVkSQ1jlEfCZfJ5c6d223+sJRS7isoKIh9+/axePHiDJ8Uk6NXjG4uuStGpZRKa1FRUTz22GPkzp2bvXv3usUbNPSKUSmlVJqZOXMmx48fZ/369W6RFJPj/keglFLKZa5cucKYMWN4+umn+c9//uPqcFKFJkallFL37MMPP+TixYtMnDjRfV4rlQxNjEoppe5JREQEU6ZMoU2bNtSqleDtOrekiVEppdQ9GTVqFFFRUYwdO9bVoaQqTYxKKaVSLCwsjLlz59KtWzfKli3r6nBSlSZGpZRSKTZgwAB8fHwYMWKEq0NJdTpcQymlVIps3ryZb775hgkTJlCoUCFXh5Pq9IpRKaWU02JiYujbty+lSpWid+/erg4nTegVo1JKKactXryY/fv3s3TpUry8vFwdTprQK0allFJOuX79OkOGDOGJJ56gTZs2rg4nzegVo1JKKad89NFHnD17lqCgoEzx6LfEZN4jU0oplWoiIiKYOHEiLVq0oEGDBq4OJ01pYlRKKZWs4cOHc+vWLSZOnOjqUNKcJkallFJJ+uWXX5g/fz49e/bMdIP5E6KJUSmlVJL69+9P7ty5GT58uKtDSRfa+UYppVSigoODCQ4O5qOPPiJ//vyuDiddGBFxdQzqPtSqVUtCQ0NdHYZSKhOKiYmhWrVqXL9+nUOHDpE9e3ZXh5RqjDE/iUiCrwTRK0allFIJWrBgAQcPHmT58uWZKikmR+8xKqWUusvVq1cZNmwY9erVo2XLlq4OJ13pFaNSSqm7jB8/nnPnzrFq1SqMMa4OJ13pFaNSSqk4jh07xkcffUT79u2pXbu2q8NJd5oYlVJKxdG/f3+yZcvG+PHjXR2KS2hiVEopZRcSEsJXX33FoEGD8PX1dXU4LqHDNdycDtdQSqWWmJgYatSoweXLlwkLCyNHjhyuDinN6HANpZRSyZo7dy4HDhzgyy+/zNRJMTnalKqUUoq///6boUOH0rBhQ1q1auXqcFxKE6NSSinGjBnDpUuXmD59+gM3PCM+TYxKKfWACw8P5+OPP+aNN96gevXqrg7H5TQxKqXUA65v377kyJGDsWPHujqUDEE73yil1AMsODiYtWvXMmnSJIoUKeLqcDIEHa7h5nS4hlLqXkVHR/P4448THR3NwYMHH6gHhetwDaWUUneZOXMmYWFhrF69+oFKisnRe4xKKfUAioyMZPjw4TRp0oTmzZu7OpwMRROjUko9gIYMGcK1a9d0eEYCNDEqpdQDZs+ePcydO5fevXtTsWJFV4eT4WhiVEqpB0hsbCw9e/akSJEiDB8+3NXhZEja+UYppR4gCxYsYM+ePSxevJjcuXO7OpwMSa8YlVLqAfHXX38xcOBAGjRoQLt27VwdToaliVEppR4Qw4cP588//+TTTz/VDjdJ0MSolFIPgJ9//pkZM2bw1ltv8fjjj7s6nAwtQyRGY4yvMeYrY0ykMSbWGDPFGPOJMWaNq2NzljFmrjFGjDFT7mMb7xpjDhhjMsTPRSmVOYgIvXr1In/+/IwePdrV4WR4GeUEHACUBzoC9YCVQDdglAtjcpoxJgfwCiBAO2PMvXZqmgkUxvoclFIqVQQGBrJt2zYmTJhAvnz5XB1OhufyxGiMKQo8A0wUkXUishN4FfhZRNL9IaDGmBPGmJEprPYykBuYhJXYnruXfYvIDWAR0O9e6iulVHxXrlyhf//+PPHEE3Tu3NnV4bgFlyZGY8xXQARggIW2psjxgD+wNF7ZssaYaGPMqHjLPzfGXDXGJPgw2HTSETgODAMigQ7xC6Qg/mVARWNMvbQOWimV+Y0ePZrz58/z2WefkSWLy6+F3IKrP6VhwFzgGlDXNq0H8gI/OBYUkSO2su8aYwoCGGOGA68DL7vi6tIWQ3GgCbBERKKxEtsLxpg47RUpiH8/cIV7vOpUSqk7DoeEsHHaNN544w2eeOIJV4fjNlyaGEXkV6wkeEBEdtqaUetg3as7kECVUUBW4H1jzBvACKC9iGy8l/0bSzbHybYqS7zlWZPYTHusz3GJbX4RkB1ocy/xi0gs1rHXceYYwi+FE7A/AIDomGj8AvxYcsAK5Xr0dfwC/Ag6GATA5ZuX8Qvw46uwrwC4eP0ifgF+rAm3+jidu3YOvwA/1h9ZD8Dpy6fxC/Bj4zErvGN/HcMvwI+tJ7Za+74Yjl+AH9tPbwfg4IWD+AX4sefMHgD2n9uPX4Af+8/tB2DPmT34Bfhx8MJBALaf3o5fgB/hF8MB2HpiK34Bfhz76xgAG49txC/Aj9OXTwOw/sh6/AL8OHftHABrwtfgF+DHxesXAfgq7Cv8Avy4fPMyAEEHg/AL8ON69HUAlhxYgl+AH9Ex0QAE7A/AL8DP/lnO+WkOTRY1sc/P2DOD/wb+1z4/fed0XvjiBfv85O2TafllS/v8hG0TeHXFq/b5MVvH4P+Vv31++JbhdF79b1PWoI2DeHPNm/b5ft/1o+fanvb5Puv70Gd9H/t8z7U96ffdv63sb655k0EbB9nnO6/uzPAt/z7JxP8rf8ZsHWOff3XFq0zYNsE+3/LLlkzePtk+/8IXLzB953T7/H8D/8uMPTPs800WNWHOT3Ps834Bfvq7l5F/9758FV54gXWxsYzTJ9ykiKuvGAGqAfsc5osDV0QkKn5BETkHTAPeBmYB74jIl45ljDFljDHbjDGHjTH7kmlibQxEx5tKYV3JOi7blMQ2OgC7ROSwLcZQ4DcSaE51Jn6bSKzPIUHGmDeNMaHGmNDo6OgkQlNKPZAklusbNnCo2FVm9f4vhUqWdHVE7kVEXDYBuYBY4A2HZbOA80nU6YR1RflDIus3AF1s/28KhGN7IXMi+68VbzoLzI63rHwi9Z+wxdIz3vLBtuXlUhq/rUwgcNSZz7BmzZqilFKObvTsKQLy4UMPSUxMjKvDyZCAUEnkvOrqK8ZqWB1vHK8YLwEJ9ic2xjyNlTh3APWNMY/HW18IqwlyIYCIbLCtqpnQ9kTkqoiEOk5AFHA23vLwROLviHVFGRRv+RKs5BfnqjG5+B3kBy4msk4ppRI3cyZen33GJ8Djy+Zy4foFV0fkdlydGKtjJZZfHZb9BngYY0o4FjTG1ABWYXVg8QNOAePibe8hrKTm2L540rY8VRljPLGGlXwrInGSmIicArYC7Y3tuUtOxn/Hw1hXukop5bzgYKRXL9YCx3r14oPwD+Lce1TOcfXbNaoDh0TklsOy723/Pgn8AdZQB+Bb4DvgbRGJtQ17mG+MaSQi35O4tHog4P8BBYDTxpiXElh/DCsB+hljTuNk/MaYvEA5YHIC21RKqYQdPIi88gqHPT3pkzs3P33wAdsvbHd1VG4pIyTG/Y4LROSEMWY30Bz4yvYAgO+AMKCdWL02wer9OQCYgPW0HLCuwoobYzwcrhpL2ZantjtPp+lpmxLzPlaicyZ+gGZYzblfp2q0SqnM69w5aNaMf4zhmRs3+GjBAnLnzs1zuXXU171waWIUkWqJrPocmG6M6SlWT85HEqgbA1SItyzSllQ7AXOMMU2xrhh/SkFMpZ0s96Kz20yk/l3x2/gDy0Xk0v1sXyn1gLh+HV54gdjISJplyUKl//yH1q1bA9iHm5TMo71SU8LV9xgTsxg4A7x1D3W7A52NMYeBD7Gu0iQ1g0srxphqwFO4yTNilVIuFhsLHTpAaChTatZk1+3bfPbZZ/ZXSrX/uj3tv27v4iDdj6ubUhMkIjHGmNeBGvdQ93fiNk26k6JAZ7GekqOUUkkbNAhWriS8Wzf6z5rFyJEjKVu2rH310EZDXRic+zJucjGlElGrVi0JDXXJ0/CUUq40Zw68+Sa333yTips3I8Avv/yCl5eXqyNzC8aYn0QkwQfAZMgrRqWUUknYuBF69IDnnmNckSL8fuQI33333V1J8c4j7h7Jd1c3DZUETYxKKeVODh2CVq2gYkWOjhvHuLp1adOmDU2bNr2r6OurXwcgpFNIOgfp3jQxKqWUuzh/Hpo1gxw5kDVr6PbGG3h6ejJlypQEi4/y035890ITo1JKuYMbN+DFF63k+P33LN66lU2bNjFjxgyKF0/4nQONSzdO5yAzB02MSimV0cXGQseOsHs3rFxJZKlSvPfcc9SrV49u3bolWu3Oa7XKFyyfXpFmCpoYlVIqoxs6FJYvhw8/hJdf5l1/f65cucLs2bPJkiXx4ejdvrGSpt5jTBlNjEoplZEtWADjx8Obb0LfvgQHBxMYGMiwYcOoVKlSklXHPZPYewpUUnQco5vTcYxKZWKbN8Ozz8JTT8HatfwTFUXlypXJnj07+/fv1zGL90HHMSqllLv57Tdo2RLKl7eaUT08GDl4MCdOnGDr1q1OJcWDFw4CULlw5bSONlPRxKiUUhlNZCQ8/zx4esI330CePOzdu5cpU6bQtWtXGjVq5NRmeq3rBeg9xpTSxKiUUhnJzZvw0ksQEQEhIVC6NLdv36ZLly4ULlyYSZMmOb2pD5t+mHZxZmKaGJVSKqOIjYXOnWH7dqv5tHZtAKZPn86+ffv48ssvyZs3r9Obe8L3ibSKNFPLqK+dUkqpB8+IEbBsGUyYYD32DTh+/DjDhw+nefPmtLItc9b+c/vZf25/8gVVHHrFqJRSGcHChTB2LLzxBgwYAICI0KNHD7JkyRLnPYvO6rO+D6D3GFNKE6NSSrlaSAh07QrPPAOffw62BLho0SKCg4P5+OOPKVmyZIo3O+25aakc6INBxzG6OR3HqJSbCw+HunWhaFHr3qLtHmJERAQVK1akUqVKfP/990k+4UalXFLjGPWTVkopV7l40XpbRrZssHatPSneaUK9efMm8+bNu+ekuOfMHvac2ZOaET8QtClVKaVc4dYtePll+OMP2LIFHn7YvurLL79k9erVTJo0ifLl7/0B4P039Af0HmNKaWJUSqn0JgKvvw7btkFQkNWUahMZGUmvXr144oknePfdd+9rN58+/+n9RvpA0sSolFLpbdQoWLoUPvgAWreOs+rtt9/m8uXLzJ8/n2zZ7u8UrY+Cuzd6j1EppdLTkiVWYuzUCQYNirPq66+/JigoiGHDhlG58v0nte2nt7P99Pb73s6DRnulujntlaqUG/nhB2jSBOrXh/XrrWeh2vz5559UqlSJIkWKsGfPHjw8PO57d34BfoDeY0yIvl1DKaVc7fffrWegPvwwrFwZJykCvPfee0RGRrJ27dpUSYoAs/5vVqps50GjiVEppdLapUvWsIwsWaxhGfnyxVn97bffsnDhQoYMGUKNGjVSbbflC957j9YHmSZGpZRKS7duQYsWcPKk9eLhMmXirL58+TJvvvkmFSpUYNiwYam6660ntgLQuHTjVN1uZqeJUSml0oqI9ai377+3eqHWr39Xkf79+3P27Fm2b99O9uzZU3X3I0JGAHqPMaU0MSqlVFoZOxYWL4bRo6Ft27tWf/vtt8yZM4f+/ftT2/aKqdQ0/8X5qb7NB4H2SnVz2itVqQzqiy/gtdegQwcICLA/GPyOP//8k8qVK5M/f35CQ0Px8vJyTZwPKO2VqpRS6enHH61xio0awezZdyVFgHfeeYfIyEjWrFmTZklx47GNADR5pEmabD+z0sSolFKp6ehRa1hGqVLw1VeQwH3DlStXEhgYyMiRI6lZs2aahTL2+7GAJsaU0qZUN6dNqUplIH/9ZT33NDISdu6ERx+9q8j58+epXLkypUqVYseOHak2ZjEhpy+fBqBknpS/yzGz06ZUpZRKa1FR1rCM48dh48YEk6KI0L17d65evcrChQvTNCmCJsR7pYlRKaXulwh06wYhIVYv1IYNEyy2ZMkSVq1axYcffkilSpXSPKz1R9YD8FzZ59J8X5mJJkallLpf48dbPU9HjAB//wSL/PHHH7z99ts0aNDgvl8n5awJ2yYAmhhTShOjUkrdj6AgGDIE2rWzEmMCRIQ33niD6OhoAgICyJo1a7qEtqzVsnTZT2ajiVEppe7Vjh3QsSM0aADz5iU4LANg1qxZfPfdd8yYMYMy8R4Jl5aK5iyabvvKTLRXqpvTXqlKucixY1CnDuTObfVALVgwwWJHjhyhWrVq1KtXj+DgYEwiyTMtrAlfA0Dz8s3TbZ/uQnulKqVUavr7b+ttGbdvw7p1iSbF27dv4+/vj4eHB/Pnz0/XpAjw0Y6PAE2MKZWixGiMqQM8B9QBigM5gItAOLAVWCUif6V2kEoplWFER0OrVtZA/g0boFy5RIuOHTuWXbt2ERQURIkSJdIxSMuK1ivSfZ+ZQRZnChljOhpjfgG2A30Ab+B3YBfwF1AbmAucMcYEGGMeTqN4lVLKdUSgRw/YtAnmzoXGib/OaceOHYwZM4b27dvTunXrdAzyXwW9C1LQO+GrWZW4ZK8YjTE/A4WBRUAHYL8kcGPSGJMH+D+gHfCrMaaziASlcrxKKeU6kyZZnWyGDrUeDp6Iq1ev4u/vz0MPPcSnn36ajgHG9VXYVwC0qNDCZTG4I2eaUhcAM0XkZlKFROQyEAgEGmMeB7Q7lFIq81ixAgYOhFdftV4jlYTevXtz4sQJtm7dSu7cudMpwLt9vOtjQBNjSiWbGEVkWko3KiI/Az/fU0RKKZXR7NoF7dtDvXqwYEGiwzLAekD4ggULGDJkCA0aNEjHIO+2+tXVLt2/u0rRcA1jTA0R2ZuG8agU0uEaSqWxEyegdm3ImdMallGoUKJFz5w5Q9WqVXnkkUfYvn17mj8LVd27pIZrONX5xsEWY8xTqRCTUkplfJcvW8MyoqJg7dokk2JsbCydO3fm5s2bBAYGZoikGHQwiKCD2tUjpVKaGJcC64wxLeOvMMY0MMZsS52wlFLKxaKj4ZVX4PBhWLkSHnssyeIff/wxGzZsYMqUKZRLYghHevo89HM+D/3c1WG4nRSNYxSRHsaYCGCZMeZtEZlpjKkCjAOaAWFpEaRSSqUrEejVyxqnOH8+PP10ksV/+eUXBg4cSPPmzXnzzTfTKcjkrWu3ztUhuKUUP/lGREYbY84Anxtj2gL1gdPA61hDOpRSyr199BHMng2DBkHnzkkWvX79Om3atCFv3rzMnTs33Z9ukxRvD29Xh+CWUpwYjTH5gXJADNAQa9C/n4jcTuXYlFIq/X39NQwYYDWjjh2bbPE+ffoQFhbGd999R+HChdMhQOctObAEAP+qCb8KSyUsRfcYjTEjgGNAT+AjrKvEWsCU1A9NKaXSWWio9fqo2rVh4ULIkvQpcvny5cyZM4eBAwfStGnTdArSeXP3zmXu3rmuDsPtpHS4RhTWo99Gich527Knga+B9YC/iESnRaAqYTpcQ6lUcuqUlRC9vKxxi8lc/R0/fpzq1atToUIFvv/++wzRCzW+6BjrdOyRNePF5mqp+XaNCiJy1HGBiGy2DeFYh5Ucn7m3MJVSykWuXIH/+z+4ccN6DmoySTE6Opq2bdsiIixdujRDJkXQhHivUtor9Wgiy/caYxoAwakSlVJKpZfbt6FNGzh0CNavh4oVk60yfPhw+1szHn44474zIWB/AACdqnVyaRzuJtXexygiR4wx9VJre0opleZEoHdvKyHOng1NmiRbZePGjUycOJEuXbq47K0ZztLEeG+SvcdojFkNjBSRfU5t0Bgv4C3guojMvP8QVVL0HqNS92HaNHj3Xejf33pzRjLOnz9PtWrVyJcvH6GhoXh763AId3W/9xhPAjuNMfuxnnzzA3DAcXiGMaY48CTQHGgBnMHqsaqUUhnT//4H770HLVrAhAnJFo+NjaVjx478/ffffPfdd5oUMzFnEmMU8BTQFhgB5AHEGHMFuAXkAzwAA+zGepHxYhGJTZOIlVLqfu3dC23bQq1asHhxssMyAKZMmUJwcDAzZsygSpUq6RDk/Zvz0xwAutbs6uJI3IszibEP8KWIvG2MuYrV87QuUAzwAi4BvwHfi8jJNItUKaVSwx9/QPPmULCgddXoxJXfjz/+yMCBA2nRogXdu3dPhyBTR9Cv1gPENTGmjDP3GC8A7UUk2BgTA9QVkd3pEp1Klt5jVCoFrl6Fhg3h2DHYvh0qV062SmRkJNWrV8fLy4uffvqJPHnypEOgKq3d7z3GbcBkY0whrOZS558IoJRSGcXt21bz6cGD1iuknEiKMTEx+Pv7c/HiRXbs2KFJ8QHhTGLsBSy0TQJsNMYcAPY5TL/qE2+UUhnae+9ZCfHzz+HZZ52qMm7cOL777jtmzZpF9erV0zjA1DdjzwwA3nriLRdH4l6SveMsImdFpCngi3XFGAREAM9hPR7uJ+CqMWavMWZeWgarlFL35JNPrOm998DJe4SbNm1ixIgRtGvXjq5d3fMe3ZrDa1hzeI2rw3A7KX1W6kpgqIiE2eZzAtWA6ncmEXG/r1VuTO8xKpWMtWvhhResDjcrV0LWrMlWOXv2LNWrV6dAgQLs3r2bnDlzpkOgKj2l2rNSRaRlvPlrWPcgt917eEoplUb277ce91a9OgQGOpUUb9++Tdu2bbl27RpbtmzRpPgASrVHwimlVIZy5oz1YPB8+axhGT4+TlUbNmwY33//PYsXL6aiE89Nzcim75wOQO86vV0ciXtJ0fsYlVLKLVy7ZjWdXr5sNaUWL+5UtW+++YYJEybQtWtX/P3d/+W+m45vYtPxTa4Ow+2k6B6jynj0HqNS8cTEwMsvWwlxzRp4/nmnqh07doxatWpRqlQpduzYgZeXVxoHqlwpNd/HqJRSGVu/flZC/PRTp5Pi9evXadGiBSLCihUrNCk+4DQxKqUyjxkzrDdm9O4NPXs6VUVE6Nq1KwcOHGDt2rWUKVMmjYNMP5O3TwagX71+Lo7EvWhiVEplDt9+C2+/bd1b/Ogjp6t98sknLF26lDFjxvDf//43DQNMfzv+2OHqENyS3mN0c3qPUSngwAGoXx8efRS+/x6cHGLx/fff88wzz/D888/z9ddfk8WJt2yozCGpe4z6W6CUcm8REdawjDx5rHuLTibFM2fO0Lp1ax5++GEWLVqkSVHZaVOqUsp9/fOP1XT655+wbRv4+jpV7datW7Rq1Ypr166xadOmTPtw8AnbrBcwD2ww0MWRuBdNjEop9xQTA/7+sG8frF4N1ao5XbVPnz7s3LmT5cuXU6lSpTQM0rX2n9vv6hDckiZGpZR7ev99WLUKpk+3mlKdNH/+fGbOnMmAAQNo1apVGgboestaLXN1CG5JG9WVUu5n1iyr52mvXvDOO05X2717N2+99RbPPPMMH3zwQRoGqNyZJkallHsJDrbGKD7/PEyd6nS1s2fP8tJLL1GsWDGWLVtGtmyZv8FszNYxjNk6xtVhuJ3M/5uhlMo8Dh6EV16BypVh2TJwMrnduHGDl156iatXrxIcHEzBggXTONCMIfxSuKtDcEuaGJVS7uHcOWjWzBqO8c03kCuXU9XuPNlmz549rFq1iipVqqRxoBnHkhZLXB2CW9LEqJTK+K5ft142fPGiNYC/RAmnq06aNInAwEDGjh3Liy++mIZBqsxCE6NSKmOLjYUOHSA0FL7+GmrWdLrqN998w6BBg2jTpg2DBw9OwyAzpuFbhgMw+qnRLo7EvWhiVEplbIMGwcqVMGUKpOCK79dff+W1116jRo0azJ8/H2NMGgaZMZ2+ctrVIbglTYxKqYxrzhyYNAl69IA+fZyudunSJV544QW8vb1ZtWoV3t7eaRhkxrXgxQWuDsEtaWJUSmVMGzdaCfG55+Djj8HJK77o6Ghat27NH3/8wdatWymRgvuRSoEmRqVURnToELRqBRUrQlCQ08MyRISePXuyefNmAgICqFOnThoHmrEN2jgIgPFNxrs4EveiiVEplbGcP28Ny8iRwxqWkTu301UnT57MnDlzGDx4MB07dkzDIN3DpRuXXB2CW9LEqJTKOG7csDrYnD9vDct46CGnq65cuZIBAwbQpk0bxozRp70AzG4+29UhuCVNjEqpjCE2Fjp2hN27rV6otRJ8h2yCdu/ejb+/P3Xr1mXBggX6bkV1XzQxKqUyhqFDYfly+PBDePllp6udOHGC5s2bU6xYMVavXk2OHDnSMEj30u+7fgBM/s9kF0fiXjQxKqVcb8ECGD8e3nwT+vZ1utrly5dp1qwZt27dIiQkhEKFCqVhkO7nRvQNV4fgljQxKqVca/NmKyE2bQqffpqiYRmvvPIKhw8fJjg4mAoVKqRxoO7ns2afuToEt6SJUSnlOr/9Bi1bQvnyVjOqh4dT1e4My9iwYQPz58/n6aefTuNA1YNE71ArpVwjMtJ6p6KnpzUsI08ep6uOGzeOOXPmMGjQIDp37pyGQbq3Puv70Ge9808MUha9YlRKpb+bN+GllyAiAkJCoHRpp6sGBAQwdOhQ/P39GTt2bJqFqB5cmhiVUukrNhY6d4bt263m09q1na767bff0qVLF5o2bcq8efN0WEYy2t3FLgAAIABJREFUpj03zdUhuCX9rVJKpa8RI2DZMpgwwXrsm5P27NlDq1atqFq1KitXrsTT0zMNg1QPMk2MSqn0s3AhjB0Lb7wBAwY4Xe3IkSM0a9aMIkWKsG7dOnLlypWGQWYePdf2pOfanq4Ow+1oU6pSKn2EhEDXrvDMM/D5504Pyzh//jzPPvsssbGxrF+/nqJFi6ZtnJlIDg992MG90MSolEp74eHQogWULQsrVjg9LOPatWs0a9aMiIgItmzZQrly5dI40MxFn3hzbzQxKqXS1sWL1tsysmWDtWshb16nqkVFRdGqVSv279/PqlWrqJ2CTjpK3Q9NjEqptHPrlvXc0z/+gC1b4OGHnaoWExNDhw4dCA4OZs6cOfzf//1fGgeaOb255k1A37KRUpoYlVJpQwRefx22bbNeNly3rpPVhLfeeougoCAmTpxIly5d0jjQzKtAjgKuDsEtaWJUSqWNUaNg6VL44ANo3drpaoMGDWL27NkMHDiQASnouaruNr7JeFeH4JZ0uIZSKvUtWWIlxk6dYNAgp6tNnDiRiRMn0r17d8aNG5d28SmVBE2MSqnU9cMP1jjFp56CWbOcHpYxa9YsBg4cyKuvvsqnn36KcbIewMaN1iiQf/6516Azp86rO9N5tT5LNqU0MSqlUs+RI9YzUB9+GFautB4Q7oRly5bRo0cPmjVrxqJFi8iaNatT9f75B3r2tN5Ydfas9ehV9a+SuUtSMndJV4fhdvQeo1Iqdfz5pzUsI0sWa1hGvnxOVVu3bh3t27enYcOGLF++HA8nxzj++CN07AjHjsG771q3MnPoePY4Rj812tUhuCW9YlRK3b+oKGsA/4kTsGoVlCnjVLWNGzfSokULqlatyv/+9z9yOJHZbt6E99+Hhg0hJsYaBTJliiZFlXr0ilEpdX9ErEe9bd1q9UKtX9+palu3buWFF16gXLlyBAcHk8eJ9zHu3QsdOsCvv1q7/Ogj0MemJs7/K38AlrRY4uJI3IsmRqXU/fngA1i0CEaPhrZtnary448/0qxZM0qXLs3GjRspWLBgkuWjo2H8eBgzBgoVslpqn38+NYLP3MoXKO/qENySERFXx6DuQ61atSQ0NNTVYagH1RdfwGuvWZdxAQFO9UDdtWsXTZs25f/bu/PwqIp0j+PfIiwDzICsElCUZUC8rMom3mEaDVd2RkEBWQRUUARFBxFFwSgio6IiCiKLQUBZNCiLBNm9CgIJhOXqBFlFZN8hCoGu+8cJmQaztgmnT/L7PE8/2qdPn34DRb+pOlX1hoeHs3LlSsLDw9M9/4cfnMvHxjp59913oWTJbIpf8ixjTJy1tn5qr+keo4gEZ/Vqp+Bw06bwwQeZSopxcXHcddddlC1bluXLl6ebFP1+595hvXqwaxfMnu2M1CopSk7TUKqIZN2OHdC+PVSsCNHRUKhQhm/ZtGkTzZs3p0SJEixfvpwKFSqkee7OnU7O/fpraNvWybuqNpV1nT/tDMDMjjNdjsRblBhFJGuOH3eWZfj9zs2+Uhnvx7l161YiIiIoWrQoy5cvp2LFiqmeZy1MnAhPPQVhYfDhh86SjCys9ZcAdcvVdTsET1JiFJHMO38eOnRwunRLl8Jf/5rhWzZt2kRERAQFCxZkxYoVVEqjwsa+ffDQQxATA3fc4STFNPKnZNKQ/x7idgiepHuMIpI51sIjjzgLB6dMce4tZiAuLo5mzZpRuHBhVq1aRdWqVVO97McfQ82azoqPsWNhyRIlRXGPEqOIZM6oUU43btgw6NYtw9PXrl3LnXfeSbFixdJMiocPw733QteuUKMGbNoE/fs7m+fIH9dhdgc6zO7gdhieo6FUEcnY7Nnw3HPO0owXX8zw9G+//ZaWLVtSpkwZVqxYkeo9xS++gD594MQJJ+cOGuTcV5Tsc9t1mauBKZdTYhSR9H33nbOQ8PbbYfLkDGfCfP3117Rq1YoKFSqkOvv0xAkYOBCmToU6dZxblbVq5eQPkHcNajLI7RA8SQMWIpK2XbugXTu47jpnD9Q//Snd05ctW0aLFi2oWLEiK1eu/F1SvJQEp0+H55+HdeuUFCX0KDGKSOpOnHCWZVy44CzLyGDbtoULF9KmTRuqVKnyux1tzp517h02bw5Fizp7A7z8cqarUkmQ2n3SjnaftHM7DM/RUKqI/F5SEnTs6NRX/OorqJ7+npuffPIJPXr0oHbt2ixevPiyvU9Xr3bWIm7f7gyhjhypShhXy52V7nQ7BE9SYhSRy1kL/frBsmXO/qc+X7qnjx8/nscee4ymTZsyb948ihUrBsC5c84E1jfecJZerFiR4aUkmz3R+Am3Q/AkDaWKyOVefx0mTYKhQ52uXhqstYwcOZJ+/frRpk0bFi1alJIUN26E+vXhtdegd2/YvFlJUbxDiVFE/uOzz5wqwJ06OWWk0mCtZfDgwQwdOpRu3brx2WefUbhwYS5ccO4dNmwIR486tyYnTlTNRLe0nNGSljNauh2G52goVUQc69Y5C/dvu80ZQk1jlf3Fixfp27cvkydPpn///owZM4Z8+fLxww9OB3P9epWHChVtq7V1OwRPUmIUEdizx1mWER7urLxPY1nGr7/+SteuXZk7dy7Dhg3jxRdfxFrDW2856/+LFnX2Arj33qscv6SqX4N+bofgSUqMInndyZPOsozffnNmyJQpk+ppR48epV27dqxZs4YxY8bw+OOPs2sX9Oyp8lCSuygxiuRlSUlw332QkACLFzsblqZi9+7dtGjRgt27dzN79mw6dOiYUh4qXz6VhwpVER9FALC0x1KXI/EWJUaRvMpaGDDAWac4ebJT6ykVGzdupFWrVpw7d44lS5ZQpcrfaN0aFi1SeahQ1+m/OrkdgicpMYrkVW++CRMmwJAhzpqKVCxevJiOHTtSsmRJli1bTnx8Ddq3d0Zdx451ljuqEkboevjWh90OwZPUpEXyos8/h6efdmbJvPJKqqdMnTo1ZYu3BQvWMnx4Dbp2dTbBiY9XeSjJvdRjFMlrYmOd8lENGzolLq7Ibn6/n+HDhzNixAgiIiJ48MEviIgowvHj8OqrTj5VeShv8EX5AFjZc6WrcXiNEqNIXvLTT8700WuvdZZlXLFpaWJiIg888ACffvopPXoMAN6iS5cw6tSBJUugdm13wpbg9Kzb0+0QPEmJUSSvOHUK2rSBxERnH9Rrr73s5V9++YX27dsTFxdH376z+fLLjuzbZxg61NnzVJUwvEeJMThKjCJ5wYULzjZv33/vTCe9+ebLXt6wYQPt2rXj+PHztGy5gwkTKlG9ulMZo1Ejl2KWPyzpYhIABcIKuByJt+jWuUhuZy088QTExMD48U5RxABz587lb3/7G0lJDShVai9fflmJgQNhwwYlRa9rPq05zac1z/hEuYx6jCK53ZgxMG6cM2vm4f9M379UHeP5518iPPwDDh7swfXXG5WHykUeuuUht0PwJCVGkdxs3jxne5p77oFRo1IOnz59mgceeIC5c3dSvPh29u+/nocegtGjIblylOQC3Wp3czsET1JiFMmtNmxwylzUrw/TpqUsy0hISOAf/+hIQsLd5Ms3h8KF8zFjhrNdquQuiUmJABQpUMTlSLxFiVEkN/r5Z2dZRunSTq+xiPPFOH/+fLp0ieTcualYewudOjnloUqVcjleyRGtZrQCtI4xq5QYRXKb06edZRmnTzvTSsuVw+/3Exn5Mi+9dBxjvuGaawrw/vvO/uGSez1a/1G3Q/AkJUaR3OTCBWf4dOtWWLgQatbk+PHj3HvvYJYt6wr4aNHiIlOmhKk8VB7QqaY2EQ+GEqNIbvLUU05CHD8e7rqLtWvX0apVNMeOvUmhQoUYN87Sq1eYykPlESd/OwlA8T8VdzkSb9E6RpHcYuxY5/HUU9i+fYmMnMhttx3h2LFR3HqrJSGhIL17GyXFPKT9zPa0n9ne7TA8Rz1Gkdxg4UIYOBDat+f4s8/xPw3fIja2J2FhRRkx4iyDBxdTJYw86PFGj7sdgicpMYp4XXy8s91bvXose/AZ2t/4HWfPPsUNN/zC4sUlqF69kNsRikvuqXGP2yF4kn6HzARjTAVjTLQx5rAxxm+MedMYM9YYMz+bP+dJY8xmY4z+XiRz9u2DNm2wJUrwaJWBRLSrxNmzzXnkkT1s316e6tU1bpqXHUk8wpHEI26H4TnqMWZOFFAeeAA4BoQBK4Am2fw57wPPJH/Oh9l8bcltzpyBtm05ftxPoyIf8OPsNhQvvot58xJp2vQGt6OTENBxdkdA6xizSokxA8aYcsCdQE9r7ZfJx8YCm6y1sdn5WdbaX40xHwGDUGKU9Fy8CPffz9L4UrTncxITK9C6dTyffVaHQoXUSxTHP2/7p9sheJKG7NJhjIkG9gMGmGqMscaYV4FuwMdXnFvVGJNkjIm84vh4Y8xpY0z9TH7sTOBmY0x290YlFznRbxAD5jenuV2Cv5Blzpx9LFhQV0lRLtO2elvaVm/rdhieo8SYvheAScAZ4LbkRwxwDfC/gSdaa7cnn/ukMaY0gDFmGNAbuDsLvct44BTQIjt+AMl95nR5h4YfPMq7DKBBg9UcOFCejh0ruh2WhKADZw5w4MwBt8PwHCXGdFhr/w8nCW621n5nrf0OaAxYYHMqb4nEuf/4jDHmQWA40N1auzQLn+lPvnbjzJyfcDSBqPgowClK6ovyMX3zdMDZQNgX5WPW1lmAs9jXF+Uj+odowLkx74vyMT/BmUN04MwBfFE+YrbHALD35F58UT6W7nTC33l8J74oH6t2r3I++0gCvigfq/euBmDroa34onys37cegPgD8fiifMQfiAdg/b71+KJ8bD20FYDVe1fji/KRcCQBgFW7V+GL8rHz+E4Alu5cii/Kx96TewGI2R6DL8qX8g99fsJ8fFG+lMkF0T9E44vypSxqnrV1Fr4oX8pGytM3T8cX5Usp3hoVH4UvypfyZzkxbiIRH0WkPB+3fhwtZ7RMeT7muzG0+6RdyvM3Vr9Bh9kdUp6P+mYUnT/tnPL85VUv0y36P9UNhq0YRq8veqU8f3bps/SZ3yfl+aCvBvHYwsdSng+MGcjAmIEpzx+K7kvF+++m88zHOGUK8dbojaxb14TixVWEVlLX+dPOl7VJyRzdY8xYXWBxwPPywClr7fkrT7TWHjDGvA38E+fP9nFr7ezAc4wxVYCpQFngLPBwKr3Jw0C1tAIyxvQB+gAUqqCp+HlBVFQ8H0afxn+uCrXrj+Dp1/9KN9/9boclIW7Ifw9xOwRPMtZat2MIWcaYvwAncZLX5ORjE4B/WGuvTeM9PXEmznxjrf1bKq8vAWZZaycZY5oD7wI32YC/CGPMDKCxtbZKRjHWr1/fxsZm6xwgCSGnTiXSqtVKvv02grCwE7z00n6ee66O22GJeJ4xJs5am+rcDw2lpq8uzsSbjQHHjgIlUjvZGHMHMAFYA9xujKlzxetlcIZIpwJYa5ckv3TrFZcqCWjxUR43bdp6rr32R779thVVq25i+/bCPPdcHfae3JsyvCySHrWV4Cgxpq8ekAT8X8CxfwMFjDHXBZ5ojLkF+BxnAo4P+AkYecX1KgK/WGuTAo7tST4eqBKQ8EeDF286fPgojRt/Qo8eNTl/viIvvLCVH39swI03/gWA7nO7031ud5ejFC9QWwmO7jGmrx7wvbX2XMCxr5P/2xD4GZylGsAi4CtggLXWn7xsY4oxpqm19mvSdtn8emPMNTj3F9/Ipp9BPMJay5tvRjNkSDkuXOhClSr/ZvHiilSpUvOy855v+rxLEYrXqK0ER4kxffVwlk+ksNbuNsasA9oC0ckbAHwF/AB0TZ5VCvARMBgYxX92yPkJKG+MKRDQa7wh+fglrYHzwNwc+HkkRCUkbKN9+3kkJPQlLCwfkZE/88ILN6VaCSOicsTvD4qkQm0lOEqM6bDW1k3jpfHAGGPMY9baA0DlVN57EahxxbHDyUm1JzAxefKNAeICTusGzLHWHs2GH0FC3NmzZxk69D3eeacW1g6ievV9LFoUTqVKRdN8z6XlLJVL/K7ZiVxGbSU4SozBmYbTG+xH1oc8H8HZRedpIBGnl2kBjDF1gWZAzXTeL7mAtZZPPvmE/v2/4fjxEYSFFSUy8hTPPlshw/JQvb/oDWj/S8mY2kpwlBiDYK29aIzpDdwSxHt/JO3Nx8sBvZJ30ZFcKi4ujn79hrFuXQ9gHDfffJq5cwtRrVrm1qRG+iIzPkkEtZVgaR2jx2kdo3ccPHiQoUOHMnnyIfLlm4QxpYiMNDzzTD7y61dUkasqvXWM+ucoksNOnz7Nm2++yeuvTyAxcSTQkxo1LjJjRhh1glirf2kLveqlq2dvoJLrqK0ER4lRJIecP3+eiRMn8tJLL3Ho0H9RpMgmjCnNs8/C8OFhFCwY3HX7LugL6L6RZExtJThKjCLZzO/3M2fOHIYOHcqOHb9QocI0oAPXXQcffQSNGv2x64+888p9I0RSp7YSHCVGkWzi9/v54osviIyMZNOmTVSp0pXy5Sewb19RHn8cXn0VihT545/T5HqV6pTMUVsJjraEE/mD/H4/0dHR1KtXj3vuuYczZ5Jo23Yzu3ZNI3/+oixfDmPGZE9SBKe816XSXSLpUVsJjhKjSJD8fj+fffYZdevWpUOHDvz222+MGDGfIkW2Mn9+LXr1MmzZAs2aZe/n9v+yP/2/7J+9F5VcSW0lOBpKFcmi3377jWnTpjF69GgSEhKoVq0aU6fOYM+ezkRG5qNUKZg/H9q0yZnPf7356zlzYcl11FaCo8QokknHjh1j/PjxjB07loMHD1KvXj0+/vhj6tS5j169wli3Djp1gvfeg1Klci6OBhUa5NzFJVdRWwmOhlJFMvDjjz/yxBNPULFiRZ5//nnq1q3L0qVLWb8+jkOHunDrrWFs3w4zZzqPnEyKAPEH4ok/EJ/xiZLnqa0ERz1GkVRcvHiRhQsXMm7cOBYvXkz+/Pnp0qULgwYNonbt2uzeDRERsHIltG4NEydCePjViW1gzEBAa9MkY2orwVFiFAlw+PBhJk+ezPvvv8+ePXsoX748kZGRPPzww4SHh2MtTJ4MTz7pnD95MvTqRarloXLK2y3evnofJp6mthIcJUbJ8y5cuEBMTAxRUVHMmzePpKQkmjVrxujRo2nXrh0FChQAYP9+ePhhWLgQfD748EO48carH2/dcmlVQxO5nNpKcJQYJc/aunUrUVFRTJ8+nYMHD1KmTBkee+wx+vTpQ40al5XSZOZM6NcPfv3VWZPYvz8ZlofKKev3rQc0sUIyprYSHCVGyVN27drFnDlzmDVrFhs2bCB//vy0bduWnj170rJly5Te4SVHjzoJcfZsZyu3qVOhusv7MT+95GlA940kY2orwVFilFzvUjKcM2cOl0p0NWzYkLfffpv777+fMmXKpPq+BQucodOjR+GVV2DwYEKiPNS7rd51OwTxCLWV4ITAP3OR7OX3+4mNjWXBggUsXLiQDRs2ANCgQQNef/11OnbsyI3p3Bw8dcqZXDNlCtSuDTExBFUeKqfULFvT7RDEI9RWgqPEKLnCiRMnWLp0KQsWLGDRokUcOnSIfPny0bhxY/71r39x7733UqlSpQyvs3y5M8v0559JLg8FhQpdhR8gC1bvXQ1og2jJmNpKcJQYxZNOnTrFN998w4oVK1ixYgUbN27E7/dTokQJWrRoQevWrWnRogWlMrnaPjERhgyBsWOhWjX49lto3DiHf4ggPbfsOUD3jSRjaivBUWKUkGetZdeuXaxdu5a1a9eyZs0a4uLiuHjxIgULFqRx48a88MILRERE0LhxY/Jn8Ubgd9/BAw/Atm1ka3monDKhzQS3QxCPUFsJjhKjhJSLFy+yY8cOtmzZwpYtW4iNjWXdunUcPnwYgMKFC3PrrbcyZMgQmjVrRpMmTShcuHBQn3XuHLz0EowaBdddB8uWwR13ZOdPkzOql3Z5Wqx4htpKcJQYxRUnTpxgx44dKY9t27axZcsWvv/+e3799VcAjDHcdNNNtG7dmkaNGtGoUSNq1qz5uyUVwdi0CXr0gM2boXdveOstKFbsD1/2qli1exUAf7/x7y5HIqFObSU4SoySLfx+P6dPn+bUqVOcPHmSkydPcuTIEfbv3/+7x65duzh27Nhl7y9Xrhw1a9bkkUceoVatWtSqVYubb76ZItk8pnnhArz2Grz4IpQsCfPmQdu22foROW74yuGA7htJxtRWgqPEKJnWuXNn4uPjSUpKuuxx/vx5zpw5g7U21fcZYyhbtizlypUjPDycBg0aUKVKlZRH5cqV+fOf/5zj8SckOPcS166F++6DceNyvhJGTpjSforbIYhHqK0ER4lRMq1SpUpYaylQoMDvHsWKFaN48eIpj2uuuYaSJUsSHh5O2bJlszwhJjv5/fDuu86s08KFne3dOnVyLZw/rHKJym6HIB6hthIcJUbJtFdffdXtELJszx5nXeKKFVe/PFROWbpzKQARlSNcjkRCndpKcJQYJVey1tm55sknnf+fNMmZZHM1y0PllBFfjwD0ZScZU1sJjhKj5Dr790OfPs5ep26Wh8op0+6e5nYI4hFqK8FRYpRcZfZsePRRZyebt9+GAQPcKw+VU64vfr3bIYhHqK0EJ5d9ZUhedfQodO7sTKqpWhXi4+GJJ3JfUgSI2R5DzPYYt8MQD1BbCY56jOJ5CxfCQw+FXnmonDLqm1EAtKjawuVIJNSprQQnF399SG4XWB6qVq3QKw+VU2Z2nOl2COIRaivBUWIUT1qxwlmGsXdv6JaHyinl/lzO7RDEI9RWgpML78BIbpaYCAMHOpt9FyzolIcaOTLvJEWA+QnzmZ8w3+0wxAPUVoKjHqN4xtq1zsbf27Y5s01HjQrt8lA5ZfSa0QC0re6xTV7lqlNbCY4So4S88+chMtJ75aFyyqf3fep2COIRaivBUWKUkLZ5M3Tv7vy3Vy+nPFTx4m5H5a7SRUq7HYJ4hNpKcHSPUULShQvw6qtQvz4cPOiUh5oyRUkRIPqHaKJ/iHY7DPEAtZXgqMcoIcnvh1mz4O67vVseKqe8s/YdAO6pcY/LkUioU1sJjkmrhp54Q/369W1sbKzbYeSIkyfVQ0zNyd9OAlD8T/rDkfSpraTNGBNnra2f2mvqMUrIUlJMnb7kJLPUVoKje4wiHjNr6yxmbZ3ldhjiAWorwVGPUcRjxseOB6BTzU4uRyKhTm0lOEqMIh7zZdcv3Q5BPEJtJThKjCIeU6RAHtzuR4KithIc3WMU8Zjpm6czffN0t8MQD1BbCY56jCIeM2nDJAC61e7mciQS6tRWgqPEKOIxS7ovcTsE8Qi1leAoMYp4TIGwAm6HIB6hthIc3WMU8Zio+Cii4qPcDkM8QG0lOEqMIh6jLzvJLLWV4GivVI8zxhwG9rgdh1x1pYEjbgchnqC2krobrLVlUntBiVHEg4wxsWltgCwSSG0l6zSUKiIiEkCJUUREJIASo4g3feB2AOIZaitZpHuMIiIiAdRjFBERCaDEKCIiEkCJUUREJIASo0gIMsZUNcYkGWMirzg+3hhz2hijdWmSI4wxk4wx1hjzptuxuEWTb0RClDFmPNAVqGytPWKMGQYMBVpba5e6G53kRsaYwsAB4C/AYaCCtfaCu1FdfeoxioSuSCAMeMYY8yAwHOiupCgZMcbsNsa8GMRb7waKAa8BZYEW2RmXVygxioQoa+0B4G1gADABeNxaO/vS68aYKsaYb4wx24wxGzW86g0hPkz+ALALeAGnx9jjyhNCPP5socQoEtp+BAoBa6y1713x2vtAlLW2GjAYmGGMMVc7QMkaa+12YBLwpDGmNEDyMHlv4G5rbawbcRljygMRwHRrbRIwE2hnjCkReF6oxp+dlBhFQpQx5g6cnuIa4HZjTJ2A18oAjYGpANbaS6Xab73acUpQsm2Y3DjyBz6SX8p3xfGwDC7VHScnTE9+/hHOL2WdcjL+UKTEKBKCjDG3AJ/j/GbuA34CRgacUhH4Jfk3+0v2JB+XEJfNw+R/B5KueNyAMxwaeGxZBmH1ANZaa7clxxgL/JtUhlMzit/r8md8iohcTcaYqsAi4CtggLXWn3w/Z4oxpqm19uu03nrVgpTscGmY/Jt0hsknGWOa4wyT32RTX0YQBzS44tg8YAGX75N6Oq1AjDENgJuB/le8NA14xRhT7VLCzGT8nqblGiIhxBhTDliN00O8y1p7Lvl4GLAVOG6tbZI8lLoTKHmp12iM2Qbcnxvu8eR2ycPki3CSWmOgnrV2U/Jrqf3dJgBdM/t3a4zZjZNYX8zk+e8CfYDy1tojAccrAruBkdba5zMTf26goVSREGKtPWCtrWyt9V1KisnHL1pra1hrmyQ/PwysA3oCJPcqDM4XlYSwUBsmN8YUBDoDiwKTIoC19idgFdD90sSuTMTveRpKFfGuR4CpxpingUScHoWGgEJYiA6TtwFKAXuNMf9I5fWdOAnQZ4zZS3Dxe4qGUkVEroKrOUyelaFUY8wXQLtMXHYxUC2j+DMTX6hTYhQRCTHGmGXATGvtxORh8nFANY0IXB1KjCIiIcYY81ecNaqlcYbJ+1hr17kbVd6hxCgiIhJAs1JFREQCKDGKiIgEUGIUEREJoMQoIiISQIlRREQkgBKjiIhIACVGERGRAEqMIiIiAZQYRUREAigxiki2MMZUNcYkJVdbCDw+3hhzOoMq9CIhQ4lRRLKFtXY7To2+J40xpQGMMcOA3sDdKqAsXqG9UkUk2ySXVtqBUw3i38AHQBdr7WxXAxPJAhUqFpFsY609YIx5G/gnzvfL40qK4jUaShWR7PYjUAhYY619z+1gRLJKiVFEso0x5g5gArAGuN0YU8flkESyTIlRRLKFMeYW4HOcCTg+4CdgpJsxiQRDiVFE/jBjTFVgEfAVMMBaex6IBFoZY5q6GpxIFmk7hyq0AAAAY0lEQVRWqoj8IckzUVfj9BDvstaeSz4eBmwFjltrm7gYokiWKDGKiIgE0FCqiIhIACVGERGRAEqMIiIiAZQYRUREAigxioiIBFBiFBERCaDEKCIiEkCJUUREJIASo4iISID/B5rPRaB+0XE0AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(6, 6))\n",
    "ax1 = plt.subplot(111)\n",
    "ax1.set_xlim(-0.1, 1.1)\n",
    "ax1.set_ylim(-0.1, 1.1)\n",
    "ax1.set_title('Forward difference example', fontsize=16)\n",
    "ax1.set_xlabel('$x$', fontsize=16)\n",
    "ax1.set_ylabel('$f(x)$', fontsize=16)\n",
    "# define our x for plotting purposes\n",
    "x = np.linspace(0, 1, 1000)\n",
    "\n",
    "# define our example function and its exact derivative\n",
    "\n",
    "def f(x):\n",
    "    return x**2\n",
    "\n",
    "def df(x):\n",
    "    return 2 * x\n",
    "\n",
    "# plot the 'exact' solution\n",
    "ax1.plot(x, f(x), 'k')\n",
    "# choose and plot two x locations to take the difference between\n",
    "x0 = 0.35\n",
    "dx = 0.5\n",
    "x1 = x0 + dx\n",
    "# plot a line representing the discrete derivative\n",
    "ax1.plot([x0, x1], [f(x0), f(x1)], 'r', label = 'Forward diff. approx. deriv.')\n",
    "# plot a line representing the exact derivative (given by function f(.)) at x=x0\n",
    "h = dx / 2\n",
    "ax1.plot([x0 - h, x0 + h], [f(x0) - h * df(x0), f(x0) + h * df(x0)], 'b', label = 'Exact derivative')\n",
    "# add some axes labels and lines etc\n",
    "ax1.set_xticks((x0, x1))\n",
    "ax1.set_xticklabels(('$x_0$', '$x_0+\\Delta x$'), fontsize=16)\n",
    "ax1.plot([x0, x0], [-0.1, f(x0)], 'g:')\n",
    "ax1.plot([x1, x1], [-0.1, f(x1)], 'g:')\n",
    "ax1.set_yticks((f(x0), f(x1)))\n",
    "ax1.set_yticklabels(('$f(x_0)$', '$f(x_0+\\Delta x)$'), fontsize=16)\n",
    "ax1.plot([-0.1, x0], [f(x0), f(x0)], 'g:')\n",
    "ax1.plot([-0.1, x1], [f(x1), f(x1)], 'g:')\n",
    "ax1.legend(loc='best', fontsize=14);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Taylor series to estimate accuracy\n",
    "\n",
    "We can use a [Taylor series expansion](http://mathworld.wolfram.com/TaylorSeries.html), or *Taylor series analysis*, to estimate the accuracy of the method. \n",
    "\n",
    "Recall that Taylor series in one dimension tells us that we can approximate the value of the function at a location in terms of its value, and value of its derivative, at a nearby point:\n",
    "\n",
    "\\begin{align*}\n",
    "f(x_0+h) & = f(x_0) + hf'(x_0) + \\frac{h^2}{2!}f''(x_0) + \\frac{h^3}{3!}f'''(x_0) + \\ldots\\\\[5pt]\n",
    "& = f(x_0)+hf'(x_0) + \\frac{h^2}{2!}f''(x_0) + \\mathcal{O}(h^3),\n",
    "\\end{align*}\n",
    " \n",
    "where $\\mathcal{O}(h^3)$ represents the collection of terms that are third-order in $h$ or higher.\n",
    "\n",
    "We call this the Taylor series expansion *about (or around) the point $x_0$* (since all the functions in the expansion on the RHS are evaluated at this point).\n",
    "\n",
    "An equivalent way of writing this expansion would of course be (just define $x$ to be $x_0+h$)\n",
    "\n",
    "$$f(x) = f(x_0) + (x - x_0) f'(x_0) + \\frac{(x - x_0)^2}{2!} f''(x_0) + \\frac{(x - x_0)^3}{3!} f'''(x_0) + \\mathcal{O}((x - x_0)^4).$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Taylor series example\n",
    "\n",
    "*Wikipedia image: The exponential function (in blue), and the sum of the first (n + 1) terms of its Taylor series expansion around the point 0 (in red).*\n",
    "\n",
    "<img src=\"https://upload.wikimedia.org/wikipedia/commons/6/62/Exp_series.gif\" width=200x>\n",
    "\n",
    "\n",
    "More terms equate with a better approximation valid a larger distance from $x_0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "If we rearrange this expression to isolate the gradient term $f'(x_0)$ on the left hand side, we find:\n",
    "\n",
    " $$ hf'(x_0)=f(x_0+h)-f(x_0) +O(h^2) $$\n",
    " \n",
    "and therefore, by dividing through by $h$,\n",
    " \n",
    " $$ f'(x_0)=\\frac{f(x_0+h)-f(x_0)}{h}+O(h) $$\n",
    "\n",
    "As we are left with $O(h)$ at the end, we know that the forward difference method is first-order (i.e. $h^1$) -- as we make the spacing $h$ smaller we expect the error in our derivative to fall linearly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Order of convergence - higher generally better\n",
    "\n",
    "For general numerical methods we generally strive for something better than this - if we halve our $h$ we would like our error to drop super-linearly: i.e. by a factor of 4 (second-order method) or 8 (third-order method) or more.\n",
    "\n",
    "But note that as we shall see when we start to solve some differential equations, reducing the size of $h$ (i.e. refining our computational mesh) does generally come at the expense of doing more work overall.  So there is a trade-off in accuracy vs cost."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## <span style=\"color:blue\">Exercise 2.1: Compute first derivative using forward differencing</span>\n",
    "\n",
    "Use the forward difference scheme to compute an approximation to $f'(2.36)$ from the following data:\n",
    "\n",
    "$f(2.36) = 0.85866$\n",
    "\n",
    "$f(2.37) = 0.86289$\n",
    "\n",
    "You should get an answer of 0.423."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## The central difference\n",
    "\n",
    "In an attempt to derive a more accurate method we use two Taylor series expansions; one in the positive $x$ direction from $x_0$, and one in the negative direction. \n",
    "\n",
    "Since we hope to achieve better than first-order accuracy, we include an extra term in the series:\n",
    "\n",
    "\\begin{align*}\n",
    "f(x_0+ \\Delta x) &= f(x_0)+\\Delta x f'(x_0)+\\frac{\\Delta x^2}{2}f''(x_0) + \\mathcal{O}(\\Delta x^3)\\\\[5pt]\n",
    "f(x_0- \\Delta x) &= f(x_0)- \\Delta x f'(x_0)+\\frac{(-\\Delta x)^2}{2}f''(x_0) + \\mathcal{O}((-\\Delta x)^3).\n",
    "\\end{align*}\n",
    "\n",
    "Using the fact that $(-\\Delta x)^2=\\Delta x^2$ and the absolute value signs from the definition of $\\mathcal{O}$ (that is we don't worry about signs when using the $\\mathcal{O}$ notation), this is equivalent to:\n",
    "\n",
    "\\begin{align*} \n",
    "f(x_0+\\Delta x) &= f(x_0)+\\Delta xf'(x_0)+\\frac{\\Delta x^2}{2}f''(x_0) + \\mathcal{O}(\\Delta x^3),\\\\[5pt]\n",
    "f(x_0-\\Delta x) &= f(x_0)-\\Delta xf'(x_0)+\\frac{\\Delta x^2}{2}f''(x_0) + \\mathcal{O}(\\Delta x^3).\n",
    "\\end{align*}\n",
    "\n",
    "Remember that we are looking for an expression for $f'(x_0)$. Noticing the sign change between the derivative terms in the two equations, we subtract the second equation from the first to give:\n",
    "\n",
    "$$ f(x_0+\\Delta x)-f(x_0-\\Delta x)=2\\Delta xf'(x_0) + \\mathcal{O}(\\Delta x^3).$$\n",
    "\n",
    "Finally, we can rearrange to get an expression for $f'(x_0)$:\n",
    "\n",
    "$$ f'(x_0)=\\frac{f(x_0+\\Delta x)-f(x_0-\\Delta x)}{2\\Delta x} + O(\\Delta x^2),$$\n",
    "\n",
    "which, contrary to the first-order forward and backward differences seen above, is an approximation to the derivative that is second-order accurate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "By considering an interval symmetric about $x_0$, we have created a second-order approximation for the derivative of $f$. \n",
    "\n",
    "This symmetry gives the scheme its name: the central difference method. \n",
    "\n",
    "The figure below illustrates this scheme. The derivative is approximated by the slope of the red line, while the true derivative is the slope of the blue line.  \n",
    "\n",
    "Even without the analysis above it's hopefully clear visually why this should in general give a lower error than both the forward and backward difference approaches. \n",
    "\n",
    "The analysis of the two methods does tell us that as we halve $h$ the error should drop by a factor 4 rather than the 2 we get for the first-order forward or backward differencing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcYAAAGTCAYAAACyIbJtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeZxN9f/A8dfHbuy77LJlX0KWynwtJb5KiGJEIrJEIdmJbEkkW2QfkaVEGIztK3u2ME3WmhgyZB3MmHn//jjX/c2MWbkz5854Px+P++Cc8znnvM+dmfu+57MdIyIopZRSypLK7gCUUkopd6KJUSmllIpAE6NSSikVgSZGpZRSKgJNjEoppVQEmhiVUkqpCDQxqmTHGFPbGPO9MeaCMSbEGHPFGLPJGNPBGJM6kc5ZxRgzwhiTMxGO3dEYI8aYYo+4vxhjRkRYHmGMkShl8htjfjLGXHWU7+NY38wY85sx5q5jffbHuBTlZowxxRw/1452x5KcaGJUyYrjA/0XICcwAGgIdAL+AGYA/02kU1cBhjvO6+7mALWjrBsG1APedWxbaoxJA3gD54GXHOtvJmGcSrmlNHYHoFR8GWNeBCYBX4vIB1E2rzbGTAIyJX1kkTnuWo2I3Lfj/CLyN/B3lNVlgSMi8sODFcaYokAW4HsR2fG45zXGGCCtiIQ87rGUspPeMark5BPgKvBxdBtF5LSIHH2wbIwpbozxNsZcNsbcM8YcNsa8HnGfB9WOxphSxpifjTG3jDF/GmOGGWNSOcp0BOY5djnpKO+s+nT8/zNjzCfGmLNACFDRGJPBGPOlMeaY47gXjTFrjDHPPMrFG2NSG2NGG2MCjTHBxphtxpjy0ZRzVqU+qEoDPIEXIsQ+Hzjn2OVbx7ptEY7Rwhizx3Gea8aY5caYIlHOc84Ys9gY08kY87vjups6tnkYY8YbY846qrvPGmMGP3hPHWU8Hed91RjztTEmyPGzWhy1StcYk8YYM8AYc8JR7XvZGLMh4ntpjMltjJlhjDnv+Hn/box5L57vbaz7GmNSOd7vc8aYbBHWVzTG3DHGfB5h3ZvGmC2OGG8ZYw4ZYzpEc05x/Dz7On7nbjt+B/M6Xt8bY64bYwKMMQOi7Pug+v1FY8yPjvNcMcZMM8ZkjMf11jPG+BpjbjrO62OMqRCf9+qJICL60pfbv4DUQDCwJJ7lCwP/AMcAL+BlYC4QDrwaodwIQBzl+mJVzU5xrHvHUSYPMMqxrhVQy/FK79guWNWR/wNaAo2BfEA2rGrNN7GqMV8HNgHXgPwRYujoOEaxOK5plCP+iVhVn4OA0459R0S9Jsf/0ztiPQIcjBB7Uce1iOO4tYByjn26OdbPBZoAbQA/4CyQJcJ5zjmu+xjwFtAAKIFVE/U/4ArQx7F+MHAX+CLC/p6O85wFpjquqRdwB1gQ5dpXAPcd194YaI5Ve/Afx/asgD/wF9DF8XP8HAgDesXxvsZrX6CQ45qWOpYzAseBA0C6COUGAd0d19MQ+BQIBbpFOa8AfwI/Y32h6ATcADZgNRcMcew/y1G2STS/M39F+H0YgvXlZH6EcsUc5TpGWNfU8V6uBl5zvHYB/wKF7f5bd4eX7QHoS1/xeWElGgHGxrP8t8BlIFeU9ZuAwxGWRxAhCUZY/xuwMcLygw+iktGcS4ALQMY4YkoNeGC1430YzbGLxbJvDuAWMDPK+gHEkhgjrNsJbIuyrmQ0H5qZgevA3Chlizk+dPtEWHcO68tK/ihl2zuO+2KU9YMdx8jrWPZ0lIuaBL/GSqLGsVzfUe6DWN6foY59SkVZPxsIAtK4Yl+sLzcCvAN84/iZlI7l2KmwvijMxqrKjvp780eU409yrB8SYV0arC9586L5nYn6+zAYK6GXjvBzi/ozPgX4Rtkvq+NaJz/O32lKeWlVqkqpGgPrgOuOarg0xups4gNUNsZkjVL+5yjLx4AixN8GEbkTdaUxprUxZq8x5hrWt/TbWMmnTAKODVARq/30+yjrlybwOHGpjfUh6R3lffsb+B14MUr5PSJyMcq6xlh3QruiHGMjkBbr7jSiqO/9b1h3uvkcyy9hfbjPjiXuxsBe4Gw0P+9cQDlX7CtWG+0srI5eXbDuKP+IeDBjVct/Z4w5j3WnGAp0Jvqf+SaJ3Bb9u+NfnwjnvI+VzApHs390vw+pgJrRXagxphTWXX3Un28wsJuHf75PJO18o5KLK1hVbEXjWT4v8LbjFZ1cWNVWD1yNsv0ekCEB8QVGXWGMaQYsAxYAI7G+kYdjJeyEHBvgKce/l6Ksj7r8uPI6/t0cw/Z/oyw/dN2OYxTFSgjRyRVlObr3Hv7/PcoFXI3ui0eUc5ZMwDkfZ98FQFesu7glETcYYzJj1UoEY7WJn8a6S34fq6o0qqjvZ0gs66P7nYnp96FgNGXh/3++3zpeUf0Vw35PFE2MKlkQkfuOziGNjDHpReReHLtcwWrnGh/D9guujA/rjiaqN4FTItLxwQpjTFoebcjHgwSUD6tdiwjLrnTF8W/HKOd5IOpwjuiu+wpWu2HrGM5xLoExBQE5jTEZY0mOV7ASVe8YtvvHcvx472uM8cBqez0GlALGAR9GKFsb60vBCyKyM8J+ifVZG9Pvw/kYyj/4+Q4k+i8/2qMYTYwqeRkHbMPqGBF1uAbGmOJYnUOOYnVgqA0cj+NOI74eJOI4e/xF4IFVfRpRe6y2xoQ6ilUN2xrYEmH9m49wrNjswkp+JUVkwSMeYwNWJ6RbIvJ7XIXjYSPW3VdnrE46MZ2zF/CXiPyTwOMnZN8pWHdjVbDGzE42xviIyAbHdg/Hv867T2NMDqwOLokhut+HcGBfDOX9sb6YlBeRcYkUU7KniVElGyKywxjzETDJGFMWmI9V9ZMDq+djZ6AtVhIZhvXhsMMY8zXWh0EOoALwtIhEV60VmxOOf3sYYxZgffAdldjH7G0AmhtjvgTWAs9iJfRrCTw3InLNcZzBxpibWMmiBtaAfZcRkRvGmP7ANGNMHmA9Vmecglg9a7eJyJLYjoE1acA7gK8x5gusHrHpsNq2XgWai0hwAmLaaoxZifVzL4yVCNJitYf9LCLbgC+xes/+z/E++WO1yT6DdfcWW2KK177GmJZYv2PtReQM8JUx5iVgvjGmkiOp7sKqop9mjBnuOM4QrLvebLheE8dQkY1Y7YrDgYVR2z0fEBExxvTAGvebDquNMgjrTrMO1peDSYkQZ7KiiVElKyIy2RizD6v6aiKQG+sO5wBWu88aR7m/jDHVsXpojsEacnEFqwoswXdCInLEWNOuvYfV6SIVUJzYqwVnY3WY6OSIbT/QDPghln1iMwIwWB/OPbE6jDQj+irPRyYis4wxAUB/rC8aabGq5nYAh+Oxf6gx5mWsu7z3sN6n21jtbT/zaNV1b2L1wO2ANQTkOtb7OcdxzuvGmDpYX4gGYCXya1hJbmUc8ca5ryMhzwa8RWRxhN3fwfoiNt8Y01RELhtrrOwXWENMLmDdZebESlqu5oU1zOh9rPd1NtAvth1EZJ2xJssYjPX+ZQQuAnuw2sSfeA+6QyullEomzP9POlFKRE7ZHE6Ko8M1lFJKqQg0MSqllFIRaFWqUkopFYHeMSqllFIRaGJUSimlItDhGslc7ty5pVixYnaHoZRSycqvv/4aJCJ5otumiTGZK1asGAcOHLA7DKWUSlaMMX/GtE2rUpVSSqkINDEqpZRSEWhiVEoppSLQxKiUUkpFoIlRKaWUikB7paZwN27c4J9//iE0NKaHkyul7JQ2bVry5s1L1qxZ7Q5FOWhiTMFu3LjBpUuXKFiwIBkzZsQYY3dISqkIRIQ7d+5w/vx5AE2ObkKrUlOwf/75h4IFC+Lh4aFJUSk3ZIzBw8ODggUL8s8//9gdjnLQxJiChYaGkjFjRrvDUErFIWPGjNrc4UY0MaZweqeolPvTv1P3oolRKaWUikATo1IuNGLECCpUqJCgfYKCgjDGsG3bNgDOnTuHMSbSHLi//PILlSpVIl26dHh6esa4TsVs6dKlpEnz+P0N79+/jzGGH3/80QVRKXekiVG5pUuXLtG7d29KlChB+vTpKViwIK+88grr1q1z6XkeJZEltsKFCxMYGEiVKlWc63r37k3lypU5ffo0q1atinGdSnxp0qQhMDCQV155xe5QVCLR4RrK7Zw7d466deuSJUsWxo4dS+XKlQkPD8fX15du3brx119/JXlMISEhpEuXLknOlTp1avLnzx9p3alTp+jRoweFCxeOdZ07S8r3MLE8uIaoPx+Vsugdo3I73bt3R0Q4cOAArVu3pkyZMpQtW5aePXty5MgRZ7nr16/z3nvvkTdvXrJkyUK9evUiVT/Onz+fzJkz4+vrS4UKFciUKRP/+c9/OHv2rHP7yJEjOX78OMYYjDHMnz8fsDpDTJs2jRYtWpApUyYGDRpEWFgY7777LsWLFydjxoyUKlWKCRMmEB4enqDr279/P88++ywZMmSgatWq7N27N9L2iFWpD/5//fp1OnXq5IwxunXxERoaSqdOnZzXULp0aSZOnIiIOMt4eXnRvHlzRo4c6XxvO3fuzN27d51lnn/+eXr27EnPnj3JkSMHOXPm5JNPPon0XhQqVIhRo0bRsWNHsmXLRocOHQA4cuQI9evXJ2PGjOTKlYtOnTpx48YNAIKDgylbtizvvfee8zh///03OXPmZPLkyQl6n+fNm0eRIkXw8PDg1VdfjXY4xOrVq6lWrRoZMmTg6aefZujQoYSEhMR6DVGrUmvUqMGAAQMiHffatWukT5+eNWvWJChm5SZERF/J+PXss89KTE6cOBHjNnd15coVMcbIZ599Fmu58PBwqVu3rjRp0kT27t0rJ0+elCFDhkiWLFnkwoULIiIyb948SZMmjTRo0ED27t0rR44ckSpVqshLL70kIiLBwcHSt29fKVOmjAQGBkpgYKAEBweLiAggefLkkdmzZ8vp06flzJkzEhISIkOHDpV9+/bJ2bNnZdmyZZItWzaZM2eOM67hw4dL+fLlY4z71q1bkidPHmnVqpX89ttvsmHDBnnmmWcEkK1bt4qIyNmzZwWQ/fv3y/379yUwMFA8PDxk8uTJEhgYKLdu3Xpo3YO443Lnzh0ZPny48xq+++47yZIli8yfP99Zpl27dpI5c2Zp06aNHDt2TNavXy9PPfWUfPjhh84ydevWlcyZM0vv3r3Fz8/PeZwpU6Y4yxQsWFCyZs0qn3/+uZw8eVL++OMPuXnzpuTLl09atGghR48ela1bt0qJEiWkdevWzv0OHjwo6dKlk5UrV0pYWJh4enrKSy+9JOHh4fG6RhGRX375RYwxMmbMGPH395fp06dLzpw5JXXq1M4yP//8s2TNmlXmzZsnp06dEl9fXylZsqQMGDAg1msIDQ0VQH744QcREZk8ebIULlw4UnzffPON5MqVS0JCQuIdc3L8e03OgAMSw+eq7R/s+krixNi7t0i9ekn76t07xhij2rt3rwCyatWqWMv5+vpKpkyZHkoIlStXlvHjx4uIlRgB+f33353bFy9eLGnTppWwsDARiTmRAdKzZ8844x0wYIA0aNDAuRxXYpw1a5Zky5ZNbt686Vy3aNGiGBPjA5kyZZJ58+ZFOlZ06x5F37595eWXX3Yut2vXTnLmzCm3b992rps3b55kyJDB+X7XrVtXypYtGykZDB8+XIoWLepcLliwoDRv3jzSuaZPny45cuSQW7duOddt2rRJADlz5oxz3eeffy45c+aU3r17S+7cuZ1fduLrjTfekMaNG0da16FDh0iJsXbt2jJmzJhIZZYvXy5Zs2aN9RqiJsZLly5JmjRpZNu2bc4y9erVk+7duycoZk2MSSu2xKhVqcqtWL+vcfv1118JDg4mT548ZM6c2fk6duwYp0+fdpZLnz49ZcqUcS4XKFCA0NBQrl27Fuc5qlev/tC6mTNnUr16ded5v/zyywS1efr5+VGpUiUyZ87sXFe7du147+8K06ZNi3QNU6dOfegaKleujIeHR6QY796966yGBqhVq1ak8Xe1a9fmzz//5Pbt2851Ud9DPz8/KleuTKZMmZzr6tat69z2QN++fSlfvjxTpkzhm2++4amnnkrQNfr5+T30vkZd/vXXX/n0008j/f68/fbb3Lhxg8uXL8d4DVHlzZuXhg0b4u3tDVhVvzt27MDLyytBMSv3oZ1vnjQJbKdJaqVKlcIYg5+fH6+//nqM5cLDw8mXLx//+9//HtoWcb7JqN3zH3yQx6ddMOKHN8CyZcvo06cPEydOpE6dOmTNmpVp06bxww8/xHmsB+Kb+BOLt7c3/fr144svvqBWrVpkzZqVr776yuW9fR+I+h6KyEOD2R8sR1z/zz//4O/vT+rUqTl16lSCzxuf91lEGDlyJC1atHhoW86cOZ3/j3oN0fHy8qJXr158/fXXLFmyhKeffjrJv/Ao19HEqNxKzpw5efnll/n666/54IMPIt1ZgdWpIXv27FSrVo1Lly6RKlUqnn766Uc+X7p06QgLC4tX2Z07d/Lcc8/Rs2dP57qId6fxUa5cORYsWMDt27edH7h79uxJ0DEex86dO6lTpw7du3d3rosu8Rw5coQ7d+44pxTcs2cP6dOnp3jx4s4yUePes2cPRYoUiTWRlCtXDm9v70jXv3PnTgDKli3rLNepUyeeeeYZevToQfv27WnYsCFVq1aN93WWK1cu2vgiqlq1Kv7+/pQsWTLex43J66+/Trdu3Vi3bh3e3t60a9fusY+p7KNVqcrtTJ8+HRGhevXqLF++HH9/f37//XdmzJhBpUqVAGjYsCF169bltddeY/369Zw9e5bdu3czfPjwaO8iY1KsWDH+/PNPDh48SFBQEPfu3YuxbOnSpTl48CDr16/n5MmTjBo1iu3btyfo2tq2bUuaNGno1KkTx48fZ9OmTXz22WcJOkZMwsLCeOaZZ5g5c2aMZUqXLs2BAwfw8fHh5MmTjBgxgl9++eWhciEhIbz77rucOHECHx8fBg0aRLdu3SLNvRsQEMBHH32Ev78/33//PZMmTeLDDz+MNcb27duTLl06OnTowLFjx9i2bRvvv/8+rVu3plixYgB8/fXX/PLLLyxatIjWrVvTrl072rVrx507dwD466+/eOaZZ2Lt8fnBBx+wYcMGJkyYwMmTJ5k5cyY//fRTpDLDhw9n4cKFjBgxguPHj/P777+zfPlyPvnkk1ivIToeHh40b96cESNGcPTo0YeqUfv27UuTJk0SfFxlD02Myu0UL16cgwcP0qhRIwYMGEClSpWoX78+P/30E7NmzQKsard169ZRv359unTpQpkyZWjdujX+/v4UKFAg3udq2bIlTZo0oUGDBuTJk4fvvvsuxrJdu3aldevWtG3blho1anDu3Dn69u2boGvLnDkza9eu5eTJk1SrVo1+/foxfvz4BB0jJiKCv78/QUFBMZbp3r07LVq0oE2bNtSsWZPz58/Tp0+fh8o1aNCAUqVKUa9ePVq2bMnLL7/M2LFjI5V5++23uXPnDjVr1qRr16507dqVDz74INYYM2fOjI+PD1evXqVGjRq0aNGCF154gdmzZwNw4sQJPv74Y6ZPn06RIkUA+OqrrwgNDaVfv36AlbT9/f25fv16jOd5/vnn+eabb5g6dSqVKlVi7dq1DB8+PFKZJk2asGbNGjZt2kSNGjWoWbMmEyZMcJ43odq3b8+RI0eoWbMmpUqVirTt/PnznDlz5pGOq5KesbvNQz2e6tWrS8SxexH5+flFqp5SKj68vLy4detWrFOePf/881SvXj3BYwtVzPTvNWkZY34VkWh7Vukdo1JKKRWBJkallFK2GjlypFvVPmivVKVUJIsXL46zzIOepEo9roCAAMaMGeNW4z71jlEppZRtPvvsM0SEoUOH2h2KkyZGpZRStjhz5gzffvstXbp0cQ7XcQeaGJVSStni008/JU2aNAwePNjuUCLRxKiUUirJ+fn5sWjRIrp3756gscdJQROjUkqpJDd06FA8PDwYOHCg3aE8RBOjUkqpJHXgwAFWrlxJ3759yZ07t93hPEQTo1IuNnHiRJd0JDh37hzGGGKa2ciVPD09I02OrlRiGjJkCDlz5uSjjz6yO5RoaWJUbqdjx44YYx561apVK8liGDFiBBUqVEiy80WncOHCBAYGUqVKFZcdc/78+Q89sQRg1apVD82FqlRi2L59Oz4+PgwcODDSI+LciQ7wV26pYcOGLFq0KNK6dOnS2RRN0gsJCSFdunTkz58/Sc4X8fmDSiUWEWHQoEEUKFCAHj162B1OjPSOUbml9OnTkz9//kivBx/e27dvJ23atGzbts1ZfubMmWTNmtX5BIMNGzbwwgsvkCNHDuczHiM+IR7gwoULtGvXjly5cuHh4UGVKlXYunUr8+fPZ+TIkRw/ftx5tzp//vwYY50wYQL58+d3PgH+1q1bD5WZN28e5cqVI0OGDJQuXZovv/wy0sOSjTFMmzaNFi1akClTJgYNGhSpKjU8PJxChQoxderUSMf9448/MMZw6NAhACZNmkSlSpXIlCkTBQsWpHPnzly7dg2Abdu28c4773D79m3ndY0YMQKIXJU6cOBAnn322YeuoU6dOvTu3Tve16RUVOvWrWPXrl0MGzYs0iPM3I6I6CsZv5599lmJyYkTJ2Lc5s46dOggTZs2jbXMwIEDpVChQnLlyhXx8/MTDw8PmT9/vnP7ihUrZMWKFfLHH3/IkSNH5I033pASJUrIvXv3RETk1q1bUrJkSalTp45s375dTp06JStXrpQtW7ZIcHCw9O3bV8qUKSOBgYESGBgowcHB0caxbNkySZs2rcycOVP8/f1l9OjRkiVLFilatKizzDfffCP58+eX5cuXy5kzZ+Snn36SfPnyydSpU51lAMmTJ4/Mnj1bTp8+LWfOnJGzZ88KIPv37xcRkX79+slzzz0X6fzDhg2TcuXKOZe//PJL8fX1lbNnz8q2bdukYsWK4uXlJSIi9+7dk8mTJ4uHh4fzum7evCkiIvXq1ZMePXqIiMjx48cFED8/P+dxz5w5I4Ds27cv3tekEia5/r3GV1hYmFSuXFlKlCghISEhdocjwAGJ4XPV9g92fSVtYuzdW6RevaR99e4dY4jR6tChg6ROnVoyZcoU6fXxxx87y4SEhEj16tXl9ddfl6pVq0rr1q1jPeatW7ckVapU8r///U9ErA/2zJkzy+XLl6MtP3z4cClfvnycsdauXVs6d+4caV2DBg0iJcbChQvLwoULI5X58ssvpWzZss5lQHr27BmpTNTEeOTIEQHk5MmTzjIlS5aUMWPGxBjf+vXrJV26dBIWFiYiIvPmzZNMmTI9VC5iYhQRqVKligwZMsS5PGrUKCldunSCrkklTEpPjEuXLhVAvL297Q5FRGJPjNrGqNzSiy++yDfffBNpXfbs2Z3/T5s2LUuWLKF8+fLkzZuXLVu2RCp7+vRphg4dyt69e7l8+TLh4eGEh4fz119/AXDo0CEqVar02F3F/fz86Ny5c6R1tWvX5tSpUwBcvnyZgIAAunbtyvvvv+8sc//+feubaQTVq0f7aDinSpUqUbFiRZYsWcKwYcPYu3cvp0+fpm3bts4yW7ZsYezYsfj5+XH9+nXCwsIICQnh4sWLCRpE7eXlxfTp0xk1ahQA3t7ezkmeE3JNSgGEhoYydOhQKlasyJtvvml3OHHSxPiEcaMnu8TKw8ODkiVLxlpmz549hIeHc+3aNS5fvhwpcTZr1oyCBQsya9YsChYsSJo0aShXrhwhISEASfYB/qDNbebMmdSpUyfWspkyZYrzeO3atWPu3LkMGzYMb29vXnjhBYoWLQrAn3/+SdOmTenSpQuffvopuXLl4uDBg7z11lvO646vtm3b8vHHH7N7927Sp0/P77//Trt27RJ8TUoBLFiwgJMnT7J69WpSpXL/ri3uH6FS0Th37hw9e/Zk2rRpNGrUiHbt2nH//n0Arly5gp+fH4MGDaJhw4aULVuWmzdvOrcDVKtWjaNHjxIUFBTt8dOlS0dYWFiccZQtW5Y9e/ZEWhdxOV++fBQsWJDTp09TsmTJh14J1a5dO06dOsWePXtYtmxZpEf1HDhwgJCQEL788ktq165N6dKluXDhwiNd11NPPUX9+vXx9vbG29ubOnXq8PTTTyfKNamU7e7du4wcOZJatWrRrFkzu8OJF71jVG7p3r17XLx4MdK61KlTkydPHsLCwvDy8qJevXp07dqVVq1aUbFiRUaOHMmoUaPIkSMHuXPnZvbs2RQuXJjz58/Tv39/0qT5/1/3tm3bMm7cOJo3b87YsWMpVKgQv/32G1myZOE///kPxYoV488//+TgwYMUKVKELFmykD59+ofi7N27N2+//TY1atTA09OTFStWsHfv3kjDH0aMGEGvXr3Inj07TZo0ITQ0lIMHD3L+/PkET4dVqFAhXnzxRbp168b169d54403nNtKlSpFeHg4kydPpkWLFuzZs+ehh78WK1aMu3fvsmnTJqpWrYqHhwceHh7RnsvLy4t+/fqRLl06hgwZEmmbK69JpWwzZ87k77//ZuHChRhj7A4nfmJqfNRX8nil1F6pwEOvggULiojIp59+Kvny5ZN//vnHuc/GjRslTZo0zs41vr6+Ur58eUmfPr2UL19eNmzYIJkyZZJ58+Y59wkICJDWrVtLtmzZJGPGjFKlShXZunWriIjcvXtXWrZsKdmzZxcg0n5RjRkzRvLkySOZMmWSt956S4YPHx6p842IyJIlS6Rq1aqSPn16yZ49u9StW1e+++4753ZAli9fHmmfqJ1vHvj2228FkBYtWjwUy5QpU6RAgQKSIUMGqV+/vixbtkwAOXv2rLNMt27dJFeuXALI8OHDReThzjciIjdv3hQPDw9JmzatBAUFPXSuuK5JJUxy/XuNzY0bNyR37tzSsGFDu0N5CLF0vjGijeXJWvXq1SWmKcP8/PwoW7ZsEkeklHoUKfHvddSoUc6OYjVr1rQ7nEiMMb+KSLQ93rSNUSmllMsFBQUxceJEmjdv7nZJMS6aGJVSSrnc6NGjuXXrFmPGjLE7lATTxKiUUsqlzp49y/Tp0+nUqaP+JqUAACAASURBVFOyrB7WxKiUUsqlhgwZQpo0aZxz8SY3mhhTOO1cpZT7S0l/p4cOHWLJkiX06dOHggUL2h3OI9HEmIKlTZuWO3fu2B2GUioOd+7cIW3atHaH4RIDBgwgZ86cDBgwwO5QHpkmxhQsb968nD9/nuDg4BT1jVSplEJECA4O5vz58+TNm9fucB7bpk2b2LRpE0OGDCFbtmx2h/PIdBxjMhfbOEaAGzdu8M8//xAaGpqEUSml4itt2rTkzZvXbZ9mH1/h4eFUr16dq1ev4u/vH+1MUe4ktnGMOiVcCpc1a9Zk/wenlHJ/y5Yt49ChQyxatMjtk2Jc9I4xmYvrjlEppRJbSEgIzzzzDFmzZuXgwYPJ4gkaeseolFIq0cycOZOzZ8+yYcOGZJEU45L8r0AppZRtbty4wahRo6hfvz4vvfSS3eG4hCZGpZRSj+zzzz8nKCiI8ePHJ5/HSsVBE6NSSqlHEhgYyKRJk2jTpg3Vq0fbXJcsaWJUSin1SEaOHElISAijR4+2OxSX0sSolFIqwfz8/JgzZw5du3alZMmSdofjUpoYlVJKJdjHH39MpkyZGD58uN2huJwO11BKKZUgW7ZsYe3atYwbN448efLYHY7L6R2jUkqpeAsLC6Nv374ULVqU3r172x1OotA7RqWUUvG2aNEiDh8+zJIlS8iQIYPd4SQKvWNUSikVL8HBwQwePJgaNWrQpk0bu8NJNHrHqJRSKl6++OILLly4wLJly1LE1G8xSblXppRSymUCAwMZP348LVq04Pnnn7c7nESliVEppVSchg0bxr179xg/frzdoSQ6TYxKKaVi9dtvvzF37lx69OiR4gbzR0cTo1JKqVj179+frFmzMmzYMLtDSRLa+UYppVSMfHx88PHx4YsvviBnzpx2h5MkjIjYHYN6DNWrV5cDBw7YHYZSKgUKCwujSpUqBAcHc+LECdKnT293SC5jjPlVRKJ9JIhbVKUaYwoaY1YZYy4bY8KNMZOMMVONMWvsji2+jDFzjDFijJn0GMf40Bhz1BjjFj8XpdSTbd68eRw7dozx48enqKQYF7e4YzTGbAIKAP2Bq0BqYCtQR0Tc/nbIGJMRuAhkAS4DBUXk/iMe5ywwUETmxWcfvWNUSiWGmzdvUrp0aZ5++ml27tyZYh5C/IBb3zEaY/IDDYDxIrJORPYAbwJH7EiKxphzxpgRCdztdSArMAHICzR+lHOLyB1gIdDvUfZXSilXGTt2LBcvXmTSpEkpLinGxdbEaIxZBQQCBljgqIocC3gBS6KULWmMCTXGjIyyfoYx5qYxxs7HR3fAutMbinXH+HbUAgmIfylQzhhTJ7GDVkqp6Jzx8+PbiRNp3749zz33nN3hJDm77xiHAnOAW0Btx2sDkB34X8SCInLKUfZDY0xuAGPMMKAT8LpdVa7GmAJAQ2CxiIRiJbZXjTE5IpZLQPyHgRvE867T39+f+fPnAxAaGoqnpyeLFy8GrHkNPT09WbZsGQDXr1/H09OTVatWARAUFISnpydr1lhNuRcvXsTT05MNGzYAEBAQgKenJ5s3bwbgzJkzeHp6sn37due5PT092bVrFwDHjh3D09OT/fv3Wxdy+DCenp4cPnwYgP379+Pp6cmxY8cA2LVrF56envj7+wOwfft2PD09OXPmDACbN2/G09OTgIAAADZs2ICnpycXL14EYM2aNXh6ehIUFATAqlWr8PT05Pr16wAsW7YMT09PgoODAVi8eDGenp6EhoYCMH/+fDw9PZ3v5ezZs2nYsKFzefr06bzyyivO5SlTpvDqq686lydOnEjLli2dy+PGjePNN990Lo8aNQovLy/n8rBhw3jnnXecywMHDuS9995zLvfr148ePXo4l/v06UOfPn2cyz169KBfv/+vTHjvvfcYOHCgc/mdd96J1J3ey8uLUaNGOZfffPNNxo0b51xu2bIlEydOdC6/+uqrTJkyxbn8yiuvMH36dOdyw4YNmT17tnPZ09NTf/dS4u/enTsMqV6ddWFhjI3w+/MksTUxishxrCR4VET2OKpRawECHI1ml5FY7Y8DjDHvAsOB9iKy+VHObyxpIr4cm1JFWZ86lsO0x3ofFzuWFwLpgehm2I0zfhEJx7r2WrHE/Z4x5oAx5sCDPzSllHpswcFcfeEFAoODmVOhAgWLFrU7IlvY3vnGGHMS8BGRno7lKUAHEckeQ/nPgL5YYzB7i8i0KNtLAAuw2vpuA11iups0xnhidfKJy3YR8YzhGMeBmyJSK8I6P+BfEXmoOjSu+B1lVgGlRaRCXIFp5xullEvcuoU0bUr4jh30z5WLzwICyJgxo91RJZrYOt/YOsDfGJMFKAEcirA6A3Avlt1OYt2R7YwuqQAzgfkiMscY0wjwNsY8I9F/A/gVqBFl3U/AWuCbCOtuxhB/DaAc0DPKpkXAZ8aY0iLyRwLjB7gDpNzfSKWUe7l+HZo0QfbswQtoMWNGik6KcbG7jbEKVsebiInxCpAjusLGmPrALGA3UNcYUznK9jxYVZALAERkk2PTs9EdT0RuisiBiC8gBLgQZb1/DPF3AEKBZVHWL8aqDo7UCSeu+CPICQTFsE0ppVzn33/hpZeQfft4N1Mmzr/wArNmzYrU5vmksTsxVsVKLMcjrPsdSGuMKRSxoDGmGvAjVgcWT+AvYEyU4xXBSmoRG97+dKx3KWNMOqxhJetFJFISE5G/gO1Ae+Po5xzP+B8oDsSUjJVSyjWuXIEGDeDwYeY1bcqCW7eYMmUKbdq0SdEPIo6L3XOlVgVOiEjEqtMdjn9rAn+DNdQBWA9sBHqJSLhj2MNcY8yLIrKDmCXWAJz/ArmAAGNM82i2n8FKgJ7GmADiGb8xJjtQGpgYzTGVUso1/vkHGjaEP/4gYNo0ur7/Pu+++y5Vq1alatWqdkdnK3dIjIcjrhCRc8aYfUAzYJVjAoCNgB/QztFrE6zenx8D44AHnVz+AgoYY9JGuGss6ljvah0c//ZwvGIyACvRxSd+gKZY1bk/uDRapZR6IDDQulM8dw7WruX9yZPJmDEjo0ePtjsyt2BrYhSRKjFsmgFMMcb0EJGLwNPR7BsGlI2y7rIjqXYEZjs63xisTjbxjalYPMu9Ft9jxrD/Q/E7eAHLReTK4xxfKaWi9fffUL8+XLgA69fjc/cuP//8MxMmTCBfvnwAzjGW27Ztsy9OG9l9xxiTRVh3U91JeJViN6xZdPoDwVh3afZPCBsPxpgqwH+AOIdpKKVUgv35p5UUL18GHx9Ca9bkw8qVKVmyJB988IGzWMeOHe2L0Q24ZWIUkTBjTCeg2iPse5LIVZPJSX7gHccsOUop5TpnzsB//gM3bsDmzVCzJjOnTsXPz4/Vq1dHenrGk54YbR/grx6PDvBXSsXpjz+sO8U7d2DTJqhWjcuXL1O6dGmqV6/Oxo0bI00U/mBGrbRp09oVcaJz2wH+SimlEtmJE1ZHm7Aw2LoVKlUCYPDgwdxyDM+I+vSMRo0aAdrGqJRSKqX57TcrKaZODdu2QblygDWp+pw5c/joo48o51gXUefOnZM4UPeiVanJnFalKqWidegQNGoEGTLAli1QujQA4eHh1KpVi4CAAPz9/cmaNavNgdpDq1KVUupJsm8fvPwyZM1qJcUSJZyb5s2bx/79+1m0aFGMSfHB47I8PDySJFx3o4lRKaVSkl27oHFjyJPHSooRHh3177//8sknn/D888/Trl27GA/RpEkTQNsYlVJKJXfbt0PTplCggJUUC0Wacpphw4Zx9epVvv7664c63ET0/vvvJ3akbk0To1JKpQS+vtCsGRQrZv3/qacibT5y5AjTp0+ne/fuVK4c04N9LE/yBOJg/9M1lFJKPa4NG+C//4WSJa3ep1GSoojQs2dPcubMyaeffhrn4a5fv87169cTKVj3p3eMSimVnK1ZA61aQfny1uD9XLkeKuLt7c3OnTuZM2cOOXJE+7jbSF57zZoKWtsYlVJKJS8rV8Kbb0LVquDjA9EkvRs3btC/f39q1KjBO++8E6/DRpw39UmkiVEppZKj776D9u3huedg3TrIli3aYp9++imXLl3ip59+IlWq+LWetWjRwpWRJjvaxqiUUsnNwoXg5QV161rtizEkxRMnTjBlyhTeffddatSoEe/DBwUFERQU5Kpokx29Y1RKqeTk22+hSxdrUvDVqyFTpmiLiQi9evUic+bMjBkzJkGnaNWqFaBtjEoppdzd9OnQo4c1gH/VKsiYMcaiS5YsYcuWLcyYMYM8efIk6DR9+/Z93EiTNZ0rNZnTuVKVekJMngwffmiNVVy+HCI8PzGqa9euUaZMGYoVK8bu3bvj3bb4JNG5UpVSKjmbMAEGDICWLWHJEkiXLtbigwcPJigoiA0bNjxSUrx48SIA+fPnf6RwkztNjEop5c5GjYJhw6xhGYsWQZrYP7b37dvHjBkz6NWrF1WrVn2kU7755puAtjEqpZRyJyJWQhw9Gt5+G+bOtZ6rGIuwsDDef/998ufPz6hRox751J988skj75sSaGJUSil3I2JVnX7+OXTuDLNmQTyqRKdPn87BgwdZunTpYz1nsXHjxo+8b0qgLbJKKeVORKxONp9/Dt27xzspBgYGMmTIEF566SVat279WCEEBAQQEBDwWMdIzvSOUSml3EV4OPTsCTNmQJ8+MGkSxPJ4qIg++ugj7t27x7Rp02J9pFR8tG/fHtA2RqWUUnYKC4OuXa0B/AMGwNix8U6KmzZtYunSpYwYMYKSJUs+dihDhgx57GMkZzqOMZnTcYxKpQD370OnTlav02HDYMSIeCfFu3fvUqlSJUSE3377jQwZMiRurCmEjmNUSil3FRpqTQa+bJk1NCOBd2vjx4/n5MmTbNy40WVJ8cyZMwA8/fTTLjlecqOJUSml7BISAm+9ZU3vNmEC9O+foN1PnTrF2LFjadOmDY0aNXJZWJ06dQK0jVEppVRSuncP3njDetDw5MnQu3eCdhcRunXrRrp06Zg0aZJLQxs5cqRLj5fcaGJUSqmkducOvP669XDhGTOgW7cEH2LRokX4+voyffp0ChQo4NLw6tWr59LjJTeaGJVSKindvg2vvgpbt1o9UB3Vlglx+fJlPvroI+rUqUPXrl1dHqK/vz8AZcqUcfmxkwNNjEoplVRu3oSmTeGXX/7/YcOP4MMPP+TGjRt88803ifLkjAfJVtsYlVJKJZ7r1+GVV2DfPusJGW3aPNJhfHx88Pb2ZujQoZQvX97FQVoS+mDjlEbHMSZzOo5RqWTg33/h5Zfh8GFYuhRatHikw9y+fZsKFSqQPn16Dh8+rGMWH4OOY1RKKbsEBUGjRnDihDUs47//feRDjRgxgnPnzrF9+/ZETYrHjh0DoEKFCol2DnemiVEppRLLpUvQsCGcOgU//WTdNT6igwcPMmnSJLp06cKLL77owiAf1rNnT0DbGJVSSrnShQvQoAH8+SesXWv9/xHdv3+fzp07kzdvXiZMmODCIKP3+eefJ/o53JkmRqWUcrWAAKhfHy5ehA0b4DHv8KZMmcKhQ4f4/vvvyZ49u4uCjFmNGjUS/RzuTBOjUkq50rlzVlK8cgU2boTatR/rcGfPnmXYsGE0a9aMVq1auSbGOBw+fBiAKlWqJMn53I0mRqWUcpXTp62keOMGbN4Mj3nnJSK8//77pEqVyiXPWYyvPn36ANrGqJRS6nH4+1tJ8d492LIFqlZ97EMuXLgQHx8fvvrqKwoXLuyCIONn8uTJSXYud6TjGJM5HceolBs4ccJKiiLWnWLFio99yMDAQMqVK0f58uXZsWNHosxw8ySLbRyjvtNKKfU4jh4FT09IlQq2bXNJUnxQhXr37l2+/fbbJE+K+/fvZ//+/Ul6TneiValKKfWoDh60Bu97eFjVp6VKueSw33//PatXr2bChAm2TOTd3/FcyCe1jVGrUpM5rUpVyiZ791oD9rNnt56UUby4Sw57+fJlypUrR/Hixdm1axdp0iT9/cuTMPONTgmnlFKutHMnNGkCefJYSbFIEZcdulevXly/fp25c+fakhQhZSfE+NA2RqWUSoht26BxY3jqKdixw6VJ8YcffmDZsmUMHTrU1uS0a9cudu3aZdv57aZVqcmcVqUqlYQ2b7YeMly8OPj6Qv78Ljv01atXKV++PPny5WP//v2kTZvWZcdOKE9PTyBltzFqVapSSj2udeusx0WVKWMlyDx5XHr4jz76iMuXL/Pzzz/bmhQBZs2aZev57aaJUSml4rJ6NbzxhjUUY+NGyJXLpYdfv349CxYsYPDgwVSrVs2lx34UdvSEdSdalZrMaVWqUolsxQp46y2oVg18fKxeqC50/fp1KlSoQJYsWTh06BDp06d36fEfxfbt2wGoV6+ezZEkHq1KVUqpR7FkCbz9NtSqZVWlZs3q8lP079+fCxcusGvXLrdIigDDhw8HUnYbY2w0MSqlVHQWLIB33oF69WDNGsic2eWnWL9+PbNnz6Z///4899xzLj/+o5o7d67dIdhKq1KTOa1KVSoRzJ4NXbtCw4bw44/WzDYudvXqVSpUqEDOnDk5cOAAGTJkcPk5VMy0KlUppeJr2jTo2dMawL9yJSRSwvrggw+4fPkya9ascbukuHnzZgAaNmxocyT20MSolFIPfPklfPQRvPYaLFsGidTmt3LlSry9vRkxYgTPPvtsopzjcYwePRp4chOjVqUmc1qVqpSLjBsHAwdCq1ZWp5tEGkt46dIlKlSoQNGiRdm9e7ftYxajExAQAJCkz4BMalqVqpRSMRGBUaNg+HBo29bqdJNIc5SKCN26dePmzZssWLDALZMipOyEGB+aGJVSTy4RGDIExoyBDh3g228hdepEO93ixYv58ccf+fzzzylfvnyinedxbdiwAYDGjRvbHIk9tCo1mdOqVKUekQj07w9ffAFdusDMmdbDhhPJ33//TYUKFahYsSLbtm0jdSIm4Melc6UqpdSTRgR694apU6FHD/jqq0RNiiLCu+++S2hoKPPnz3frpAiwdOlSu0OwlSZGpdSTJTwcuneHWbOsHqgTJ4IxiXrKWbNmsXHjRqZPn06JEiUS9VyukN+FTw1JjvR5jEqpJ0dYGHTubCXFgQOTJCmeOnWKfv360ahRI7p165ao53KVNWvWsGbNGrvDsI3eMSqlngz370PHjuDtbfVAHT480ZPi/fv38fLyIm3atMydOxeTyOdzlS+++AKAZs2a2RyJPRKUGI0xtYDGQC2gAJARCAL8ge3AjyLyr6uDVEqpxxIaCl5e8P338NlnMGhQkpx29OjR7N27l2XLllGoUKEkOacrrFixwu4QbBWvqlRjTAdjzG/ALqAP4AGcBPYC/wLPAXOA88aY+caY4okUr1JKJUxICLRpYyXFiROTLCnu3r2bUaNG0b59e1q3bp0k53SV3Llzkzt3brvDsE2cd4zGmCNAXmAh8DZwWKIZ42GMyQb8F2gHHDfGvCMiy1wcr1JKxd/du9ZMNj//bPU87dUrSU578+ZNvLy8KFKkCF9//XWSnNOVVq1aBUCLFi1sjsQe8alKnQfMFJG7sRUSkeuAN+BtjKkMPNndmpRS9goOhtdfh40brc42772XZKfu3bs3586dY/v27WRNhGc4JravvvoK0MQYIxGZnNCDisgR4MgjRaSUUo/r9m1o1gy2bYO5c63nKiaRlStXMm/ePAYPHszzzz+fZOd1pdWrV9sdgq0SNPONMaaaiBxMxHhUAunMN0pFcfOm9cioXbuseU+9vJLs1OfPn6dSpUo8/fTT7Nq1y23nQlWxz3yT0HGMW40x/3FBTEop5XrXrsFLL8Hu3fDdd0maFMPDw3nnnXe4e/cu3t7eyTopLlu2jGXLntwuIglNjEuAdcaYllE3GGOeN8bsdE1YSimVQFevQsOG8OuvsGIFJHFP0K+++opNmzYxadIkSpcunaTndrUZM2YwY8YMu8OwTYLGMYrI+8aYQGCpMaaXiMw0xlQExgBNAb/ECFIppWJ1+TI0agR+fvDDD9C0aZKe/rfffuOTTz6hWbNmvJeEnXwSy7p16+wOwVYJnvlGRD41xpwHZhhj3gLqAgFAJ6whHUoplXQuXrTuFE+fhjVrrKrUJBQcHEybNm3Inj07c+bMSTaz28TGw8PD7hBsleDEaIzJCZQGwoAXsAb9e4rIfRfHppRSsbtwAerXh4AAa6xi/fpJHkKfPn3w8/Nj48aN5M2bN8nPnxgWL14MgFcSttG6kwS1MRpjhgNngB7AF1h3idWBSa4PTSmlYhEQAPXqwfnzsGGDLUlx+fLlzJ49m08++YRGjRol+fkTy5w5c5gzZ47dYdgmocM1QrCmfhspIpcc6+oDPwAbAC8RCU2MQFX0dLiGeiKdPWslwn//tZJirVo2hHCWqlWrUrZsWXbs2JGse6FGFRpqfYynpGuKypUPKi4rIqcjrhCRLY4hHOuwkmODRwtTKaXi4dQpKyneugWbN0P1aD/bElVoaChvvfUWIsKSJUtSXAJJadeTUAntlXo6hvUHjTHPAz4uiUoppaLz++9WUgwNhS1boEoVW8IYNmyY86kZxYunvGcmzJ8/H4COHTvaGoddXPY8RhE5ZYyp46rjKaVUJMeOWb1PAbZuhQoVbAlj8+bNjB8/ns6dOye7p2bE15OeGONsYzTGrAZGiMiheB3QmAxAdyBYRGY+fogqNtrGqJ4IR45YSTFtWutO8ZlnbAnj0qVLVKlShRw5cnDgwIEnflhDcva4bYx/AnuMMYexZr75H3A04vAMY0wBoCbQDGgBnMfqsaqUUo/nwAFrbGLmzFZSLFnSljDCw8Pp0KED165dY+PGjZoUU7D4JMYQ4D/AW8BwIBsgxpgbwD0gB5AWMMA+rAcZLxKR8ESJWCn15NizB15+GXLmtJKije15kyZNwsfHh+nTp1OxYkXb4kgKs2fPBqBLly42R2KP+FSl3gfqiMg+Y8wYrJ6ntYGngAzAFeB3YIeI/JnI8aootCpVpVg7d8Irr0C+fFZSLFLEtlB++eUX6tWrx2uvvcaKFStSxOw2sWnoaMvdvHmzzZEkntiqUuOTGP8B2ouIjzEmDKgtIvsSIU71CDQxqhRp61b473+hcGErKRYoYFsoly9fpmrVqmTIkIFff/2VbNmy2RaLcp3HbWPcCUw0xuTBqi6N/4wASimVUBs3wmuvQYkS4Otr3THaJCwsDC8vL4KCgti9e7cmxSdEfBJjT2CB4yXAZmPMUeBQhNdxnfFGKfXYfv4ZWrSAsmVh0ybIk8fWcMaMGcPGjRuZNWsWVatWtTWWpDR9+nQAunfvbnMk9ohzrlQRuSAijYCCWHeMy4BAoDHW9HC/AjeNMQeNMd8mZrBKqRTsxx/h9dehYkWr+tTmpOjr68vw4cNp167dE9cJZc2aNaxZs8buMGyT0LlSVwJDRMTPsZwZqAJUffASkSfna5Ub0DZGlSIsXw5t21rTu61fD9mz2xrOhQsXqFq1Krly5WLfvn1kzpw50vbQUDh+3LaJd5QLuGyuVBFpGWX5FlYb5M5HD08p9UTz9oa334Y6dWDdOsiSxdZw7t+/z1tvvcWtW7fYunXrQ0kxJATefBN8fOCPP6BgQZsCVYnGZVPCKaVUgs2bB+++C56e1kOGM2WyOyKGDh3Kjh07WLRoEeXKlYu07e5deOMNWLsWpkxJuUlxypQpAPTu3dvmSOyRoOcxKqWUy8yaBZ06QaNGVqZxg6S4du1axo0bR5cuXR56SO+dO1Zn2bVrYeZM+OADm4JMAr6+vvj6+todhm0S1Mao3I+2MapkaepUK7M0bQorVkCGDHZHxJkzZ6hevTpFixZl9+7dZIgQ0+3b0KwZbNsGc+ZY+Vwlb7G1Meodo1IqaX3xhZUUmzeHVavcIikGBwfTokULRIQVK1ZESoo3b1oT8GzfDgsXalJ8Emgbo1Iq6YwZA4MHQ+vWsHix9bQMm4kIXbp04ejRo/z888+UKFHCue36dSsp7tsHS5ZAmzY2BpqEJk6cCEC/fv1sjsQemhiVUolPBEaOtF7t2sH8+ZDGPT5+pk6dypIlSxg1ahSvvPKKc/3Vq9b85UeOWKNJXn/dxiCT2O7du+0OwVbaxpjMaRujcnsiMGgQjBsHHTtajXSpU9sdFQA7duygQYMGNGnShB9++IFUqazWpaAgq0/QiROwcqU1batKWVw2jlEppRJEBPr1g0mToGtXmD4dUrlH14bz58/TunVrihcvzsKFC51J8dIlaNAATp+Gn36y7hrVk0UTo1IqcYSHQ+/e8PXX0KuXNfDPTR7XdO/ePVq1asWtW7fw9fV1Tg5+4YKVFP/6y5q2tX59mwO1ybhx4wD45JNPbI7EHpoYlVKuFx4O3brB7NnWHeOECW6TFAH69OnDnj17WL58OeXLlwcgIMBKhBcvwoYN8MILNgdpo8OHD9sdgq00MSqlXCsszJrNZsECq21x9Gi3Sopz585l5syZfPzxx7Rq1QqAc+espHjlivXUq9q17Y3RbkuXLrU7BFtpYlRKuc79+9ChgzW2YeRIGDrUrZLivn376N69Ow0aNOCzzz4D4NQpKynevGk9/rF6tN0x1JNEE6NSyjVCQ62hGMuXw9ix4GbtUxcuXKB58+Y89dRTLF26lDRp0uDvbyXFe/dg61Z9WsYDo0aNAqx5Y59EmhiVUo/v3j1r9Pvq1VYP1A8/tDuiSO7cuUPz5s25efMmPj4+5M6dm+PHrY42ItZUbxUq2B2l+/D397c7BFtpYlRKPZ67d6FlS+uRUV9/DT162B1RJA9mttm/fz8//vgjFStW5MgRaNjQmnhnyxZ45hm7o3QvixcvtjsEW2liVEo9uuBg65ETvr7W0zLee8/uiB4yYcIEvL29GT16NK+99hq//moN3s+UyUqKpUrZHaFyN5oYlVKP5tYt65ET27fD3LnWrDZuFgDTYQAAIABJREFUZu3atQwcOJA2bdowaNAg9uyBxo0he3arTbF4cbsjdE/Dhg0D4NNPP7U5EntoYlRKJdyNG9CkCezZY00G3rat3RE95Pjx47Rt25Zq1aoxd+5cfvnF0KQJ5M1r3SkWKWJ3hO4rICDA7hBspXOlJnM6V6pKcteuWbddv/4K330HjrGA7uTKlSvUrFmT27dvc+DAAU6dKsR//wuFClm1vgUL2h2hspvOlaqUco0rV+Cll+C336wHDL/2mt0RPSQ0NJTWrVvz999/s337dvz8CvHaa1a1qa8v5M9vd4TK3WliVErFz+XLVldOf3/48UerKtXNiAg9evRgy5YtzJ8/n6tXa9GiBZQpA5s3Q548dkeYPAwcOBCAsWPH2hyJPTQxKqXidvGiNejv7FlYs8bq1umGJk6cyOzZsxk0aBDZs3egeXOoWNGa5i1XLrujSz6uXLlidwi20sSolIrd+fPW9DDnz1tjFT097Y4oWitXruTjjz+mTZs2VKo0ilat4NlnrQnBs2e3O7rk5ZtvvrE7BFu5x4PRlFLu6a+/oF49CAwEHx+3TYr79u3Dy8uL2rVr07jxQtq2TUWtWtadoiZFlVB6x6iUit6ZM9ad4rVrsGkTPPec3RFF69y5czRr1oynnnqKt97aQKdO6ahXz6rxzZzZ7uiSp379+gFW1fSTSBOjUuphJ09aSTE42OrK+eyzdkcUrevXr9O0aVPu3btH9+6H+eCDrDRqZPUN8vCwO7rk686dO3aHYCtNjEqpyPz8rI42oaHWSPjKle2OKFqhoaG88cYb/PHHH7z//jFGjHiKJk1g5UrIkMHu6JK3adOm2R2CrbSNUSn1/44ds9oRw8OtR064aVJ8MCxj06ZNtGz5C1OnluG112DVKk2K6vFpYlRKWQ4dspJimjTW/Kfly9sdUYzGjBnD7NmzqVdvA8uW1eSNN6zHQKZPb3dkKUOfPn3o06eP3WHYRhOjUgr277faFD08rKRYpozdEcVo/vz5DBkyhEqVVrJ9+8u0bQtLlliPkFLKFbSNUakn3e7d1tynuXJZbYrFitkdUYzWr1/Pu+92pnhxb44ebUHHjjBnDqRObXdkKcvkyZPtDsFWeseo1JNsxw5r7tN8+aw7RTdOivv376dly1bkzj2Xs2fb8t578O23mhSV62liVOpJ5esL/9fencfpWK8PHP98s6WUrIccfrbU6dDBQcjPYxllyciShoiRJUspRTqaRNJxRCqFME3NVCZZx5rdQScHk3BsJdFizTLWGTPX74/7cX5PzJj1ee7nfu7r/Xrdr5r72a5xzT3XfJf7+23VytpyYt06KF/e7ogy9N1339G6dRsKFHiPY8eeYNAgmDoVbtLfYH4xcOBABg4caHcYttGuVKXcaPlyeOQRqFrVWl37D3+wO6IMHT16lAcfbElS0nguX+7BkCHw5ptgjN2Rha7ChQvbHYKttDAq5TaLFkHHjnDvvdaKNiVL2h1Rhs6dO0fr1m05dGgkqandeekleP11LYr+5tYVb67Sjgil3GTePOjQAe67z+pKDeKimJycTIcOnUlMfJbU1O68+qoWRRUYWhiVcov4eHj0UahTx+o+LV7c7ogylJqaSrdukaxY0RORrowdCyNHalEMlL59+9K3b1+7w7CNdqUq5QaxsdCzJzRqZHWl3nab3RFlSETo1+8ZZs/uBLRnwgQYMsTuqNylhMs3r9TCqFSoi46G3r2haVNYuBBuvdXuiG5o6NAoZs5sDbTh3Xdh0CC7I3KfN954w+4QbKVdqUqFsqlT4cknrXsVFy0K+qL42msTmDChMdCGqVNFi6KyhRZGpULVO+9A//7w8MPWPkxBPgX/nXdm8sortYAwZsxIo18/HVC0S2RkJJGRkXaHYRvtSlUqFI0fD8OGWTNQP/sMCha0O6Ibio7+gsGDqwENiYlJo0cP/dVkp/JBvNhDIOhPn1KhZswYiIqCxx6zJt0E+eran3/+JU8+WR74K3FxV3j8cd0iw26jR4+2OwRbaVeqUqFCBF55xSqK3btDXFzQF8W5c9cSEVECY2oTF3dZi6IKCtpiVCoUiMBLL8G4cdCrF3zwQdCvrr1gwSY6dSoO3E1c3AW6di1qd0jKq1u3bgDExcXZHIk9tDAq5XQi1o1+kyZZk20mTw761bUXLtxM+/Z3AJWZNesinTvfYXdIysfdQbwfZyBoYVTKydLS4Omn4f33YfBgeOutoF8eJiFhm7co/pHPPz9Pp07uvpk8GEVFRdkdgq20MCrlVGlp0K+ftVPv0KFWN2qQF8XFi7/lkUfuQKQ0c+acp337UnaHpNR1gru/RSmVvtRUiIy0iuLLLzuiKC5Zspvw8KKIlGTu3HNaFINYREQEERERdodhG20xKuU0V67AE09Y9yeOHm3NQg1yS5bso23b2xApwvz5SYSHl7M7JHUDNWvWtDsEW2lhVMpJkpOha1eYM8dqJQ4bZndEmZo/fw8dO94BFGDBgrO0bVvB7pBUJoYPH253CLbSwqiUU1y+bG0blZBgTbJ59lm7I8pUfPwuunQphTGwaNF5WrXSoqiCnxZGpZzg4kVrebdly+C992DAALsjylRs7A569LiTm25KYfnyKzRvrkXRKTp27AjAnDlzbI7EHloYlQp2Fy5AeDisXg3Tp1tbSAW5mTO307t3BfLlu8jq1dC48R/tDkllQ4MGDewOwVZaGJUKZufOWbtj/POfEBNjTboJclOmJDJgQGXy50/in//MT/36ZewOSWXTCy+8YHcIttLCqFSwOnMGWreGr7+GTz4BB0yfnzRpK889V40CBX5j06abqVPnD3aHpFS26X2MSgWjU6egRQvYvBni4x1RFMeO/RfPPXcPBQueoGHDfowe3cfukFQOhYeHEx4ebncYttEWo1LB5uRJqyju2mXdluGAX1AvvbSGv/+9Pjff/DNbtxZjxYpWdoekcqF58+Z2h2ArIyJ2x6ByoU6dOrJlyxa7w1B55dgxCAuDfftg/nxo2dLuiDI1YMBipkwJ49ZbD/Ptt3+gcuXb7A5JqUwZY7aKSJ30HtMWo1LB4tdfoXlzOHgQFi+2/j+IiQjdus3h00/DKVr0R3bt+iPlyhW2Oyylck0Lo1LB4KefoFkz+OUXWLoUPB67I7ohESE8PI5Fi7pQsuQP7N5dkZIl/39T5FatrK7UpUuX2hWiygW3508Lo1J2+/FHqygePw5ffgkNG9od0Q2lpqbSvPmHrFsXSdmyB9m9uwpFi/5+Hl/btm1tik7lBbfnT8cYHU7HGB3uwAFo2hTOnoXly6FePbsjuqGLFy/SqNFMtm0bQMWKB9mxoxJFigT3rh5KpUfHGJUKRvv2WS3Fixdh1SqoXdvuiG7o5MmT1Ks3kwMHhnHPPT+ydWtlbrnF7qiUynt6H6NSdvjPf6xxxORkWLs26IviwYMHuffeKRw4MIzatX8lMfF/blgUw8LCCAsLC1yAKk+5PX/aYlQq0L791rolI18+qyjee6/dEd1QYmIiHs9CkpJG8r//e4KVK8tSsOCNX/PYY48FJjjlF27Pn44xOpyOMTrMtm3WzfuFC1uLglerZndEN7R8+XLCw/9FcvJIWrY8y8KFt1OgQOavUyrY3WiMUbtSlQqUzZutexOLFIF164K+KMbEfESrVl+TnDySDh0ukJCgRVG5gxZGpQJh0yar+7R4cVi/HqpUsTuiDKWlpfHyy1FERv6CyCt065bM55/fQv5sDLw0adKEJk2a+C1G5V9uz5+OMSrlb+vWQZs2UK6cNfv0j8G7N+GFCxd44okezJnTABhCnz6pTJ1akJuy+Sd0z549/RGeChC350/HGB1OxxiD3MqV1iLgFStaRbFsWbsjytAvv/xCePgjbN3aAxjI008Lb79tMHqbogpBOsaolB2WLbM2Ga5a1Zp9GsRFcdu2bdStez/bt/cHBvLCC+SqKKakpJCSkpKnMarAcXv+tDAq5Q8JCdCunXUrxpo1ULq03RFlaN68eTRq5OH06YlcuRLJiBHwj3+Qq5ZiixYtaNGiRd4FqQLK7fnTMUal8tqcOdbGwrVrW63GYsXsjihdIsLYsWN5+eWRlCixiJMnWzJ6NERF5f69e/funfs3UbZxe/50jNHhdIwxyHz2GXTvDvffD0uWQNGidkeUrqSkJHr06MG8eQmUL7+Bw4fv5403YPhwuyNTKjB0rVSlAuGjj6BXL2jUCBYtgtuCc8PevXv30r59e/buPUj16rvZubMqEyfCc8/l3WdcuHABgFt0MVVHcnv+dIxRqbwwYwZERlo7ZSxZErRFMSEhgXr16nHs2Fnq1j3Ezp1VmTw5b4siQOvWrWndunXevqkKGLfnT1uMSuXW++/DwIHQsiXMnWst9xZk0tLSeO2113j11VepWbMhRYqsZOPGwnzwAfTpk/ef179//7x/UxUwbs+fjjE6nI4x2mzSJKu5FR4On38OhQrZHdF1Tp06RY8ePUhISKBLlz78/PMUNmzIR3Q09Ohhd3RK2UPvY1TKH8aNs4pix44we3ZQFsXNmzdTq1Ytli1bxrhxUzh0aBobN+YjNta/RfHMmTOcOXPGfx+g/Mrt+dPCqFROvPaaNYUzIgJmzSLTfZgCTESYNGkSjRo1AmDp0q+YM+cpvv7aMGsWdO3q389v164d7dq18++HKL9xe/50jFGp7BCBV16BMWPgiScgOtraVzGInD59ml69ejFv3jzCw8OZMCGGxx4rxo4d8MUX1roD/vbMM8/4/0OU37g9fzrG6HA6xhhAIvDiizB+PPTuDdOmke3Vtf1sy5YtdO7cmcOHDzNu3Dgef/w5HnzQsHcvzJsHrVrZHaFSwUHHGJXKLRFrPHH8eBgwIOiKYlpaGm+++SYNGzbkypUrrF+/ni5dhtCsmWH/fuu2ykAWxRMnTnDixInAfaDKU27Pn3alKpWZtDTrdoypU+HZZ2HixNwtJJrHfvrpJ3r06MHq1atp374906dP59KlEjRpAj//bN1WGeit9Tp16gTA2rVrA/vBKk+4PX9aGJW6kdRU6NvXGkt88UV4442gKopz5syhT58+XL58mRkzZtCrVy8OHzY0awbHjsHy5fDAA4GP6/nnnw/8h6o84/b86Rijw+kYox9duWIt8RYba024efXVoCmK586dY/DgwURHR1O3bl0++eQT7rrrLg4cgGbN4PRpqyjef7/dkSoVnHSMUansSkmBbt2sojhmDIwaFTRFcePGjdSqVYsPP/yQESNGsHHjRu666y727wePB5KSYPVqe4vikSNHOHLkiH0BqFxxe/60K1WpayUnW/cnzptnbUw4dKjdEQHWws5RUVG89dZbVKhQgbVr19K4cWMAdu+2WopXrljbP953n72xRkREAO4do3I6t+dPC6NSvi5dgkcftaZxTpoEgwfbHRFgtRIjIyPZv38//fv3Z9y4cdzmXah8xw5o3tyaJLt2Lfz5z/bGCjBc969yNLfnTwujUlddvAjt21uDc1OmwFNP2R3Rda3EVatW0axZs/8+npgILVpYq9GtXg13321jsD5atmxpdwgqF9yePy2MSgGcP28tBL5mDcycaU26sdn69evp06cP+/btu66VCPDvf8ODD8Ltt1tFsUoVG4O9xuHDhwEoX768zZGonHB7/rQwKpWUBG3awMaN8PHH1qQbG508eZJhw4YRHR1NpUqVrmslAmzaZN2wX6KEVRQrVrQn1ox0794dcO8YldO5PX9aGJW7nTljVZjNm+HTT+Gxx2wLRUSIjY3l+eef5/Tp0wwfPpyoqKjrdlFfv96q42XLwqpVEIx/1L/88st2h6Bywe3508Ko3OvUKasvcvt2ay/FDh1sC2Xfvn089dRTrFmzhgYNGjBt2jRq1Khx3fNWrbJ6fCtUsFqKZcvaEGwWhIWF2R2CygW350/vY1TudOKEdX/Dt9/C3Lm2FcXz588TFRVFjRo12LZtG1OnTmXDhg3pFsXly+Hhh6FyZWv2abAWRYADBw5w4MABu8NQOeT2/GmLUbnP0aMQFgbffQcLF8JDDwU8BBHhs88+Y9iwYfz888907dqVCRMmUKZMmXSfn5AAnTrBvffCihVQsmSAA86mXt7JS24do3I6t+dPC6Nyl19+sW76O3QIFi+2Wo0BtnXrVgYPHszGjRupXbs28fHxPHCDBU3nzrWGPmvVslqNxYoFMNgcGjVqlN0hqFxwe/60MCr3OHzYKoRHjsDSpeBdNSZQjh49yogRI4iOjqZUqVLMmDGDnj17ku8GGx3Hx8Pjj0O9elbIRYsGMOBc8Hg8doegcsHt+dPCqNzh4EGrKJ48CV9+CQ0aBOyjk5KSmDhxIm+++SaXLl1iyJAhREVFUTSTKhcbCz17QqNG1kI8PrcwBr29e/cCcHewrDigssXt+dPCqELf999bRfHsWVi5EurWDcjHJicnM336dEaPHs2xY8fo1KkTr7/+OtWqVcv0tdHR0Ls3NG1qDYPeemsAAs5D/fr1A9w7RuV0bs+fFkYV2vbutYri5cvWqjY1a/r9I9PS0pg9ezYjRozg+++/x+PxsHDhQu7P4nYXU6bAgAHWnKB586BwYT8H7Adjx461OwSVC27PnxZGFbp27bIm2ohY9zdUr+7Xj0tLS2PBggWMGjWK7du3U6NGDRYvXkyrVq0wWdyy6u234dlnrdsyZs+Gm2/2a8h+07BhQ7tDULng9vzpfYwqNG3fDk2a/P+WE34simlpacydO5datWrRoUMHLly4wMcff0xiYiKtW7fOclEcP94qih06wJw5zi2KADt37mTnzp12h6FyyO350xajCj3btllbTtxyi7U8zF13+eVj0tLSmDdvHqNGjWLHjh1Uq1aN2NhYIiIiyJ8/e5fWmDEQFWVtA/nxx1CggF9CDphBgwYB7h2jcjq3508LowotX39tDc7dcYc1plipUp5/xKVLl4iNjWXChAns3buXatWqERcXR0RExA1vvUiPCIwcCa+9Bt27w4cfQjbfIiiNHz/e7hBULrg9f1oYVejYsAFat4bSpa2WYoUKefr2v/32G1OmTOHdd9/l6NGj1KpVi08//ZTOnTtnuyCCVRSHD4d//AOefBKmTQuNoghQN0Azf5V/uD1/WhhVaFi71pqxUq6cVRTLlcuzt96/fz+TJ09m5syZnD9/noceeoihQ4fSrFmzLI8fXksEhgyBSZOgf3+YPNkaDg0V33zzDQA1AzALWOU9t+dPC6NyvhUroF07q9t01SrIYL3R7EhNTWXx4sW8//77LF++nPz589OlSxdeeOEF7rvvvly//969MHUqDB4Mb70FOayvQevZZ58F3DtG5XRuz58WRuVsS5ZY0zjvvtu6eb9UqVy93fHjx5k5cyZTp07lxx9/5M4772TUqFH06dOHsnm4ncU998A330C1aqFXFAEmTZpkdwgqF9yePy2MyrkWLIBHH4UaNaxl3kqUyNHbXLlyhWXLlhETE8PChQtJSUmhadOmTJgwgfDwcAr4aYpoKK+25dYuuFDh9vxpYVTONHs2dO0Kf/0rLFtmzULNpp07dxITE0NcXBxHjx6lVKlSDBw4kL59+/KnP/3JD0G7x7///W9AJ3E4ldvzp4VROc+nn1r3NjRoYHWl3n57ll/6ww8/MHv2bOLj49m2bRv58+enbdu29OzZk1atWvmtdeg2Q4cOBdw7RuV0bs+fFkblLB99BJGR4PFYu/cWKZLpS64Ww9mzZ7NlyxYA6tWrx6RJk+jatSulcjkuqa43efJku0NQueD2/AVFYTTGlAPeBf4XKAFMAgoAFUWkrZ2xZZUxZgbwJPCWiAzJ4Xs8B0QCNUUkLS/jCwnTp0O/fhAWBvPnWyvbpCMtLY0tW7awaNEiFi9ezLZt2wCrW2j8+PF06tSJihUrBjBw96nu53VplX+5PX9GROyOAWPMCuBOYCjwG5APWAM0FJEtdsaWFcaYwsAR4DbgOFBORK7k8H1+AF4SkQ+z8po6derI1VZQSHvvPRg0yLqBP52FRE+fPs3KlStZtGgRS5cu5dixY9x0003Ur1+fdu3a8eijj1LJD6vgqPRt2rQJ0MWoncoN+TPGbBWROuk9ZnuL0RhTBmgO9BSRJd5z7wLbnVAUvdoDtwPjgBeBlsCi7L6JiFw0xnwMvABkqTC6wsSJ8Pzz1r2K8fFQqBBnz55lw4YNrFmzhjVr1pCYmEhaWhrFihWjZcuWtGnThpYtW1IihzNVVe787W9/A9w7RuV0bs+frS1GY8xcrKLi6+/AU8BoEXnL57lVgd3AWBEZ6XN+CtANaGpXITXGLAfuAu4GfgbWikjna56TpfiNMbWBrcADIrIps88O+Rbj3/8OL73EuVatWNSlC//aupWvvvqKrVu3kpqaSsGCBalfvz5NmzYlLCyM+vXrZ3sBb5X33L4DvNO5IX83ajEiIrYdwJ+B6UASUN97eAAB6qTz/CnAWaCk9+tXgMtAmI3fw51AKlYhB3gHuAQUy0n8WFuBnbn6fpkdRYoUkQ8//FBERJKTk8Xj8UhsbKyIiJw/f148Ho/MmjVLREROnz4tHo9H5syZIyIix48fF4/HIwsXLhQRkV9//VU8Ho8sXbpUREQOHTokHo9HVqxYISIi33//vXg8Hlm7dq2IiOzZs0c8Ho9s3LhRRER27NghHo9HNm/eLCIiiYmJ4vF4JDExUURENm/eLB6PR3bs2CEiIhs3bhSPxyN79uwREZG1a9dK48aNZeXKlfLF7Nky/r77REC+KFRI8lk/E1K4cGFp1KiRjBgxQlauXCkXLlwQpZTKLmCLZPB71dY/rUVklzHmDuBbEfkXgDHmRaxfgt+m85JRwBPAi8aYPcBIoIuIrAxUzOnojlXM4rxffww8DTwGTL3muZnGLyJpxphvsf5ISJcxpi/QF6BQoUJ59G0E1unTp9mzZw/Hjx9n2rRpnDp1ik2bNvHdd98RFhZGYWAAMObWWznQsSOTGzTg/vvvp3r16npLhQOsW7cOAI/HY3MkKifcnj/bJ98YY/YDy0VkkPfrt4EeIpLuHdvGmNeB57HGRweLyHvXPF4F+AgoDZwH+kgGXazGmDBgRRbCXCciTTJ4j11AkojU9zm3GzglIteNXGcWv/c5c4FqIpLp1LBg6UpNS0sjKSmJs2fPcubMGc6cOcOJEyf49ddfrzt++OEHfvvtt9+9vkyZMlSvXp0aNWr89xjyzDPcVLCga8c5nKxJkyaAe8eonM4N+QvayTfGmNuAKkCiz+mbsboXM7IfKARsSK+oYLXSYkRkhjGmBfCJMeYeSf8vgE1AVpY4uZBB/HWBe4FB1zwUC7xujKkmIvuyGT/ARaBwFuIKqIiICL755htSUlJ+dyQnJ3Pu3Dky+iPLGEPp0qUpU6YMZcuWpW7dulSpUuW/R+XKlSmSzv2IMXFx6bybcoLo6Gi7Q1C54Pb82T1LoSZg+H1hPAkUS+/JxphmwDTgK+ABY8xfRGS7z+OlsLogWwOIyArvtkB/Ba5rVonIBWBPLuLvAaQA8decjwPGYHWbvpzV+H0UB07kIi6/qFSpEiJCgQIFrjtuv/12ihYt+t/jjjvuoHjx4pQtW5bSpUvnaEJM5cqV/fBdqEDQ3Dmb2/Nnd2GshVVYdvmc2wMUMMb8UUR+unrSO1tzPjADeA7YB4wF2vi8tgLwi4ik+Jz70Xs+T/sbjTEFgQhgqYj8roiJyCFjzDqguzEmSkQki/FfVQnYnJfx5oU33ngjoJ+3cqU19BoWFhbQz1W5p7lzNrfnLxgK439ExLfrdL33v/WAn+C/tzosBb4EnvZOUBkFRBtjGovIejLmr019HsZapeewMeaRdB4/ADQBmhhjDpPF+L2TkaoBb/opbscYM2YM4N6L08k0d87m9vzZfR/jN8A3ItLzmvNfYxXMSO8CAJuAQ8BDV4uoMSYfsBOfSS7ertQDQPGrrUZjzD6ga0YTcHIR+wIgPAtPXY5V6DKN33v+cWAm1uo5JzN782CZfOMPhw8fBqB8+fI2R6KyS3PnbG7I340m39g+KzU9xpiewNtAWe84YHZeuwqYJSLTvZNv3sea4Rl832g6jDFLgRMi0j0rzw/lwqiUUv5yo8J4U6CDyaJYrBVkBuTgtU8Bkd6W4njgcQcVxZpAU6z7HV1v2bJlLFu2zO4wVA5o7pzN7fkLyhYjgDGmPlBbRN63O5ZAMca0xFox57OsviaUW4xuuJcqVGnunM0N+XNcV6rKulAujEeOHAGsm/+Vs2junM0N+QvaG/yVupFQvihDnebO2dyev2AdY1SKhIQEEhIS7A5D5YDmztncnj/tSnW4UO5KdcM4R6jS3DmbG/KnXanKkb744gu7Q1A5pLlzNrfnTwujClolS5a0OwSVQ5o7Z3N7/nSMUQWtuXPnMnfuXLvDUDmguXM2t+dPxxgdTscYVTDS3DmbG/Kn9zGGsFAujGfOnAGgaNGiNkeisktz52xuyJ9OvlGOFMoXZajT3Dmb2/OnY4wqaMXHxxMff+0e0MoJNHfO5vb8aVeqw4VyV6obxjlClebO2dyQPx1jDGGhXBgvXLB2HLvllltsjkRll+bO2dyQPx1jVI4UyhdlqNPcOZvb86djjCpoxcXFERcXZ3cYKgc0d87m9vxpV6rDhXJXqhvGOUKV5s7Z3JA/HWMMYaFcGFNSUgAoUKCAzZGo7NLcOZsb8qdjjMqRQvmiDHWaO2dze/50jFEFrZiYGGJiYuwOQ+WA5s7Z3J4/7Up1uFDuSnXDOEeo0tw5mxvyp2OMIcwYcxz40e44/KgkcMLuIFSOaO6cLdTz9z8iUiq9B7QwqqBmjNmS0V91Krhp7pzNzfnTMUallFLKhxZGpZRSyocWRhXsPrA7AJVjmjtnc23+dIxRKaWU8qEtRqWUUsqHFkallFLKhxZGpZRSyocWRhVSjDEzjDFijJlodyxuYIypaoxJMcaMuuZHAvmqAAAE8ElEQVT8FGNMkjHGlffBuVEoXXs6+UaFDGNMYeAIcBtwHCgnIlfsjSr0GWOmAI8DlUXkhDHmFWAE0EZEVtobnQqEULv2tMWoQkl74HbgH0BpoKW94bjGKCAf8KIx5klgJNBdi6LzGGMOGmNezcFLQ+ra08KosiXIu856AD8AUVh/tT5x7ROCPH5HEpEjwCTgaWAa8IyIfH71cWNMFWPMBmPMPmNMov4b50yQ/+yG1rUnInroka0DmAKcBUp6v34FuAyE2RjTnUAqMNr79TvAJaCYE+J3+gH0BAT4ZzqPrQB6e/+/BbAX7zCOHtn+d/brzy5wEHg1m68JuWvP9gD0cN4BlAHOA+OBJ70XRWebY3rR+4u5mvfrOt6vn3JC/E4+gGbeX26bgDTgLz6PlQKSgAI+5/YCdeyO24lHXv7sAgZrs3rf4yAw+ppz+TJ5n5C79mwPQA9nHsDr3r8KrwADr3msCrAB2Ack3uiXIBDmvYgyO9ZmEs8u4F/XnNsNbMpu/Hpk6+egtrcF8B5Q0PuLdbHP438F9l7zmi+BDnbH7tQjD6+9JnrtpX/kR6mc2Q8UAjaIyHvXPDYViBGRGcaYFsAnxph7xHtVXGMT8KcsfN6FjB4wxtQF7gUGXfNQLPC6MaaaiOzLRvwqC4wxVYGlWIXuaRFJ844fRRtjGovI+oxeGrAgQ1NeXXtbgbrXnFsILOL366QmZRRIyF57dldmPZx3EGRdZ8BkIBnvuIXP+Qre+MZkNX49svxvXgY4AKwFCvmcz4dPayGDn4d9/vx5COXD39ce2RxjDNVrT2elqmwxxtQG5gMzsLpiDgFjfZ5SAfhFRFJ8zv3oPe+PeAoCEcBSEfndbuMicghYB3Q3xpgsxq+yQESOiEhlEWkiIpd9zqeKyJ9EpKH36+PAZqzJOXhbMQartaKyQa+9wNGuVJVlQdp19jBQAjhsjHkknccPYF2ETYwxh8lZ/Cp3ngI+MsYMxeoSf1y8zQeVNXrtBZaufKOyxBhTBqv74xDw0NVWgjEmH7ATOCUiDY0xpbAuiOJX/3I1xuwDuorIFj/EtQAIz8JTlwPVMos/r+NTKrcCee0ZYw5ijVG+moXnhuy1p4VR5TljzCpglohM93advY81lVt/2JTyI7328oYWRpXnjDF3AR8BJbG6zvqKyGZ7o1Iq9Om1lze0MCqllFI+dFaqUkop5UMLo1JKKeVDC6NSSinlQwujUkop5UMLo1JKKeVDC6NSSinlQwujUkop5UMLo1JKKeVDC6NSSinlQwujUipPGGOqGmNSvLsm+J6fYoxJMsbUsSs2pbJDC6NSKk+IyHdYe+09Z4wpCWCMeQXoBbT3x+4qSvmDrpWqlMoz3i2Svsfa1WEP8AHQRUQ+tzUwpbJBNypWSuUZETlijJkEPI/1++UZLYrKabQrVSmV1/YDhYCvROQ9u4NRKru0MCql8owxphkwDfgKeMAY8xebQ1Iq27QwKqXyhDGmNjAfawJOE+AQMNbOmJTKCS2MSqlcM8ZUBZYCXwJPi0gyMApobYxpbGtwSmWTzkpVSuWKdybqJqwW4kMictl7Ph+wEzglIg1tDFGpbNHCqJRSSvnQrlSllFLKhxZGpZRSyocWRqWUUsqHFkallFLKhxZGpZRSyocWRqWUUsqHFkallFLKhxZGpZRSyocWRqWUUsrH/wGrSAsmQVwWaAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(6, 6))\n",
    "ax1 = plt.subplot(111)\n",
    "ax1.set_xlim(-0.1, 1.1)\n",
    "ax1.set_ylim(-0.1, 1.1)\n",
    "# ax.grid(True)\n",
    "ax1.set_title('Central difference example', fontsize=16)\n",
    "ax1.set_xlabel('$x$', fontsize=16)\n",
    "ax1.set_ylabel('$f(x)$', fontsize=16)\n",
    "# define our x for plotting purposes\n",
    "x = np.linspace(0, 1, 1000)\n",
    "\n",
    "# define our example function and its exact derivative\n",
    "\n",
    "def f(x):\n",
    "    return x**2\n",
    "\n",
    "def df(x):\n",
    "    return 2 * x\n",
    "\n",
    "# plot the exact solution\n",
    "ax1.plot(x, f(x), 'k')\n",
    "# choose and plot two x locations to take the difference between\n",
    "dx = 0.4\n",
    "x0 = 0.5\n",
    "xl = x0 - dx\n",
    "xr = x0 + dx\n",
    "# plot a line representing the discrete derivative\n",
    "ax1.plot([xl, xr], [f(xl), f(xr)], 'r', label = 'Central diff. approx. deriv.')\n",
    "# plot a line representing the exact derivative at x=x0\n",
    "h = dx / 2\n",
    "ax1.plot([x0 - h, x0 + h], [f(x0) - h * df(x0), f(x0) + h * df(x0)], 'b', label = 'Exact derivative')\n",
    "# add some axes labels and lines etc\n",
    "ax1.set_xticks((xl, x0, xr))\n",
    "ax1.set_xticklabels(('$x_0-\\Delta x$', '$x_0$', '$x_0+\\Delta x$'), fontsize=16)\n",
    "ax1.plot([xl, xl], [-0.1, f(xl)], 'k:')\n",
    "ax1.plot([xr, xr], [-0.1, f(xr)], 'k:')\n",
    "ax1.plot([x0, x0], [-0.1, f(x0)], 'k:')\n",
    "ax1.set_yticks((f(xl), f(xr)))\n",
    "ax1.set_yticklabels(('$f(x_0-\\Delta x)$', '$f(x_0+\\Delta x)$'), fontsize=16)\n",
    "ax1.plot([-0.1, xl], [f(xl), f(xl)], 'k:')\n",
    "ax1.plot([-0.1, xr], [f(xr), f(xr)], 'k:')\n",
    "ax1.legend(loc='best', fontsize=14);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## <span style=\"color:blue\">Exercise 2.2: Compute first derivative using central differencing</span>\n",
    "\n",
    "Use the data below to compute $f'(0.2)$ using central differencing:\n",
    "\n",
    "$$f(0.1) = 0.078348$$\n",
    "$$f(0.2) = 0.138910$$\n",
    "$$f(0.3) = 0.192916$$\n",
    "\n",
    "You should get 0.57284"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## <span style=\"color:blue\">Example: Write a function to perform numerical differentiation</span>\n",
    "\n",
    "As covered above, the expression\n",
    "\n",
    "$$\\frac{f(x+\\Delta x) - f(x-\\Delta x)}{2\\Delta x},$$\n",
    "\n",
    "can be used to find an approximate derivative of the function $f(x)$ provided that $\\Delta x$ is appropriately small. \n",
    "\n",
    "Let's write a function `diff(f, x, dx = 1.0e-6)` that returns the approximation of the derivative of a mathematical function represented by a Python function `f(x)`.\n",
    "\n",
    "Let's apply the above formula to differentiate $\\,f(x) = e^x\\,$ at $\\,x = 0$, $\\,f(x) = e^{−2x}\\,$ at $\\,x = 0$, $\\,f(x) = \\cos(x)\\,$ at $\\,x = 2\\pi$, and $\\,f(x) = \\ln(x)\\,$ at $\\,x = 1\\,$, i.e. functions we know the exact derivative of.\n",
    "\n",
    "In each case, using $\\,\\Delta x = 0.01$, let's write out the error, i.e. the difference between the exact derivative and the result of the formula above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The approximate derivative of exp(x) at x = 0 is: 1.000017. The error is 0.000017.\n",
      "The approximate derivative of exp(-2*x) at x = 0 is: -2.00013.  The error is 0.00013.\n",
      "The approximate derivative of cos(x) at x = 2*pi is: 0.00000.  The error is 0.00000.\n",
      "The approximate derivative of ln(x) at x = 0 is: 1.00003.  The error is 0.00003.\n"
     ]
    }
   ],
   "source": [
    "def diff(f, x, dx=1e-6):\n",
    "    numerator = f(x + dx) - f(x - dx)\n",
    "    derivative = numerator / ( 2.0 * dx )\n",
    "    return derivative\n",
    "\n",
    "dx = 0.01\n",
    "x = 0\n",
    "f = np.exp\n",
    "derivative = diff(f, x, dx)\n",
    "print(\"The approximate derivative of exp(x) at x = 0 is: %f. The error is %f.\"\n",
    "      % (derivative, abs(derivative - 1)))\n",
    "x = 0\n",
    "\n",
    "def g(x):\n",
    "    return np.exp(-2*x)\n",
    "\n",
    "derivative = diff(g, x, dx)\n",
    "print('The approximate derivative of exp(-2*x) at x = 0 is: {0:.5f}.  The error is {1:.5f}.'\n",
    "        .format(derivative, abs(derivative - (-2.0))))\n",
    "\n",
    "x = 2*np.pi\n",
    "f = np.cos\n",
    "derivative = diff(f, x, dx)\n",
    "print('The approximate derivative of cos(x) at x = 2*pi is: {0:.5f}.  The error is {1:.5f}.'\n",
    "        .format(derivative, abs(derivative - 0)))\n",
    "\n",
    "x = 1\n",
    "f = np.log\n",
    "derivative = diff(f, x, dx)\n",
    "print('The approximate derivative of ln(x) at x = 0 is: {0:.5f}.  The error is {1:.5f}.'\n",
    "        .format(derivative, abs(derivative - 1)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## <span style=\"color:blue\">Exercise 2.3: Compute the derivative of $\\sin(x)$</span>\n",
    "\n",
    "Compute \n",
    "\n",
    "$$\\frac{d(\\sin x)}{dx}\\qquad\\textrm{at}\\qquad x = 0.8$$\n",
    "\n",
    "using (a) forward differencing and (b) central differencing. \n",
    "\n",
    "Write some code that evaluates these derivatives for decreasing values of $h$ (start with $h=1.0$ and keep halving) and compare the values against the exact solution.\n",
    "\n",
    "Plot the convergence of your two methods."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "You should get something that looks like this\n",
    "\n",
    "![\"Convergence plot\"](images/fd_cd_convergence.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Approximating second derivatives\n",
    "\n",
    "Numerical differentiation may be extended to the second derivative by noting that the second derivative is the derivative of the first derivative. \n",
    "\n",
    "That is, if we define a new function $g$ where:\n",
    "\n",
    "$$ g(x) = f'(x), $$\n",
    "\n",
    "then\n",
    "\n",
    "$$ f''(x) = g'(x), $$\n",
    "\n",
    "a consequence of this (obvious) observation is that we can just apply our differencing formula twice in order to achieve a second derivative.\n",
    "\n",
    "And so on for even higher  derivatives!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We have noted above that the central difference method, being second-order accurate, is superior to the forward difference method so we will choose to extend that.\n",
    "\n",
    "In order to calculate $f''(x_0)$ using a central difference method, we first calculate $f'(x)$ for each of two half intervals, one to the left of $x_0$ and one to the right:\n",
    "\n",
    "\\begin{align*} \n",
    "f'\\left(x_0+\\frac{\\Delta x}{2}\\right) & \\approx \\frac{f(x_0+\\Delta x)-f(x_0)}{\\Delta x},\\\\[10pt]\n",
    "f'\\left(x_0-\\frac{\\Delta x}{2}\\right) & \\approx \\frac{f(x_0)-f(x_0-\\Delta x)}{\\Delta x}.\n",
    "\\end{align*}\n",
    "\n",
    "[NB. We make this choice since we have a vision to ultimately end up with an approximation that utilises a three-point stencil].\n",
    "\n",
    "Of course the formula on the RHS of these two equations we recognise from above as first-order forward and backward differences, if we were to consider the derivatives on the LHS to be evaluated at $x_0$. \n",
    "\n",
    "However, by considering the LHS at $x_0\\pm \\Delta x/2$ they are in actual fact second-order *central* differences where the denominator of the RHS is $2\\times (\\Delta x/2)$.\n",
    "\n",
    "[This may seem a bit weird, but is just showing that there are sometimes more than one way to interpret a numerical approximation].\n",
    "\n",
    "We can now calculate the second derivative making use of these two values within another central difference. \n",
    "\n",
    "But we must note again that these two values are telling us $f'(x)$ at the points $x_0\\pm{\\Delta x}/{2}$, which are only $\\Delta x$ rather than $2\\Delta x$ apart. \n",
    "We must remember this in the denominator of the central difference formula to yield\n",
    "\n",
    "\\begin{align}\n",
    "    f''(x_0)&\\approx\\frac{f'(x_0+\\frac{\\Delta x}{2})-f'(x_0-\\frac{\\Delta x}{2})}{\\Delta x}\\\\[5pt]\n",
    "    &\\approx\\frac{\\frac{f(x_0+\\Delta x)-f(x_0)}{\\Delta x}-\\frac{f(x_0)-f(x_0-h)}{\\Delta x}}{\\Delta x}\\\\[5pt]\n",
    "    &\\approx\\frac{f(x_0+\\Delta x)-2f(x_0)+f(x_0-\\Delta x)}{\\Delta x^2}.\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## <span style=\"color:blue\">Exercise 2.4: Compute second derivative</span>\n",
    "\n",
    "Calculate the second derivative $f''$ at $x = 1$ using the data below:\n",
    "\n",
    "$f(0.84) = 0.431711$\n",
    "\n",
    "$f(0.92) = 0.398519$\n",
    "\n",
    "$f(1.00) = 0.367879$\n",
    "\n",
    "$f(1.08) = 0.339596$\n",
    "\n",
    "$f(1.16) = 0.313486$\n",
    "\n",
    "You should get 0.36828"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Aside: Non-central differencing and differentiation by polynomial fit\n",
    "\n",
    "In this particular case we were given more data than we actually used. An alternative approach would be to use *non-centred differencing*, e.g. the following is also a valid approximation to the second derivative\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "    f''(x_0)\\approx\\frac{f(x_0+2\\Delta x)-2f(x_0+\\Delta x)+f(x_0)}{\\Delta x^2}\n",
    "\\end{align}$$\n",
    "\n",
    "This can come in handy if we need to approximate the value of derivatives at or near to a boundary where we don't have data beyond that boundary.\n",
    "\n",
    "If we wanted to use all of this data, an alternative would be to fit a polynomial to this data, and then differentiate this analytical expression exactly to approximate the derivative at any point between 0.84 and 1.16 (recalling that extrapolation is dangerous)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Numerical methods for ODEs\n",
    "\n",
    "One of the most important applications of numerical mathematics in the sciences is the numerical solution of ordinary differential equations (ODEs). This is a vast topic which rapidly becomes somewhat advanced, so we will restrict ourselves here to a very brief introduction to the solution of first order ODEs. A much more comprehensive treatment of this subject is to be found in the Numerical Methods 2 module.\n",
    "\n",
    "Suppose we have the general first-order ODE:\n",
    "\n",
    "\\begin{align}\n",
    "u'(t)&=f(u(t),t) \\\\\n",
    "u(t_0)&=u_0\n",
    "\\end{align}\n",
    "\n",
    "[Notation: For $u=u(t)$, $\\frac{du}{dt}=u'=\\dot{u}$.]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "That is, the derivative of $u$ with respect to $t$ is some known function of $u$ and $t$, and we also know the initial condition of $u$ at some initial time $t_0$.\n",
    "\n",
    "If we manage to solve this equation analytically, the solution will be a function $u(t)$ which is defined for every $t>t_0$. In common with all of the numerical methods we will encounter in this module, our objective is to find an approximate solution to the ODE at a finite set of points. In this case, we will attempt to find approximate solutions at $t=t_0,\\,t_0+\\Delta t,\\,t_0+2\\Delta t\\,,t_0+3\\Delta t,\\,\\ldots$.\n",
    "\n",
    "It is frequently useful to think of the independent variable, $t$, as representing time. A numerical method steps forward in time units of $\\Delta t$, attempting to calculate $u(t+\\Delta t)$ in using the previously calculated value $u(t)$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Euler's method\n",
    "\n",
    "To derive a numerical method, we can first turn once again to the Taylor\n",
    "series. In this case, we could write:\n",
    "\n",
    "$$ u(t+\\Delta t)=u(t)+\\Delta t u'(t) + O(\\Delta t^2) $$\n",
    "\n",
    "Using the definition of our ODE above, we can substitute in for $u'(t)$:\n",
    "\n",
    "$$ u(t+\\Delta t)=u(t)+ \\Delta t f(u(t),t)+ O(\\Delta t^2).$$\n",
    "\n",
    "Notice that the value of $u$ used in the evaluation of $f$ is that at time $t$. This simple scheme is named **Euler's method** after the 18th century Swiss mathematician, Leonard Euler."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Dropping the $O(\\Delta t^2)$ terms, we have a scheme we can easily implement in code:\n",
    "\n",
    "$$ u(t+\\Delta t)=u(t)+ \\Delta t f(u(t),t)$$\n",
    "\n",
    "\n",
    "This is what is known as an explicit method, because the function $f$ in this relation is evaluated at the old time level $t$\n",
    "-- i.e. we have all the information required at time $t$ to explicitly compute the right-hand-side,\n",
    "and hence easily find the new value for $u(t+\\Delta t)$.\n",
    "\n",
    "This form of the method is therefore more correctly called either Explicit Euler or Forward Euler.  We could also evaluate the RHS at some time between $t$ and $t+\\Delta t$ (in the case of $t+h$ this method is called Implicit or Backward Euler) this is more complex to solve for the new $u(t+\\Delta t)$ but can have advantageous accuracy and stability properties."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The formula given is used to calculate the value of $u(t)$ one time step forward from the last known value. The error is therefore the local truncation error. If we actually wish to know the value at some fixed time $T$ then we will have to calculate $(T-t_0)/h$ steps of the method. This sum over $O(1/\\Delta t)$ steps results in a global truncation error for Euler's method of $O(\\Delta t)$.\n",
    "\n",
    "In other words, Euler's method is only first-order accurate -- if we halve $\\Delta t$ we will need to do double the amount of work and the error should correspondingly halve; if we had a second-order method we would expect the error to reduce by a factor of 4 for every doubling in effort!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "To illustrate Euler's method, and convey the fundamental idea of all time stepping methods, we'll use Euler's method to solve one of the simplest of all ODEs:\n",
    "\n",
    "$$ u'(t)=u(t),$$\n",
    "$$ u(0)=1.$$\n",
    "\n",
    "We know, of course, that the solution to this equation is $u(t)=e^t$, but let's ignore that for one moment and evaluate $u(0.1)$ using Euler's method with steps of $0.05$. The first step is:\n",
    "\n",
    "$$\\begin{align}\n",
    "  u(0.05)&\\approx u(0)+0.05u'(0)\\\\\n",
    "  &\\approx1+.05\\times1\\\\\n",
    "  &\\approx 1.05\n",
    "\\end{align}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Now that we know $u(0.05)$, we can calculate the second step:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "  u(0.1)&\\approx u(0.05)+0.05u'(0.05)\\\\\n",
    "  &\\approx 1.05+.05\\times1.05\\\\\n",
    "  &\\approx 1.1025\n",
    "\\end{align}$$\n",
    "\n",
    "Now the actual value of $e^{0.1}$ is around $1.1051$ so we're a couple of percent off even over a very short time interval and only a couple of steps of the algorithm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## <span style=\"color:blue\">Exercise 2.5: Implementing Forward Euler's method</span>\n",
    "\n",
    "Write a function *euler*( *f*, *u0*, *t0*, *t_max*, *h*) that takes as arguments the function $f(u,t)$ on the RHS of our ODE,\n",
    "an initial value for $u$, the start and end time of the integration, and the time step.\n",
    "\n",
    "Use it to integrate the following ODE problems up to time $t=10$\n",
    "\n",
    "$$u'(t)=u(t),\\quad u(0)=1$$\n",
    "\n",
    "and \n",
    "\n",
    "$$u'(t)=\\cos(t),\\quad u(0)=0$$\n",
    "\n",
    "and plot the results. A template to get you started is below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (<ipython-input-5-fba0b6c9dcea>, line 8)",
     "output_type": "error",
     "traceback": [
      "\u001b[1;36m  File \u001b[1;32m\"<ipython-input-5-fba0b6c9dcea>\"\u001b[1;36m, line \u001b[1;32m8\u001b[0m\n\u001b[1;33m    while ... add your code here\u001b[0m\n\u001b[1;37m                ^\u001b[0m\n\u001b[1;31mSyntaxError\u001b[0m\u001b[1;31m:\u001b[0m invalid syntax\n"
     ]
    }
   ],
   "source": [
    "\n",
    "def euler(f,u0,t0,t_max,dt):\n",
    "    u=u0; t=t0\n",
    "    # these lists will store all solution values \n",
    "    # and associated time levels for later plotting\n",
    "    u_all=[u0]; t_all=[t0]\n",
    "    \n",
    "    \n",
    "    while ... add your code here\n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    return(u_all,t_all)\n",
    "\n",
    "\n",
    "def f(u,t):\n",
    "    val = u\n",
    "    return val\n",
    "\n",
    "(u_all,t_all) = euler(f,1.0,0.0,10.0,0.1)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Heun's method\n",
    "\n",
    "Euler's method is first-order accurate because it calculates the derivative using only the information available at the beginning of the time step. As we observed previously, higher-order convergence can be obtained if we also employ information from other points in the interval. Heun's method may be derived by attempting to use derivative information at both the start and the end of the interval:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "  u(t+\\Delta t)&\\approx u(t)+\\frac{\\Delta t}{2}\\left(u'(t)+u'(t+\\Delta t)\\right)\\\\\n",
    "  &\\approx u(t)+\\frac{\\Delta t}{2}\\big(f(u(t),t)+f(u(t+\\Delta t),t+\\Delta t)\\big)\n",
    "\\end{align}$$\n",
    "\n",
    "The difficulty with this approach is that we now require $u(t+\\Delta t)$ in order to calculate the final term in the equation, and that's what we set out to calculate so we don't know it yet! So at this point we have an example of an implicit algorithm and at this stage the above ODE solver would be referred to as the trapezoidal method if we could solve it exactly for $u(t+\\Delta t)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Perhaps the simplest solution to this dilemma, the one adopted in Heun's method, is to use a first guess at $x(t+\\Delta t)$ calculated using Euler's method:\n",
    "\n",
    "$$ \\tilde{u}(t+\\Delta t)=u(t)+\\Delta tf(u(t),t) $$\n",
    "\n",
    "This first guess is then used to solve for $u(t+\\Delta t)$ using:\n",
    "\n",
    "$$ u(t+\\Delta t)\\approx u(t)+\\frac{\\Delta t}{2}\\big(f(u(t),t)+f(\\tilde{u}(t+\\Delta t),t+\\Delta t)\\big)$$\n",
    "\n",
    "The generic term for schemes of this type is **predictor-corrector**. The initial calculation of $\\tilde{u}(t+\\Delta t)$ is used to predict the new value of $u$ and then this is used in a more accurate calculation to produce a more correct value. \n",
    "\n",
    "Note that Heun's method is $O(\\Delta t^2)$, i.e. a second-order method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## <span style=\"color:blue\">Exercise 2.6: Implementing Heun's method</span>\n",
    "\n",
    "Repeat the previous exercise for this method.\n",
    "\n",
    "For some ODEs you know the exact solution to compare the errors between Euler's and Heun's method, and how they vary with time step."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "You should be able to get a plot that looks like this for the case $u'=u$.\n",
    "\n",
    "![\"Comparison between the Euler and Heun method for the solution of a simple ODE.\"](images/euler_vs_heun.png)"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": false,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "974.226px",
    "left": "29px",
    "top": "110.73px",
    "width": "427.513px"
   },
   "toc_section_display": false,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}