{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Demonstrates accuracy of one- and two-sided finite-difference derivatives\n", "\n", "\n", "**Randall Romero Aguilar, PhD**\n", "\n", "This demo is based on the original Matlab demo accompanying the Computational Economics and Finance 2001 textbook by Mario Miranda and Paul Fackler.\n", "\n", "Original (Matlab) CompEcon file: **demdif02.m**\n", "\n", "Running this file requires the Python version of CompEcon. This can be installed with pip by running\n", "\n", " !pip install compecon --upgrade\n", "\n", "Last updated: 2022-Oct-22\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## About\n", "\n", "Demonstrates accuracy of one- and two-sided finite-difference derivatives of $e^x$ at $x=1$ as a function of step size $h$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Initial tasks" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "plt.style.use('seaborn-dark')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Setting parameters" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "n, x = 18, 1.0\n", "c = np.linspace(-15,0,n)\n", "h = 10 ** c" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "exp = np.exp\n", "eps = np.finfo(float).eps\n", "\n", "def deriv_error(l, u):\n", " dd = (exp(u) - exp(l)) / (u-l)\n", " return np.log10(np.abs(dd - exp(x)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## One-sided finite difference derivative" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "d1 = deriv_error(x, x+h)\n", "e1 = np.log10(eps**(1/2))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Two-sided finite difference derivative" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "d2 = deriv_error(x-h, x+h)\n", "e2 = np.log10(eps**(1/3))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plot finite difference derivatives" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEUCAYAAAAvLpGtAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABTb0lEQVR4nO3dd1xT1//H8VcSRoAwBBeioKIg7r33qHtLtbVW6+r42vHT2l2t2w47bK21rbbWWrFqa9Xa1m3rwr2wOHChCMqGQPb9/cFXvqIoBBJugPN8PHxokjvegOSTc8+55ygkSZIQBEEQhP9Syh1AEARBcCyiMAiCIAh5iMIgCIIg5CEKgyAIgpCHKAyCIAhCHqIwCIIgCHmIwiA8VGhoKAMHDmTw4MF5/ty4ccNu5zxz5gwvvfSSVft8/vnntGvXjjt37uR5fsCAAURGRtoyXqGsWbOGr7/+ukj7RkZGMmDAgAeev3HjBmFhYbk/g4EDBzJq1Ci2bt1apPNMmjSJS5cuFWnf06dPM2PGDKBoPy/B8TnJHUBwbCtXrsTX17fEzteoUSMWL15s9X6ZmZm8/vrrLF++HIVCYYdkhffEE0/Y5bhqtZrffvst9/HNmzcZN24cKpWK3r17W3Wsb775psg5Ll26REJCAlD0n5fg2ERhEIokMjKSefPm4e7ujlar5bXXXuPDDz/MfbxhwwZ+/fVXVq1ahVKppGLFirz77rvUqlWLN954g9TUVGJjY+natSvTp0/Pc9w5c+awZcsW3njjDTQaDefPnyc+Pp7Q0FDef/99PDw8HsgzaNAgTp06xYoVK5gwYcIDr4eGhnLw4MHcInf38cWLF/n444/x9/fnypUruLm5MXnyZFatWsWVK1d47LHHeOuttwDYtWsXS5cuxWg0olaref3112nWrBmff/45J0+e5Pbt24SGhhIUFERKSgozZszgypUrzJgxg+TkZJRKJc8//zz9+vVj9+7dLFu2DIPBQHJyMkOGDOGVV16x6mcQEBDASy+9xPLly+nduzcGg4GPPvqII0eOYDabqV+/Pu+88w4ajYbu3bvTuHFjzp8/z9SpU1mwYAGfffYZ33//PQ0aNGD8+PEA/PTTTxw+fJiPP/6Y+fPnc+rUKbRaLZIkMXfuXKpVq8bixYvJyMjgzTffZMiQIcyZM4c1a9bQpUsX/vrrLypVqgRAeHg4U6ZMoV27dg/N9dNPPxEREYGzszOurq7Mnj2bOnXqWPV9EOxAEoSHCAkJkQYMGCANGjQo988LL7wgSZIkHTp0SKpXr55048aNfB8fOHBA6tmzp5SUlCRJkiRt2LBB6tu3r2SxWKTXX39dGjt2bL7nPHTokNS/f39JkiTp9ddfl0aOHCnp9XrJYDBIQ4YMkdavX//APosXL5ZmzZolRUdHS82bN5fOnj0rSZIk9e/fXzp06FDu13I3y72PDx06JIWFhUlRUVGSJEnShAkTcs+ZlJQkNWjQQIqPj5euXLkiDRgwQEpOTpYkSZIuXLggdejQQdJqtdLixYul3r17S0ajMU8eSZKkIUOGSD/++KMkSZIUFxcn9ejRQ0pPT5eeeuop6cqVK5IkSVJ8fLwUFhaWm+fu13+v2NhYqWnTpg88f+HCBalJkyaSJEnS559/Li1cuFCyWCySJEnSokWLpJkzZ0qSJEndunWTvvjii9z9unXrJp0+fVo6ePCgNGDAgNznR4wYIe3fv186fvy49OKLL0pms1mSJElatmyZ9Oyzz0qSlPOznDx58gM/r9dee0369ttvJUmSpEuXLkldu3aVzGbzQ3OZTCapQYMGUkJCgiRJkvTrr79KERERD3yNQskTLQbhkR51Kcnf35+AgIB8H//zzz/069cvd99hw4Yxb9683P6JFi1aFOr8nTp1wsXFBYCQkBDS0tIeum1oaCivvPIK06ZN45dffinU8QGqV69O/fr1AQgMDMTT0xMXFxd8fX3x8PAgLS2NI0eOcPv2bcaNG5e7n0Kh4Pr16wA0bdoUJ6e8v06pqalER0cTHh4O5Hx/duzYAcBXX33Fnj172LJlCzExMUiSRHZ2dqEz35tBrVYDsGfPHjIyMjhw4AAARqMRPz+/3G1btmz5wP5t2rRBr9dz5swZ3NzcSE5Opl27digUCry9vYmIiCA2NpbIyMh8W2r3Cg8PZ9asWUyYMIENGzYwfPhwlErlQ3OpVCr69OnDqFGj6Nq1Kx07dqRLly5Wfw8E2xOFQSgyd3f3hz62WCwPbC9JEiaTKd99H+bumx7kvAlKBUztNWbMGPbt28e8efMeuo3BYMjz+G7huev+N3jI+XratWvHp59+mvvcrVu3qFy5Mtu3b8/367l7nHv7PC5fvkzVqlUZOnQoPXv2pGXLlgwfPpwdO3YU+LXl58yZM4SEhORmfOutt3LfXLVaLXq9Pnfb/DIqFApGjBjBb7/9hrOzMyNGjEChULBnzx7mzZvHM888Q48ePahduzabNm16ZJaWLVtiMpk4ffo0W7ZsYe3atQXm+uijj7hw4QIHDhzg66+/5rfffuOzzz6z+vsg2JYYlSTYRadOndi6dSvJyckAbNiwAR8fH4KCgux+7gULFrB3716uXbuW+5yvry9nzpwBYMuWLVYfs127duzfv5+YmBgA9u7dy6BBg9DpdA/dR6PR0KBBAzZu3AjkFJInnniCCxcukJmZySuvvEL37t2JjIzEYDDkW0wf5cqVK3z55Ze5/QMdO3Zk9erVucd69913+fjjjws8ztChQ9m1axd//fUXw4YNA2D//v1069aNJ598koYNG7Jjxw7MZjMAKpUqt8DfLzw8nDlz5hAaGoq/v/8jcyUnJ9OlSxd8fHwYN24cr7zySu7PSJCXaDEIjzR27FiUyryfH6ZOnZrnk3x+OnTowLhx4xg7diwWiwVfX1+WLVv2wLHswdfXl4ULFzJx4sTc59555x1mz56Nl5cX7du3z+0gLaw6deowe/Zspk6diiRJODk5sXTp0gIvryxatIhZs2axatUqFAoF8+bNo3HjxnTt2pW+ffvi4uJCSEgIderU4dq1aw+0Xu6l0+kYPHgwAEqlEldXV6ZOnUrXrl0BeOGFF3j//fcZOnQoZrOZsLAw3njjjQK/tkqVKlG/fn1MJhNVqlQBYNSoUUybNo2BAwdiMpno0KED27Ztw2Kx0LRpU5YsWcKUKVMYM2ZMnmMNGTKEjz/+OE9BelgujUbD888/z7hx41Cr1ahUKubOnVtgXsH+FFJR2q+CIAhCmSUuJQmCIAh5iMIgCIIg5OGwfQxDhgzB09MTyBlOuGDBApkTCYIglA8OWRjuDmVbtWqVzEkEQRDKH4e8lBQdHU12djbjx4/n6aef5uTJk3JHEgRBKDccclTS+fPnOXXqFOHh4Vy9epVJkybx559/5rnx6M6dDBkTCuVZfHrOvQtVvR49ZLe0U2bEAWDxrCZzEsGWKlXyLHAbh7yUVKtWLYKCglAoFNSqVQsfHx/u3LmTe8OMIMhp5h/nAVg2sonMSezLc0fOdNppQ9fLnEQoaQ5ZGNavX8+FCxd47733SEhIIDMz0+obkgTBXsa3DZQ7QonIavmy3BEEmTjkpSSDwcCbb75JXFwcCoWCV199lebNm+fZRlxKEgRBsF5hLiU5ZGEoDFEYBLncSM2ZBbW6j5vMSexLmZYz15TF2/7zWwklp9T2MQiCI5vz1wWgHPQx7JoGiD6G8kgUBkGw0uT25eMTdFbraXJHEGQiLiUJgiCUI4W5lOSQN7gJgiO7mpzF1eQsuWPYnSolBlVKjNwxBBmIS0mCYKUF2y8CZb+PQbPndUD0MZRHojAIgpVe6FhT7gglQtu24EV+hLJJ9DEIgiCUI6KPQRDs4FKilkuJWrlj2J0qKRpVUrTcMQQZiEtJgmClD3deAspBH8Pf7wCij6E8EoVBEKz0UpfackcoEdr278gdQZCJ6GMQBEEoR0QfgyDYwfnbmZy/nSl3DLtT3YlCdSdK7hiCDERhEAQrfbw7ho93l/0bvzT7ZqLZN1PuGIIMxKUkQbDS3dZCaGWNzEns625rwVypgcxJBFsS024LgiAIeYg+BkGwg6j4DKLiy/4HE6eEkzglnJQ7hiADURgEwUqL915m8d7LcsewO48Dc/E4MFfuGIIMxKUkQbDS3bue61T0kDmJfd2969nsV0/mJIItiT4GQRAEIQ+xtKcg2MGpm2kANAnwLvQ+726N5s9/bxdq2yPTOhcpl6053ToKgMm/Zf6vxx/HOe4QFreK6MMex3P7i+jqP4nPxvBCHf/Of27YLKtgW6IwCIKVvtx3FSj8XEmJmXpq+ro5zBt+YXkcWgg8fK4kU9XmuF75E4tbRZTaBMwV6mAMaCfe8MsAcSlJEKx0d/W2mr7uhdr+p2M36BzsR3UfNwAOXEnmi3+uAOCiUrLiyaYoFQr7hC2Gu6u3mSsEP3I79ekVKCxm9LUew+L9v/Wwna/tRnNwAQCSyoXUEZtAIca7yE1cShIEOyioIKRmG/k3IYN2NX0BuJWuzy0KAB/tusTXo5pS0cPFrjmL61EFQaFLwTnuME6JZ7G4V0GZdjFPUQDQ/PMuqUN/QfKobO+ogo2J8i0IVjoWm8qx2NSHvn4rXcdney+ToTNxOUlL3ftGL7Wv5csTK4+xyMGn1XC+eRDnmwfzfU2VHotH5AdkN5mIsVprzH5hD2xjCOqOb0RPPP4R02qUNqLFIAhW+vrANQCWjfTJ9/WwKp70CKnEjgt3SNIaGNU8IPe1ux3XfzzXFiel410+upf74UVA/n0MpsqN0dcZgOulLSiz7pDdZEKe151uHQVJIumZ46AUbzOljfiJCYKV3u0dUuA2fcMqM/vP8zQJ8Ebj+r9fs50XEgms4IaTUoEkSWgN5jyvO5KM7ose+bouZBieO6di8m+F5JL3urVrzBbMPrVzioIkoTBmPrCN4Lgc83+kIDiwe/sLHrWNWYIWNfIOae1drxJztl3g19PxuDopeb1nHcKqOOYb5v19Bvm9rpDMGALaPfCavu5gPHe9itu51UgqNZld5mOq3NheUQUbE6OSBMFKkddSAGgTVOGR2528kUbjAC+HHHFUGM6x/wBgrNHpods4xR3Ouc9BjDYqNcSoJEGwgxWHrgMFF4am1Qt/A5wjcj/6GQBpjygMpmqtSyqOUIJEi0EQrBSfrgOgqpda5iT2pcyIA8DiWU3mJIItiRaDINhBWS8Id4mCUH6JC4OCYKUDV5I5cCVZ7hh253xtN87XdssdQ5CBaDEIgpVWHo4Fcm5UK8vcjy8BIC2om8xJhJIm+hgEwUqJWgOAw09pUVwKbc5ssGJKi7JF9DEIgh2U9YJwlygI5ZfoYxAEK/0dk8TfMUlyx7A7lyvbcbmyXe4YggxEi0EQrLT6aM56A52D/WROYl9uJ5cBYKjVS+YkQkkTfQyCYKXULCMAPu7OMiexL0V2zsgrya1sd7KXN6KPQRDsoKwXhLtEQSi/RB+DIFhp18VEdl1MlDuG3bnEbMUlZqvcMQQZiBaDIFhp7fGbAHSvW1HmJPbldnoFAIbgfjInEUqa6GMQBCtl6k0ADruOgq0o9OkASK5eMicRbEn0MQiCHZT1gnCXKAjll+hjEAQrbYu+zbbo23LHsDvXi5twvbhJ7hiCDMrHRx9BsKENp24B8Fi9sn1nsPrsDwDo6w6SOYlQ0kQfgyBYSWc0A6B2VsmcxM6M2Tl/Oxe8lKng+C7eyWTj6Xg+GNWswG0dssVgsVh47733OH/+PC4uLsydO5egoEevPysIJaXMF4S7REEoExK1Br7af5VNZ+LxUjvxQSH2saqP4dChQ0WMZp0dO3ZgMBhYu3Yt06ZNY+HChSVyXkEojK3nEth6LkHuGHbnen4Druc3yB1DKCKd0cyKQ9cZtvwwv0cl8ESLAH6Z0KpQ+1rVYvj8889p27ZtkUJa49ixY3TqlLPObNOmTTl79qzdzykIhfXbmXgA+tWvInMS+1KfWwOAPnS4zEkEa1gkib+ib7Pkn6skZOjpWsePFzvXJrBC4VuAVhUGhULBf/7zH2rVqoVSmdPYmDp1qnWpCyEzMxONRpP7WKVSYTKZcHJyyCtfQjmzZEQjuSOUiLRBa+SOIFjp1M00Ptlzmaj4DOpV1jCrbygtavhYfRyr3mmHDx9OamoqZrMZhUKBn599ZpfUaDRotdrcxxaLRRQFwWE4qcrJKG9V+ZgTqiy4kZrNkn+usONCIpU0LszsE0K/+lVQKhRFOp5V77aurq5EREQQHBzMhQsXmDJlSpFOWpDmzZuze/du+vXrx8mTJwkJCbHLeQShKDafzbmUNLBhVZmT2Jfrvz8DoA97XOYkwsNk6k2sOHSdiBM3USkUTG4XxFOtquNWzAESVhWGlStX8ssvv+Dh4UFmZiZjx45lyJAhxQqQn169erF//35GjRqFJEnMnz/f5ucQhKLaEpXT8VzWC4M6WhQGR2WySGw8fYtlB66Rmm2kf4MqvNChJpU9XR+6j3PsPjx3vgzTLxR4fKv7GDw8PICcyz2urg8PURxKpZLZs2fb5diCUFzLRjaRO0KJSBu6Xu4IQj72X0nms72XuZKURfPq3rzStTZhVR49/5Hz9b14bx2P2adWoc5hVWEIDAxk4cKFtGzZkqNHjxIYGGjN7oIgCEIRXUrU8tneyxy6mkINHzUfDqpPlzp+KAroR3C5tguvPyZh9gkmdXAEhZkT2Ko7n/ft28e1a9eIiYkhODiYxx9/HGdneTqoxJ3Pglx+PZ0zJcbQxv4yJ7EvddRqAHQNRsucpHxL0hr4+sA1Np65hYeLExPbBRLetBrOhRgE4XJ1B15/TMbkG0La4DVI6gq2n1116dKlrF692ppdBKHM2X7+DlD2C4Prpc2AKAxy0RpM/HwijpWHY9GZLIQ3rcbEdkH4uBXuw7jLlW14/fkspor1SRu4GkntU+hzW9VieOqpp/D29rb7fQyFIVoMgiCURZn6nILw07EbpOlMdKrty0tdalPT173Qx3C5/Adefz2PqWJD0gatRnL1zn3N5i2G4cPFHZCCIAj2kKk3EXH8JmuO3yRdZ6JjbV8mtg2kgb9162K4XNqC1/YpmCo1Jm3gj0VaV8OqwuDv718iU2IIgiNbdzIOgPCm1WROYl/qMysB0DUaK3OSsi1D97+CkKHPaSFMbBdE/aoFf7K/n+vFTXhufxFTlWakDVyF5GL9McBB50oqjFl/nifAW011Hzeq+6gJ8Fbj4+ZcYA+9IBTXPzFJQNkvDK5XtwOiMNhLus6YWxAy9Wa6BPsxsV0g9QoYevowrhd+xXPHyxirtiJ9wEokF03BOz1Eqe1jaD13O7czDXme83BREeCtJsDHjere6pyC4eNGgLeaql5qnJSiaAiCIK+0bCNrjt8k4vhNtAYzXev4MbFdEKGVi/5G7np+PZ47p2Ks1oa0/ivB+eH9EXbvY9Dr9dbsblO/P9sWndFMXLqOG6k6bqbpuJmazY1UHZcTtey7nITR/L+ap1Iq8PdyzW1l3P+3u0s5mWNfEARZpGYbWXPsBmtPxKE1mOletyIT2gYSUoyCADlTl3jumoYxoD1p/b+3yToahWoxvPLKK3z66acArFixgvHjxwPw9NNP88MPPxQ7RFEUNCrJIkncztBzM03HjdTs//79v3+n60y52zopFbzVq26Zn+JAsI01x28C8ETzAJmT2JfbqW8ByG4yUeYkpVtqlpHVx27w84k4so1meoRUZELbIOpU8ij2sdXn1qDZ/RrGGp1I67u8UEXBZi2GpKSk3H/v2bMntzA48qqgSoWCql45l5Dym3Y2XWfMLRa/nL7FnL8uIAGDRHEQCnDkWgpQ9guD8439gCgMRZWSZeDHozdZd/ImOqOFnqGVmNA2kOCKxS8IAOqoH/Hc8waGwK6k9f0GnGy34p7Vc1nfWwxKbUev2YiPOZkKqiQauyfRq7meV/Bj7l85k0uJ4iA8ysdDG8odoUSk9/9O7gilUnKWgR+P3GD9qTh0RguP1avE+LaB1PazTUGAnBFjnn+/jT6oO+l9vgYntc2ODYUsDPcWAIcsBmYjSl0SiuxklNlJ//2TeM/jRJS6ZBR3X9OnPXCIb6t35oUaL+UUBwkGNRLFQRCEwkvU5hSEDafiMJgtPFavMhPaBFLTr/A3phWG+vQKPP+Zgb5mL9L7fAUq209mWqjCcOnSJaZNm4YkSXn+HRMTY/NAheXzy1AUWTlv+Pm90QNICiWS2heLmx8WN1/MFRsgufliUfv99zk/JDc/VMnn0ex7j2XeCfxfjbeYu+2/LQdRHIR8rDoSC8CYVjVkTmJfbie+AiC72XMyJ3FsCRl6fjgcy8YztzBbJHqHVeaZNoFW3alcWG4nv0Gzfxb6Wr1J770UVC42PwcUsjDc7XgGGDVqVL7/LmmS0gVzpUY5b/RufljcKmJx80Vy88t543evmHMbuKLgiaaMAe0w+wTj9edkFju9ilfAu8zdlvOaKA7C/c7cKh/TsTjHHwMgW+YcjiouTcfKw7FsjorHIsGA+lUY16YG1X1sd63/Xm4nvkJzYC764H6k91pi1xX2rLqPwZHYY64kVdJ5vH8fiyI7iY810/kivh7vPFaXwY3K9mRpgiAU3vWUbL6PvM7Wf2+jVOT0SY5tXQN/L9te57+X2/ElaA4uQFdnIBk9FxerKNj8PoayzuwXSsrwTXhvfYZpt+dQqdJkZm6TkCQYUsZn0hQE4dEuJ2lZceg628/fwVmlJLxpNca0rP7IVdNswf3oYjwiP0BXdzAZPT8Dpf3ftkVhuI/kUZnUIevx2vEiYy8vo6pvHC9sz1naUBQHAeD7yOsAjGtTtheqcjv2BQDZLeyztntpceF2Jisir7PrQiJqZyWjW1RndMvq+HnY5/r+vdyPfILH4UXoQoaR0ePjEikKIApD/pzdSO+9DI+D8+l9chnrvRMYvf05JMr+HPxCwS7c0codoUQ4JZ6TO4KszsVnsPzQdf6OScLDRcUzbWrwRPPq+LiXwOJkkoR75Id4HFuMrl44Gd0+AmXJzc5gVR/Dxo0bWbZsGQaDAUmSUCgU7Ny50575Hqqk1mNQn/0Bzd/vcE1Vi5GZ/8f4Xm0YJoqDIJRZp26msfzQdQ5eTcFL7cSo5gGMahaAp7qEPkdLEh4H5uJ+chnZ9Z8gs+v7hRpEU1iF6WOwqjD079+fL7/8En///70xurjYvzmVn5JcqMfl2i48/3qeVLMbo7OnMaRnT1EcBKEMkSSJ4zfS+PbQdY5eT8XHzZnRLQIY0bQaGtcSvLAiSXjsm4n76RVkNxxLZuc5Ni0KYIfO5xo1ahAUFFTkQKWVIag7qcN+xXvL0/wizeb5nckghTOsSdmedlnI37cHrwEwsV3Z/l1wP/IpAFmtXpE1hz1JkkTktRSWH7rOyZvp+Hm48EqX2gxr4o+bcwlPrClZ0Ox9G7eoVWQ1mYS2wwyQ6YZiqwqDWq1m4sSJhIWF5d4BLde02yXNXLE+qSM24/n7M6xIXMSM3Yls4HmGi+JQ7lxLKR8j+1Wp8t3Aam+SJLHvcjIrIq9z9lYGlTUuTO8ezKCGVVGXdEEAsJjR7HkNt3/XktX8P2jbviFbUQArLyX9+uuvpKamYjabUSgU+Pn5MWTIEDvGezjZ1nw2aPH863nU13fxjakfli7vMqxp2b4DVhDKCosksediIssPXefCHS3VvFwZ2yaQAfWr4OJk20s2hQ9lwnPnVNQXfkHb8hWyWk+za1Gw+aUkV1dXIiIiCA4O5sKFC0yZUg6Hsbl4kNH/O0x/z2RS1Pf89fdtNpoXMaRFsNzJBEF4CLNFYsf5O6yIvM7lpCwCK7gxs08IfepVxkklU0EAMBvx3PEy6kub0LZ5jayWL8mX5R5WtRhGjhzJihUr8PDwIDMzk7Fjx7JhwwZ75nso2VoM93A+8S1eB2ZxxlKL422+YEDrxnJHEkrAV/uvAvBch5qy5rA398gPAchqM13mJEVnMlv449/bfH84lusp2dT2c2d8m0B6hlZCJfeKjmYDXttewPXyn2S2f6fE5qSyeYtBoVDg4ZEzdaxGo8HV1b53/Dk6Y7OJpHhVp95f/8Hv8NPsNH5Ojw6d5I4l2FlChnwrF5YkVeYtuSMUmcFkYcu5BFZGXicuXU9IJQ/eH1SfrnX8UDrCDNEmHV5/PYfr1R1kdpxFdpMJcifKw6oWw2uvvYavry8tW7bk6NGjpKamsnDhQnvmeyhHaDHkij+J869PozRns7vBB7TvNlTuRIJQLumMZn47E88PR2K5nWmgQVVPJrQNpGNtX8dZMsCUjfcfE3G5vpeMLgvQNRxToqe3+X0MJpOJtWvXEhMTQ3BwMI8//jjOziVwF2A+HKowAJbUWAxrn8TfeI1tNf4PZbNx+Hu7UdnTFSe5m6yCUMZlGcz8cvoWPx69QZLWQLMALya0DaJ1kI/jFAQAYxbevz+D880DZHb7EF39kp+h2uaFwZE4WmEAMGWnkbj6aRrpj3HYEsos41j+pSaVNK74e7lS1Uud529/TzVVvVzlGR4nFNkX/1wBYEqnWjInsS+PgwsA0LZ7U+YkD5epN7HuZByrj94gTWeiVaAPE9oG5rucr9wUhky8tozFOf4IGT0+QR86XJYcNutjWLduHeHh4SxatOiB6lte7mMoDCc3b6o8s4Grx1fT5OQitijf5rjfQNZqnuaiVsGpm2lsj9Zjvq8U+7o7/69oeN5XPLzUJXcrvlAoadlGuSOUCIUuRe4ID5WuMxJx/CYRx+PI0JvoUMuX8W0DaVzNS+5o+VLo0/HeMganhJNk9PoCfd1Bckd6pEK1GPbt20fHjh1ZsmQJAQE5C6BLksSiRYvYt2+f3UPmxxFbDPdS6NNwP/IJbqe/Q3LRkNVqKtkNn8akcCIxU8+tdD230nXEp+uJS9cRn67jVrqehAw9epMlz7Fq+bozoEEV+tWvTEVN+e7wF8q3lCwDPx27ybqTcWgNZrrW8WN820DCqhT8KVguCl0q3ptH45R4jvTHlmAI7idrHptdSlq3bh3r168nJiaGOnXqIEkSFosFo9HIxo0bbZHVao5eGO5SJV9As+89XGL/xlQhhMxOszDWePjIJUmSSM4y5haKm2k6/o5J4nRcOkoFtKvpy8CGVehU20++G3IEoYQlZupZdfQGv5y6hd5koWdoJca3CaROJQ+5oz2SIjsZ701P4JR8kfS+X2Oo2VPuSLYrDAaDgdu3b7Ns2TKeey5nrK1SqcTPz69cTKJXbJKEy5VtaPbPRpV+DX3tPmR2mIHFq/Dz+V9NzuL3qAS2nkvgdqYBL7UTvetVZmDDKtSrrHGsDrYy7tM9lwF4pWttmZPYl8f+OQBoO7wrW4b4dB0/HLnBb/eup9w6kJp+tl9P2dYUWXfw+W0UqrSrpPVbjjGwq9yRgBLofDYajWJUkjVMOtxPfoP7scUgWchq9hxZzf8DzoX/T262SBy+nsKWswnsuZSIwSwRXNGdAQ2q0jescoksHlLevb/jIgCv96wrcxL70ux9G4DMLvNK/Nw3UrP5/nAsv0clANC/QRXGtbbfesq2ptQm4P3bSFQZN0nr/z3G6h3kjpTL5oUhIiKC7777DpPJhCRJODk5sW3btmKFLKpSWRj+S5l5C4+D81Ff+BWzxh9tu7fR1x1s9fwo6Toj28/fYUtUAmdvZaBSQPtavgxoWJVOtX1xlvNWf0EogqvJWXwfeZ0//72NSqnIXU+5qh3XU7Y1ZWYc3htHosy6TfqAlRirtZU7Uh42LwzDhg3jq6++YunSpfTp04eVK1fy5ZdfFitkUZXmwnCX060jaP6ZgfOdMxj9W5PZaTamSg2LdKzLSVp+j0rg93O3SdIa8FY70SesMgMbViW0ssbGyQXBti4l5qynvOP8HVyclAxv4s9TLatTqZQNtlCmx+Lz20gUuhTSBv6IqWoLuSM9wOZTYlSoUIHKlSuj1Wpp06YNixcvLnI4AUz+rUgdsQV19Fo8Di7E5+e+6BqMRtvmNSQ3X6uOVdvPgxc71+b5jrWIvJrClqh4fjl9i7Un4qhbyYMBDarQN6wyFdzFpabiWrQ7Zzrqad3K9sSJHv/MBEDbaZbdzhGdkLN85p5LSbg7q3i6dQ2ebBGAbyn8f6pMu4rPxpEojJmkDVqDqUpTuSMVmVWFwdPTkx07dqBQKIiIiCA5OdleucoPpQpd/SfRB/fPHd7qemlz7vBWVNb14TgpFXSo7UuH2r6kZhvZFn2HLVHxfLLnMov/vkKn2r4826EmdSo69mgOoWw7E5fO8kPX2X8lGY2rioltAxnVPABvN3n6LItLlRKD92+PozAbSBu8tsgtf0dh1aWkzMxMYmNj8fPz47vvvqNbt260bt3anvkeqixcSsqPtcNbC+tSopYtZxPYEhWP1mBmYrtAxraqIe+Uw0K5cyw2leWHrnPkeireaidGt6xOeEkvn2ljquQL+GwcCUikDl6D2S9M7kiPZPM+hl27drFhwwYMBgOSJKFQKPjmm2+KFbKoymphAB4Y3qoLHkBm5zlI7pWKfeiULAMf7oph+/k7hFbWMLNPCHUriT4IwX4kSeLwtVSWH7rGiZvp+Lo781TL6gxvUg13l9I9HYwq8Rw+m55AUjiRNjgCs6/jj1SzeWHo3bs3s2fPxtvbO/e5evXqFS1dMZXpwnCXSYf7yWW4H12M5KQms+OsnPlVbHDPwq6Liby/4yJpOhMT2gQyrk0NMYqpkMRw1cLJb/nMp1vVYHAjmZbPtDGnO2fw/u0JJGc30gavxexTOu5rsXnnc926dWnTpk2RAwlWclKT1fJl9MH98dz1Kl47X8FwcSMZXd/H4hlQrEN3r1uR5tW9+WjXJb4+eI3dlxKZ2TuU0Cqi9VAQV6fS/6ZWGJJT0YaI5rd85ps96zCgQdUyc7e+U/xxvDc/heTqReqQn626WbU0sHrN54iICGrX/l9lXLBggV2CFaRctBjuJVlQn1mJ5uACJIUCbbs30TV8GhTF/0XbeymRBTsukZptZFzrGkxoGyhaD4LV8ls+c1zrGvQNk3n5TBtzunUE781jkNz8SB28FotXdbkjWcUu9zFMnDgRT8//HbhTJ3lWLCt3heG/lOmxeO55A5fYvRj9W5PR7UPMFYo/bDIt28gne2L4/dxtgiu6M6N3KPWrOu7EZILjuH/5zFp+7kxwlOUzbcz55gG8t4zDrKlK2uC1WDT+ckeyms0Lw+TJk/n666+LFcpWymthAECScD2/Hs2+91CYdGhb/R/ZTZ+1emhrfvZdTmL+9oskaw2MaVWDSe2Cykzz31bmbbsAwNuPhcicxL40u18DILPbB/m+nrt85uFY4tJ0hFTyYELbQLrWregYy2famHPsP3hvfQazZyCpgyOQPCrLHalIbN7HoFarmTBhAvXr18+dtM3W6zFIkkTnzp2pWbMmAE2bNmXatGk2PUepp1CgrxeOoUYXPP95B82hhbhe2kJm94+KPX66Y20/1o715pM9MXx/OJa9MUnM7B1CA3/HnOdeDqV1rL21JHWFfJ/Pb/nMV7sFO9bymTbmcnUnXn9OxuxTO6couPnJHcmurO5jyLOzQsGQIUNsGujatWssWLCAr7766pHblesWw31cYn7Hc+87KHTJZDV/gayWL0MROw7vdeBKMvO2XSBRa+CpltWZ3L4mrqL1UG7dv3xm0wAvJrQNpE1QhTJbEABcLv+F11/PYfKrR9qgnx5aMEsLm6/gdvnyZbv/B4iKiiIhIYExY8agVqt5880383R2Cw8yBPcnOaA9mv2z8Tj2Oa4xW8no/hEm/1bFOm77Wr6sHdeST/de5ocjN/g7Jol3e4c67CpZgn3kt3zmvP71HHL5TFtzubQFr+1TMFVsSNqg1Uiu3gXvVAYUqsXwzz//0KlTJ5u3GNatW8fKlSvzPDdjxgySkpLo27cvR48eZcGCBWzYsOGBfUWLIX/O1/fguft1lJlxZDcah7btG+BS/OkvDl1NZt62iyRk6HmyRXWe6xBUJsaiF8WsP88DMLNPqMxJ7Mt12ytcScriqaRxpWL5TFtzvfArnjtewVS1OWkDfkByKRuDMWzWYrg78kin0/HEE08AOWsxfPzxx8WIB+Hh4YSHh+d5Ljs7G5Uq5w2nZcuWJCQk5N5lLRTMGNiVlCd24nFoIW5nvsf16nYyur6PMbBLsY7btqYva8a24It/rrD62A3+uZzEu4+F0LR6+fgEda8qnqVrxk9rpWQZWH3sJlUvSBjMrrSo6e3wy2fammv0Ojx3TsVYrQ1p/Vfa5MNVaWJVH8Mbb7yB0Whk4sSJzJw5k86dOzNlyhSbBvrwww/x8fFh0qRJREdHM2PGDH7++ecHthMthoI5xR3Gc/d0nFJj0NV7nMwO79rk+ujhaynM23aBW+l6hjXxZ0ijnKm9RfEu3Urr8pm2po5ajWbPGxirdySt3wpwLh2LAxWWXVZwe/PNN9m4cSOzZ89+4NO+LaSlpTF9+nSysrJQqVTMmDGD4OAHx+mLwlBIJh0eRz7F7cRSJLUvGZ3nYAjqbtWqcfnJMpj54p8r/HIqDrME1bxc6Va3Et1DKtLQ37NMDlcsq+5fPrNPWGXGlZLlM21NfeZ7PP9+B31gN9L7fmOTQRyOxuaFYdq0aej1eiZNmsT8+fPp3r07zz77bLFCFpUoDNZxunMWza5pOCdGASA5uWFxq4jFzfe/f1dEyv23HxY3PyT3iljUfljcfB/6C5KaZWRvTCK7LiZy+FoqJotEZY0LXetUpHtIRZoGeJe5m5ze3RoNwJx+8swTZisFLZ/puf1FADJ6fS5bxpLkdvIbNPtnoa/Vm/TeX4KqbF4ytPl9DM2bN2f06NEArFq1ikWLFhUtmVDiTJUakjpiC66Xt6LMuIEyKwlldiJKXRJKbQJOSedQZiWhsBjy3d/i4omk9sVyt1i4+2Fxq4jaM4AhtR9jcKNGZOhM/HM5id0XE/ntbDw/n4yjgpszXer40T2kIq1q+JSJqRGCKpTuSwv3L585tLE/T7eq/sDymWafsr0Q0b3cjn2B5tBC9MH9Se/1hU1uFi3NrL6UlJqayvXr16levTq+vtatMmZLosVgB5KEwpCBIjsJZe6fRJTZySiyE3P/rcxORJGdjFKXhMJiQkKBMaAd+rqD0Nfuh+TmS5bBzIEryey6mMj+y8lkGc14ujrROdiXbnUr0bZmBXFPRAm7dEfLisic5TNdnZQMK6XLZ9qUJOF+9FM8Di9CV3cIGT0/BWXpXRuiMGx+KWnr1q189tlnBAcHc/HiRaZMmcLgwYOLFbKoRGFwAJKEKuUirhc34XppE06pl5EUKow1OqGrMwhD7d5Irt7oTRYOXU1h98U7/B2TTIbehLuzig61feletyLta/mW+nn5Hdn9y2eGN6vG6BYBYplXScI98kM8ji1GVy+cjG4fgbLs/z+0eWEYOXIkK1aswMPDg8zMTMaOHZvvPQYlQRQGByNJqBLPob60CdeLm1BlxCIpXTAEds1pSdTsBS4eGM0WjsamsutCInsvJZGSbcTVSUm7mhXoVrcinYP9HH41r7e2/AvA/AGOvVLXvctnero6Map5NUY2K/zymZ5/vQBARu8v7RlTHpKEx4G5uJ9cRnb9J8nsutAmMxWXBjbvY1AoFHh45Axd02g0uLqW4yaokJdCgblSA7SVGqBt+wZOt0/mtiRcr25DclKjD+qJvu4g2gV1o13NEF7vKXHyRhq7Liay+2Ji7ifagQ2r8HizAAId9Fp+iIMP37x/+cwXOtYs0vKZpor17ZRQZpKEx76ZuJ9eQXajsWR2mlNuikJhWdVieO211/D19aVly5YcPXqU1NRUFi5caM98DyVaDKWEZMH51pGcIhHzO8rsRCzOHhhq9UZfdxCGGp1B5YJFkjgTl84vp2+xLfoOZotEh9q+jGoWQOsgH3GPRAEkSSLyWgorDl3PXT5zTKsaDGvsLy7T3UuyoNn7Fm5RP5LVZDLaDu/aZEXE0sTml5L27dvHtWvXiImJITg4mMcffxxnZ3l670VhKIUsJpxvHsxpRcRsRalPw+Lqjb52H/R1B2MMaA9KJxIz9Ww4dYtfTt8iOctIbT93RjYPoF9Y5XI7DcfDlPXlM23KYkaz+zXcoteS1fw/OdPFlLOiAHYoDKNHj2b16tXFCmUrojCUcmYDLrH/4HppEy6X/0JpzMTi5oe+zkCyw57AXKkBepOF7edvE3E8jvO3M/FWOzG4kT/hTf0fGFpZkl7bdA6ADwbJd6klv+Uzx7YJZED9KjZbP8Prj0kAOTd6lXYWE547/w/1hV/Rtvo/slpNLZdFAexQGJ566im8vb2pVasWSmXOfz5br8dQWKIwlCEmHS7Xd+dcbrqyDYVZj7FSY3T1R6GvOxiLixcnbqYRcTyOvZcSUQDd6lZiVPNqNK7mVeKXmVYdiQVgTKsaJXpeyH/5zGfa1KBPPdsvn+l2Imfq++xmz9n0uCXObMRzx0uoL21G2+Z1slq+KHciWdm8MNw/uyrA0KFDrUtlI6IwlE0KXQquF37F7VwETknncjqtg/ujCxuFsVpb4tL1rDsZx29n4snQmwiromFU8wB6hVYq0+tU3798Zm0/d8aX0eUzbcpswGvbC7he/pPM9u+S3UyemRociV3mSrJYLCQnJ+Pn5ydrh6AoDGWcJOF05wzqc2twvbgRpSEDk3ctdGEj0dcLJ9O5IlvPJRBx/CbXUrLx83BhRBN/hjXxx7cMjc/PXT4z8jpx6fqc5TPbBdG1jp+Yj6ogJh1efz6L67WdZHSaja7xeLkTOQSbF4bt27ezYMECvLy80Gq1vPfee3To0KFYIYtKFIZyxJiNa8zvqP9dg0tcJJJChSGoB7r6o9AFduPQ9Qwijt/k4NUUnFUKeterzKjmAYRW1tglztRfzwLw8dDiLaP6KPktnzmhbWCJLp/p9fszAKT3/65EzmdTxmy8/5iIS+xeMrouRNfgKbkTOQybF4YhQ4awfPly/Pz8SExM5LnnnmP9+vXFCllUojCUT6rUy6j/jUD97zqU2Xcwu1dBX28EurCRxJirsvbETbZEJaAzWWhW3ZuxrWvQvqZtl55cc/wmAE80D7DZMe+6f/nMZgFeTGgbJMuQXbdT3wKQ3WRiiZ632IxZeP8+DuebB8no/hH6sJFyJ3IoNi8M48aN4/vvv3/o45IkCkM5Zzbicm0X6n8jcLm2E4VkwVCtLbr6T5AY8Bgb/03l5xNxxGfoaeTvyaT2QbR14LWJ718+s3WgD+PbBpaL5TNtSWHIwHvLWJzij5LR4xP0ocPljuRwbF4YpkyZQnZ2Nq1ateLs2bMkJibSunVroORHJ4nCINyl1MbjGr0et3NrUKVfw+LihT5kCJmhI/nldmW+O3Sd+Aw9jat5MbmdPJ++HyZdZyTi+E0ijseVy+UzbUmhT8N78xicbp8io9cX6OsOlDuSQ7LLqCSz2Zzb+Xx3yCqU/OgkURiEB0gWnOMO5XRYx2xFYdajr9mLlA5z2HhNyYpD17mdaaBpgBeT2gXRKrBoBeKlDWcAWDy8UZGjJmbqiTgRx/qTcWgNZrrW8XO45TO9N+dcl08b+KPMSQqm0KXgvWk0Tkn/kt57KYbafeSO5LBsXhi2bdvGwoULczufZ86cSceOHYsVsqhEYRAeRaFLRR31Ix5HP0VSqMhqM520sLFsjLrD94evcyfTQLPq3jzbPsjqyzXrTsYBEN60mlX7JWoN7LqQyM4LdzhxIw3AoZfPVJ9ZCYCu0ViZkzyaIjsJn9+eQJUaQ3qfrzHU7CF3JIcmOp+Fck+Zfh3PvW/hcn0PxspNyOj6AVkVwth4+hbfH44lUWugRQ1vJrcPonl1H5uf/24x2HHhDidvpCEBtfzc6RlSkT5hVRx2osDSQqG9jc+mJ1ClXSWt3wqMgV3kjuTwbD67qo+PD35+fgBUrFgRjcY+wwEFwVYsXoGkDViF68Xf0OybSYV1/VA3ncTIVtMY3Kgqv56JZ+XhWJ5de5qWgT482y6IptW9i3XOxEw9uy4msuNCYp5iMKldED1CK1Lbz/FaB6WRMvMW3r+NQpUZR9qAHzBWl2fofFlU5M7nqKgo7ty5IzqfhVJDoUvB48A83P6NwOwVSEaX+RgDu6Iz5gwRXXk4luQsI60DfZjcPogmAfkXiBfWnQbgy/DGuc/lVwxq+7nTM6RSqS0G3r+NAiBtcITMSR6kzLiJz8bHUWQnkj7gB4zV2sgdqdQokSkx7hKdz0Jp4XzzIJo9r+OUehld3SFkdnwPyb0iOqOZ9adu8cPhWFKyjbQNqsDk9kE0um+E0K+nbwHQqbZvTjE4f4eTN9PLRDG4lzoqZ8JMXYPRMifJS5l+HZ+NI1Ho00gbuApT1RZyRypV7DIlxl0pKSmsW7eOyZMnF2X3YhOFQSgWkw73Y1/gfnwJkrM72vbvogsbCQoF2UYz60/G8cORG6RmG2lXswLPtg+igb8XdzL1uR3IZbEYODpl6hV8fhuJwqglbdAaTJUbF7yTkIddCsPp06dZvXo1+/fv57HHHmPGjBlFDlgcojAItqBKvoDnnjdwvnUYQ0A7Mru+j9mnNpBzF/K6k3GsOhJLms5ELV93riZn/a8YhFaiZ0glavm5y/tFlBOqlEt4bxyJwmIkddAazJUayB2pVLJZYTAYDPz++++sXr0aFxcXMjMz+fnnn1Gr5ZsTXxQGwWYkC+pza/A4MA+FWU9WixfJav4CqHIm49MaTPx8Io5DV1NoGejDvpgk1M4qlo1sInNw+/L+dQQAaUPlGXl4L1VSND6/jQIUpA6OwOwXKnekUstmhaFjx44MGDCAUaNGUbNmTSZOnMi3335rk5BFJQqDYGtKbQIe+95DfWkzpgohZHR7H5N/qwe223w2HoCBDauWdMQS5frvzwDowx6XNYfqThQ+m0YhKV1IG7IWc4U6suYp7Ww2XPXpp59my5Yt3Lx5kxEjRlDEbglBcGgWjypk9F6KPnQEmr1vUeGXoWQ3eAptuzeRXP83QqmsF4S75C4IAE63T+G96UkkZw9SB6/F4lNL7kjlglV9DIcPH2bdunX8/fffjBgxgsGDBxMSEmLPfA8lWgyCXRm0eBz+CLfTy7G4VSSz02wMwf1BocBktgDYfMU0h2M25vytkmddd6f4Y3hvfgrJ1YfUIWuxeAXKkqOssduopPT0dH777Tc2bNjAxo0bi5Kt2ERhEEqC0+3TaHa/hnPiWfRBPcjsNItJf+ZMZyH6GOzHKe4w3lvGILlVJHXIz1g8bT/FeXll1+GqchOFQSgxFhNup1fgEfkRWEysDXgLQ81e9G0cJHcyu3I9vwGgxKeudr6xH+/fx2HWVCNtyFosHuXj0l1JEYVBEGxIqY3H4+AC1Oc3YHavgrbdm+hDh4GijF9SKkHO1/fivXU8Zu+apA6OQHKvJHekMkcUBkGwA+ONo2j2z8Er8RjGKs3I7DgLU9XmcseyPWN2zt/OJTPRn8vVnXj9MQlzhTqkDl6D5OZXIuctbwpTGMRHHUGw0pT9zoxjFuk9PkGZcZMKGwbhufP/UGoT5I5mU95bxuC9ZUyJnMvl8p94/TERk189UoesFUVBZlbNrgqwefNmdu3albtIT7du3RgwYIDNgwmCoxrexB8Afb2mGGr3xf3YYtxOfotLzFayWr6Us0ayylXmlMWna/h0iZzH5dIWvLZPwVSpEWkDf8wzNFiQh9WXkmbMmMHs2bNzH8+aNYuZM2faPFhBxKUkwZEoU6+g2T8H16vbMHsFkdlxJoaavcBBlhB1VK7nf8Fz5yuYqrYkbcBKJBfHWcGurLLLpSSDwcCePXuIjo5m79696HS6IoUThNIqU28iU2/K85zFpxbp/VeQOvBHJJUL3lvH4735KVTJF2RKWXwKfToKfbrdju/678947ngZY7W2pA5YJYqCA7G6xZCdnc22bdtISEigSpUq9O7dW5Y5k0SLQZDLs2tPAY+4j8FsxO3sStwPf4zCqCW70TiyWv0fktqn5ELagD3vY1BH/Yjnnjcw1OhMWt/lJdbBLZTQqKSvv/5alqm3RWEQ5LLrYiIA3etWfOR2iuwkPCI/RB21Gkntg7bNa+jqPwlKVUnELDaXmK0AGIL72fS46tPf4fnPu+iDupPe52twkm8yzvLILoXh5Zdfzv23JElER0ezbds269MVkygMQmmhuhOFZt8MXOIiMfnVJ7PTLIwB7eSOJQu3E8vQHJiDvlZv0nt/WSY66Usbm6/5DKDRaJg3b17uYzk6ngVBTqlZOXMI+bgXbg4hc6UGpA1Zj+ulLXgcmIPPxnB0wQPQtn8Hi1d1e0YtFkV2MgCSm69Njud27As0hxaiCx5ARq/PZZuDSSiY1S2G2NhYatSokfs4NTUVHx8fW+cqkGgxCHIpsI/hUYzZuJ9YivuJL8FiIbvRWLJavGizN19bslkfgyThfuQTPI58jK7uEDJ6fgpKqz+TCjYi7nwWBDv4OyYJgM7BRb8JS5lxE/fDH6M+vw7JyZ3sZs+R1WQSuDjO0qAuV7YDYKjVq1jHcY/8CI+jn6Kr9zgZ3T4sNX0sZZVNC8N33333wHMajYaGDRsSFhZmfbpiEoVBKAtUyRfwOPQ+rlf+wuJWEW3Ll9E1GJ27elxp53bqWzT73iO73kgyu38o5pVyADYtDNOmTePs2bN069YNgD179tCoUSMuX75Mnz59mDRpUvHSWkkUBkEuiVoDABU9bPfm7RR/DI+DC3CJO4TZKxBt61fRhwyR9Y1Uob0NgORRuUj7u0avx2vnK+hr9yW991Jx+chB2LQwTJgwgcWLF+PhkdPU1Wq1vPTSSyxZsoRhw4axdevW4qW1kigMglyK1cfwKJKE8/U9eBxaiHNiFCa/MLRt38AQ1F2WO6iL08fgcmUbXn9MwhjQjrQBK8XoIwdi01FJcXFxODv/bxSBs7MzcXFxqNVqXFzKRrNXEApjbOsaBW9UFAoFxqBupAZ2wfXiJjwiP8T797EY/NugbfcmJv+W9jnvQ2Q1/0+R9nO+eQCvv57HVKkh6X2/FUWhFCp0YRgwYAAjR46kR48eSJLE7t276d+/P1lZWQQHB9szoyA4lPa17DyCSKFEHzIEfXA/1OfW4HHkUyr8MgR9zcfQtn0ds1+ofc//X8agblbv43T7NF6/j8fsFZgzIZ6Lxg7JBHuzalTS2bNnOXbsGJIk0aJFCxo1amTPbI8kLiUJcolPz5kfrKpXCd2xa8zC/dRy3E58icKQib7eCLStptn9HghlRhwAFs9qhdpelXIJn1+GITm7kzrsFyyawu0nlCybT6Ln5OSEUqnEyckpz2Wl4tq+fTvTpk3LfXzy5EnCw8MZNWoUX3zxhc3OIwi2MPOP88z843zJndDZnayWL5I85gDZTSfjenETvqs747HvPRTZSXY7reeOl/Dc8VKhtlVm3MR705OgUJI26CdRFEq5QheGlStX8uqrr5KSkkJSUhLTp09n1apVxQ4wd+5cFi1ahMViyX1u5syZLFq0iDVr1nDq1CmioqKKfR5BsJXxbQMZ3zawxM8rqSug7fAuyaP/QRc6DLfTK/Bd1QH3I5+gMGTa/HxZLV8mq+XLBW6nyE7Ce9OTKAwZpA38EbNPbZtnEUpWoS8lDRw4kLVr1+Lu7g5AVlYWI0eOZPPmzcUKsHXrVnx9fVm7di2ffPIJmZmZhIeH88cffwA5BcloNDJx4sQ8+4lLSUJ5p0q+iEfkB7he/gOLmx/aNtNzJukrwSGuCkMG3htH4pR8nrRBP2Gs1qbEzi0Ujc3nSlKpVPn+uzDWrVvHypUr8zw3f/58+vXrR2RkZO5zmZmZaDT/67Dy8PAgNjbWqnMJgj3dSM1ZC7m6j7xTRZt965Le9xucEk7gcWAunnveQB29jowuCzBXrF/s4yvTrgFg8Q7KfwOTDq+t43FKjCK933JRFMqQQheGYcOGER4eTq9evZAkiR07djB8+PBCnyg8PJzw8PACt9NoNGi12tzHWq0WLy+vQp9HEOxtzl85i+/Y/D6GIjJVaZYzSd+FDWj2z6HCz33JbjyerNbTijUqyHNXTr9fvvcxWEx4/fUCLjcPkt5zMYaaPYt8HsHxFLowPPPMM7Ru3Zrjx48jSRLvv/++XabC0Gg0ODs7c/36dWrUqMG+ffuYMmWKzc8jCEU1uf1DPkHLSaFAHzoCQ1APPA69j/upb3C9tJnMTrMw1O5XpBvkslpPy/8FyYLn7um4Xt1GRqc56EOHFTO84GgKLAzNmjVDcc9/qnu7JBQKBcePH7d5qFmzZvHqq69iNpvp2LEjTZo4xiczQQBoUcNH7ggPJakrkNl1Ibp64XjueRPvP59FH9SdzM5zsXhZ12Ge75oRkoTH/jmoo9ehbTUVXeNnbJRccCRidlVBsNLV5CwAavq6y5ykABYTbme+xz3yQxQWE1ktXyGr2bOFnqBPlRIDgLnC/25gdT+6GI/ID8hq9AzaTrNlmapDKB4x7bYg2IHd5kqyE2VmHJp9s3CN+R1ThTpkdpmPMaB9gfvdP1eS+uwPeO59C13IsJw1FcRMqaWSKAyCYAenbqYB0CTAW+Yk1nG5uhPNP++iSr+OLnQ4me3fRXJ/+LrVTreOAmDyb4nrhY14bn8RQ82eOes0i9XXSi1RGARByMuYjfuxz3E/sRTJ2R1t2zfRNXj0vQ8u13bhtXU8xqotSBv4IzjJO0xXKB5RGATBDi4l5gynrlPRcVZbs5Yq+SKav9/C5eZBjFWakdFlIeZKDfJukxSN6s5ZvPa+gcmnDmlDfkZyFUPHSztRGATBDkpbH8NDSVLuvQ8KXcp/7314NffeB5+f++GUGIXZK5DUYb8+8rKTUHqIwiAIdhAVn/N/r0HVgn/BSgOFLhWPQwtRR63G4lGFzI7vYarYgArrByIpVKSG/273mVyFkiMKgyAIheYUfxzN3jdxToxCcnJDclKTOvQXzL515Y4m2JAoDIJgB+dv58xkGlq5DC5C8997H1yj15Hd9FlMvqEP9D0IpZsoDIJgB2Wmj6EAxVnzWXBcojAIgh2U6RbDPVR3ctZBES2GskUUBkEQBCEPmy/tKQhCzqikuyOTyjKnhJM4JZyUO4YgA1EYBMFKi/deZvHey3LHsDuPA3PxODBX7hiCDMSlJEGwUlm487kwVEnRAJj96smcRLAl0ccgCIIg5CH6GATBDk7dTMudYbUsc7p1NHeGVaF8EYVBEKz05b6rfLnvqtwx7M7j0EI8Di2UO4YgA3EpSRCsVGpWcCum/FZwE0o/0ccgg7i4myxZ8ilpaWmYzSaCg0N44YUXcXcvekflwYP7iYj4EYVCgcViYcCAwTz2WF+2bt2Ml5cXHTt2ybP9oEG92bTpr0Ide/LkccyaNR9//2pFzicIQulRmMLgVAI5yg29Xscbb0zl9dffpUGDhgD88ccW3nvvbT744NMiH/ejjxbw/fdr8PT0JCtLy9ixT9KqVRv69Rtoo+SCNY7FpgLQooaPrDnszfnmQQCMAe1kTiKUtDJbGH6PSmDT2XibHnNQw6r0b1Dloa8fOLCPpk2b5xYFgL59B/Drr+uZM2cGrq6uxMffIikpkbfeeo/Q0Hrs2rWDtWtXo1Qqady4Kc8//+IDx/X19WXdujV07dqDWrVqs3r1OlxcXFi+fBl+fn4MHDiUDz6Yx5UrlwkIqI7BYAAgISGeDz6Yj8Ggx8XFlddee4sqVaqybNkSIiMPUqVKFdLSUm36PSoPvj5wDYBlI33kDWJn7ocXAWKupPJIdD7bUFzcTQICHpy33t+/GqdOnaBqVX8+/vgLhg8fyaZNv5CensaKFcv47LOlLF26nMTE2xw5cuiB/Rcu/BidTsesWW8zeHAfVq36jnuvAB46dACDwcDXX3/Ps89OQa/XAbBkyWeMGDGSzz9fxhNPPMVXX33B5cuXOHXqBN9++wPvvDOLrKws+31Dyqh3e4fwbu8QuWPYXUb3RWR0XyR3DEEGZbbF0L9BlUd+ureHSpUqc+5c1APP37gRS5MmzahbNxSAypWrcObMKW7ciCU1NYVXX30JgKysLG7evMnOnXO4cSMWH58KvPba28THx/PCCy/xwgsvcefObd5++zVCQ8Nyj3/lSgxhYTkTnVWtWpXKlXO+7suXL7Fq1XesXr0SACcnJ65cuUy9emEolUo8PDTUrl3Hrt+Tsqi6T/lY89jiHSR3BEEmZbYwyKFjxy788MMKzp07S/36OZeTNm/eiI9PBZRKJQqFIs/2/v4BVK5chU8//RInJye2bt1M3bohDBkyPHebpKREZsx4gy+//JYqVari51cRPz8/XFxccrcJCqrJjh1/AU+QmHiHO3fuABAYWJMnnniKRo2acO3aVU6cOEZgYBDr10dgsVjQ6/VcvVr2p3awtchrKQC0CaogcxL7co79BwBjjU4yJxFKmigMNuTu7s7773/C4sWLSE9Pw2QyU6dOXd57bx6LFz/YJK9QoQIjR45mypTJmM1m/P2r0b17rzzb+PlV5P/+bzpvv/0aKpUKi8VM+/adaN26LWfO5KwL0KlTV06fPsWkSWOpWtUfHx8fAP7zn5dZtGghBoMBvV7Hyy+/St26oXTr1pOJE5+mYsVKVKjga/fvS1mz4tB1oOwXBvejnwGQJgpDuSOGqwqCleLTc/pwqnqpZU5iX8qMOAAsnmIoc1kihqsKgh2U9YJwlygI5ZcYlSQIVjpwJZkDV5LljmF3ztd243xtt9wxBBmIFoMgWGnl4VgA2tcq2/0z7seXAJAW1E3mJEJJE30MgmClRG3ODYQVPVwK2LJ0U2hvAyB5VJY5iWBLoo9BEOygrBeEu0RBKL9EH4MgWOnvmCT+jkmSO4bduVzZjsuV7XLHEGQgWgw29Pnnn3D+/L8kJyeh0+moVi0AH58KzJ37fpGPabFYWLLkU2JiLqFUKnFycubll6cREFCdmTPf5J13ZuPs7Jy7/aFDB9i5cxtvv/1egce+du0qH344ny+++LrI+cqj1UdvANA52E/mJPbldnIZAIZavQrYUihrRGGwoRdf/D8Atm7dzLVrV/OdEM9akZEHSEy8w6effgnA33/v4fPPP2bhwo+ZNWtBsY8vWO/9gfXljlAi0vuIDwzlVZktDK7R61H/G2HTY+rCRqGvN6LQ21+8eJ5vvlnKBx98yvbtf/LjjytZuXINp06d5M8/f+eFF15izpx30Wq1mM1mJk16nhYtWuU5RuXKVYmO/pedO7fRokVrOnXqQrt2HQAYMWIgq1ev59atOBYsmI1a7YabmxpPTy+AfGduTUxMZPbsd5AkCV/fsv2J11583J0L3qgMkNzK9qgr4eFEH4Md1a0bSnz8LfR6PZGRB1EoFCQnJ7F//166dOnGypXLadmyDUuWfMOcOQtZuHAOFoslzzGCg+vw+uvv8Pffexgz5nEmTBjD2bOn82zz7bdLmTjxWT777EsaNmwM8NCZWyMifqRnz958/vkyOnfuWlLfijJl18VEdl1MlDuG3bnEbMUlZqvcMQQZlNkWg77eCKs+3dtL69btOHHiGLdvJ/DYY304evQwJ0+eYPLk/7Bhw1oee6wPkDMzq7u7B//+e46lSxcD0KdPP+rVa0BgYBCzZs1HkiSOHIlkxow386zQduXKZcLCcibta9SoKdeuXX3ozK1Xrlymd+9+/922Cb/+Kubat9ba4zcB6F63osxJ7Mvt9AoADMH9ZE4ilLQyWxgcRefOXfn66y+pWzeU1q3b8eGH86lRowZOTk4EBdXi1KmThITU486d22RkpBMaWi9PZ3BExI9cunSRN9+cgUqlolat2qjVbnlmag0MrMnZs6dp27Y90dE5034/bObW69evEhV1mrp1Q/j333Ml/v0oCxYNaSB3hBKR3m+F3BEEmYjCYGeNGjUhNvYao0c/TZ06dYmPv8WTTz4NwNNPP8OCBbPZs2cner2e1157GyenvD+SESNGsWTJZ4wfPxp3dw+USiXvvjsrzzbTpr3BzJlvsmbNKnx8fHBxcX3ozK0TJz7PzJlvsmPHNqpVCyix70NZonEtH782kquX3BEEmYg7nwXBStuic+4Ifqxe2b4BzPXiJgD0dQfJnESwJXHnsyDYwYZTt4CyXxjUZ38ARGEoj0SLQRCspDOaAVA7q2ROYmfG7Jy/ncvHUqblhWgxCIIdlPmCcJcoCOWWuI9BEKy09VwCW88lyB3D7lzPb8D1/Aa5YwgyEC0GQbDSb2fiAehXv4rMSexLfW4NAPrQ4TInEUqa6GMQBCuZzDl3pzupyniD22zM+VtVPqYAKS8K08fgEP+zt2/fzrRp03Ifb9u2jZ49ezJmzBjGjBnD4cOHZUwnCHk5qZRlvyhATkEQRaFckv1S0ty5c9m3bx9hYWG5z0VFRTF9+nR69+4tYzJByN/mszmXkgY2rCpzEvty/fdnAPRhj8ucRChpsn/sad68Oe+9916e56KiotiwYQNPPvkkCxcuxGQyyRNOEPKxJSqBLVFlv/NZHf0z6uif5Y4hyKDE+hjWrVvHypUr8zw3f/58GjduTGRkJBEREXzyyScAfPfdd/Ts2ZPq1aszc+ZMQkJCeOqpp0oipiAIQrlXYpeSwsPDCQ8PL9S2w4cPx8srZ56WHj168NdffxWwhyAIgmArsl9Kup8kSQwaNIj4+JzruAcPHqRBg/Ixm6UgCIIjkL3z+X4KhYK5c+cyZcoU1Go1wcHBPP646PwSBEEoKaXyPobt27fz559/smjRIiBneOsHH3yAv78/AC+++CKtW7eWM6LwEPf/7E6ePMm8efNQqVR07NiRKVOmyJxQKIgkSXTu3JmaNWsC0LRp0zzDzQXHY7FYeO+99zh//jwuLi7MnTuXoKCgh27vcC2GgojhraVXfj+7mTNn8vnnn1OjRg0mT55MVFSUuHTo4K5fv06DBg346quv5I4iFNKOHTswGAysXbuWkydPsnDhQpYuXfrQ7R2uj6EgYnhr6XX/zy4zMxODwUBgYCAKhYKOHTty8OBB+QIKhRIVFUVCQgJjxoxh0qRJXL58We5IQgGOHTtGp06dgJwW3tmzZx+5vcO2GB42vLVfv35ERkbmeb5Dhw55hrdGRESI4a0yKuzPLjMzE41Gk/vYw8OD2NjYEsspFCy/n+WMGTOYPHkyffv25ejRo0yfPp0NG8Rke47s/t81lUqFyWR6YMXIuxy2MIjhraVXYX92Go0GrVab+1ir1eb+HAXHkN/PMjs7G5UqZ+rxli1bkpCQgCRJedYhFxzL/b9rFovloUUBSuGlpPuJ4a2ll0ajwdnZmevXryNJEvv27aNly5ZyxxIK8MUXX+S2IqKjo6lWrZooCg6uefPm/P3330DOgI+QkJBHbu+wLYbCEsNbS7dZs2bx6quvYjab6dixI02aNJE7klCAyZMnM336dPbu3YtKpWLBggVyRxIK0KtXL/bv38+oUaOQJIn58+c/cvtSOVxVEARBsJ9SfylJEARBsC1RGARBEIQ8RGEQBEEQ8hCFQRAEQchDFAZBEAQhD1EYBEEQhDxEYRAEQRDyEIVBEP6rWbNmdj2+Tqfjqaeewmw2A3DgwAGmT5+eZxuDwcDo0aPFRJCCrERhEIQSsmHDBnr16pU7z1B0dDT169fPs42Liwvt2rVj69atckQUBEAUBkF4wHfffceAAQMYMGAA33//fe7zS5YsoU+fPjzzzDNMnTqV5cuXW3XczZs306NHj9zH0dHR3LlzhyeffJIOHTpw4MABAHr27MnmzZtt8rUIQlGU+rmSBMGWzp49yy+//MLPP/+MJEk8/vjjtG7dGrPZzLZt29i4cSMmk4lhw4ZZNVmjwWAgNjaW6tWr5z4XHR3NwIED+emnn9i2bRubN2+mffv21K1blzNnztjjyxOEQhGFQRDucezYMXr27Im7uzuQM/nY0aNHsVgs9OjRA7VaDUC3bt1y94mNjWXp0qVkZmayePFisrKymDVrFs7OzrRu3ZpBgwaRkpKCp6dn7j5Go5G0tDQmTJgAgMlkyn1dpVLh7Oz8wBz6glBSxKUkQbhHUeaUrFGjRp7ZKrdt20bv3r2ZO3cuu3btAkCtVmMwGHK3iYmJoV69eiiVOb+C58+fp27durmvGwwGXF1di/plCEKxiMIgCPdo1aoVO3bsIDs7m6ysLHbs2EHLli1p3rw5u3fvRq/Xo9Vq2bNnz0OPkZCQgL+/P0BuR7O3tzdmsxm9Xg/kXEaqV69e7j7nz58nNDQUgJSUFHx9fXF2drbTVykIjyYuJQnCPRo0aMCwYcNyVy0bMWJE7sih7t27M2jQIAICAmjYsGGeS0P3qlKlCvHx8YSFhWGxWHKf79ChA8eOHaN9+/ZER0fTuHHj3NcuXryYu3hKZGQkXbp0sdeXKAgFEusxCEIhabVaPDw8yM7OZvTo0cyZM4cGDRqQkpLCJ598woEDBwgPD2fMmDHMmTMHFxcXWrRowaBBgwA4d+4c3333HR9++OEjzzNlyhSmTp1K7dq1S+LLEoQHiMIgCIU0bdo0Ll26hF6vZ+jQoTz77LNWH2P9+vUMHTo09xLT/QwGA1u3bmXIkCHFTCsIRScKgyAIgpCH6HwWBEEQ8hCFQRAEQchDFAZBEAQhD1EYBEEQhDxEYRAEQRDyEIVBEARByEMUBkEQBCEPURgEQRCEPP4foFewsH8I3vIAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots()\n", "ax.plot(c,d1, label='One-Sided')\n", "ax.plot(c,d2, label='Two-Sided')\n", "ax.axvline(e1, color='C0', linestyle=':')\n", "ax.axvline(e2, color='C1',linestyle=':')\n", "\n", "ax.set(title='Error in Numerical Derivatives',\n", " xlabel='$\\log_{10}(h)$',\n", " ylabel='$\\log_{10}$ Approximation Error',\n", " xlim=[-15, 0], xticks=np.arange(-15,5,5),\n", " ylim=[-15, 5], yticks=np.arange(-15,10,5)\n", " )\n", "\n", "ax.annotate('$\\sqrt{\\epsilon}$', (e1+.25, 2), color='C0')\n", "ax.annotate('$\\sqrt[3]{\\epsilon}$', (e2 +.25, 2),color='C1')\n", "ax.legend(loc='lower left');" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3.9.7 ('base')", "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.9.7" }, "vscode": { "interpreter": { "hash": "ad2bdc8ecc057115af97d19610ffacc2b4e99fae6737bb82f5d7fb13d2f2c186" } } }, "nbformat": 4, "nbformat_minor": 4 }