CourseContent/Notebooks/SolveMathProblems/1-IntroductionToNumericalSolutions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical models\n",
    "\n",
    "We start with the numerical solution of a very simple differential\n",
    "equation. In fact we choose something simple enough that we already \n",
    "know the answer.\n",
    "\n",
    "\\\\[\n",
    "        \\frac{d\\theta}{dt} = - k \\theta\n",
    "\\\\]\n",
    "\n",
    "This is the equation which governs radioactive decay, in which case\n",
    "\\\\(\\theta \\\\) is the amount of the radioactive isotope remaining and \\\\(d\\theta /  dt\\\\)\n",
    "is the activity that we can measure. \\\\(k \\\\) is closely related to the half life.\n",
    "\n",
    "The solution to this equation is\n",
    "\\\\[\n",
    "        \\theta(t) = \\theta_0 e^{-kt}\n",
    "\\\\]\n",
    "\n",
    "where \\\\(\\theta_0\\\\) is the amount of the radioactive material remaining.\n",
    "The same equation also describes the cooling of, say, a cup of coffee. In this\n",
    "case we interpret \\\\( \\theta \\\\) as the excess temperature (above room temperature). \n",
    "\n",
    "<img src=\"images/theta_t1.png\" width=\"66%\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xd4VNW+xvHv2lPSeyMVEkhIQofQexVQQBQVFREFywFsWI7KUc9REctFxd4LoCIKAipVegm9d0IgEEoSShJCSJvZ94+dQIBAEjKTyUzW53nmCZns2ft37pWXlVWFqqpIkiRJjkWxdQGSJEmS5clwlyRJckAy3CVJkhyQDHdJkiQHJMNdkiTJAclwlyRJckAy3CVJkhyQDHdJkiQHJMNdkiTJAelt9WB/f3+1Xr16tnq8JEmSXdq8efNpVVUDyrvOZuFer149Nm3aZKvHS5Ik2SUhREpFrpPdMpIkSQ5IhrskSZIDkuEuSZLkgGS4S5IkOSAZ7pIkSQ6o3HAXQnwnhEgXQuy6zs+FEOIjIUSSEGKHEKKl5cuUJEmSKqMiLfcfgL43+Hk/ILr49SjwedXLkiRJkqqi3HBXVXUlcPYGlwwCpqiadYC3ECLYUgVe7e9PPua7Ua+Qk3PRWo+QJEmye5bocw8FjpX6PrX4vWsIIR4VQmwSQmzKyMi4qYel7Unjor47r73/HuuSz9zUPSRJkhxdtQ6oqqr6laqqCaqqJgQElLt6tkyBwZEAtCtYwNCvEnlp1g6yLhZaskxJkiS7Z4lwPw6El/o+rPg9q/Dw9ANAmDN5p2k6v248Rq/3VzBv50lUVbXWYyVJkuyKJcJ9LjC8eNZMOyBLVdWTFrhvmdzdvQBI03lzT85U5o7pSKCHE6N/2sIjUzZzMkv2xUuSJFVkKuQvQCLQUAiRKoQYKYR4XAjxePEl84BkIAn4GhhttWoBnYszAJlOgXBiC40vJDJnTEde7h/L6qQMer+/kqmJRzCbZStekqTaq9xdIVVVvbecn6vAGItVVA4t3C9q/eyBkbBsAvqYvjzapT59GwXz8h87eWXObmZvO8HbdzQhOsijukqTJEmqMexuhWpJy92Ue5HznZ+GUzth318ARPi5MnVkGybd1YxDGTn0/2gV7y8+QH6RyZYlS5IkVTu7C3e9mxbuzgUKySFNwS8alr0FZjMAQgjubBXGknFdubVJMB8tOUi/D1exJum0LcuWJEmqVnYX7iUtd+dCHcnnU6Dbi5C+B/bMvuI6P3cnPhzagikPt8Gkqtz/zXqe/GUr6efzbFG2JElStbK/cHcyAOBaqCc5KxkaDYaAOFg+EUxF11zfJSaAhU934ame0SzYdYqe/7eCKYlHMMkBV0mSHJjdhbuiCAB8cNfCXdFB95fh9AHY/kuZn3E26HimdwwLnu5Ms3BvXp2zm8GfrWFHamZ1li5JklRt7C/cdVq4e6quJGcma2/GDYDQBK31Xnj9ee5RAe5MHdmGj+5twcmsPAZ9uoZX5+wiO0+ucJUkybHYbbi7mowczzlOXlEeCAG9/gvZx2HD1zf8vBCCgc1CWPJsV4a3q8u0dSn0nLSCOduOyxWukiQ5DDsMd61kZ5MBFZUj2Ue0H0R2hga9YNUkuFh+d4uns4H/DWrMnDGdCPZy5qnp2xj27XoOZeRYsXpJkqTqYYfhrrXcDcU9KZe6ZgB6vgZ5mbBmcoXv1yTMiz9Gd+SNQY3YkZpF3w9XMnH+XnLyrx2clSRJshd2F+6ieEBVFIIiFG1QtURwU2hyF6z7HLIrvr2NThE80L4eS5/txu3NQ/lyRTI9Jy2XXTWSJNktuwv3kpa7ml9AuEf4leEO0H08mItgxTuVvneAhxPv3dWMWaM7EOihddXc89U69p7MtkTpkiRJ1cb+wr245W7KLyTKK4pDmYeuvMA3EhIegi1T4HTSTT2jZYQPs8d0ZOIdTTiYdp5bP1rFa3N2kZUrZ9VIkmQf7C/ciwdUzSYzMR5RpGSnUGAquPKiLs+D3hmWvn7Tz9EpgnvbRLDsuW4Ma1eXqetS6D5pOdM3HJU7TkqSVOPZYbgXd8sIHTFO4ZhU07VdM+6B0PFJ2DMHjq6v0vO8XY28Pqgxfz3RmfoBbrw4aye3f7aGrUfPVem+kiRJ1mS/4a7oqK/XzuE+eO7gtRd2eAI8gmHhS5c2FauK+BBPZjzWnslDm3MqK4/Bn63l+d+2y71qJEmqkewv3JXLLfc6eGFUjGWHu9ENer4KxzfDrpkWebYQgkHNQ1n6XDce6xrF7G3H6f7ecj5bnkReodxWWJKkmsP+wr245W4WCiL3IlHeURw4d6Dsi5sOheBm8M9/b7gtQWW5O+l5qV8ci57pSocG/ry7YL88x1WSpBrFDsNdK1kVOswXLhDtHV12yx1AUaDPBMhOhcRPLV5LpL8bXw9P4KdRbXF30jP6py3c8+U6dqZmWfxZkiRJlWF/4V6qW8aUk0O0TzTpF9PJyr9OoEZ2htjbYPUHcD7NKjV1bODP30925q3BTTiUkcPAT1fz3G/bScuW/fGSJNmG/YX7pdkyCuacC0T7RANcv2sGoPfrUJQPyyZYrS6dIrivbQTLnu/Go12imLvtBN3/bzmfLD0o++MlSap2dhfuQhEgirtlcnKI9tbC/bpdMwB+9aHNo7B1KpzaZdX6PJ0NvNQvjsXjutAlOoD/W3SAnpNWMHf7CdkfL0lStbG7cAet9a4KBfOFCwS6BuJp9ORg5g3CHaDr8+DsBYvGQzWEbF0/N754oBXTH22Hl4uBJ3/Zyp2fr2XTkbNWf7YkSZJ9hrsiUJ2cMeWcRwhBtM8NBlVLuPhAt5cgeTns+7ta6gRoF+XHn0904t07m5J67iJDvkjksamb5NbCkiRZlX2Gu04BozPm81pAxvjEcPDcQcxqOYuVEkZCYDwseMmiUyPLo1MEd7cOZ/nz3Xi2dwyrD56mzwcreWX2LjLO51dbHZIk1R72Ge6KAKMzpvPabo3RPtHkFuVyIufEjT+o00P/9yDrKKz+sBoqvZKrUc8TPaNZ8UJ37m8bwS8bjtLtvWV8tOQguQVy/3hJkizHPsNdJ1CNzpizzwNUbFC1RL1O0HiINjXy7GFrlnld/u5OvD6oMYue6ULn6ADeX3yAbu9pm5IVmaq+VYIkSZLdhjtGI6bzxeFePB2y3EHVEn3eAEUPC8dbq8QKiQpw54sHWjHzX+0J83HhxVk76f/RKpbuS5MzayRJqhK7DXfVYMScrXXLuBncCHMPY9/ZfRW7gWcIdH0B9v8NBxdbsdKKaVXXl5n/6sAXw1pSaFJ5+IdN3Pv1OrYfK/8sWEmSpLLYabgroHe61HIHiPOLY++ZvRW/SbvR4NcA5v9bW+BkY0II+jYOZtEzXXh9UCMOpuUw6NM1PD51M0np58u/gSRJUil2Ge5CEaA3YD5/HrV4O9843zhSc1I5X1DBINQbod+7cPaQVfaduVkGncLw9vVY8UJ3nukVw+okbWbNc79t59jZXFuXJ0mSnahQuAsh+goh9gshkoQQL5bx8wghxDIhxFYhxA4hRH/Ll3qZohOoegOoKuYLFwCI9Y0FqHjXDECDntq+Myvfg8xj1ij1prk76XmqVzQrX+jOyE6RzN1+gh6TlvPfubvl9ElJkspVbrgLIXTAp0A/IB64VwgRf9Vl/wFmqKraAhgKfGbpQktTFIGqMwBc6neP84sDqFzXDEDfidrXec9Xy8rVyvJ1MzL+1nhWPN+NIa3CmLouha7vLeP/Fu4nO0+e6SpJUtkq0nJvAySpqpqsqmoBMB0YdNU1KuBZ/GcvoJwJ51Wj6ASqogO41O/u7+JPgEtA5VruAN4R2srVA/Nh31+WLtVigr1cmHhHUxY/04UesYF8siyJzu8s44sVh7hYIDcmkyTpShUJ91CgdJ9FavF7pf0XGCaESAXmAU9YpLrr0MJdD4CpuOUOxYOqZyvZcgdo9y8IagLzXoC87PKvt6GoAHc+ua8lfz3RiRYR3rw9fx9d31vGtHUpFMo58pIkFbPUgOq9wA+qqoYB/YGpQohr7i2EeFQIsUkIsSkjI+OmH6ZtHKa13M2lZszE+sZyOOsweUWV3EddZ4ABH8L5k1bdFtiSGod68cNDbZjxWHsifF35z+xd9Ji0nBmbjsmFUJIkVSjcjwPhpb4PK36vtJHADABVVRMBZ8D/6hupqvqVqqoJqqomBAQE3FzFaFMhL3XLZJeaDukbh0k1VWyl6tXCEqD1KFj/pXbuqp1oE+nLb4+35/sRrfF2MfLC7zvo9f4KZm5OlSEvSbVYRcJ9IxAthIgUQhjRBkznXnXNUaAngBAiDi3cb75pXg5FEajFpZvPX9ktA9xc1wxAz1fAPQj+fApM9rPXixCC7rGBzB3bka+HJ+Bq1PPsb9vp88FK5mw7jslc8waKJUmyrnLDXVXVImAssBDYizYrZrcQ4nUhxMDiy54FHhFCbAd+AUaoVlw/r+gEZrQTmUq33EPcQvA0et58uDt7Qf934dROWP+FJUqtVkIIescH8feTnfhiWCuMeoWnpm/jlg9X8uf2E5hlyEtSraGvyEWqqs5DGygt/d6rpf68B+ho2dKuT1EEZjMobm5XtNyFEMT5xrHvTCVnzJQWNxBi+mp97/EDtdk0dkZb7VqHPvFBLNh9ig8WH+CJX7byydIknuoVTd9GdS6dRStJkmOyyxWqik6gmlUUT88rWu6gDaoeOHeAQvNNzgEXQtsWGKF1z9TAue8VpSiC/k2CWfB0Fz66twWFZjOjf9pC/49WsXD3Kbk5mSQ5MLsMd6ETmE1mdB4el/Z0LxHrF0uBuYDDWVXYztc7Anr/Dw4tha3Tqlit7ekUwcBmISx+pisf3tOc/CIzj03dzK0frWb+zpOyu0aSHJBdhruiUzCbVBRPj0t7upeI99MWz+45s6dqD0kYCXU7adsCZ1t1TVa10SmC21uEsviZLvzfXc3IKzTxr5+2cMuHcuBVkhyNfYa7IjCbVHQenlfsDAlQz7Me7gZ3dp3eVdWHwMCPwFQAfz5t190zV9PrFIa0CmPxuK5MHtocIeCp6dvo9f4KZmw6JhdDSZIDsM9w1xWHu6cH5qysK38mFBr5N2Ln6Z1Vf5Bffej5KhxcCDt+rfr9ahidIhjUPJQFT3Xhi2EtcTXqeOH3HXT/v+X8tD6F/CK5rYEk2Sv7DXezilJGyx2giX8TDpw9QL7JArsntn0Mwttq+76fP1X1+9VAiqLtJf/XE534bkQC/u5OjP9jF13fXc73aw7LvWskyQ7ZZ7grxQOqnh6Yc3Iu7eleorF/Y4rUospvIlbmw3Qw6FMoyoO/n3Wo7pmrCSHoERvEH6M7MG1kWyL8XPnfn3vo/O5SvlxxiAv59rOwS5JqO/sMd72CuUhruaOqV+wvA1rLHah6v3sJ/2jo/rK2a+SumZa5Zw0mhKBTtD8zHmvPr4+2Iy7Yk4nz99HxnaV8sPgAZy8U2LpESZLKYZfhrtMrl7plgGu6ZgJdAwl0DbRMv3uJ9mMhrDX8Pc5hZs9URNsoP6aObMsfozuQUNeHyUsO0vHtpfx37m6OZ160dXmSJF2HnYa7trpSeHkDYDp37UHSTfybWK7lDlr3zOAvwVQIs0eDuXbNKGkR4cM3D7Zm0TNd6N8kmGnrUuj67jLGzdjGgTR5xqsk1TR2Ge6KTitbuBe33DOvDffG/o1JyU4hKz/rmp/dNL/6cMsESF4GG7+23H3tSEyQB5PubsaKF7ozvH095u88RZ8PVjLqx41sOnLW1uVJklTMLsO9pOWOx/XDvaTffffp3ZZ9eKuHILoPLH4VMvZb9t52JNTbhVcHxLP2xR480yuGzSnnGPJFIkM+X8uSvWly1ask2ZhdhntJyx234nDPurZ1Hu8Xj0BYtt8dtL1nBn4CBleY9QgU1e7BRR83I0/1imbNiz3474B4TmblMfLHTfSbvIpZW1LlgihJshG7DPdLLXdXd6DslruH0YNIr0jL9rtfunkQDJgMJ7fDynctf3875GrUM6JjJMuf78YH9zQDYNyM7XR9dxlfrTwkD/OWpGpml+Fe0nJXUVA8PMoMd9D63Xec3mGd3Q/jB0Lz+2HVJDi63vL3t1MGncLgFmEseLoz341IIMLPlbfm7aP9W0t4/c89HDuba+sSJalWsMtw1+m1sk1FZnReXmV2ywA09W/K2byzpOakWqeQvm+DV5jWPZNnwYFbB1CyIGr6o+3564lO9I4PYkriEbq+t4wxP29h27Gy/0GWJMky7DTctW4Zs0lF5+193ZZ788DmAGxL32adQpw94Y5vICvV4TYXs6TGoV58OLQFq/7dnUc6R7HyQAa3f7qGu75Yy4Jdp+RulJJkBXYZ7krplvsNwr2BdwPcDe5sTd9qvWIi2kKP8bB7Fmydar3nOIBgLxde6h9H4ks9efU2bfD18Wmb6TlpOVMSj5BbILc3kCRLsctw1+mKz08tCffrdMvoFB3NAppZN9wBOj4DkV1h3guQboH9bBycu5OehztFsvy5bnx2f0t83Iy8Omc37Scu5d0F+0jLzrN1iZJk9+wy3Eta7uYiVetzv07LHbSumUOZh8guyL7uNVUvSIE7vgKjG/z+EBTKZfkVodcp9G8SzB+jOzLzX+3pUN+PL1YcouPbS3nil61sTjknjwKUpJtkl+GuK54tYzJpLXdzdjZqUdm/0rcIbIGKyvb07dYtyqOOtj1B+h5Y+LJ1n+WAWtX15fNhrVj+XHce6liP5fvTufPztQz6dA2ztqTKveUlqZLsM9wNxQOqRdqAKly7eViJJv5N0Amd9btmAKJ7QYcnYdN3sGeO9Z/ngCL8XBl/azzrXurJm7c3JrfAxLgZ2+n49lLeX3yAdNllI0kVYpfhrlzVcoeyNw8DcDW4Eusby7YMK82YuVqPVyC0Fcx5As4cqp5nOiA3Jz3D2tVl8TNdmDayLc3Dvfl46UE6vL2Up6ZvZevRc7YuUZJqNLsM90tTIYvM6Ly9gLJXqZZoEdiCnRk7KTRXwypJvRGGfK9tUzDjQdn/XkUle8t/82Brlj/XjQc71GPp3nQGf7aWQZ+s5o+tsstGkspil+F+ueVeqlsm68aDqnmmPPadqaaZLD514c5vIG2Xw5/eVJ3q+rnxym3xrHu5J28MasT5/CKe+XU7Hd9exvuL9nMyS/5DKkkl7DLcL61QLSzVLZN5/RWiLQJbAFRPv3uJ6N7Q9QXY9hNsmVJ9z60F3Jz0PNC+Hv8805UpD7ehaZgXHy9LouPbS3lkyiZWHsiQu1JKtZ7e1gXcDEV35QpVuHG3TKBrIKHuoWzL2MZwhldLjQB0/TekboR5z0NwUwhpUX3PrgUURdAlJoAuMQEcO5vLzxuOMmPjMRbvSaOunyv3t43grlbh+LgZbV2qJFU7+2y5Gy6vUFXc3UGnu2G3DECroFZsTttcvfOmFZ22PYFbAMwYDrnyMAtrCfd15d99Y1n7Ug8mD21OkIczb83bR9uJSxj36za2HJVz5qXaxT7D/VLL3YwQotyFTAAJQQmczTtLUmZSdZR4mZsf3P0jZJ+EPx6rdcfzVTcnvY5BzUOZ8Xh7Fj7dhXsSwlm0J407PlvLrR+t5uf1R7mQL7c5kByfXYa7UAQIMBVpLTGdrw+mszeeGtcmuA0AG09ttHp91whLgL4T4eAiWDah+p9fSzWs48Ebtzdm3cs9mTC4MSrw8h87affWEl6ds0ue/So5NLvscxdCoNMpmItP+dH7+lF09swNPxPqHkqIWwgbT23kvrj7qqPMK7UepR3user/IKgRNL6j+muopdyd9Nzfti73tYlgy9FMpq1LYfqGY0xJTKFlhDdD20RwW9NgXI12+ddBkspUoZa7EKKvEGK/ECJJCPHida65WwixRwixWwjxs2XLvJaiF5gKi1vufr6YzpTfn926Tms2pm3ErNqga0QIuHUShLeFOWPg5I7qr6GWE0LQqq4PH9zTnHUv9+Q/t8aRdbGQF37fQZsJS3j5j53sTJX78kuOodxwF0LogE+BfkA8cK8QIv6qa6KBl4COqqo2Ap62Qq1X0OkUTFe03MsP9zbBbcjKz+LguYPWLq9seie4Zxq4+MD0++DCadvUIeHrZmRU5yj+GdeV3x5vzy2N6jBrSyoDPlnNrR+tYmriEXk0oGTXKtJybwMkqaqarKpqATAdGHTVNY8An6qqeg5AVdV0y5Z5LZ1eYC7Swl3n54s5Kwu18MZ/GVsHtQZs1O9ewj0Qhv4EFzK0GTQmGSC2JISgdT1fJt3djPUv9+KNQY1QVXhlzm7aTPiHZ2dsZ+ORs3KmjWR3KhLuocCxUt+nFr9XWgwQI4RYI4RYJ4ToW9aNhBCPCiE2CSE2ZWRk3FzFxRS9gsmk/YXT+/oBUHTuxoOqwe7BhLmHseHUhio9u8pCWsDATyBlDcz/t21rkS7xcjHwQPt6zHuqM3+O7cQdLcNYuPsUd32RSO8PVvLNqmTO5OTbukxJqhBLzZbRA9FAN+Be4GshhPfVF6mq+pWqqgmqqiYEBARU6YE6vXK55e7rA4Cpgl0zm9M226bfvbSmd0HHp2DTt7D+K9vWIl2jSZgXbw1uwobxPXl3SFM8nfW8+fde2k1cwuifNrNsXzpFJjmtVaq5KjI94DgQXur7sOL3SksF1quqWggcFkIcQAt7q/V/KDpxueXuV9xyP3PjGTOgDarOOjiL/Wf3E+cXZ63yKqbna3A6CRb8G7wjoGGZv/BINuRq1HN3Qjh3J4RzIO080zccY/a248zbeYoADyfuaBHKkFZhRAd52LpUSbpCRVruG4FoIUSkEMIIDAXmXnXNbLRWO0IIf7RummQL1nkNnV7BdKnl7gtUrOVe0u9u864Z0Faw3vk11GkKvz+sTZWUaqyYIA9eHaDtNf/lA61oHu7Nt6sP0/uDlQz6dA1T16WQlSvHUKSaodxwV1W1CBgLLAT2AjNUVd0thHhdCDGw+LKFwBkhxB5gGfC8qqrlN6OrQNFdHlCtTMs9yC2ISK9IEk8kWrO8ijO6wX2/ajNofr4Hsq7+pUiqaYx6hVsa1eHr4QmXplTmF5p4ZfYuWr/1D2N+3sLy/emY5OZlkg1VaNWGqqrzgHlXvfdqqT+rwLjiV7XQWu7aXx7FwwMMhgrNdQfoGNKR3w78Rl5RHs56Z2uWWTEedeD+GfDtLVrAPzwfnOSv+fbA392JUZ2jGNkpkt0nsvl9cypzth3n7x0nCfJ0YnCLMIa0CqNBoLutS5VqGbvcfgCKp0IWD2gJIdD7+pa7SrVEh5AO5Jvy2ZK2xZolVk5QI7j7B+0M1t8eApPc/8SeCCFoHOrFfwc2Yt3LPfliWEsah3jx9apker2/gts/XcPUxCOcvVBg61KlWsJuw10p1XIHrd+9vP1lSrQKaoVBMbD2xFprlXdzGvTSVrEmLYZ5z8lDPuyUk15H38bBfDuiNYkv9WB8/zguFpguzZ0f9eNG/tpxgrxCeYKUZD12u5lG6b1lgEq13F0NrrQMasmaE2t4juesVeLNSXgIMlNg9Qdad023Mnd7kOxEoIczj3SJYlTnSPaePM+cbceZve04/+xNx91JT7/GdRjcIpS2UX7oFGHrciUHYrfhrujFlS13P18KUlIq/PmOIR15f/P7pF1II8gtyBol3ryer0FOOiyfqO0F33qkrSuSqkgIQXyIJ/EhnrzQN5b1yWf4Y+tx5u86xW+bU6nj6cyg5iHc3iKUuGBPW5crOQC77ZbR6RRMhaVb7hXbX6ZEh5AOADWvawa0TcYGfAQxfbUzWHfPtnVFkgXpFEGHBv68d1czNv2nF5/c14LGoZ58u/ow/Savou+HK/lixSFOZMozYaWbZ7/hbrg8zx20lruam4s5N7dCn4/xicHfxb9mhjuATg9Dvtd2kZz1CBxeaeuKJCtwNui4rWkI3zzYmg3jtb1tXI063p6/j47vLGXoV4n8suEo5+RArFRJdhvueoNC0VUtd6DCrXchBB1COpB4MhGTuYYObBld4b7p4FsffrlPLnJycL5uRh5oX49Zozuy4vluPN0zhvTsfF6atZPWE/5hxPcbmLk5Ve5WKVWI3Ya7znBVt0yAPwBFldiQrENIB7Lys9hzZo/F67MYFx94YBa4eMO0IXDmkK0rkqpBXT83nuoVzZJnu/LXE50Y2TmSg2k5PPvbdhLe/IdHp2ziz+0nyC2QU2alstntgGpJt4yqqto898BAAIrSKx7uHUM6ogiF5anLaRLQxFqlVp1nCAybBd/3gx8HwkPzwKeurauSqkHJ/PnGoV682DeWrccy+XP7Cf7ecZJFe9JwMejoGRfIgGYhdI0JwNmgs3XJUg1hty13vUErvaTfXV+8y2RlWu7ezt40D2jOimMrLF+gpQXEwPDZUHAepgyE7BO2rkiqZkIIWkb48NqARiS+1JPpj7bjjpahrD10hsembqb1m9r+88v2p1Mod6ys9ew43LUWSlFB8eZhPj6g11OUXrlzQrqHd2f/uf2cyLGDsKzTBIb9ARfOaC34HKufiSLVUDpF0C7KjwmDm7D+5Z78+HAb+jauw6I9p3jo+420mfAPL83awYoDGTLoaym7DXfdVS13oSjo/f0rHe5dw7sCsCLVDlrvAGGt4P7fIPs4TLkdcis+/VNyTAadQteYgEtTK78ZnkDn6ADmbjvBg99tIOHNf3jut+0s2ZtGflENnTwgWZzd9rlf6pYpPagaGFipbhmASK9I6nnWY8WxFdwbe69Fa7Sauu1h6M/aJmNTB8PwOdqAq1TrOel19IoPold8EHmFJlYdPM38XSdZuPsUv29OxcNJT8+4QPo1CZZ99A7ObsO9pOVe0i0DWr974dGjlb5X17Cu/LzvZy4UXsDN4GaxGq2qfne4ZypMvx+m3anNqHH2snVVUg3ibNDROz6I3vFBFBSZWXPoNPN3agOxs7edwNWoo3tsIP0bB9OtYQBuTnYbB1IZ7LdbRn9ltwyAPjCg0i13gG7h3Sg0F9bcBU3XE3ML3PWDNv99yiC4WLGN06Tax6hX6N4wkHeHNGPj+F5MG9mW21uEsj75DGN+3kLLNxbz2NRNzNl2nPNyHr1DsNt/qvXG4pZ74ZUtd1NoGhqtAAAgAElEQVRmJuaCAhSjscL3ah7YHE+jJ8uPLad33d4Wr9Wq4m7TWvAzhmuDrMPngKuvrauSajCDTqFTtD+dov15Y1BjNh45y/ydJ1mw+xQLd6dh1Cl0aOBHn/g69IoPJNCjBpx5IFWa/YZ7SbdMqW1TDcVz3U0ZGSihoRW/l6KnS1gXVqaupMhchF6xs/+zNOwH9/wEvw4rDvjZ4OZv66okO1Ay66ZdlB+vDWjE1mPnmLfzFIv3pPHyHzsZPxtahHvTO74OfRoFUT9AHjpiL+y4W0YbCDJd1XIHKKzkjBmAnhE9yczPZFPaJssUWN1i+sC9P8OZg/DjAMipfPeUVLspiqBVXV9euS2eFc93Y8HTnXmmVwyFJpV3Fuyj56QV9Jy0nLfn72PL0XOY5TGCNZqdNVEvK+mWuXq2DFRuIVOJTqGdcNG7sOjIItoFt7NMkdWtQS/tPNafh8IPt2pdNJ7Btq5KskNCCGLreBJbx5Mne0ZzPPMi/+xJY/GeNL5ZlcwXKw4R4OFEr7gg+jQKokN9P5z0cuZNTWK34V4yoHp1nztUbguCEs56Z7qEdWHJ0SWMbzsenWKn/6FGddPmwf8yFL67Reui8Y2ydVWSnQv1duHBDvV4sEM9snILWbY/nUV7TjF323F+2XAUN6OObg0D6dMoiG4xgXi5Gmxdcq1nv+Fexjx3na8v6HSVXshUok/dPiw8spDNaZtpE9zGInXaRGRnGD4XfroTvuur7UtTp7Gtq5IchJergdtbhHJ7i1DyCk0kHjrDoj2nWLwnnb93nkSnCFrV9aFnbCA94wKpH+COEPKUqepmt+Fe1mwZoSjoAwJuOtw7h3XWumZSFtl3uIO2kvWhBdoipx/6w32/QURbW1clORhngzZXvntsIBNuV9l6LJNl+9JZsi+difP3MXH+PiJ8XekRG0iP2EDaRvnK7ptqYrfhfrlb5srl1PqgQIrS027qni56FzqFduKflH94qc1L9ts1UyIwFh5eAFNv1173TNX65SXJCpTiFnuruj48d0tDTmReZOm+dJbtS+eXDUf5Ye0RXI06Okf70zM2iG6xAXKapRXZbbiXtf0AgCE4hPx9+276vn3q9WFxymK2pG+hdZ3WVaqxRvCpCw8vhGl3aAOtd3wFje+wdVVSLRDi7cKwdnUZ1q4uFwtMJCafZsnedJbuS2fhbq0B1izMix6xQfSIDaRRiCeKPCTcYuw23BWdgqKIMsI9mJxlyy7t815ZXUK74KxzZtGRRY4R7gDugfDgX9og6+8Pw/lT0H60rauSahEXo644xINQVZW9J8+zdF8aS/el8+GSA3zwzwECPZzo3jCQrg0D6NjAHy8XOShbFXYb7qANqhaVEe5qfj6mc+fQ+1Z+paarwZXOYZ1ZlLKIF9q8gEFxkP/AXLzhgT+081gXvgSZR+GWCWDvXU+S3RFCEB/iSXyIJ2N7RHMmJ5/l+zNYui+debtO8uumY+gUQYtwb7rGBNC1YQCNQ7xkq76S7D7cr26564PrAFB48uRNhTvAgKgBLE5ZTOKJRLqEdalynTWGwQXu+hEW/QfWfQZZx+DOb7T3JclG/NyduLNVGHe2CqPIZGbbsUxWHMhgxYEMJi0+wKTFB/BzM9I52p+uDQPoHB2Av7uTrcuu8ew63PUGhaKia/vcAYpOnoRGjW7qvp1CO+Ht5M2fh/50rHAHraXedyJ4hcPCl7XVrPdOl9sVSDWCXqeQUM+XhHq+PNunIadz8ll98DQrDmSw8kAGs7dph+o0CfW61KpvEe6NXme3i+2txq7DXWdQMBVcOVvGEKKtyCw8cfKm72vQGbil3i3MTppNTkEO7kYH3E+j/WjwCtO6ab7pBcNmgl99W1clSVfwd3e6NKfebFbZfSKbFQfSWXEgg89XHOKTZUl4OOvp1MCfrjEBdI4JINRb/iYKdh7u+jL63HU+PggnJwpP3ny4AwyoP4Bf9//K4pTFDI4eXKV71VjxA8Gjjnboxzc9tS6bqK62rkqSyqQogiZhXjQJ82Jsj2iyLhayNun0pS6c+btOARDl70bHBtqul+2i/GrtwGyFwl0I0ReYDOiAb1RVffs6190J/A60VlXV6jtw6fTKFfu5F9eAITi4yuHe1L8pER4R/JX8l+OGO0B4G3hkCfxyr7bgqd870OYRW1clSeXycjHQr0kw/ZoEo6oqB9NzWHXwNGuSTjNzSypT16WgCGga5k2n4rBvEeFdaxZRlRvuQggd8CnQG0gFNgoh5qqquueq6zyAp4D11ii0LHqj7oqTmC69H1xH63OvAiEEt0XdxufbP+fUhVPUcatTpfvVaL5RMHKx1kUz7zlI36uFvK52tngk+yOEICbIg5ggD0Z2iqSgSBuYXZ10mtUHL3fhuBh0tIn0pXO0Px0b+NMwyMNhZ+FUpOXeBkhSVTUZQAgxHRgE7LnqujeAd4DnLVrhDegMCvm5Rde8bwgO4cLq1VW+/21Rt/HZ9s/489CfPNLUwVuzzp7auaxL/gdrJsPpA3D3FHnwh2SXjHqFNpG+tIn0ZVzvGLLzClmffJbVBzNYnXSaN//eC4C/u5EO9bVWfacG/oQ4UH99RcI9FDhW6vtU4IpNSoQQLYFwVVX/FkJUW7gbjDouZOZf+35wMEUZGagFBYhKnMh0tXDPcBKCEph1cBYjm4xEEQ4+Iq/ooPfrENgI5j4BX3eHob9AULytK5OkKvF0Nlw6TxbgZNZFVhd34axOOsPc7dosnEh/N9pF+dG+vh/tonztenuEKg+oCiEU4H1gRAWufRR4FCAiIqKqj0bvpFB01WwZKJ4xo6oUpqdjDAur0jOGxAzhxVUvsuHUBvvd572ymt2jzZyZfp820DrwY2gyxNZVSZLFBHu5cFdCOHclhKOqKvvTzrP64GkSD53hr+0n+GXDUQDqB7jRvr4f7aP8aRfli58dza+vSLgfB8JLfR9W/F4JD6AxsLx4uX8dYK4QYuDVg6qqqn4FfAWQkJBQ5WNcDEYdhfllhbs2170wNbXK4d6rbi+8Nnjx+4Hfa0+4A4QlwGMr4bcRMHMkHNsAfd4E/c3/JiRJNVHpg0lGdY6iyGRm94lsEpPPkHjoDLO2HGfaOi3sY4LcaV/csm8b6YePW839+1CRcN8IRAshItFCfShwX8kPVVXNAi6tgBFCLAeeq47ZMnqn64R7uPZbQcHRo7i1q1ogO+mcGBA1gOn7p3M27yy+zrWoD9qjDjz4Jyx+DdZ9Cie2wt0/gmeIrSuTJKvR6xSahXvTLNybx7vWp9BkZufxLBIPnWFd8hlmbErlx8QUAGLreBS37LWwr0mHlJQb7qqqFgkhxgIL0aZCfqeq6m4hxOvAJlVV51q7yOsxOGmzZVSziig14m0IrgMGA4XHjt3g0xU3JGYI0/ZOY27SXEY0HmGRe9oNnQH6vgXhrWHOWPiyCwz5DiIdbOWuJF2HQafQMsKHlhE+jOnegIIiMztSM0k8dIbE5DP8vP4o3685ghAQH+xJ20g/2kb50rqeL742bNlXqM9dVdV5wLyr3nv1Otd2q3pZFWMwavNViwrNGJwuz10VOh3GkBAKjqVa5Dn1vevTIrAFMw/O5MFGD9bOU2UaDYbAePj1AZgyCHr8Bzo+A4qDDzJL0lWM+stbJDzRM5r8IhPbjmaSmKy17H9an8J3aw4DWjeONmvHj7aRvgR5Vt8ArV2vUC0J9MJ80xXhDmCIiKDw6FGLPWtIzBDGrx7PhlMbaBtcS080CmgIjyyFP5+EJa9D8gptf3gPB14DIEnlcNLraBvlR9soPwDyi0zsTM1i/eGzbDh8lj9K9dnX9XOlTT1fhraJoFVdH6vWZdfNLr3xcrhfzRgeTsGxY6hqlcdtAbil3i34OPkwbe80i9zPbjm5w53fajNojm2AzzvCwcW2rkqSagwnvY6Eer6M6d6AHx9uw/bX+vDn2E7859Y4YoI8WLw3jdRzuVavwyFa7mVOh4wIx3z+PKbMTPQ+Vf8X0knnxJCYIXyz8xuOZR8j3DO8/A85KiGg5XAIb6sd/vHTEGg/Fnq+JmfTSNJV9Drl0p44ozpHYTarmCzU6LwRu265l+6WuZoxXAtfSw2qAgyNHYpO6Ph5388Wu6ddC2gIo5ZA60cg8RP4tjecOWTrqiSpRlMUgaEatii283DXyi8sq+VeHO4FRy0X7oGugfSu15vZSbO5UHjBYve1awZnuPX/4J6f4NwR+KITbPwGqqFlIknS9dl1uF/qc8+7UcvdcoOqAMPihpFTmMOcpDkWva/di7sN/rVW66r5+1mYdidkn7B1VZJUa9l1uN+oz11xcUEfEGDRljtA04CmNPVvys/7fsZkvva5tZpXqHZO662T4GgifNYOdsyQrXhJsgGHCPey+twBjHXrUpCSYvHnDm80nJTsFJYcXWLxe9s9IaD1KHh8Nfg31LYR/u1BuHDG1pVJUq1i1+Fe0i1T1p7uAMaoKAqSky3+3F4RvajnWY9vdn5jsamWDsevPjy8QJtBs28efNYWds2UrXhJqiZ2He6XW+7X7ukO4FQ/ClNmJkVnz1r0uTpFx8ONH2bv2b2sPbHWovd2KIoOOo+DR5dr57X+/rB24lPW8fI+KUlSFdl1uOv0CooiKMy/fssdoOCQ5afn3RZ1G0GuQXy982uL39vh1GkMI//RdpVMXg6ftoWN34K57P+/SZJUdXYd7gAGZ12ZUyEBnIrDPT/5sOWfqzMwotEINqdtZlv6Novf3+Ho9NDhCRi9FkJbwt/j4Idb4fRBW1cmSQ7J7sNdb9RRdJ0BVX1wMMLFhYJk6yysuSP6DnycfPhixxdWub9D8o2C4XNg0KeQvlvbvmDZW1B40daVSZJDsftwN1xnT3cAoSgYI+uRf8jyg6oArgZXRjQewZrja9iavtUqz3BIQkCLYTBmI8QNgBXvaF01+xfYujJJchgOHe4ATlH1rTJjpsS9sffi7+LPR1s+kjNnKssjCIZ8qx0IoneGX+6Bn4dqK10lSaoSuw93o7OOgryyZ8uANmOm8MQJzLnW2YXNRe/CI00eYVPaJhJPJlrlGQ4vsos2L77363B4pdaKX/EuFObZujJJslv2H+4uevJzrx/uxsjiQdXDlh9ULTEkZgjBbsF8vOVj2Xq/WXojdHwKxm6EmL6wbAJ82gZ2z5Zz4yXpJth9uDu56Cm4eIOWe0wMAPn7D1itBqPOyL+a/YtdZ3ax9OhSqz2nVvAK1c5pfWA2GN211a3f94Pjm21dmSTZFbsPd2M54W6sG4FwcSF//z6r1jGg/gCivKJ4f/P7FJoKrfqsWqF+d3h8FQyYDGeS4OseMOtRuQBKkirIMcI934RqLvtXd6HT4RQTTd5e64a7XtHzXMJzHD1/lF/2/WLVZ9Uaig5ajYAntkCncVoXzcetYOkEyD9v6+okqUZziHBHhYIbzJhxjo0jb/9+q/eHdw7rTMfQjnyx/QvO5Z2z6rNqFWdP6PWa1h8f2x9WvguTm0PiZ3LQVZKuw+7D3clFOykwP/f6XSHOsQ0xZ2VRdPKk1et5PuF5coty+WzbZ1Z/Vq3jUxeGfAejlkJQI1j4ktaS3zIVTNfvmpOk2sjuw91YHO4FF28w1z02FoC8ffutXk997/oMiRnCbwd+48A56w3i1mphreDBudpKV48gmDtW2ztezqyRpEvsPtydLoX79VtuzjExIAR5+/ZWS01jm4/F0+jJ64mvY1bl5lhWE9VNO8P1nmla//xvD8JXXbUthmXIS7Wc3Ye7sQLhrri5YYyIIN/Kg6olvJ29eTbhWbZnbGfmwZnV8sxaSwhtC4N/rYXbv4C8bJh+L3zRGfbMkTtPSrWWA4S7tqd7/g3CHcC5UTwXd++qjpIAGFh/IK3rtOaDzR9w+uLpanturaXooPm9MHYTDP4Sii7CjOHweQftkBB5JKJUyzhAuJffcgdwbtqUohMnKUxPr46yEELwn3b/Ia8oj/c2vlctz5TQthZuNhTGbIA7vwXVrB0S8lk72D4d5BoEqZaw+3C/NFumnHB3adYMgLwdO6xeU4koryhGNRnFvMPz5MrV6qbooMkQGL0O7voBFAP88VjxFMpP5Tx5yeHZfbjrDAqKTpTfco+PB4OBi9urL9wBHmnyCHG+cfwv8X+czbPscX9SBSgKNBqsbUx23wzwqQcLX4b3G8E//4Xzp2xdoSRZhd2HuxACJ9cbb0EAoDg54Rwby8Xt26upMo1BZ+DNTm9yvuA8b657U24sZiuKAjG3wEN/a/Pk63eDNZPhwyYwZwykV89guyRVF7sPdwAnVwN5F8rvS3Vp2pSLu3ahmqp3cC3GJ4YxzcewOGUxfyX/Va3PlsoQ1grungJPbIaWD8LOmfBZW5g6WDswRA6+Sg6gQuEuhOgrhNgvhEgSQrxYxs/HCSH2CCF2CCGWCCHqWr7U63NxN5CXU4Fwb94MNTeX/IPVf27niEYjaBHYggnrJ3A0+2i1P18qg28U3Pp/8Mxu6P4fSN+rHRjycUtY+zFclFtISPar3HAXQuiAT4F+QDxwrxAi/qrLtgIJqqo2BX4H3rV0oTfi7G7gYkXCvUULAHI3Vf/2sTpFxzud30EndDy34jnyTfnVXoN0HW5+0PV5eHqnNvjqEQKL/gOT4mDuk5C229YVSlKlVaTl3gZIUlU1WVXVAmA6MKj0BaqqLlNVteSoo3VAmGXLvDEXD2OFwt0YFoYhJITcDRuqoaprBbsHM6HTBPae3SunR9ZEOoM2+PrwfG0AtuldsGOGNlf+u76w7RcosM6JXpJkaRUJ91DgWKnvU4vfu56RwPyyfiCEeFQIsUkIsSkjI6PiVZbDubhbpiKDla5t2pC7YQOqjVYudgvvxoPxD/Lr/l9ZcEQeCF1j1WkCAz+GcXug9xtwIQNmPw6TYuHvZ+Fk9Q7MS1JlWXRAVQgxDEgAymyWqqr6laqqCaqqJgQEBFjsuS7uBlSzesPj9kq4tm2LKTPTJv3uJZ5q9RRNA5ry6ppX2X/W+puZSVXg6gsdn9RWvo6YBw37wtZp8GUX+LIrbPxW2/JAkmqYioT7cSC81Pdhxe9dQQjRCxgPDFRVtVo7lF3cDQAVGlR1a9sGgNz1661a040YFAMfdvsQD6MHTyx9gjMXz9isFqmChIB6HeGOr+DZfdDvPTAXwd/jYFJDmPkIJP0jtx6WaoyKhPtGIFoIESmEMAJDgbmlLxBCtAC+RAv26lnfX4qzhxGgQv3uhpAQDOHhXFhvm373EgGuAXzU4yPO5Z1j3PJxFJgKbFqPVAkuPtD2Ua1fftRSaHo3HFwI0+6ED+Jhwctat41c0yDZULnhrqpqETAWWAjsBWaoqrpbCPG6EGJg8WXvAe7Ab0KIbUKIude5nVVcbrlXLCDd2rcnd9061ALbBmojv0a80fENtqRv4X+J/5MLnOyNENqc+QGT4dkD2tz50ATY8JXWbfNZe1j9gTz3VbIJfUUuUlV1HjDvqvdeLfXnXhauq1Kci8O9Ii13APduXcmcMYPczZtxa9/emqWVq29kXw5nH+azbZ/h5+LHuFbjbFqPdJMMzhA/SHvlntV2otzxq7bFwT//g4j22kyc+IHgUcfW1Uq1QIXCvaZzcde6ZSrS5w7g1q4dwmgkZ/kKm4c7wONNH+fMxTN8v+t7fJ18GdF4hK1LkqrC1RfaPKK9zhyCnb/Dntkw/3mY/wLU7QDxt8ugl6zKIbYfMDjp0BsVcrMr1s2iuLri2rYtOcuXW7ewChJC8FKbl7il3i1M2jyJ2UmzbV2SZCl+9aHbv2F0IoxeD91e1Fr285/XplV+fyts+BqyT9i6UsnBOES4A7h5O3Ehs+KTdNy7daUgJYX85MNWrKridIqOtzq9Rfvg9ry29jX+PPSnrUuSLC0wVgv3Meu0rYi7/lubPz/vOXg/Dr7qBive01bEyvEXqYocJtzdfSoX7h7duwNw/p9/rFVSpRl1Rib3mEzrOq0Zv3q8bME7ssA46P4SjN2gBX3PV0HoYNmb2orYyc1g/otweKWcXindFMcJd29ncs5VPNwNISG4NG9O9t9/W7GqynPRu/BJj09oF9yOV9e8yqyDs2xdkmRtgXHQ+Vl4ZAk8u1+bfRMQC5u+gx8HwHv1YeYo2P4r5FhuZbfk2BxiQBUud8uoZhWhiAp9xrN/f9Leeov8pCScGjSwcoUV56x35qMeH/H0sqd5be1rnM07y8jGIxGiYv+7JDvmUQdajdBeBRfg0FLYNw8OLoKdv2nXBDeDBr20V1hrbU8cSbqK47TcfZwwm9UKT4cE8OzXFxSF7HllboVjU856Zz7u8TH9I/szectkJqyfgEnuM167GN0gbgAM/hyeOwiProAer4DBDVZ/CN/3g3ej4NdhsPkHyJRbSUuXOVTLHSDnXB6unsYKfUYfEIBrmzZk/fUX/mPHIJSa9W+dQWdgYueJBLkG8f3u78nIzWBi54m4GlxtXZpU3RQFQpprry7PQV4WJK/QtjxIWgJ7iwfgfepBvc4Q2RUiO8uplrWYw4S7u48W7hcy86ESR4V43zGYEy/8m9z162vEnPerKUJhXMI4Al0DeXfjuzww/wEmd59MmEe17qos1TTOXto8+fiB2syajH1a2B9ZBXvnwtap2nX+MRDZRXvV66zNwZdqBWGrJe8JCQnqpk2bLHa/C1n5/PDvNXQZGkOTbhUPPnN+PklduuLarh1hkz+0WD3WsPr4al5Y+QKKUHi3y7t0COlg65KkmshsglM74PAqbbZNyloovKD9LKixtlo2op328pKNBHsjhNisqmpCedc5TMvd1cOI3qCQlXGxUp9TnJzwGjyYs9OmUZSRgd6CWxFbWqfQTvx66688uexJ/vXPvxjdbDSjmoxCp+hsXZpUkyg6CGmhvTo+CaZCOLEVDq/QAn/bz7Dxa+1ar3AIb3s57APjtc9Lds9hWu4A09/YgIevE7eOaVapz+UnHya5f3/8x4wh4ImxFq3JGnILc/lv4n+Zf3g+rYJaMbHTRILdg21dlmQvTEWQtguOroNj67Sv509qP3Py1GbghLeF0FYQ2lJ25dQwFW25O1S4L/hyJ2dOXOD+/7Wr9GePjR7Dxc2babB0CYqbm0XrsgZVVfkz+U8mrJuATtHxSrtX6Fuvr5wuKVWeqmozbY6th6OJcHQ9pO8BirPBJ1IL+dBW2qtOUzDKQX1bqXXdMgBeQa4c3n4as8mMoqvczBf/Rx/hyNClZP7+O74PPmilCi1HCMHA+gNpEdCCF1e9yAsrX2De4XmMbzueOm5yhoRUCUKAT13t1fRu7b28bG1P+uObtdexDdpOl6CtpA2MLw78ltq8+4A4bWdMqcZwqJb73rUnWDplH/e/3g7vwMq3LFIeGE5BSgr1Fy5AcXGxaG3WVGQu4qe9P/HJ1k/QKTqeavkU9zS8B0XUrKmdkp07nwYntsDxLZdDPy9T+5miB/+G2tmzpV+yS8fiamW3zMmkTGb93xZuG9uMuo39Kv353E2bSBn2AAFPP43/449ZtLbqcOz8Md5IfIPEk4nE+8XzfMLzJNQp978BSbo5qgrnDsOpnXByh/b11E44X2qHS6/wK8M+IA58I+WgbRXUynDPu1DIt8+uov3g+rS8pRKT3Us5NmYsuevWUX/RQvR+lf8HwtZUVWXe4Xl8sPkD0nLT6BXRi3GtxhHuGV7+hyXJEnIyIG3nlaF/5iCoZu3nOidt/n1grLaHTmCc9tWnngz9CqiV4Q4wZfxagiI9uWVU45v6fH5yMskDBuI1cCAhE9+ycHXV52LRRabsnsK3u76l0FzIHQ3uYFSTUXJWjWQbBbmQsRfS913+mr4XslMvX6N3Af/oy2HvHwN+DbSWvt7JdrXXMLU23Bd8uZPTqTkMe+PmV5umv/8BZ776iojvvsWtg30vFErPTefL7V8yK0nbXVKGvFSj5GVDxv5rg790145QwDtCC/qrX56h2tYMtUitDfdN84+wfk4yoz7ogpPLzU0GMuflcXjQ7ahFRUT+MQudp6eFq6x+J3NO8u2ub5l5cCaocEvkLTwQ9wCN/BvZujRJutbFTDh7SDum8PRBOJNU/Dp0ebUtgN4ZfOuDfwNtyqZPveJXXa2/3wF3zKy14Z6y+wx/fbydQU83Jyz25kfqc7duJeWB4Xh070boRx85zPzxkzknmbJnCn8k/cGFwgu0CGzB/XH30yO8BwYH/IsgORhV1RZclYT96aTLf848CuZSu8IKRdteoSTwvesW/7n4HwFXX20aqJ2pteFecLGIb55dRctbImg3qH6V7nXm+x9If+cdAl94Ab+HH7JQhTVDTkEOs5Nm89Pen0jNScXHyYfb6t/G4AaDifaJtnV5klR5ZpN2Fu25I5CZon299EqBC+lXXm9018LfK0zr3rn6z56hNXLufq0Nd4CZ725CVWHIv6s2DVBVVY4/9TTnFy8m9P1JePbrZ6EKaw6T2cSaE2uYnTSbZceWUWQuopFfI/pH9qdPvT5yQZTkOAouaCFfOvizUrVX9nHtPNuruQWUEfyh4F5H207Zo4627341qtXhvn5uMpvnH2HkpM44uVatq8Gcl8fRUaO4uH0H4Z9+gnuXLhaqsuY5l3eOv5P/Zs6hOew7uw+ApgFN6VO3D73q9iLUPdTGFUqSFRXmaSGffbw49I9D1rErvy84f+3nnDzBPehy2HvUAY/g4veCte/dg8DJ3SJl1upwP5Wcxcx3N9PzwThi21d9VogpO5uUESPIP5hEyMSJeN12qwWqrNlSslNYnLKYRUcWsffsXgCivKLoFNqJjqEdSQhKwKir2KEokuQw8rK0kM85pa3YPX8Scoq/lv6+KO/azxo9wD1Qe3V4AmJvLkdqdbirqsq0VxLx9Hdh0NMtLHJP0/nzpI4eQ+7Gjfg/MRb/xx9H6GrHgouj2UdZfmw5q4+vZlPaJgrNhbjoXWgZ1JKEoARaBj81Pl4AAApPSURBVLaksX9jGfaSBNqgb15mGeF/CnLStVf7MRDb/6ZuX6vDHS53zQx7oz2e/pbZJ8acn8+pV18la85c3Dq0J/jNNzGEhFjk3vYitzCXTWmbWJW6io2nNnIo6xAARsVIk4Aml4I+3i+eINcgh5llJEk1Ra0P95xzeUx9JZG4DiF0u6+hxe6rqipZM2dy6s0JoCgEjBmD7wPDEMba2Wo9l3eOLelb2JKmvfae3YtJ1Q7y9nX2Jd4vnni/eOJ846jvXZ9wj3D0ikNtRipJ1arWhzvA8p/2sXftSe4Z3wbfEMuOaBekHiftzTfJWb4cfUgw/o8+htcdg1FqaciXuFh0kf1n97PnzB7tdXYPyZnJlwJfr+ip51mPKK8ooryjqO9Vn3DPcMLcw/By8rJx9ZJU88lwB3KzC/jlf+vx9Hdm8LMt0Rst20euqioXVq/h9CefcHH7dnReXnjdPgivO+7AKSZGdkkUyyvKIykziUOZh0jOSiY5M5nkrGRSc1Ixl2wmBXgYPAjzCCPUPfTS1zpudQh0DSTQNRAfJx95pKBU61k03IUQfYHJgA74RlXVt6/6uRPw/+2dW2wcZxWAvzOzN9vr2E4cW2lIYhKlanNBpVSQAAKiUqjykAhBaSoqQIqo0gpeEA+gSigqD4gHeECqBJGooCBogAdkiVYVKS6RoiQ0kkOTpiRNQ6Abt7ZzceKNHXu9e3iY8Xp9293Yu7O380nj+S9nZs7x/8+ZmX/+nfMi8DHgGvC4ql7Ot88gnDvApf5hXjl0hp5tq3hk/1YisdIPCagqYydPcuPwYUaPvAapFOG1a4nv2kXLpz5J0wMPEOroKPlxa52J9ASXb14mMZogkUyQGE1wJXmFRDLBQHKAifTELHlXXFY1raKrqYvVzavpau5iZWwl7dF2b4m10xHtoCPWQVu0jaZQ7XyT3zCKpWTOXURc4ALwCJAA3gCeUNVzOTLPAB9R1QMisg/4kqo+nm+/QTl3gLP/SHD0pQu0dzfz6cc2s27LyrLdVU9dv87o346Q7Ovj9vHj6ITnoCI9PcS2biWyaSPRjZuIbtpIaM09uPHqD+lXCTKa4er4VYbGhhgaG2J4bJjBsUGGx4cZHhtmaNwrG5kYWXQfMTdGe6yd1kgrreFWWsItxMNx4pE48XDcy/vpeDhOS6SFplATMTdGLBTLrptCTYSdsD2JGVVBKZ37TuCgqn7Rz/8AQFV/nCPzqi9zXERCwAfAas2z8yCdO8B7b1+n73f/ZvTaHVZ0xtiwvZOuDa20dzXT3BahuTWCG3ZKegJnxse5c/YsY/2nGe/vZ+L8eVIDA7NknHicUHc34e4u3M5O3BVtuCtacVpX+OtWnOYWnGgEiUaRaBTHX0s0ikSiSDiEOA64rrduIKYyU9yavMXInRFuTNxgZGJkJu2vRydHuZ26TTKVJDmZJJlKcjt1e96TQT4ccYi60VnOP+pGiYViRJwIITdE2AnPXtyF8xEnkk2HnBCuuDji4IqL67jZfMgJ4YjjpcVLT9dnt8mVlxCO42TzgiAi2bXDTP+erp+O1iUiM9vkbOdIac8JY/mUMobqWuC9nHwC+MRiMqo6JSI3gVXA1eLULT/r7l/J1w7u4J1Tg7xzapC3jw1wpi8zW0ggFHYIRVzckON9U0i8jj/dv0VkXllhtkP3dugGVNHJSTSVQqemsgvpNDo8BYMZNDOt1x1/WeBn0fnIVUwEmTZOFqgvvLNZq7tUZBnbLo8QbXTSRieFgrYoM7cgmvM3p0xzy/KncjaZXZBDyl9yU0Z5qNbLUte2JHu+V95ob4HOSRORp4CnANavXx/koQFwww737VzDfTvXkMkoIx+McevaOGM3JxkbnSSdypCaTJOezJCeyqA6c2Jn017GK8sUOOCSUc/hp6bQqRSaTkM6g2bSkPGdfzrtrzOeIjrtpHwlswvovLJFFNd5idnV+Z7yCmy7WHFjMX3x8PrPdD+adVHJlpGVzaZ1RnZ+3dxLzSJtOP03bzMt3ljz6udf1/JuvxwKaVXGA+QRWdqxm9rLPzOsGOd+BciN0fYhv2whmYQ/LNOG92J1Fqp6CDgE3rDMUhQuFY4jrLynpeRTJA3DMKqBYgZo3wA2i8iHRSQC7AN658j0At/w018B/p5vvN0wDMMoLwXv3P0x9G8Dr+JNhXxBVd8SkeeAU6raC/wK+K2IXASu410ADMMwjApR1Ji7qr4MvDyn7Ic56TvAY6VVzTAMw1gqjTVvzjAMo0Ew524YhlGHmHM3DMOoQ8y5G4Zh1CHm3A3DMOqQin3yV0SGgf8ucfNOqujTBgFhNjcGZnNjsBybN6jq6kJCFXPuy0FEThXz4Zx6wmxuDMzmxiAIm21YxjAMow4x524YhlGH1KpzP1RpBSqA2dwYmM2NQdltrskxd8MwDCM/tXrnbhiGYeShqp27iDwqIudF5KKIfH+B+qiIHPbrT4pIT/BalpYibP6uiJwTkTdF5DURKRRqqOopZHOO3JdFREWk5mdWFGOziHzVb+u3ROT3QetYaoro2+tFpE9E+v3+vbsSepYKEXlBRIZE5Owi9SIiP/f/H2+KyIMlVUBVq3LB+7zwu8BGIAL8C9gyR+YZ4Bd+eh9wuNJ6B2DzLqDZTz/dCDb7cq3AUeAE8FCl9Q6gnTcD/UCHn++qtN4B2HwIeNpPbwEuV1rvZdr8GeBB4Owi9buBV/CiAe4ATpby+NV85/5x4KKqXlLVSeAlYO8cmb3Ab/z0n4GHpbaj+Ra0WVX7VHXMz57Ai4xVyxTTzgA/An6CF1S21inG5m8Bz6vqDQBVHQpYx1JTjM0KrPDTbcAANYyqHsWLb7EYe4EX1eME0C4ia0p1/Gp27gsF5l67mIyqTgHTgblrlWJszmU/3pW/lilos/+4uk5V/xqkYmWkmHa+F7hXRI6JyAkReTQw7cpDMTYfBJ4UkQRe/IjvBKNaxbjb8/2uCDRAtlE6RORJ4CHgs5XWpZyIiAP8DPhmhVUJmhDe0Mzn8J7OjorIdlUdqahW5eUJ4Neq+lMR2YkX3W2bavlC0dcz1XznfjeBuckXmLuGKMZmROTzwLPAHlWdCEi3clHI5lZgG/C6iFzGG5vsrfGXqsW0cwLoVdWUqv4HuIDn7GuVYmzeD/wRQFWPAzG8b7DUK0Wd70ulmp17IwbmLmiziHwU+CWeY6/1cVgoYLOq3lTVTlXtUdUevPcMe1T1VGXULQnF9O2/4N21IyKdeMM0l4JUssQUY/P/gIcBROR+POc+HKiWwdILfN2fNbMDuKmq75ds75V+o1zgbfNuvDuWd4Fn/bLn8E5u8Br/T8BF4J/AxkrrHIDNR4BB4LS/9FZa53LbPEf2dWp8tkyR7Sx4w1HngDPAvkrrHIDNW4BjeDNpTgNfqLTOy7T3D8D7QArvSWw/cAA4kNPGz/v/jzOl7tf2C1XDMIw6pJqHZQzDMIwlYs7dMAyjDjHnbhiGUYeYczcMw6hDzLkbhmHUIebcDcMw6hBz7oZhGHWIOXfDMIw65P83vdkhvWFyrQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%pylab inline\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "theta_0 = 1.0\n",
    "\n",
    "time_values = np.linspace(0,1.0,1000)\n",
    "\n",
    "for const_k in [1.0, 3.1, 10.0, 31, 100.0]:\n",
    "\n",
    "    exact_theta_values = theta_0 * np.exp(-const_k * time_values)\n",
    "    plot(time_values, exact_theta_values)\n",
    "\n",
    "\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to be able to march forward in time from our starting point\n",
    "(just like the picture above)\n",
    "where $\\theta = \\theta_0$ to obtain the value of $\\theta$ at\n",
    "later times. To do this, we need to approximate the original\n",
    "differential equation, and, in particular, the value of the time\n",
    "derivative at each time. There are a number of ways to do this.\n",
    "\n",
    "### First order numerical approximation\n",
    "\n",
    "Assume that the variation in \\\\(\\theta(t) \\\\) is linear, i.e.\n",
    "\\\\[\n",
    "    \\theta(t') = \\theta_n + \\beta t'\n",
    "\\\\]\n",
    "\n",
    "where we use a local time coordinate \\\\(t' = t - n\\Delta t\\\\), so that when we differentiate\n",
    "\n",
    "\\\\[\n",
    "        \\frac{d \\theta}{dt} = \\beta\n",
    "\\\\]\n",
    "\n",
    "To determine the approximation for the derivative therefore\n",
    "becomes the solution to the following equation:\n",
    "\n",
    "\\\\[\n",
    "        \\begin{split}\n",
    "            & \\theta_{n+1} = \\theta_n + \\beta \\Delta t \\\\\n",
    "            & \\Rightarrow\t\\beta = \\frac{d \\theta}{dt} = \\frac{\\theta_{n+1} - \\theta_n}{\\Delta t}\n",
    "        \\end{split}\n",
    "\\\\]\n",
    "\n",
    "This is a first order difference expression for the derivative which we\n",
    "substitute into the original differential equation for radioactive decay at\n",
    "the current timestep\n",
    "\n",
    "\\\\[\n",
    "        \\frac{\\theta_{n+1} - \\theta_n}{\\Delta t} = - k \\theta_n\n",
    "\\\\]\n",
    "\n",
    "This rearranges to give us a time-marching algorithm:\n",
    "\n",
    "\\\\[\n",
    "        \\theta_{n+1} = \\theta_n (1-k \\Delta t)\n",
    "\\\\]\n",
    "\n",
    "It is an indication of the fact that this problem is really not all that difficult\n",
    "that this difference equation can be written recursively\n",
    "to give:\n",
    "\\\\[\n",
    "        \\theta_{n+1} = \\theta_0 (1-k \\Delta t)^n\n",
    "\\\\]\n",
    "\n",
    "In a moment we will compute some values for this expression to see how\n",
    "accurate it is. First we consider whether we can improve the accuracy of the\n",
    "approximation by doing a bit more work."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f49fabc60f0>]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl0XOWd5vHvT6XNsiRvkrGNjTdJXsEYPMYEyNIkaUg6pruz2Rx6JmcYOEkPyeTASQ45nZA06enJcsLJApMJ3cmkk5xACOmk3QlpJpMmAzgYkDEGWwZblpEtr1osy5JtSaV6549bJdUmqWTXeuv5HHR0672vqn6+yE9dv/fW+5pzDhER8ZeSXBcgIiLpp3AXEfEhhbuIiA8p3EVEfEjhLiLiQwp3EREfUriLiPiQwl1ExIcU7iIiPlSaqxeuq6tzS5YsydXLi4gUpJ07d3Y55+on65ezcF+yZAnNzc25enkRkYJkZu2p9NOwjIiIDyncRUR8SOEuIuJDCncRER9SuIuI+NCk4W5mPzCzU2a2Z5z9ZmbfNrNWM3vNzK5Jf5lj2rsH+MmOdj7+4520dw9k8qVERApWKrdC/hB4GPjROPtvBRrDX9cB3w1/T6t/fK6Nn2xvg7Y2VnS101dRzbMNc/ir66en+6VERArepOHunHvWzJZM0OU24EfOW69vh5nNNLP5zrnjaaoRgLMXgly9/d/45q+/AcD/adzEkwfew19dP1FpIiLFKR1j7pcDR6Ied4TbEpjZ3WbWbGbNnZ2dU3qRmxrr2F+3ePRxU2c7LxzsJjgSuoiSRUT8LasXVJ1zjzrnNjjnNtTXT/rp2RjrFs3kxIIljJhX8hW9Jxju72d3x5lMlCoiUtDSEe5HgUVRjxeG29KqLFDCNU3zeWvWAgBKcDR0d/Dcgan9C0BEpBikI9y3Af8xfNfMJuBMusfbI97eVMebdVeMPl7R2c7zB7oy8VIiIgUtlVshHwNeAFaYWYeZ3WlmHzezj4e7PAW0Aa3APwB/nalib2yo40D0uHtXO7uO9HL2wnCmXlJEpCClcrfM1kn2O+C/pq2iCSytm07nkkb4o/d4RWc7IyHHCwe7ee+aedkoQUSkIBTUJ1TNjFkb148+buw6DMDzrRqaERGJVlDhDrDyxvUMBrx/cFx+tpOawQGNu4uIxCm4cH/byvm0zV44+rix8zBtXQN0nD6Xw6pERPJLwYX7rOnlnLqiYfTxii5vURKdvYuIjCm4cAewtWtHt5vC4f6cxt1FREYVZLjXX3/t6HYk3Le3djEScrkqSUQkrxRkuC9719ikk02d3h0zveeG2XtMUxGIiECBhntFw3IGyysBqD/Xy5yBXgCe07i7iAhQoOFOSQl9yxpHHzZF7ndXuIuIAIUa7kD5uqtGtyPj7jvbT3NuKJirkkRE8kbBhnvthqtHtyPhPjQS4sVDPbkqSUQkbxRsuNuVV45uRy6qgoZmRESggMOdqHvdV3S1g/Nug1S4i4gUcrgvWEBoxgwAagcHmHe2G4A3T57lZN+FXFYmIpJzhRvuZpTEn72H6exdRIpd4YY7xAzNNEaHu6YiEJEi55twXxF1UfW5A104p6kIRKR4+Sbcm6LO3Lv6B3njxNlcVCQikhcKO9zXrBndbOo6jLnQ6GONu4tIMSvscK+vh7lzAZgWHGRR78nRXZoCWESKWWGHO8QNzYyNu790qJsLwyO5qEhEJOd8Fu5j4+4XhkPsbD+di4pERHLOV+G+orM9ZpemABaRYlX44R5zUTU23J9v7cx2NSIiecFX4b6sp4PSkbEpf/cc7aO7fzAXVYmI5FThh/uMGbBoEQAVI0EWnz4es3v7we5cVCUiklOFH+6QOENklOcPaGhGRIqP/8K9Mz7cNRWBiBQff4R71Lj7iu7YcD925gIHOweyXZGISE75I9yjztyv7O1I2K2hGREpNimFu5ndYmZvmlmrmd2fZP8VZvaMme0ys9fM7H3pL3UCq1aBGQDzOzuoCA7F7NYUwCJSbCYNdzMLAI8AtwKrga1mtjqu2+eBJ5xz64EtwP9Md6ETqqqC5csBKAmFWN4de/b+wsFuhkdCyX5SRMSXUjlz3wi0OufanHNDwOPAbXF9HFAb3p4BHEtfiSmKGppZd+ZIzK6BoRF2He7NdkUiIjmTSrhfDkSnZUe4LdqXgDvMrAN4CvhkWqqbiqiLqu8YOpWwW+PuIlJM0nVBdSvwQ+fcQuB9wI/NLOG5zexuM2s2s+bOzjSHbdSZ+9okF1U1BbCIFJNUwv0osCjq8cJwW7Q7gScAnHMvAJVAXfwTOecedc5tcM5tqK+vv7iKxxMV7vOOtCbs3n2klzPnhtP7miIieSqVcH8ZaDSzpWZWjnfBdFtcn8PAzQBmtgov3LM7DtLUBKWlAJQebmfF9NjdIQcvtOnsXUSKw6Th7pwLAvcATwP78O6K2WtmD5rZ5nC3+4C7zGw38BjwMZftj4WWl3sBH3ZbeeIFVE0BLCLFojSVTs65p/AulEa3PRC13QLckN7SLsLatdDSAsBNF07wNWbG7Fa4i0ix8McnVCNiVmU6QonF7j7cc47D3eeyXJSISPb5Ntwr3mzhqoUzE7o8pwU8RKQI+Dbc2bOHmxoTbtjheQ3NiEgR8Fe4L1sGlZXe9vHjvKsukNBle2sXIyFNASwi/uavcA8EvEnEwq4608H08tiA77sQ5LUOTUUgIv7mr3CHmKGZ0jf2sWnZnIQuGpoREb/zdbizZw83Jhl311QEIuJ3vg/3mxoTpzl4pf00/YPBLBYlIpJd/gv3qNkh2bOH5XVVzJ9RGdMlGHK82Nad5cJERLLHf+F+xRVQXe1t9/RgJ09yY0OSoRmNu4uIj/kv3M1ih2b27uWmpsShmec0v7uI+Jj/wh0Sxt1vWJ54x8zBzgGOnzmfxaJERLKnKMJ9TnUFaxbUJnTT0IyI+JU/wz3uoiqQ9JZI3e8uIn7lz3CPG3PHOd6e5JbI51u7CGkqAhHxIX+G+2WXwZzwOHt/Pxw+zLWLZ1FRGvvH7RkYouV4Xw4KFBHJLH+Ge/wdM3v2UFkWYOPS2Qldn9enVUXEh/wZ7pAQ7kDSoRndEikifuTfcE/xourLb53mwvBItqoSEckK/4Z7/EVVYOW8GuqqK2K6DQVDvHSoJ5uViYhknH/DPfrMvaUFRkYws+SrM2ncXUR8xr/hPns2LFjgbQ8OwsGDAEnnmXl2v8bdRcRf/BvukPSiarJx9zdOnKXz7GC2qhIRyTh/h3uSi6qX1VbSdFl1QtftGpoRER/xd7gnuagKJF3A41ndEikiPlI84R4+c4fx55lxTlMRiIg/+DvcV68e296/37uwCly3dDblgdg/+qmzgxw41Z/N6kREMsbf4V5dDUuXetvBoBfwQFV5KdcunpXQXVMAi4hf+DvcIelFVUg+NKOpCETEL/wf7uNeVE0M9xfbehgMaioCESl8xRXuUWfuaxbMYGZVWUzX88MjvNLem63KREQyJqVwN7NbzOxNM2s1s/vH6fMRM2sxs71m9tP0lnkJxgn3QIlxQ5JPq2poRkT8YNJwN7MA8AhwK7Aa2Gpmq+P6NAKfA25wzq0BPp2BWi/OihUQCHjbbW0wMDC666Yk4a55ZkTED1I5c98ItDrn2pxzQ8DjwG1xfe4CHnHOnQZwzp1Kb5mXoLISGhq8bedg377RXckuqr5+9AynB4ayVZ2ISEakEu6XA0eiHneE26I1AU1mtt3MdpjZLcmeyMzuNrNmM2vu7Mzi8Mc4F1UXzqpiWd30mK7OwR8PdmerMhGRjEjXBdVSoBF4J7AV+AczmxnfyTn3qHNug3NuQ3194hQAGTPOuDvolkgR8adUwv0osCjq8cJwW7QOYJtzbtg5dwjYjxf2+WGicE96UVVTEYhIYUsl3F8GGs1sqZmVA1uAbXF9foV31o6Z1eEN07Slsc5LM0G4b1o+h0CJxbQd7T3PW93nslGZiEhGTBruzrkgcA/wNLAPeMI5t9fMHjSzzeFuTwPdZtYCPAN8xjmXPwPXDQ1QXu5td3RA79i97LWVZaxflDCCxPMamhGRApbSmLtz7innXJNzbrlz7r+H2x5wzm0Lbzvn3L3OudXOuSudc49nsugpKy2FlSvHHre0xOxONu7+rOaZEZEC5v9PqEZMMDSTbCqCHQe7CY6EMl2ViEhGKNyBdQtnUlNRGtN2djDI7g5NRSAihUnhDpQGSrh++ZyEH3l2v4ZmRKQwFU+4jzP1b0SyoRlNRSAihap4wn3JEqiq8rY7O+FU7AwJNyZZV/XVI730XRjOQnEiIulVPOFeUhJ79h41DQHAkjlVLJw1LaZtJOTYoakIRKQAFU+4w4Tj7maWdGhGS++JSCFSuEe5sSFxaEbj7iJSiIor3CcYlgG4oWEOFjsTAYe6BjjSo6kIRKSwFFe4x5+5x00ONrOqnKsun5HwYzp7F5FCU1zhvmABzAzPI3PmDByNn9wy+VQEz2vcXUQKTHGFu9lFjbtvP9jFSEhTAItI4SiucIdJP8x0zeKZVJUHYtp6zw2z99iZTFcmIpI2xRfu4yy5F1FRGuC6pbMT2nVLpIgUkuIO9yRn7pD806paek9ECknxhXv87ZChxGl9357kourO9tOcGwpmsjIRkbQpvnCvr4fLLvO2z5+HQ4cSujTMreay2oqYtuERx4ttPdmoUETkkhVfuMOkF1XNLOldMxp3F5FCUZzhPslFVRhvCmCNu4tIYVC4j3NR9YaGxHDff7Kfk30XMlWViEjaKNzHCff6mgpWza9NaNfQjIgUguIM9+gx9zfegOHkC3IkHZrRLZEiUgCKM9xra2HRIm97eBgOHEjaLfm4ezchTUUgInmuOMMdUrqo+h+WzKa8NPYQdfUP8saJs5msTETkkincYdxx98qyABuXJE5FoLtmRCTfKdxh3HCH5FMA66KqiOQ7hTtMGO7Jxt1fOtTDheGRTFQlIpIWxRvuK1cyuqZeaytcSH7/+qp5tcyZXh7TNhgM0fzW6UxXKCJy0Yo33KuqYPlybzsU8m6JTKKkxJJ+oOk5jbuLSB4r3nCHSxqa0dJ7IpLPUgp3M7vFzN40s1Yzu3+Cfh80M2dmG9JXYgalHO6Jk4jtPdZHd/9gJqoSEblkk4a7mQWAR4BbgdXAVjNbnaRfDfDfgBfTXWTGpBju82ZU0jC3OqH9+VadvYtIfkrlzH0j0Oqca3PODQGPA7cl6fdl4KtA4cysNcnUv9FuTDLurqEZEclXqYT75cCRqMcd4bZRZnYNsMg595s01pZ5TU1QWuptt7fD2fE/efr2pmRTEXThnKYiEJH8c8kXVM2sBHgIuC+FvnebWbOZNXd25sHdJuXlsGLF2OOWlnG7Xrd0DmUBi2k7fuYCBzv7M1WdiMhFSyXcjwKLoh4vDLdF1ABrgT+Y2VvAJmBbsouqzrlHnXMbnHMb6usTL1LmRIrj7tMrSll/xayEdn1aVUTyUSrh/jLQaGZLzawc2AJsi+x0zp1xztU555Y455YAO4DNzrnmjFScbimGOyRfOFvj7iKSjyYNd+dcELgHeBrYBzzhnNtrZg+a2eZMF5hx0RdVx5kdMuLGJLdE7mjrZigYSndVIiKXpDSVTs65p4Cn4toeGKfvOy+9rCyawpn7lZfPYMa0Ms6cH1vcY2BohF2HT3PdsjmZqlBEZMqK+xOqAMuWQWWlt338OHR3j9s1UGK8bXliiH/v2TbdNSMieUXhHgjA6qjPZE0yNJPs06r//sYpfrD9rTQXJiJy8RTuMKWhmT9bN59ZVWUJ7V/57T5e7ziT7spERC6Kwh2mdFG1trKMr31oXUL78Ijjnsde4eyF5Itti4hkk8IdpnTmDvCe1ZfxsbctSWhv7z7H53+1R+PvIpJzCndIDPcUwvlz71vJmgW1Ce3/8uoxfr6zI53ViYhMmcIdYNEiqKnxtnt64MSJSX+kojTAw7dfw/TyQMK+L/7LXlpPjT9PjYhIpincwVtub4pDMwBL66bzd3+xNqH9/PAI9/x0l9ZZFZGcUbhHTOGiarS/WL+QD127MKH9jRNn+bvfjD8RmYhIJincIy7izD3ibzevYVn99IT2n+w4zG9fP36plYmITJnCPeISwn16RSkPb72G8tLEw/nZX7zGkZ5zl1qdiMiUKNwjosN9714ITW0ysNULavnC+1cltJ+9EORTj+9ieESTi4lI9ijcI+bOhbrwlL79/XD48JSf4o5Ni/nTNZcltO863MtDv9t/qRWKiKRM4R5hdtEXVceewvjaB9dx+cxpCfu++4eDPLs/D1afEpGioHCPdgnj7hEzqsr49tarCZRYwr57n3iVU2cLZ/1wESlcCvdoaQh3gGsXz+be9zQltHf1D3Hvz3YTCml6AhHJLIV7tDSFO8An3rGcGxuSLMvX2sV3/9/BS3puEZHJKNyjRY+579sHIxf/CdOSEuOhj66jrro8Yd9Dv9vPzvaei35uEZHJKNyjzZoFCxZ424ODcPDSzrDn1lTy0EeuTmgfCTk+9dirnDmn6YFFJDMU7vHSODQD8Pamej7+juUJ7Ud7z/PZX+zW9MAikhEK93hpDneA+97bxPorZia0P733JD/Z0Z6W1xARiaZwj5eBcC8LlPDtLeupqSxN2Pfl3+yj5VhfWl5HRCRC4R4vA+EOsGh2FV/74FUJ7UPBEPc89goDg8G0vZaIiMI93qqo+WEOHPAurKbJrVfO545NVyS0t3UO8MVtU/9ErIjIeBTu8aqrYelSbzsYhP3pnRPm8+9fzcp5NQntT+7s4Je7tDyfiKSHwj2ZDA3NAFSWBXj49vVMK0tcnu/zv9zDoa6BtL6eiBQnhXsyGQx3gIa5NfztbWsS2geGRrjnp68wGNTyfCJyaRTuycTP7Z4BH752IbddvSChfe+xPr7y2zcy8poiUjwU7slET0OQgTN38KYH/rs/X8viOVUJ+/739rf4XcvJjLyuiBQHhXsyK1ZAIDwm3tYGA5kZB6+pLOPhrddQFkicHvgzT+7mWO/5jLyuiPifwj2ZykpobPS2nfMmEcuQKxfO4P5bE5fn6z03zKcff5WglucTkYuQUrib2S1m9qaZtZrZ/Un232tmLWb2mpn93swWp7/ULMvwRdVo//mGJdy8cm5C+0tv9fDt3x/I6GuLiD9NGu5mFgAeAW4FVgNbzWx1XLddwAbn3FXAk8DX0l1o1l3ikntTYWZ8/cPrmFdbmbDvO8+08seDXRl9fRHxn1TO3DcCrc65NufcEPA4cFt0B+fcM865c+GHO4CF6S0zB7J45g4we3o539xyNfGr8zkHn378Vbr70/dJWRHxv1TC/XLgSNTjjnDbeO4Efptsh5ndbWbNZtbc2Znni0VnOdwBNi2bw6dubkxoP3V2kPt+ruX5RCR1ab2gamZ3ABuAryfb75x71Dm3wTm3ob6+Pp0vnX4NDVAeXkWpowN6e7Pysp/8k0auWzo7of0Pb3by/ecPZaUGESl8qYT7UWBR1OOF4bYYZvZu4G+Azc65wh9DKC2NnUQsw+PuEYES41tb1jOrqixh31f/7Q12H8nOm4yIFLZUwv1loNHMlppZObAF2BbdwczWA9/DC/ZT6S8zR7J4UTXavBmVfOMj6xLagyHHPY+9Qt8FLc8nIhObNNydc0HgHuBpYB/whHNur5k9aGabw92+DlQDPzezV81s2zhPV1hyMO4e8ScrL+O/3Lg0of1Iz3k+98+va3k+EZlQ4tJASTjnngKeimt7IGr73WmuKz/kMNwBPnvLSl481MPrR8/EtP/mtePc1FDHlo2Jc8OLiIA+oTqxHId7eWkJD9++nuqKxPfgL/3rXvafPJv1mkSkMCjcJ7J4MUyf7m13dsKp7F9OWDxnOn//l1cmtF8YDnHPT1/h/JCmBxaRRAr3iZSUwOqoD+Nm8aJqtM3rFvDRDYsS2vef7OfBX7fkoCIRyXcK98nkeGgm4oubV9Mwtzqh/bGXDvPr147loCIRyWcK98nkSbhXlZfy8O3rqShN/F/2uV+8zuHuc0l+SkSKlcJ9MnkS7gAr59XywAfi52yDs4NB7vj+i2zbfYwRTVEgIijcJxcf7jm+v/z2jVfwvivnJbQf7jnHpx7bxZ9+81mFvIgo3Cc1fz7MnOlt9/XB0YSZF7LKzPgff3kVC2dNS7q/9VS/Ql5EFO6TMsuroRmAGdPK+M7W9VSVB8bto5AXKW4K91TkWbgDrL9iFr/4xNu4ftmcCfsp5EWKk8I9FXkY7gCr5tfy2N2bePzuTQp5EYmhcE9FnoZ7xKZlcxTyIhLDcjW74IYNG1xzc3NOXnvKurogsrjItGnQ3+99ejVP7Wjr5lv/9wAvtHVP2rdhbjWfurmR9185n0D8Gn8iknfMbKdzbsOk/RTuKZo3D06e9LZbW2H58tzWkwKFvIj/pBru+Xv6mW/yfGgmGQ3XiBQvhXuqCjDcIxTyIsVH4Z6qAg73iEjIP3bXJjYtS1yEO5pCXqSwKdxTlaP1VDPh+uVzePzu6xXyIj6mC6qp6uuDGTO87bIyGBjwvvvACwe7+dbv97OjrWfSvrrwKpJbulsmExYvhsOHve29e2MX8vABhbxI/tPdMpngg3H3iWi4RsQ/FO5TER3uBT7uPpGLDfknmo/Q1tlPSEEvknOluS6goERfVPXhmXu865fP4frl16c0XNN6qp/PPvkaANPKAqyYV8Oq+bWsnu99Xzm/luoK/bqJZIv+tk2Fz4dlxjOVkAc4PzzCq0d6efVIb0z74jlVrJpXy6r5tawKh/7CWdMw05i9SLrpgupUnD8P06d7qzGVlHhzzExLvmiGn03lwutkaipLw4Ffw+oFXvA3XVZDZdn4c9WLFDPdLZMpTU1w4IC3/corsH59buvJoXSGfLQSg2X11TFn+Kvn1zK3pkJn+VL0Ug13DctM1dq1Y+G+d29Rh3tkuOalQz089fpxWo71se9EH2cvBC/peUPOG8NvPdXPv+4ea589vdw7w58fGdqpZXl9NeWlui9AJJ7CfarWrIFf/tLbLqJx94lsXDqbjUu9u2qcc3ScPs++433sO37W+36ij/buc5f8Oj0DQ2xv7WZ769gsl2UBo2FuTUzoz59RyYxpZcyYVkZpQMEvxUnhPlVFelE1VWbGotlVLJpdxXvXzBtt7x8M8uaJPlrCgd9yrI83T5zl/PDIJb3e8IgLv5H08c8kLl5eXVE6GvTRXzOryqiN2h5tn1bOjGll1FSWUqIPZ0kBU7hPlcL9olRXlHLt4tlcu3jsvvmRkKO9e2DsDD/8dezMhbS9bv9gkP7BIEd7z0/p58ygtjLuTSHmTWCc9qpyppcHdG1Aci6lC6pmdgvwLSAA/KNz7itx+yuAHwHXAt3AR51zb030nAV7QXVoCKqrYXjYe9zXBzU1ua3JZ3rPDcUG/ok+9p/sZygYynVpKSktMWoqS6koDVBWapQFSigPlFBe6n0vC5RQVhppG9s/1lZCWcAoD3g/P9bmfXnPY0naxrbLIvtLSigxo6QE73v8tqE3ogKTtguqZhYAHgHeA3QAL5vZNudcS1S3O4HTzrkGM9sCfBX46MWVnufKy2HFirGz9pYWuO663NbkMzOrysMXa8fmnh8eCdHWOTAa+C3hMf2u/sEcVppcMOQ4fW4YGM51KSkx88I+YDa2XRK7HXkTCES9IUTavTeMqO1x3kAs/FqGEf5v9HHk/SXxcfTPJbYR7hvz3JHthNcy4t/Got/XLKY9tqeN8yD6GePfI2OfL/nP3PfeJubWVpIJqQzLbARanXNtAGb2OHAbEB3utwFfCm8/CTxsZuZydZ9lpq1ZMxbuX/gCNDbmtp4iUAasCH/9eVT7+aERTp8b4vS5IXoGhui7EGQoGGIwGGJ4JIRPfwPFB77zti3c9fZlzK3NzPOnEu6XA0eiHncA8aeqo32cc0EzOwPMAbrSUWTeWbsWfvYzb/t3v/O+JCemhb8W5LoQkSn6p2v/LKPPn9X7xMzsbjNrNrPmzs7ObL50en3gA4n/BhMRmaJMxkgqZ+5HgUVRjxeG25L16TCzUmAG3oXVGM65R4FHwbugejEF54V162DnTvjjH3NdieSh4RFvWCg4EmLEOUIhRzDkGAl/BeO+j4RC4ccQCoWS9h1JaPN+ZsQ5RkYgGHkccgSdY2QkvC/kcM7h8GbNCDkH3n+Ewu2SG53TZ2X0+VMJ95eBRjNbihfiW4Db4/psA/4T8ALwIeDffTveHrF+fVF/OlXGVxb+KgTOOZzDexOKbIe87ZCDUPS2S94eeQMJRf1s7HN6+5wj5o3GMfZGQ1Rb+P1ntC+RvpE+cftduFMkcdw4z5X4Z4/ajnqbi+8a2y/22CVrj28Y77nXAHXVFQl1pcuk4R4eQ78HeBrvVsgfOOf2mtmDQLNzbhvwfeDHZtYK9OC9AYhInrPIHTIJ95FIoUvpQ0zOuaeAp+LaHojavgB8OL2liYjIxdLEGyIiPqRwFxHxIYW7iIgPKdxFRHxI4S4i4kM5W2bPzDqB9ov88Tr8OrXBxdHxiKXjMUbHIpYfjsdi51z9ZJ1yFu6XwsyaU5nysljoeMTS8RijYxGrmI6HhmVERHxI4S4i4kOFGu6P5rqAPKPjEUvHY4yORayiOR4FOeYuIiITK9QzdxERmUBeh7uZ3WJmb5pZq5ndn2R/hZn9LLz/RTNbkv0qsyeF43GvmbWY2Wtm9nszW5yLOrNhsmMR1e+DZubMzNd3SKRyPMzsI+Hfj71m9tNs15hNKfxducLMnjGzXeG/L+/LRZ0Z5cLzLefbF970wgeBZUA5sBtYHdfnr4H/Fd7eAvws13Xn+Hi8C6gKb3/Cr8cjlWMR7lcDPAvsADbkuu4c/240AruAWeHHc3Ndd46Px6PAJ8Lbq4G3cl13ur/y+cx9dGFu59wQEFmYO9ptwD+Ft58Ebrb4Zcv9Y9Lj4Zx7xjl3LvxwB96qWX6Uyu8GwJeBrwIXsllcDqRyPO4CHnHOnQZwzp3Kco3ZlMrxcEBkaeoZwLEs1pcV+RzuyRbmvny8Ps65IBBZmNuPUjke0e4EfpvRinJn0mNhZtcAi5xzv8lmYTmSyu9GE9BkZtvNbIeZ3ZK16rIvlePxJeCmKlOGAAABl0lEQVQOM+vAW6vik9kpLXtSWqxDCouZ3QFsAN6R61pywcxKgIeAj+W4lHxSijc08068f9E9a2ZXOud6c1pV7mwFfuic+4aZXY+3ktxa51wo14WlSz6fuU9lYW4mWpjbJ1I5HpjZu4G/ATY75wazVFu2TXYsaoC1wB/M7C1gE7DNxxdVU/nd6AC2OeeGnXOHgP14Ye9HqRyPO4EnAJxzLwCVePPO+EY+h/vowtxmVo53wXRbXJ/Iwtzg/4W5Jz0eZrYe+B5esPt5THXCY+GcO+Ocq3POLXHOLcG7/rDZOdecm3IzLpW/K7/CO2vHzOrwhmnasllkFqVyPA4DNwOY2Sq8cO/MapUZlrfhHh5DjyzMvQ94woUX5jazzeFu3wfmhBfmvhcY95a4Qpfi8fg6UA383MxeNbP4X2hfSPFYFI0Uj8fTQLeZtQDPAJ9xzvnyX7kpHo/7gLvMbDfwGPAxv50Y6hOqIiI+lLdn7iIicvEU7iIiPqRwFxHxIYW7iIgPKdxFRHxI4S4i4kMKdxERH1K4i4j40P8HPv8p/gl53A8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "steps = 10\n",
    "theta_0 = 1.0\n",
    "const_k = 10.0\n",
    "delta_t = 1.0 / steps\n",
    "\n",
    "\n",
    "theta_values = np.zeros(steps)\n",
    "time_values  = np.zeros(steps)\n",
    "\n",
    "\n",
    "theta_values[0] = theta_0\n",
    "time_values[0] = 0.0\n",
    "\n",
    "for i in range(1, steps):\n",
    "    theta_values[i] = theta_values[i-1] * (1 - const_k * delta_t)\n",
    "    time_values[i] = time_values[i-1] + delta_t\n",
    "\n",
    "exact_theta_values = theta_0 * np.exp(-const_k * time_values)\n",
    "    \n",
    "plot(time_values, exact_theta_values, linewidth=5.0)\n",
    "plot(time_values, theta_values, linewidth=3.0, color=\"red\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Higher order expansion\n",
    "\n",
    "First we try fitting the local expansion for \\\\(\\theta\\\\) through an\n",
    "additional point.\t \n",
    "This time we assume that the variation in \\\\(\\theta(t)\\\\) is quadratic, i.e.\n",
    "$$\n",
    "    \\theta(t') = \\theta_{n-1} + \\beta t' + \\gamma {t'}^2\n",
    "$$\n",
    "\n",
    "The local time coordinate is $t' = t - (n-1)\\Delta t$, and when we differentiate\n",
    "\n",
    "$$       \n",
    "    \\frac{d \\theta}{dt} = \\beta + 2 \\gamma t'\n",
    "$$  \n",
    "\n",
    "To solve for \\\\(\\beta\\\\) and \\\\(\\gamma\\\\) we fit the curve through the sample points:\n",
    "$$\n",
    "    \\begin{split}\n",
    "            \\theta_n &= \\theta_{n-1} + \\beta \\Delta t + \\gamma (\\Delta t)^2 \\\\\n",
    "            \\theta_{n+1} &= \\theta_{n-1} + 2 \\beta \\Delta t + 4 \\gamma (\\Delta t)^2\n",
    "     \\end{split}\n",
    "$$\n",
    "\n",
    "Which solve to give\n",
    "$$\n",
    "     \\begin{split}\n",
    "     \\beta &= \\left( 4 \\theta_n - \\theta_{n+1} - 3\\theta_{n-1} \\right) \\frac{1}{2\\Delta t} \\\\\n",
    "     \\gamma &= \\left( \\theta_{n+1} + \\theta_{n-1} -2 \\theta_n \\right) \\frac{1}{2\\Delta t^2} \n",
    "     \\end{split}\n",
    "$$\n",
    "\n",
    "We can subsitute this back into the equation above and then into the original differential equation and we obtain the following\n",
    "\n",
    "$$\t \t\n",
    "     \\left. \\frac{d\\theta}{dt} \\right|_{t=n\\Delta t} = \\beta + 2\\gamma \\Delta t =\n",
    "     \\frac{1}{2\\Delta t} \\left( \\theta_{n+1} - \\theta_{n-1} \\right)  = -k \\theta_n \n",
    "$$\n",
    "\n",
    "The difference approximation to the derivative turns out to be the average of the expressions for the previous derivative and the new derivative. We have now included information about the current timestep and the previous timestep in our expression for the value of \\\\(\\theta\\\\) at the forthcoming timestep:\t\n",
    "\n",
    "$$\n",
    "     \\theta_{n+1} = \\theta_{n-1} -2k \\theta_n \\Delta t\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f49fa978c18>]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XmYXGWZ9/HvU1uv6WydfesQEkIQCRiWAURE0IBKYEDEQUVgYBS3d2RcGZXxHRTFwX2L4AAK4oiiEdEhLLIHCUsCIUAaQvaEztbppLurq6ue+eOp6nOq6lTSla5equv3ua6+6qm7Tp+ll/vcdZ+ljLUWEREZ/kKDvQIiIjIwlPBFRCqEEr6ISIVQwhcRqRBK+CIiFUIJX0SkQijhi4hUCCV8EZEKoYQvIlIhIoO9An6NjY22qalpsFdDRKSsPP3009utteMONN2QSvhNTU0sX758sFdDRKSsGGPW9WY6tXRERCqEEr6ISIVQwhcRqRBK+CIiFUIJX0SkQijhi4hUCCV8EZHB0NwMXV0DukglfBGRgXb99TB7NhxxBHR3D9hilfBFRAbab3/rHpub4bnnBmyxSvgiIgMtHvfGA9jWUcIXERloiYQ3VktHRGQY8yd5JXwRkWGs3Ct8Y0zYGPOsMebu9POZxpgnjTHNxpjfGGNipVqWiEhZGwYV/qeB1b7n3wS+Y609FNgFXFbCZYmIlK9yTvjGmKnAu4Eb088NcBpwZ3qSW4BzSrEsEZGy52/pJJMDtthSVfjfBT4HpNLPxwK7rbWZXddGYEqJliUiUt7KtcI3xrwHeMNa+/RBfv8VxpjlxpjlLS0tfV0dEZGhr4wP2p4EnG2MeR24A9fK+R4wyhiT+QjFqcCmoG+21i621i6w1i4YN+6AH8koIlL+yrXCt9Z+0Vo71VrbBFwIPGCtvQh4EDg/PdnFwB/7uiwRkWGhXBP+fnwe+IwxphnX07+pH5clIlIeUin3lTGACT9y4El6z1r7N+Bv6fFrwHGlnL+ISNnLTfDDpMIXEZFcSvgiIsOQtfDSS9ktHP8ZOqCELyIyLFx+ORx+OJx5phdThS8iMgxlPujk3nuhtdWNVeGLiAxDnZ3eOPNBJ6rwRUSGIX81nxkr4YuIDDPJpDtom5FJ+GrpiIgMM7mfVasKX0RkmMqt5FXhi4gMU6rwRUQqRKFKXglfRGSYUUtHRKRCqKUjIlIhVOGLiFSIQglfFb6ISBlrb4dHHnEXW2WopSMiMsxYCyeeCKecApde6sXV0hERGWa2boUVK9z4nnu8eG6Fr9MyRUTKXNAdMUEVvojIsONP8r1J+KrwRUTKVDzujQslf1DCFxEpe/7Enkp5Z+qopSMiMsz4K3wonNhV4YuIlLnchJ+p+Au1dFThi4iUqdzEnnmuu2WKiAwzhSp8tXRERIaZQhW+WjoiImXKWrjySjjqKHffnAxV+CIiw8yKFfCTn8DKlXDddV5cFb6IyDDT0hI8VoUvIjLM+O+ZU+jqWv9zJXwRkTJVKOH39jz8TGJXS0dEZIjzJ/xCyR9U4YuIlL1StXTKucI3xkwzxjxojHnRGLPKGPPpdHyMMWapMWZN+nF031dXRGSQ9LWlM0wq/G7gKmvtPOAE4OPGmHnAF4D7rbWzgfvTz0VEypMO2oK1dou19pn0uA1YDUwBFgG3pCe7BTinr8sSEel31sLtt8MNN8C+fV682B7+EDwPP1LKmRljmoCjgSeBCdbaLemXtgITCnzPFcAVANOnTy/l6oiIFO/xx+Gii9w4lYJ/+zc39if5VMol6kik94m9UIWfSrmvUP8fUi3ZEowx9cDvgP9nrd3jf81aawEb9H3W2sXW2gXW2gXjxo0r1eqIiByczAeS5479CR+8yr7Ys3Ry4+B9aEo/K0nCN8ZEccn+Nmvt79PhbcaYSenXJwFvlGJZIiL9qr3dG3d0eONCCb+3t1YodHvkQrF+UIqzdAxwE7DaWnuD76UlwMXp8cXAH/u6LBGRfudP+IX69v7nfT0Pv1CsH5Sih38S8CHgeWPMc+nYl4DrgP8xxlwGrAMuKMGyRET6l7+q702FX4qWTrkkfGvto4Ap8PI7+jp/EZEB1V8tnSFQ4etKWxERv2ITfhlV+Er4IlKZNm6Ek06C006DXbu8eG8TfuZ5Xy+8KhTrB0r4IlKZfvELd879gw/Cr3/txUvVw1dLR0RkiNi61Rtv2+aNe3uWzoF6+IWuqFVLR0RkgO3d6439t1AoVQ+/mAq/nC68EhEpO31N+Afbw1eFLyIywPwJv1CS7+72knGpztJRD19EpJ9s2waf+QwsXpwd91f1hSp88HYAw+CgbUnvlikiMuR8/evw/e+78bHHwtFHu3FvWjrgEv6IEX0/aKuWjohIP3vxRW+8erU3LibhWxvcw7e2uIRvrVo6IiL9prXVG7e1eeOghG9tfsLv7HQJOZXKjsfjwdV6JrEHvZZKqcIXEek3e3wfz3GghJ9I5Cf2jo786h5cws/t34Or8Asl8ERCFb6ISJ9t2wbXXw9//3t23F/hZ5J/MpldyWfGudU97D/h57ZzwMWC4jDoCV8HbUVkePjEJ+DOO90B1g0bYORIFw+q8HMTe6bCLybhd3YWrvCD2jaZ14IuslKFLyJShExl39YGa9a4cSKRncQzCd/fzoGDS/j7q/ALJXz/Of5+SvgiIkXYudMb797tHvfsyZ6mUMKPx13lHZSQD6aHX6ilo4QvIlKEjg64777saryrKzuJF0r4mee5CR9clR9U4Xd2Fp/wC1X4QfMHJXwRkUCLFsEZZ8DChV7Mfz978BK+/4AtFK7wwSXjYnv4xbZ0guYDSvgiInm6umDpUjd+5BEvgfvbOVB8SwcKV/hq6YiIDIKWluznO3a4x1JU+Pv2Fd/DL7bCV0tHRCTA174GRx0Ff/qTF3vjjexpMgk/t8LP7AAKJXz/rRQySlnh6ywdEZFe2rQJrrkGVq6EL3zBi+cm/O3b3WNvWzoHc9C22B5+Mlm4V68KX0QqmrX5lfJrr7k4wNq13ji3pXOghJ9b4Wf66/3Zw8/MK4gqfBGpWPE4HHccjB0LS5Z48U2bvHFHh9eK6W1Lp1CFD25exZyls7/TMgsdnA2aPyjhi0gFu/deWL7cVcSZe9YDbN6cPV3mA8cLtXR6e9AWCif8gzloW6jCL9S6UUtHRIa9ri53r5sPf9irysG1bjLWrfPG/gofDpzwS1HhH0wPv9iWziCfh6+bp4lI//v1r+FHP3LjmTPhP/7Djdeu9abZsMH16o0pXOEXOi2ztz18cDuBUvXw1dIRkYq1YQOceiqcc0528nz6aW/83HPe+PXXvXE87iX0vlb47e0uGRfb0unvg7Zq6YhI2Ukm4aGHYOPG7Pj117v4H/8It97qxV96yRs3N3tjf4UPbocBxffwcxM+uGRfbEunFBdeZeYVRBW+iAxpQVXs17/uKvl587xkDPDUU9742We98csve+NXX3WfKmVtdoUPXluntxV+oSttwbV1iqnwiz1LJ5EonMALtXRU4YvIkBD0wRw33AD19fC+93nnwgP88pfusa0N/vIX7/tXrvSmeeEF97hvH6xf78XjcZfQd+/Or8DXr3fx3AS7dWtwy2X7drfzKJTw+/OgbWZeQVThi8igSySyk2/GVVe5xP6Vr2RPe801LhndeadXsb/xhvcBI+B98MiaNdkJc9Uqt5N45ZX85TU357dzwFX4udU9uISfe8AW3M5j8+b8z6EFtxPoz4O2UDjhD9GDtmV/ls41S1YRCRlG1kQZWRtlZE2UhpooDdURGqqjjKiOUl8doTYaJhQyg726IgPjkUfcKY8XXQSR9L95Rwccfzw8/zx8+cvuXjXgprvhBje+9lr41KegsRGeeCI7oT34IBxzDDz+ePayMm0c/8FYcMl206bsdk5Gc7O72CrX/hJ+bjsnw7/z8du8OfhdS6l6+BD8DiIzryDDPeEbYxYC3wPCwI3W2utKNW9rLb984nU+f/9NtFbXs626ntb0V1tVHW1VteyN1dJWVUt7VQ11VVHqqsLUVUWoi0WojYXdV5XbIdTEwtREw1RF3WN1NER1+rEqEqYq4h5jkZD7CrvHqkiIaDhENGyIpuPRcIiwdjBSrK4uiEbdqYkZqZSrsCdMgLe9zYtv2QIf/KD7nttug+nTXfzhh+G001yye+YZ+N73XPz2212yB/iv/4LPftZ9/utvfpO9rLvvho98BP761+x1e/BB907gscey4ytXuoSYm/DBtXX8B2wzXn01OFlu2JB/wBZcsvcfK/ALegcB2ef1+7W1lebmaaCE72eMCQM/As4ANgJPGWOWWGtfLMX898a7qe5s54qn7jrgtEkTYm+shn2xGtqj1bTHqmmPVrM3VkN7rIaOSBWd0RjxcIyOaBW70693RmIkwlG60l8d0So6IzHikRjJUIhEKEIiHOmJxcMxUqEQSRMiFQ4RCUeIhA2RkCEaDqXH7jEccvFwKEQ4hHs0EAmFCIUgHDKEjJsubAyh9GM4ZDCGnrgxhpCBkDGEQmQ/N27akDEYIBRyj6YnDgY3Nibzmi+WeW5MT5yeabKn9f3e3SNe3jL+uC+Gb1rv+31j3yuZZe5vmiwmcJg9iW9GB9o9h7rihLriJEc0ZMWrdrRQt+F1dh95NDbi/Us1rFnNuL8/xua3L6Rj8tSe+IRHH2DGH+5g69veyYb3vq8nPuu2Gznie9fSetg8Hv/Br0iMGg3AUf/5eWb+/jYA/v7Nn7L5jPeAtZzw6Y8w8dEHANj2gYt54oe/Amt566euYmy6sk3+5Kcsfc/FxBvH87brv8OozMLa23n22z9j/bkf4NT//iUjfduz5dY7+Pv8M3jbXX/ypgcSf3uIvyxfx8n3PsgY/w8gkeDh3y5l7sPLGJ/zM3vh3scZ9eJKpubENy9/ns7G8RySE+94dS2vL1/N4TlxkklW3b+MI3LjQPNjz3JoQHzds6uZERDf1LyBKQEJtnX7bqr27KU6J57q6GDXjj0EvB9h19btjA6It7bsyvqZZuzZuYeGgPj6N/bw7HPunc0Jh4xlQkPuWpRGf1f4xwHN1trXAIwxdwCLgJIk/NaOBKM6C/TKcoRtipHxfYyMF9jz9pN4OEI8HKMrEsXiMlbKhIiHo8QjMbrCUawv6SRCEeKRaE/cWEvIpugORXp2SNZAJJUilEqSCoXpjMTojMSwxhBJdhNNJUkZQ0e0io5IFalQiFh3gliym5BN9ezwusJRqhJxarrjRFJJtwOsqqUzEqMmEac20Ul1Is6+WA17qurYF6uhpjtOQ+c+6rvaaY/WsKtmBLtrRlCdiDOmYw+jOtrojMTYXjeK7XWjiSYTjNu3i3H7dpMMhdgyopGtIxoBmLynhcl7WginkmwYNZH1oybSGYnRtGsLTbs2U9/VztpRk1nTOI1dNQ0cumMjh7W8zrh9u3htzBSen3go60dN4rCW15m/5RUO2bmJtaMn8/dpb+L5ibOY27KOE9et4MitzWxpaOSBQ47lkZlHc8jOTSx8+TFOXL+SPVV13D33rdwz92TGtLfy/pVLOfPlx0iGQtx1xNu5ff6ZAFz21B/40DN/pi7Rye/edBrfPfmf2FnTwMeW3clHn/wdVckEz0+YxefO+jQvN87gY8vu5F8fvY2ITdH0g+v5zHuu4r5Dj+Ofn7qLcx/8BQBTHvgLj/7vk3zn5Iu46Lm/cM69PwZg9KoVjPnIP/HhC77Ge1c/zDl/vq3n72POV/+NT6yr4ZjNL3HOo/f3xCc8/jfu+sZNAJyzYnlPPJzo4rWvfIP7Zh3PopdXZf1tdt/03/xgUx2L1qzOio965G9c/9O/suilF7Li0X17+dV37uBdq1bk/Z3fd+ufmbciv8JfvfQJ5r3xWl7C3/38ajaPaMlL+LE3tvHsYyvzEz6w5r7HAxP+2ieeC0z4G557KTDhb1n9GlMC4ju372bsvvyEH0omeWPz9sCEv2NzcMLf6Uv4KQwh3AHv3TtaexJ+IhQmmnI75qdfbeFf73A/v5svObbfEj7W2n77As7HtXEyzz8E/LDQ9G95y1tsMZ7fuNse9anb7bWnXmJ/eML77K/mL7R/Ouxk+/CM+fbZSXNs85ipdlvdaLsvWmWtO3ykL30FfsVDkcB4Zzhq2yP5fz+d4ajdUj8mL94VCtsXxzUFzuvx6UcGxu899DibxATEjw9c9iMzjgpc9itjp9nnJs7Oi7dW1dl7Dz0ucNl/PPyUwPg9c04MjC+ddWxg/KGmowPjKyYeGrgNe6PV9uWx0wO/Z+WEWYHx1Y0zAuNrxkztGft/jq+NntQz7jahnvG6kRN6xglffGv9GNsVCgf+TWwcMS7rd58Z+38P/jyz1RffG63uGW+rG90zbovV9IyXzH2rnfH5u+2Mz99tH3xpW1F50Lrkurw3OdmkE3G/MMacDyy01v5z+vmHgOOttZ/wTXMFcAXA9OnT37KuUN8twPa9cZa+uI3WjkTW156OBHs6u2nrTNDW2c3ezm4SnXHqu9qpTXRS29VJXVcHtYlO6ro6qetqpyYRpyqZoLo7Tk0i3lPhVnV3EU12E0smqOruoqY7TnWii6pkgnAq6b2WTFCdiKer6CThVIqIDThzQCpKVyhCLJXfPtgTq6WhK79f3FI3inH7dufFN40Yx6S27T2Von/6mkSc+q7s3nBnJEZL3WimtW7Lm9crY6czZ0f+mTqF4mvGTmP2jg158eYxUzl058a8+KtjpjBrp2tP+KvbnTUNRFJJGtLvspMmRDj9P7Ju1ERm7HY9+m4T6vnf2TByQs82+CvijQ3jmLrHnbXj/xlvHtHI5DZ3MVY8HKEq6eJb68cwce/Onp9Ndbe7YGpb/RgmpOPt0SpqE/GedWitrmdsh+vR743V9PyM/b+jrHjtKMa1u3hbrIYRAXH/73177Uga291ZRPfMOZErz/0S4Cr8Uw/LbY7tnzHmaWvtggNN198tnU3ANN/zqelYD2vtYmAxwIIFC4ra+zTWV/GB46b3atruZIp98ST7urrZF+9mb7ybjq4k7V0u1plIuueJJJ2JFDsTSTYlksQTKTq7k3QmknR1p+hKpogn3GNXt/uKd6dIJN1XV3eKRMqSSKawKUss2U11d5xYMoGxAK5FE0t2U9XdRVW3d6WegfSOJUEs6T4xJ2UMFkMk1U11dxfV3V1YXFsoGQoTTiWp7o5T3d2FsZZEOEp3KIyx1u28uuOEUyni4WhPWymzo4umEnRGXNsnGQpTk+hkRFc71Yk4HdFq9sWqiYdj1CU6eto4+2I1tMVq2VdVQ21XJ6M72hjdsYeOaBU7axporRlBdaKLxvZdNO7bTSIUoaV+NC11owmnkkxq28GkNvePurlhHJsaxmONYdrubUzfvYXq7i7Wj5rI66MnszdWyyE7NzJ7xwZGt+9h7ZjJvDSuiW0jxjJ7+3qO3NrM1NZtvDZmKs9NnsOasdM5rOV1jt+wirktr7N+1AQemzGf5VPnMXv7ek579SmO2rqGnTUNLD30eJbOPoFprVs5Z9XfOGrrGpImxAOzFvCbN7+Luq52Llv+R9681V0V+uL4mfzgH97P7poRfP6hW5i/xR0o3F47km+dcjFPTzmc6/76A47d5LqVSRPi+ydeyG1Hn8l3//RtTl7ntUHuPuxkvnjmJ7nh7hs4o/nJnvgrY6dzwUXX8eUHbuS8Fx7oibdHq7j4gv/g/Svv5fKn/pD1d33NO/6FGbu38LmHb82K/2r+maxpnM43//qDrPjTk+fy0xPO5+e//8+s+N5YDV8946P8+o4v5f3vXH/Kh1l817V58VuPeTdfuf/neYXNo03zaejcx7j23Vk7qFfHTKUq2dXzM80k+7ZYDavHz+xJ+P75vTBhVk/CzyR7gFfHTutJ+P4d6oZRE3sSfibZA2wd0diT8Kt9/3M7akf2JPxMsgfojFbRGanqee7foe6pqutJ+P54e6wa0vtw/7zikVjP2L/seNiLD1Rx2N8J/ylgtjFmJi7RXwj8Uz8vM1AkHGJkbYiRtdEBW2YynfgTyRTdSUt3ytKdcuNkynueTHnPU+nHZMqSsv5H8mIpa0ml8MbWjW3OYzLl/uky01gLFvd6zFoivli7hXZrSe+biGHpstACbMu8G7Rggb3W0mYh854sM49OYJeFNT1x2/Nvv9vCarzlp2dHK/BCeh4eyypgVWbePdNbVpOej2/6CPBq+ss/mxrrzhi4FdfTToUj2FCIBiytwC1cSf2uHSQjUTrqG3p6rDdfdBFT164mkuji9TlHUWMM1cAti87kyeceY0zLZp456UwSdSM40lpuP+1Ymh9aQtOalTx56iLWzTmKk4Hfn/Bz9t51I/OXLeXZf3gX95/7z5wWCnHPW77PmJ98lbc8/le2TZrBr/79Z7xtzHgeW/B15l/7UWa97Hq6v7v83znirSfw0ilHs+Xq55m08VUAXpx/EpELL2BboosdLy5l7PYtAMSrqnntsk+SqKmjddmvGbl7e8/P4oVzP0Td8aez574f07DHux3B6mNPZfx73knrX2/Imr51VCPV5/8j2x+9mcaW7NMlY+96J1vXPMjUddlnycTecgy7Orcx7sXlWfHu2XMIdbbD1uas+J6JU6k+pAleeYJcyaPmwyuP58VTcw+Htc/kxe3MmbDhhbw4U6fClvyzeSKNY0nu2EA4mf0uLFVdQ6SuFgJOs4+MGgU7808dDY9sgPROK+xL4OERdZA+mce/c4rUe/OfXB/hvUdNBmD8iH7q3wMH7Pn09Qs4C3gF93949f6mLbaHLzIsNDdb29GRHWtrs/bb37Z2yZLs+KpV1r75zdYef7y1GzZ48Tvv9PraX/mKF7/uOi8+fry1nZ0uftVV2b3wP/3Jxf/lX7LjF1/s4pdemh0fO9baVMrayy/P76svX27tJz+ZH//Wt6y9+ur8+Nlnu23NjY8bZ+0tt+THwdrFi4PjX/tacPxTnwqOv/e91o4cmR9vanI/56DvOf304PhJJwXHFywIjvvn/6539elPiF728Pv9Sltr7T3W2jnW2lnW2vz3hSKVbtYsqM6p6urr3Tnv731vdnzePFixApYtc1VrxnnnwdKl7jbE11zjxa+80n2ilDHw3e9CVbpNcckl3jSjR8M73+nGixZlL2/hQvf49rdnx0880c3z2GOz4+EwHHEEvOlN+ds5d67b1lxNTTBtWn58yhSYODE/DnB40Hk8wCG55/345hWkvh7q6vLjtbVQUxP8PSODTrjEXdMQpLY2OO7/nQ+H8/BFZACdfnp+bMQId2Wstd4Vt+CS8rXXuouuvvxliKX7yaed5hJaa6ub/owzXDw34Z90knvMTfiHH+4SWVDCP+wwGDUqPz5zpnfRmN/kycEJv64ueHoonPBHj3aJN/fCqmITfiQSPH1mXkEKJXz//HUvHREpiXA4O9lnfOlL7t3C+ed7saoquPlmeOtb4Wc/825/MGWKS9gZJ5/sHo84IjtxHX20F/eLRl1iPzTgjPliK/zx44NvywDuXU/QtjY0BFfg9fXBCbmmJjjhV1d775JyFarwC+0gBqHCV8IXkWznnONuz3Dppdnxb33LJdSLL3YtHXCJPJPkAebPd48jR2a3nGbNctNOnJifYJuaXDw3UU+e7BJ7OJwdHz/ezSO3DQYwZkzwu4iRI13Sz1Vshb+/hK8KX0SGjbPPdve6ufnm7HtbvP/97jEazT7m4G/rzJ3rHo3Jr/KbmlxSz+2zT5ni4uNzzkkfN87Np7ExOx6LueQalPD3V+GXKuGXQQ9fCV9E+ubjH4f773c3UZs924sfeaQ3ziR8yD5wO2qUl6Bz2zqT3WmKeW2dzA4gt60zZozbERSq8EuV8IPeWWTmFUQJX0SGjXDYHez1J3WAK65wd/icOBEuv9yL+yv8mTO9cW7Cz1T8hRJ+boU/On1Xm2Ir/KCEXKhlVMoKfxBaOjpLR0T6x6GHep956+/P+xN+U5M3zj3zplCFP26ce8xN+GPS9/AcHXA7s2Ir/Jqa4MR+MD18HbQVkYoQieQfjD3vPFf5R6Nw2WVe3F/hR6NeQp80Kfv7C1X4mYSfW+Eb45Jxfx60jUQKt3pU4YtIxRo71n3MYltbdh/en/AnTYJQuh4tpocP+Ql/xAg3r2J7+NGA27AU6uFXVXnXMgTNK4guvBKRihCL5SfsI4901bi12Qd8+9rDz1wZW2zCDzqfv1CFH4sVTvhDqKWjhC8iQ8PMmXDjjfDQQ/DFL3rxYnv4uQk/08op9qBt7vn/UDjhV1UFvyPIzCuIWjoiUtEuvTT/gq9CCb+3LZ2DqfBraryWkt/BVPhDqKWjg7YiMrRNnuxVw9OmeYm12Aq/VAdtS9XD15W2IiI56urgRz9y9/f54Q+9eF97+Ma4pFuKs3TKpMJXS0dEhr5LLsm+pTP0vqVTqIdfX++SfqGEn/1pPM7+evhlcNBWFb6IlKfcq2ELXXhVqMLPJOJCB22LqfDL5LRMJXwRKU/GeHfnHD3auxVDb3v4mStjCx20LaaHv7+WTqEPUvHPJ5kMfkdRYmrpiEj5uukmdyrnokVeYq2udsm3K/2B4YUq/P0l/Nra4Kq72Ao/Gi18umYs5s4ESqU//zaVCj4VtISU8EWkfM2bBzfckB3L3DHzjTfc84NN+IlEfrzYg7axWPAFXODddiKzY+ru7veEr5aOiAw//rZOppWTe3HUgRJ+KXr4+6vwc+8zNAB9fCV8ERl+/Adu/R867q/yMwm/2IO2xfTw95fwo1ElfBGRPvvHf3SPEybACSd4cf+B20zCD4ezq3ZjXPIuRYUfi7n5BbVqVOGLiJTAZz/rPoFrzZrsCj6owofstk5trUvS0Wh+ot5fDz+oV5+p7oOqfFX4IiIlYIy742bugdreJvyM3PZNdbU7syY3uVdVee8M/PaX8FXhi4j0o2ITfm5bJ7MDyN0RZKr+3MSe2QEo4YuIDLBCCd+f5HuT8HPbOpnEXqjCL9TuUcIXEekn/o9LnDDBG/srfH+S723CzzzPTfhDrMLXhVciUjmuvBIefdR9ataiRV682JZOsRX+EDloq4QvIpVj7lx45pn8eF8TfqEKXwdtRUSGmGJ7+LkHbQtV+EO8jBT0AAAJc0lEQVSspaOELyJSzGmZ0PcKPxx2p3Eq4YuIDLDeHrTNJPa+9vAziV4JX0RkgPWmh+9vwRR7lk7uaZlBp2sq4YuIDIDMh6dA9qmb/oTvb+8UuvBqOFf4xpjrjTEvGWNWGmPuMsaM8r32RWNMszHmZWPMu/q+qiIi/eQDH4DzzoN3vxsuu8yLF0r4ZdrS6etpmUuBL1pru40x3wS+CHzeGDMPuBA4ApgM3GeMmWOtTfZxeSIipTdiBNx5Z368twm/2AuvyrGlY62911qbWctlwNT0eBFwh7U2bq1dCzQDx/VlWSIiA65QG6dMK/xS9vAvBf6SHk8BNvhe25iOiYiUj7728IdYhX/Alo4x5j5gYsBLV1tr/5ie5mqgG7it2BUwxlwBXAEwffr0Yr9dRKT/VFoP31p7+v5eN8Z8BHgP8A5rrU2HNwHTfJNNTceC5r8YWAywYMECGzSNiMigKLaHX6iSHyIJv69n6SwEPgecba1t9720BLjQGFNljJkJzAb+3pdliYgMuFIdtB0i5+H39SydHwJVwFJjDMAya+1HrbWrjDH/A7yIa/V8XGfoiEjZKXSbhd7eS2eIVfh9SvjW2kP389q1wLV9mb+IyKA65RSX9Nvb3Tn6GX29l06ZVvgiIsPXuHGwdi1s3AhHH+3Fe3vQttBZOuVY4YuIDHvjx7svv75W+OV40FZEpCL19X745XilrYhIRfJX+JEIhNKpVBW+iMgw40/4/iRfKOEPkdMylfBFRIrlT/iFkr//uSp8EZEyVegiLLV0RESGmWJbOjpoKyJSptTSERGpEDNmePfZmTvXiw/xCl8XXomIFGv0aPj97+H+++HKK734EK/wlfBFRA7GwoXuy683lbz/uVo6IiJlaoi3dJTwRURKZYi3dJTwRURKRRW+iEiF0IVXIiIVQi0dEZEKoZaOiEiF6O2HmKvCFxEpc6rwRUQqRLEHbcNhL6aELyJSRvwJPxLxEroO2oqIDDOxGJx7rhtfcIEXHyItHd1LR0SklH73O3jlFZgzx4sNkQpfCV9EpJSMgcMOy44NkQpfLR0Rkf6m0zJFRCrEEGnpKOGLiPS33rR0ksl+Xw0lfBGR/qYKX0SkQuigrYhIhfBfUQuq8EVEhi1jsqt8JXwRkWHMn9zV0hERGcaGS4VvjLnKGGONMY3p58YY831jTLMxZqUx5phSLEdEpGz5E365VvjGmGnAO4H1vvCZwOz01xXAT/q6HBGRsjZMKvzvAJ8DrC+2CLjVOsuAUcaYSSVYlohIeSr3hG+MWQRsstauyHlpCrDB93xjOhY0jyuMMcuNMctbWlr6sjoiIkPXvHnusaEBJk9246F2t0xjzH3AxICXrga+hGvnHDRr7WJgMcCCBQvsASYXESlPP/4x/PznsHAh1NW5WCjkTtm01n2lUi7WTw6Y8K21pwfFjTFHAjOBFcYYgKnAM8aY44BNwDTf5FPTMRGRynTIIfCNb+THIxFIJNy4uzv/YxJL6KB3Jdba56214621TdbaJlzb5hhr7VZgCfDh9Nk6JwCt1totpVllEZFhZADbOv31ASj3AGcBzUA7cEk/LUdEpLyVY8JPV/mZsQU+Xqp5i4gMWwOY8HWlrYjIYFLCFxGpEEr4IiIVQglfRKRCKOGLiFQIJXwRkQqhhC8iUiGU8EVEKoQSvohIhVDCFxGpEEr4IiIVQglfRKRCKOGLiFQIJXwRkQpRjrdHFhGRg/DRj8JZZ7nEf8QR/booJXwRkcF09tkDtii1dEREKoQSvohIhVDCFxGpEEr4IiIVQglfRKRCKOGLiFQIJXwRkQphrLWDvQ49jDEtwLqD/PZGYHsJV6dcVOJ2V+I2Q2VudyVuMxS/3TOsteMONNGQSvh9YYxZbq1dMNjrMdAqcbsrcZuhMre7ErcZ+m+71dIREakQSvgiIhViOCX8xYO9AoOkEre7ErcZKnO7K3GboZ+2e9j08EVEZP+GU4UvIiL7UXYJ3xiz0BjzsjGm2RjzhYDXq4wxv0m//qQxpmng17L0erHdnzHGvGiMWWmMud8YM2Mw1rOUDrTNvunOM8ZYY8ywOJujN9ttjLkg/fteZYy5faDXsdR68fc93RjzoDHm2fTf+FmDsZ6lZIz5hTHmDWPMCwVeN8aY76d/JiuNMcf0eaHW2rL5AsLAq8AhQAxYAczLmeZK4Kfp8YXAbwZ7vQdou98O1KbHHyv37e7NNqenGwE8DCwDFgz2eg/Q73o28CwwOv18/GCv9wBs82LgY+nxPOD1wV7vEmz3KcAxwAsFXj8L+AtggBOAJ/u6zHKr8I8Dmq21r1lru4A7gEU50ywCbkmP7wTeYYwxA7iO/eGA222tfdBa255+ugyYOsDrWGq9+V0D/H/gm0DnQK5cP+rNdl8O/MhauwvAWvvGAK9jqfVmmy3QkB6PBDYP4Pr1C2vtw8DO/UyyCLjVOsuAUcaYSX1ZZrkl/CnABt/zjelY4DTW2m6gFRg7IGvXf3qz3X6X4SqDcnbAbU6/xZ1mrf3zQK5YP+vN73oOMMcY85gxZpkxZuGArV3/6M02XwN80BizEbgH+OTArNqgKvb//oD0EYfDjDHmg8AC4G2DvS79yRgTAm4APjLIqzIYIri2zqm4d3IPG2OOtNbuHtS16l8fAG621v6XMeYfgF8aY95krU0N9oqVk3Kr8DcB03zPp6ZjgdMYYyK4t387BmTt+k9vthtjzOnA1cDZ1tr4AK1bfznQNo8A3gT8zRjzOq7HuWQYHLjtze96I7DEWpuw1q4FXsHtAMpVb7b5MuB/AKy1TwDVuPvNDGe9+r8vRrkl/KeA2caYmcaYGO6g7JKcaZYAF6fH5wMP2PQRkDJ2wO02xhwN/AyX7Mu9pwsH2GZrbau1ttFa22StbcIdtzjbWrt8cFa3ZHrzN/4HXHWPMaYR1+J5bSBXssR6s83rgXcAGGMOxyX8lgFdy4G3BPhw+mydE4BWa+2WvsywrFo61tpuY8wngP/FHdn/hbV2lTHma8Bya+0S4Cbc271m3AGRCwdvjUujl9t9PVAP/DZ9jHq9tfbsQVvpPurlNg87vdzu/wXeaYx5EUgCn7XWlu272F5u81XAz40x/4o7gPuRci/kjDG/xu24G9PHJr4KRAGstT/FHas4C2gG2oFL+rzMMv+ZiYhIL5VbS0dERA6SEr6ISIVQwhcRqRBK+CIiFUIJX0SkQijhi4hUCCV8EZEKoYQvIlIh/g/zfssfnDV0NQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "steps = 100\n",
    "theta_0 = 1.0\n",
    "const_k = 10.0\n",
    "delta_t = 1.0 / steps\n",
    "\n",
    "theta_values = np.zeros(steps)\n",
    "time_values  = np.zeros(steps)\n",
    "\n",
    "theta_values[0] = theta_0\n",
    "time_values[0] = 0.0\n",
    "\n",
    "theta_values[1] = theta_values[0] * (1 - const_k * delta_t)\n",
    "time_values[1] = delta_t\n",
    "\n",
    "for i in range(2, steps):\n",
    "    theta_values[i] = theta_values[i-2] - 2.0 * theta_values[i-1] * const_k * delta_t\n",
    "    time_values[i] = time_values[i-1] + delta_t\n",
    "\n",
    "exact_theta_values = theta_0 * np.exp(-const_k * time_values)\n",
    "    \n",
    "plot(time_values, exact_theta_values, linewidth=5.0)\n",
    "plot(time_values, theta_values, linewidth=3.0, color=\"red\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results are more accurate when a smaller timestep is used although it\n",
    "requires more computation to achieve the greater accuracy. Higher order expansion\n",
    "also increases the accuracy and may be more efficient in terms of the number of computations\n",
    "required for a given level of accuracy.\n",
    "\n",
    "Note, however, that the supposedly better quadratic expansion produces an error which\n",
    "oscillates as time increases. Does this error grow ? Does this make second order\n",
    "expansions useless ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Second Order Runge-Kutta\n",
    "\n",
    "<img src=\"images/theta_rk2-1.png\" width=\"66%\">\n",
    "\n",
    "The Runge-Kutta approach to higher order integration methods is\n",
    "illustrated in the figure above. The idea is to estimate the \n",
    "gradient \\\\(d \\theta / d t\\\\) at the half way point between two\n",
    "timestep values.  This is done in two stages. Initially a \n",
    "first order estimate, \\\\( \\hat{\\theta} \\\\) is made for the value of the function\n",
    "\\\\( \\theta\\\\) at \\\\(t=t+\\Delta t /2\\\\) in the future. This value is then\n",
    "subsituted into the differential equation to obtain the\n",
    "estimate for the gradient at this time. The revised gradient is\n",
    "then used to update the original \\\\(\\theta(t)\\\\) by an entire timestep.\n",
    "\t\n",
    "The first order step is\n",
    "$$\n",
    "\t\t\\begin{split}\n",
    "\t\t\\hat{\\theta}(t+\\Delta t /2) & = \\theta(t) + \\left.  \\frac{d \\theta}{d t} \\right|_t \\frac{\\Delta t}{2} \\\\\n",
    "         &= \\theta(t) \\left[ 1-\\frac{k\\Delta t}{2} \\right]\n",
    "\t\t\\end{split}\n",
    "$$\n",
    "\n",
    "Substitute to estimate the gradient at the mid-point\n",
    "$$\n",
    "\t\\left. \\frac{d \\theta}{d t} \\right|_{t+\\Delta t /2} \\approx -k \\theta(t)  \\left[ 1-\\frac{k\\Delta t}{2} \\right]\n",
    "$$\n",
    "\n",
    "Use this value as the average gradient over the interval \\\\( t\\rightarrow t+\\Delta t\\\\) to update \\\\(\\theta\\\\)\n",
    "\n",
    "$$\n",
    "    \\begin{split}\n",
    "        \\theta(t+\\Delta t) & \\approx \\theta(t) + \\delta t \\left(  -k \\theta(t)  \\left[ 1-\\frac{k\\Delta t}{2} \\right]  \\right) \\\\\n",
    "            & \\approx \\theta(t) \\left( 1 - k \\Delta t + k^2 \\frac{\\Delta t^2}{2} \\right)\n",
    "    \\end{split}\n",
    "$$\n",
    "\n",
    "It's worth noting that the Taylor expansion of the solution should look like\n",
    "$$        \n",
    "    e^{-kt} = 1 - kt + \\frac{k^2 t^2}{2!} - \\frac{k^3 t^3}{3!} + \\ldots\n",
    "$$\n",
    "\t\n",
    "The Runge Kutta method can be extended by repeating the estimates on smaller regions of the interval. The usual choice is fourth order RK. This is largely because, obviously, it's accurate to fourth order, but also because the number of operations to go higher than fourth order is disproportionately large. See Numerical Recipes for a discussion on this and better methods for ODE's.\n",
    "\t\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f49fa8ad6a0>]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8VPW9//HXZ2aykA0IWYBAwpIgIChCBBegKrairUtrr6J2sbXavb23/fVeu97W3v5uW+/113prF+1mrVaprRYV67WA4oYQQNmXsCasCUtCSMg2398fMxkmYckAk0xm5v18POYx53zPd2Y+h8A7h+/5zjnmnENERBKLJ9YFiIhI9CncRUQSkMJdRCQBKdxFRBKQwl1EJAEp3EVEEpDCXUQkASncRUQSkMJdRCQB+WL1wXl5eW7EiBGx+ngRkbi0fPnyWudcfnf9YhbuI0aMoKKiIlYfLyISl8xsRyT9NCwjIpKAFO4iIglI4S4ikoAU7iIiCUjhLiKSgLoNdzP7rZntN7M1p9huZvagmVWa2Sozmxz9MgNa2/1U7j/Ci6v38OCCzSzZeqCnPkpEJK5FMhXy98DPgD+cYvu1QFnwMQ34RfA5qh5ZvJX/mr+WogO7KT1QxZG0TI5+9CYuGTUo2h8lIhL3ug1359xiMxtxmi43An9wgfv1LTGzAWY2xDm3J0o1AjAoK5Vr1i3mwef+C4CXS6fxp2uujuZHiIgkjGiMuRcBVWHr1cG2E5jZPWZWYWYVNTU1Z/QhpQVZVA4qPr5+YCeV+xvOolwRkcTXqydUnXMPO+fKnXPl+fndfnu2k9H5WWzJLcKPAVB8eB/7ag7T1NLeE6WKiMS1aIT7LmB42PqwYFtUZab5yMsfwM4BgwHwOj8jD+xiS42O3kVEuopGuM8DPhacNXMJUBft8fYOZYVZbM47/nukrFZDMyIiJxPJVMg/AW8B55lZtZndZWafMbPPBLvMB7YClcAjwOd6qtiyE8bdq9i8/0hPfZyISNyKZLbMbd1sd8Dno1bRaZQVZLOky5H7s/t05C4i0lVcfUN1dEEWm8OO3MsOVGlYRkTkJOIq3EsLstgyaFhofcSh3eyuqaO5TTNmRETCxVW49++XQnbeAKpzCgBI8bcz7MButtc2xrgyEZG+Ja7CHQLj7l1nzOikqohIZ3EX7oFvqoaF+4EqNuukqohIJ3EX7mWFXU6qaq67iMgJ4i7cS/OzqAwbltFcdxGRE8VduJcVZncalhl1sJqd++tpbffHsCoRkb4l7sI9NzOVlEG57M3KBSCtvY0hB/aw44BmzIiIdIi7cIfASdXNXS5DoHF3EZHj4jLcywo7j7sHTqpq3F1EpEN8hntBNpvzul5ATEfuIiId4jLcA8Myx4/cx9TuZJPmuouIhMRluJcVZrEpryS0Xnqgiu176zRjRkQkKC7DPT8rDV/eIPZkDQIgva2FoTXVbKs9GuPKRET6hrgMdzNj7JBsNhSMCLWNq9nG+j31sStKRKQPictwBxg7OIf1BSOPr9dsZ8NezZgREYG4DvdsNuSHhft+HbmLiHSI23AfNySH9fkjjq/v386GPTpyFxGBOA730oIsduYNo9kbuA1s0ZEaGvfXcuhoS4wrExGJvbgN9/QUL8WF/dkcNiVS4+4iIgFxG+4AY4fksKHT0IzG3UVEIN7DfXB2p3H3sfu3sWGvwl1EJK7DffyQztMhx2lYRkQEiPNwD3yR6Xi4j6ndwebdh2n3uxhWJSISe3Ed7oNz0mkflMf+zIEAZLQ2M7h2ty5DICJJL67D3cyCX2YaEWrTuLuISJyHOwS+zLSu62UI9GUmEUlyCRDuncfdAydVdeQuIskt7sN97OCTzXXXkbuIJLeIwt3MZpvZRjOrNLN7T7K92MwWmdlKM1tlZtdFv9STG1OYzbawyxAMr9tHw579HNRlCEQkiXUb7mbmBR4CrgXGA7eZ2fgu3b4FzHXOXQTMAX4e7UJPpV+ql5IhA9kYdvQ+Yd8WVu+q660SRET6nEiO3KcClc65rc65FuBJ4MYufRyQE1zuD+yOXondm1DUn9WDS0PrE/dWskbhLiJJLJJwLwKqwtarg23hvgt8xMyqgfnAF6NSXYQmFvVndWF4uG9mdbXCXUSSV7ROqN4G/N45Nwy4DnjMzE54bzO7x8wqzKyipqYmSh8dDPcuR+4alhGRZBZJuO8ChoetDwu2hbsLmAvgnHsLSAfyur6Rc+5h51y5c648Pz//7Co+ifFDc6gsKAmdVC0OnlTVtd1FJFlFEu7LgDIzG2lmqQROmM7r0mcnMAvAzMYRCPfoHZp3IyPVp5OqIiJhug1351wb8AXgJWA9gVkxa83sPjO7Idjtq8DdZvYu8CfgTudcr169a0JRf9YUamhGRATAF0kn59x8AidKw9u+E7a8Drg8uqWdmdC4+7vB9b2beV7hLiJJKu6/odrhZCdVV2nGjIgkqYQJ9/FDc6jM10lVERFIoHDPSPUxXCdVRUSABAp3CAzN6KSqiEiihfuwruPum3UZAhFJSokV7vqmqogIkGDhPn5oDpu7fFP16K691DY0x7gyEZHelVDhnpHqY+TQXNYVjA61TdqziXd2Ho5hVSIivS+hwh3gouIBrBx6Xmh98q4NrNh5KIYViYj0vgQM94GsKBp7fH33BlbqyF1EkkzChfvk4gGdwn3Snk2s3nmAtnZ/DKsSEeldCRfuo/KyaMgfwr6sXACyWpoYtnsbm/Y1xLgyEZHek3Dh7vEYk0pyWTE0fGhmo8bdRSSpJFy4Q2BoputJVY27i0gySchwP/lJVR25i0jySMhwnzR8AKsLS2n1eAEoPVjNgaq9ukKkiCSNhAz3/v1SGF40iHUFo0Jtk/Zs4p0qDc2ISHJIyHCHU427a2hGRJJDwob7ScfddeQuIkkiYcN9cvHATtMhJ+3ZxLs7DtLu79X7douIxETChntpQRaHC4qoyRwAQE7zUYbs2srGvUdiXJmISM9L2HD3eozJI3JZOuz8UNu0natZuu1ADKsSEekdCRvuANNG5rJ0+ITQ+tSqtby97WAMKxIR6R0JHe5TR+bydli4T6taw9KtB3BO4+4iktgSOtwvGNaf7UNGcjg9C4D8xsP0r9rGlhpdRExEEltCh3uaz8ukklyWhY+7V63R0IyIJLyEDneAaSMH8fbw4+E+tWoNSxXuIpLgkiDcO59UnVa1hre3aNxdRBJbwof7RcUD2TS0lIbUfgAMPVKLr2oHVQebYlyZiEjPSfhw75fq5fziQSwvGhdqu6RqDW9rvruIJLCED3c4cUrkVJ1UFZEEF1G4m9lsM9toZpVmdu8p+txiZuvMbK2ZPRHdMs/NSee7K9xFJIH5uutgZl7gIeC9QDWwzMzmOefWhfUpA74OXO6cO2RmBT1V8NkoLxnImqFlHPOlkt7WQsnhvbRu38Huw00MHdAv1uWJiERdJEfuU4FK59xW51wL8CRwY5c+dwMPOecOATjn9ke3zHOTnZ7CmOF5VISNu1++411er6yNYVUiIj0nknAvAqrC1quDbeHGAGPM7A0zW2Jms0/2RmZ2j5lVmFlFTU3N2VV8li4vzeONEZNC69O3r+T1zQp3EUlM0Tqh6gPKgCuA24BHzGxA107OuYedc+XOufL8/PwofXRkZpTlsXjERaH1y7e/y5ub9uPX9d1FJAFFEu67gOFh68OCbeGqgXnOuVbn3DZgE4Gw7zOmlAxkS1EpB/rlAIHrzOTv2MT6vfUxrkxEJPoiCfdlQJmZjTSzVGAOMK9Ln2cJHLVjZnkEhmm2RrHOc5ae4uXiUXm8WXJhqG36tnc0NCMiCanbcHfOtQFfAF4C1gNznXNrzew+M7sh2O0l4ICZrQMWAV9zzvW5bwnNKMvjtbChmRnbV+qkqogkpG6nQgI45+YD87u0fSds2QFfCT76rOml+fx+5PGTqlOr1/L5zXs51tpOeoo3hpWJiERXUnxDtcPYwdm0DB3GltzAZJ/0thYmbl9DxfZDMa5MRCS6kircPR7j8tI8Xg+bEjlz+0peq+zdaZkiIj0tqcIdYHppHq+HjbtP366TqiKSeJIu3GeU5bOkeCJtFtj1Cfu2sKeyitqG5hhXJiISPUkX7oP7p1M4vJCVQ8eG2mZuW8HiTRqaEZHEkXThDnDFmHwWjS4Prc+qXMqC9X3qcjgiIuckKcP9qnEFLCidGlp/z7YVvLl+Dy1t/hhWJSISPUkZ7hePyGV30WiqcwLXt8lpPsp5W1dRsV3XeBeRxJCU4Z7i9TBzbAELRx8/ep9VuZQFGzQ0IyKJISnDHWDW2AIWjr44tH7VlmUsWL+PwJdtRUTiW9KG+xXnFbBkxAU0pqQBMPrgLmzzZrbWHo1xZSIi5y5pwz03M5UJowo7faHpqsqlLNSsGRFJAEkb7hCcNRM2NDNryzIWbNgXw4pERKIjqcN91tjCTuPuF1evZcOGKuoaW2NYlYjIuUvqcB9TmEXa8CLeHRy4aVSKv52ZlTp6F5H4l9ThbmZcPa6Ql8umhdqu3fgmL67ZG8OqRETOXVKHO8DsCYN58bzLQ+tXbq1g2ZqdNDS3xbAqEZFzk/ThfvGIXOpKRrMhrwQI3MDj8o1LWagvNIlIHEv6cPd6jGvO73z0fu3GN5i/ak8MqxIROTdJH+4A7584hPlh4X7V1mW8tWYnRzU0IyJxSuEOTB2Zy4GSUipzhwGQ0drMJZsqeGWjrvEuIvFJ4Q74vB6umTCYF8ZOD7Vdt/EN5q/R0IyIxCeFe9C1E4Z0GneftWUpb66uoqmlPYZViYicHYV70KWjB7G3pIytA4cCkNXSxFTNmhGROKVwD0rxenjf+YOZHzY0c9PaV3hm5a4YViUicnYU7mFuuLCIZ8ZfGVq/astSVryzhYNHW2JYlYjImVO4h7l09CAaRpWGrjWT1t7G7PWv8fyq3TGuTETkzCjcw3g9xk2Tinjm/ONH7x9cs4i/rtDQjIjEF4V7Fx+cXMRz42bSZoE/mot3rePA6g1sqWmIcWUiIpFTuHcxdnAOBaXFvDpqSqjtg2sX8axOrIpIHIko3M1stpltNLNKM7v3NP1uNjNnZuXRK7H3feiiIv56/lXH19cs5JkV1fj9unm2iMSHbsPdzLzAQ8C1wHjgNjMbf5J+2cCXgbejXWRvu3HSUBaWTaU+NQOAEYf3ULB2Jcu2H4xxZSIikYnkyH0qUOmc2+qcawGeBG48Sb/vAz8CjkWxvpgoyEnn4vHDOs15v2XVyzy1rCqGVYmIRC6ScC8CwlOtOtgWYmaTgeHOuReiWFtM3Ty5iLkXvDe0fsP6V3llWSWHGzXnXUT6vnM+oWpmHuAB4KsR9L3HzCrMrKKmpm9fcfGa8wezrXRi6CYeGa3NXLdqIX/RtEgRiQORhPsuYHjY+rBgW4dsYALwipltBy4B5p3spKpz7mHnXLlzrjw/P//sq+4F6Slebp4ynCcmzQ613fHOizyxZDvO6cSqiPRtkYT7MqDMzEaaWSowB5jXsdE5V+ecy3POjXDOjQCWADc45yp6pOJedNu0Yp49/0qafGkAjKvZTs67K1i6TSdWRaRv6zbcnXNtwBeAl4D1wFzn3Fozu8/MbujpAmNpdH4W48cX89y4GaG22979O08s3RnDqkREuhfRmLtzbr5zboxzbrRz7gfBtu845+adpO8ViXDU3uH2aSU8Mena0Pr161/j9WWVupiYiPRp+oZqN645v5CdpRNYVzASgH5tzVz/7j+YW6FpkSLSdyncu5Hm8/JP5cN5POzo/c7lz/HY61tobffHsDIRkVNTuEfg9mnFPDPhKg6nZwGBb6yOX/Ea81frHqsi0jcp3CNQMiiTGZNKOk2L/NSyZ/nN69s0LVJE+iSFe4Q+NWMUj07+AK0eLwDTqtbgr1hOxY5DMa5MRORECvcIlZcMpHDcaJ4bNzPUdteyZ/n1a1tjWJWIyMkp3CNkZtw1fSS/KT9+zbQPbHiNVUvWsuPA0RhWJiJyIoX7Gbhu4hAOnjeBJcMnAJDib+fOink8vFhH7yLStyjcz0CK18Odl43g1xd/MNT2kZXz+cfiteypa4phZSIinSncz9CcqcUsPf9S1uePACCz9RgfffsZfvWqjt5FpO9QuJ+h/v1SuHP6KP7nsjmhto8vf475r65lf33c36dERBKEwv0sfHL6SBZPnMGmQcUAZLc0cfvbz2rsXUT6DIX7WRiQkcpHLx/Fzy67JdT2yYp5PPvqOmobmmNYmYhIgML9LH1q+kgWTLyCLbmBOw7mNB/l9iXP8stXtsS4MhERhftZG5SVxh2Xj+Lnlxw/er976V957h+rqD7UGMPKREQU7ufk7hmjeGnSLCpzhwGBsffPvP4ED/zvphhXJiLJTuF+DvKz0/jEe0r50RV3htruWPkiFa8sZ+3uutgVJiJJT+F+ju6ZOYrlF85gWdF4AFL9bXz11cf44YsbYlyZiCQzhfs5yk5P4YuzyvjPKz4Rartx/ascem0Jr2+ujWFlIpLMFO5RcMe0EmovmMLfx1waavv2wl9z33NrdLcmEYkJhXsUpPo8/J9rzuPHMz/e6XrvY195gUff3B7b4kQkKSnco+QDE4eQM2kCv5tyQ6jtWwt/wyMvvMM+XZZARHqZwj1KPB7j+zdO4MHpt7EvKxeAgqOH+NTCx/i/89fHuDoRSTYK9yiaOKw/N84Yyw+uvCvU9omKeax/+U3e2nIghpWJSLJRuEfZ1645j9emzArd0MPn/PzgpZ/zzb+8y7HW9hhXJyLJQuEeZQMyUvm3a8fx7fd+NnRy9eJd65j5v0/x3/+7McbViUiyULj3gFvKh5M15UJ+Me3DobZ/XfwoLz/3Jst3HIphZSKSLBTuPcDjMe7/8AX8aubtbMgrASCjtZkfzv8p/zp3pYZnRKTHKdx7SGlBNl+YfT5fu+6fabPAH/MlVWu47OU/818vaXhGRHqWwr0H3T1jJJ6Ly/nVtJtDbd9Y9DteefZVXtm4P4aViUiiU7j3IJ/Xw48/fCG/mHkHG/MCt+Tr19bMz/72I77x+FLdc1VEekxE4W5ms81so5lVmtm9J9n+FTNbZ2arzGyBmZVEv9T4dN7gbL5w3QS+eMO/csyXCsDY2h189vlf8C9z38HvdzGuUEQSUbfhbmZe4CHgWmA8cJuZje/SbSVQ7py7AHga+HG0C41n98wYReFl5Xz/qk+F2j66cj5ZLzzHQ4sqY1iZiCSqSI7cpwKVzrmtzrkW4EngxvAOzrlFzrmOe8stAYZFt8z45vEYD9wyiZem38T8MZeF2u+f/xOefWohCzfsi2F1IpKIIgn3IqAqbL062HYqdwEvnmyDmd1jZhVmVlFTUxN5lQkgPzuN/zdnEvde+yWqcwoAyGlp5OG//Aff/N0bVO5viHGFIpJIonpC1cw+ApQD959su3PuYedcuXOuPD8/P5ofHRdmlOXz0Wsv5NMf/AZNvjQARh+s5r6//ohPP7qU+mOtMa5QRBJFJOG+Cxgetj4s2NaJmV0NfBO4wTnXHJ3yEs9X3nse+e+5lH+79kuhtvdWLuWmZx7m84+v0M09RCQqIgn3ZUCZmY00s1RgDjAvvIOZXQT8ikCwawL3aXg9xk/nXMTqmdfxy6kfCrV/8a2nGPr04/zbX1bhnGbQiMi56TbcnXNtwBeAl4D1wFzn3Fozu8/MOu5McT+QBfzZzN4xs3mneDsB+vdL4ZGPTeHn77uLV0ZOCbX/4KWHODz3GX6sb7CKyDmyWB0llpeXu4qKiph8dl+xYP0+/vmRxTzxxNeZuG8LAE2+NG6f8wOu//SH+OT0kTGuUET6GjNb7pwr766fvqEaQ7PGFfL1W6fxyQ9/l539C4HAN1h/+/T3+PNvn+exJTtiXKGIxCtfrAtIdrdPK2Zv/VQ+3nIfT//xawxqqmfgsSM8/uS3uM08eO16bp9WHOsyRSTOKNz7gH+5uoz99ce4s/l7PP7Ut8hpPkpuUz2PP/lNbsNwvJ87pumKDiISOQ3L9AFmxn/cNIGRs9/DR2+5j/rUDADyGuv405Pf4Ilf/o2HFlVqFo2IREzh3kf4vB4euOVChs++gjtv+R4Nqf2AQMA/+cS9vPbIn/n+8+t1oTERiYjCvQ/xeT385NZJDJ59JR+75T7q0jIByG5p4tG5/86u3z3Ol596R3dyEpFuKdz7GJ/Xw0/nXETx+2dxyx0/Yl9WLgBp7a38/NkfUvDIQ9z6yzfZp2vBi8hpKNz7oBSvhwdumcRlN13BzXf8mG0DhwDgdX6+veg33PHIfdz8k4W8W3U4xpWKSF+l2TJ9lMdjfOcD48nLSuPDKffzq2d+QPmu9QDcsvofjDy4m8/V3svdt07n45eNwMxiXLGI9CU6cu/DzIzPX1nKNz/xHj5+xw95esKs0LaLd63jb7/5Igse/COf/eMK6pp0RUkROU5H7nHgQ5OHMSIvk3uy+rExr4R7X/09Xucnr7GOR+f+Ow9VreX6nZ/kP2+dzOWlebEuV0T6AF1bJo7sPtzEpx9bTuZbr/HgvPspOHootG3V4FK+8v6vcNn1M7j32rFkpOr3tkgi0rVlEtDQAf14+rOXct6t13PdJx7ktZJJoW0X7K3khd9/mfSfPMDs+xexYL1u3SeSzHTkHqdeWLWHr/95JXNef5qvvvYYae1toW1rCkfzzfd9jsKrZ/Kd68czbGBGDCsVkWiK9Mhd4R7Hqg428rWn3+XA2yt54IUHQpcNBvBjPDFpNj+78mN86H2T+OwVo8lOT4lhtSISDQr3JOH3Ox59azv//cIaPvbG03zpzSdJb2sJba9Py+R/Lr2V52d8iE9fM545U4tJT/HGrmAROScK9ySzrfYo3352DTsq1vC9l3/JVVs7/9lW9S/kwcvm8Oal1/Hpq8dwS/lwhbxIHFK4JyHnHM+t2sP3n1vLpJWL+fqi3zLq0O5OfbYPGML/XDaHN6e9j4/OLOWOqSX0z9BwjUi8ULgnsfpjrfx80Rb+sHgz/1TxAl9+40/kNtV36lOdk8/vym/kb+XX8oHLz+MjlxRTWpAdo4pFJFIKd6H6UCP3v7SRBW9Xcufyedy99Bn6Nx/t1Kc+LZO5E6/m8Yuuo2DKRG6fVsw15w/WkI1IH6Vwl5ANe+v56T828/ryLXx8+XPcufw58hrrTuj3Wskknrrwfbw1YTrvnVzCBy8q4uIRuXg8um6NSF+hcJcTrN9Tz88WVbJwxXY+uGYhn1r6zAlj8gB1aZk8P24Gz46/gqpxFzH7giKumziEKSUD8SroRWJK4S6nVHWwkd+8vo25S3dw8aYKPrLyBa7aUoHX+U/ouzcrl/nnTefv513GtrILmTl+CLPGFXD56DydiBWJAYW7dKuuqZW/rqjmj0t2cKxyGzevWcDNaxZQcnjvSfsf7JfDotHlLBpVzpsjJzFiTDHTS/O4ZPQgJhcP1Di9SC9QuEvEnHO8ve0gTy+vZv6q3Zy/dRU3rXuF2RvfZFCXWTYd/BirhpTx+ohJLBk+kdXF5zNmdCGTSwYyuTjwyM9O6+U9EUl8Cnc5K40tbfx9zV6eX7WHNzfupXzru1y76Q1mVS5lcMPBU76u1eNl9eBSlheNY3nRON4Zch6eYUVMHD6AC4YNYNyQbMYP6U9hTppuLCJyDhTucs7qGlt5ad1eXl63j9c37aO0ajOztixl5rYVXLhn80nH6MPtzxzIqsGlrBlcyrqCUawrHEXD4CLGFOZQVphFWUEWo/KzGJWfydD+/TQrRyQCCneJqmOt7by+uZZXNu1n8aZa6nbt49Kdq7hk52ou2bmasbU7Inqf+tQMKvOGs3lQMZvyitk6aBjbBg5lf95QivKyKRmUwfDcDIYPzGDYwH4UDezH0P79GJCRoiN+ERTu0sN2HDjK65W1vL31IEu2HqB1334u2r2RKbvWM2XXeibs20JWS1PE79fq8VLdv4AdA4ayY+BgqnMK2dW/gF05+ezOyadhQC6DB2RSmJNGQXY6hTlp5GcHHnlZaQzKTGNQViq5mamkeHWbAklcCnfpNc45dhxoZMXOQyzfcYgVOw9TuecwJbXVXLB3M+P3bWX8/m2M37+VAccazuozWjw+9mUPYl9WLnuzBlGTNZD9WbnUZA6kNmMAtZkDOJiRw4F+/UnNymRgZioDM1Lon5FK/34pDOiXQk4/HznpKWSnp5Cd7gs9stJSyEj1kpXmIzPNR6pPvxyk74pquJvZbOCngBf4tXPuh122pwF/AKYAB4BbnXPbT/eeCvfEdqy1nbW761ldfZj1e46wbk89G/fWM+BwLWNqdzKmdiejD1Yx8uBuRh2sPu3J2jPVmJLGofQcDvfL5lC/bOrSs6hLz6Y+PZP6tEzq07OoT8vgSFomDan9OJqaEXzux9HUdNpS08lI85GR6iMj1Ut6ipd+qV76pXhJT/GQluIl3RdYTk/xkubzkObzkurzhB5p3sBzitdDitdI8XlI9QbWfV4jxRN89hq+0LIHr8fweSz4HFj3egyPoWEpASIP925vtGlmXuAh4L1ANbDMzOY559aFdbsLOOScKzWzOcCPgFvPrnRJBOkpXqaUDGRKycBQW1u7n50HG9m0r4HN+45QUXuUubVH2VrTQFt9A8V1eyk5tIfhh/cyrH4/w+r2MaxuP4OPHGDgsSMRf3ZGazMZrTUUHak5q9r9GEdT02lKSacpJY3G4HNTShrHfGk0+dJo9qXQHHpOpdmbyhFfCi1eHy3eQFtrcDn07PHR6vXR6vHR5vXS5vHR6vHS5g0+e3y0e7y0eTzBZ+/xZ/Pg9XrwmuHxEHwOBH9o2QK/BDwew2OBbWYElu34sscTeDYCvzBC7QaGgRFa7tjW8XulYzn0WghuC28//lpC24+/d2D5+C8r69Snc3tIR21hfbv26dxuJ28/5e/HU//iPNXnRdK/O1+aVUZBdnrkLzgDkdxFeSpQ6ZzbCmBmTwI3AuHhfiPw3eDy08DPzMxcrMZ8pE/yeT3B2TFZzJ4wONTunONQYys7Dhxl58FGqg42sulwEwsONrH7cBO765qwo40MOVJLwdGDFDQcZPCRA+QdPUxe42HyGw4xqKmOQY115DbWkeJvP6c6PTiyW5rIPoMah+EXAAAHo0lEQVRzBr2h3Ty0mwe/x0Obx4vfPPjNQm2B7V6cWaDd48FvHhxhy2ah1zkMF3z2B9udgeu0HfwEXgeE+jsj9LrAto73onO/4DIE+hNa7tzugoHYUVNHn5CwxHSc+Lrw9o7PP95+cq5TCp86kd0pNrlTvMZFmO4PXnYbd142kp66GGsk4V4EVIWtVwPTTtXHOddmZnXAIKA2GkVKYjMzcjMDJ0MvKh54wnbnHHVNreypO8b+I83sqz/G/vpj7G5o4d2GZmqONHPwaAsHj7Zw6Ggzmc1NDGiqZ2DTEQYcO0L/Yw0MONZAzrEGspsb6X+sgayWRrKaG8lubgwstzSRGXyktbfG4E+he17nD0w/Pf0MVIkTv59yfY++fyThHjVmdg9wD0BxcXFvfrTEMTNjQEYqAzJSGTfk9H3b/YFfBIcaWzjc2EJdUyt1Ta0cbmylvqmNA8da2drUSkNzGw3NbdQfa+Noc+DR0NxGY0s7tLWR2dJEemszGa3H6NfWTL+WZtLbmgPLrc2ktbWS3tZMelsLqe2tpLW1kNbWSmp78NHWSoq/jdS2VtLaW/H520hpDz78bfj87aS0B54D29rx+tvxuXZ8wWWv8+Pzt+H1+/F1850Cka4iCfddwPCw9WHBtpP1qTYzH9CfwInVTpxzDwMPQ+CE6tkULHI6Xs/x/wWcDeccLe1+jja309TaTlNL8NHazrHW48/NrX6a29o51uqnpd1Pc2s7h9v8NLcF1lva/LQGn1va/LT6Ha3Btla/o609sNzmd7S1B9bb/I52v6O13U+734XW253D+V3gyN3vx+Pa8fn9eFzg4fWHPePH4/fjcS7Q5vxYp2WwsHYDPM4PjsBrnQv2d5jzB8bAg6/vWO4YkAj1wwXfK/hP2hHWFhC+PbydYD8g7P05/l6AhSVF+HvQTZ/j/U7UeUDm1FEUyes7OYOR6JrME/+XGk2RhPsyoMzMRhII8TnA7V36zAM+DrwFfBhYqPF2iUdmRprPS5qvb10EzR8M+Xa/w+8Cwe/3O/yOUFvHs99P4Dn0ILTduUD+tDuHcy4w9u0C7R39nAPH8b6OwHt07R/YFt4OdLwuWHdHCnS8X2A58BpO1ycsPcL7dwqVTn1OfL8T+oe/9BSvPX2/CJxB7H0NyM/quesvdRvuwTH0LwAvEZgK+Vvn3Fozuw+ocM7NA34DPGZmlcBBAr8ARCRKPB7Dg6ELb0qkIhpzd87NB+Z3aftO2PIx4J+iW5qIiJwtfRVPRCQBKdxFRBKQwl1EJAEp3EVEEpDCXUQkAcXskr9mVgNEdoeHE+WRnJc2SMb9TsZ9huTc72TcZzjz/S5xzuV31ylm4X4uzKwikkteJppk3O9k3GdIzv1Oxn2GnttvDcuIiCQghbuISAKK13B/ONYFxEgy7ncy7jMk534n4z5DD+13XI65i4jI6cXrkbuIiJxGnw53M5ttZhvNrNLM7j3J9jQzeyq4/W0zG9H7VUZXBPv8FTNbZ2arzGyBmZXEos5o626/w/rdbGbOzOJ+VkUk+2xmtwR/3mvN7InerrEnRPB3vNjMFpnZyuDf8+tiUWc0mdlvzWy/ma05xXYzsweDfyarzGzyOX+o67iucx97ELi88BZgFJAKvAuM79Lnc8Avg8tzgKdiXXcv7POVQEZw+bPxvs+R7newXzawGFgClMe67l74WZcBK4GBwfWCWNfdS/v9MPDZ4PJ4YHus647Cfs8EJgNrTrH9OuBFAvcBuQR4+1w/sy8fuYduzO2cawE6bswd7kbg0eDy08AsszO593if0+0+O+cWOecag6tLCNwZK95F8rMG+D7wI+BYbxbXQyLZ57uBh5xzhwCcc/t7ucaeEMl+OyAnuNwf2N2L9fUI59xiAve6OJUbgT+4gCXAADPr5qaSp9eXw/1kN+YuOlUf51wb0HFj7ngVyT6Hu4vAb/t41+1+B/+bOtw590JvFtaDIvlZjwHGmNkbZrbEzGb3WnU9J5L9/i7wETOrJnAfiS/2Tmkxdab/9rvVqzfIlugxs48A5cB7Yl1LTzMzD/AAcGeMS+ltPgJDM1cQ+B/aYjOb6Jw7HNOqet5twO+dc/9tZpcSuMvbBOd0l/Az0ZeP3M/kxtyc7sbccSSSfcbMrga+CdzgnGvupdp6Unf7nQ1MAF4xs+0ExiTnxflJ1Uh+1tXAPOdcq3NuG7CJQNjHs0j2+y5gLoBz7i0gncD1VxJZRP/2z0RfDvfQjbnNLJXACdN5Xfp03JgbEuPG3N3us5ldBPyKQLAnwhgsdLPfzrk651yec26Ec24EgXMNNzjnKmJTblRE8vf7WQJH7ZhZHoFhmq29WWQPiGS/dwKzAMxsHIFwr+nVKnvfPOBjwVkzlwB1zrk95/SOsT6L3M0Z5usIHK1sAb4ZbLuPwD9sCPzQ/wxUAkuBUbGuuRf2+R/APuCd4GNerGvujf3u0vcV4ny2TIQ/ayMwHLUOWA3MiXXNvbTf44E3CMykeQd4X6xrjsI+/wnYA7QS+B/ZXcBngM+E/awfCv6ZrI7G3299Q1VEJAH15WEZERE5Swp3EZEEpHAXEUlACncRkQSkcBcRSUAKdxGRBKRwFxFJQAp3EZEE9P8BKPVHD6HTzegAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "steps = 100\n",
    "theta_0 = 1.0\n",
    "const_k = 10.0\n",
    "delta_t = 1.0 / steps\n",
    "\n",
    "\n",
    "theta_values = np.zeros(steps)\n",
    "time_values  = np.zeros(steps)\n",
    "\n",
    "\n",
    "theta_values[0] = theta_0\n",
    "time_values[0] = 0.0\n",
    "\n",
    "for i in range(1, steps):\n",
    "    theta_values[i] = theta_values[i-1] * (1 - const_k * delta_t + const_k**2 * delta_t**2 / 2.0)\n",
    "    time_values[i] = time_values[i-1] + delta_t\n",
    "\n",
    "exact_theta_values = theta_0 * np.exp(-const_k * time_values)\n",
    "    \n",
    "plot(time_values, exact_theta_values, linewidth=5.0)\n",
    "plot(time_values, theta_values, linewidth=3.0, color=\"red\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}