{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Numerically solving differential equations with python\n", "\n", "*This is a brief description of what numerical integration is and a practical tutorial on how to do it in Python.*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Software required\n", "\n", "*In order to run this notebook in your own computer, you need to install the following software:*\n", "\n", "* [python](http://python.org)\n", "* [numpy](http://numpy.org) and [scipy](http://scipy.org) - python scientific libraries\n", "* [matplotlib](http://matplotlib.org) - a library for plotting\n", "* the [ipython notebook](http://ipython.org/notebook.html) (now renamed to [*Jupyter*](https://jupyter.org/)\n", "\n", "On Windows and Mac, we recommend installing the [Anaconda distribution](https://www.anaconda.com/distribution/), which includes all of the above in a single package (among several other libraries), available at https://www.anaconda.com/distribution/\n", "\n", "On Linux, you can install everything using your distribution's prefered way, e.g.:\n", "\n", "* Debian/Ubuntu: sudo apt-get install python-numpy python-scipy python-matplotlib python-ipython-notebook\n", "* Fedora: sudo yum install python-numpy python-scipy python-matplotlib python-ipython-notebook\n", "* Arch: `sudo pacman -S python-numpy python-scipy python-matplotlib jupyter\n", "\n", "Code snippets shown here can also be copied into a pure text file with .py extension and ran outside the notebook (e.g., in an python or ipython shell).\n", "\n", "### From the web\n", "Alternatively, you can use a service that runs notebooks on the cloud, e.g. [SageMathCloud](https://cloud.sagemath.com/) or [wakari](https://www.wakari.io/). It is possible to visualize publicly-available notebooks using http://nbviewer.ipython.org, but no computation can be performed (it just shows saved pre-calculated results)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## How numerical integration works\n", "\n", "Let's say we have a differential equation that we don't know how (or don't want) to derive its (analytical) solution. We can still find out what the solutions are through **numerical integration**. So, how dows that work?\n", "\n", "The idea is to approximate the solution at successive small time intervals, extrapolating the value of the derivative over each interval. For example, let's take the differential equation\n", "\n", "$$\\frac{dx}{dt} = f(x) = x (1 - x)$$\n", "\n", "with an initial value $x_0 = 0.1$ at an initial time $t=0$ (that is, $x(0) = 0.1$). At $t=0$, the derivative $\\frac{dx}{dt}$ values $f(0.1) = 0.1 \\times (1-0.1) = 0.09$. We pick a small interval step, say, $\\Delta t = 0.5$, and assume that that value of the derivative is a good approximation over the whole interval from $t=0$ up to $t=0.5$. This means that in this time $x$ is going to increase by $\\frac{dx}{dt} \\times \\Delta t = 0.09 \\times 0.5 = 0.045$. So our approximate solution for $x$ at $t=0.5$ is $x(0) + 0.045 = 0.145$. We can then use this value of $x(0.5)$ to calculate the next point in time, $t=1$. We calculate the derivative at each step, multiply by the time step and add to the previous value of the solution, as in the table below:\n", "\n", "| $t$ | $x$ | $\\frac{dx}{dt}$ |\n", "| ---:|---------:|----------:|\n", "| 0 | 0.1 | 0.09 |\n", "| 0.5 | 0.145 | 0.123975 |\n", "| 1.0 | 0.206987 | 0.164144 |\n", "| 1.5 | 0.289059 | 0.205504 |\n", "| 2.0 | 0.391811 | 0.238295 |\n", "\n", "Of course, this is terribly tedious to do by hand, so we can write a simple program to do it and plot the solution. Below we compare it to the known analytical solution of this differential equation (the *logistic equation*). **Don't worry about the code just yet**: there are better and simpler ways to do it!" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAECCAYAAADelD2uAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deVzU1f7H8dcZNsUF3DfAXVNTU6zQFrPMyjaz5WrZpqlltv7MttsN267d23bbM8sWy/bMPfc9NxLTRFRQcBcRBESWmfn+/kBJBEtzhu8MvJ+PB4+Z+Z4Z5jPgvD2c+Z5zjGVZiIiI/3DYXYCIiJweBbeIiJ9RcIuI+BkFt4iIn1Fwi4j4GQW3iIifCbS7gOPdc889VkREhN1liIjYbsyYMR9ZlnVPWW0+FdwRERHExsbaXYaIiO3GjBmz82RtGioREfEzCm4RET/jU0MlZcnMzGTPnj12l1EpNWrUiPDwcLvLEJET+HxwHzhwgGbNmlG1alW7S6lUjhw5wq5duxTcIj7I54dKCgsLqVKlit1lVDpVqlShsLDQ7jJEpAw+H9wAxhi7S6h09DMX8V1eCW5jTENjzBRjTM+jtxsYY14wxvzbGBPqjeeMS8ngnQVbiUvJ8Ma395gLLriAuLi4035cYmIiAKNHj+bNN9/0dFlSSf3d982ZvN/87bG+mC1eGeO2LGuvMaYmcKzb1g9wAhHARcDPnny+uJQMbhu/ggKnm+BAB1/cE0N001qefAqPmTBhAi1atDitx2zevJkHHniA2bNn88gjjxASEuKl6sQfxaVksCI5nZgWdU7r333c9oPcNn4lBS43QQEO/jfgHNo2rEmhy02B002By03hsUuXmwKnRYHLzZZ92by/KAmnyyLAYRgUE0Xj8Kq43OC2LJwuC5dl4XYfvXS5cbsKMa58MrKyWZa4B4fl5HvcxDSvSZ2qDozbhbGcOKyiS+N24bAKMZYLh+XC4XaSm5dHyoFsAnGyEzdxtUMIDTI4LBdg4cANloXBjbHcGKziy4JCJweyj2CwWDXPYlv1IEICDAY3Do49xjp6+4/rLpeLrCMFNMTNrgXQrFUd6oQGARZY1p9c8sft3mOgXhuP/s7L68PJ4OOul/gb3BgzDBgGcM011/zpN2n2xPS/fKK8Qjc3vrf8T++zfezVZR5fuXIlc+fOZcmSJfTo0YNnn32WsWPHMmnSJObNm8djjz1GfHw8Xbt2pXr16owdO5bo6GhiYmL45ZdfmDlzJk899RQdOnSgfv36XHnllXTv3p1hw4aRk5PD/fffT48ePRg/fjzr16/ngw8+YPDgwcyfP58hQ4bw0ksvsXbtWiZNmoTb7Wb8+PH88MMPjB8/noSEBKZPn84bb7xBhw4dGDVqFE899RSNGzfm4osvpmPHjkRFRfHiiy/y0UcfsWDBAqKiov7y5yX+yelyk364gMWJaTw1eT1Ol0WgwzDgvChqVQvmcL6T3AInOfkucvOd5OQ7OVzgJDffRU6+k9wCF4fznVhYVKGAMOdh3vtuBg2ruKnuyKeGI59qJp/qR79CzRGqkk9N8uiYm83rjhxCHAWEUEid9RbhwW4CrUIC3QUEWgUEugsIsAoIcBddt0wATkcw+VYguUEOCgnEZTkISQuhakgIliMAtwnEchz9MoFYjoCjl4EQEMBut5NURwFOAnASQGRADaJqVQdHAJYJAOMAY7CMAwgEhwMoOpa4/zB7sgtxWUW3m4aF075xOBhH0f2PPc4cfYwjAMs4WJt6iJXbM3BbBoyhanBjLm/TECi6DRy9NCe/rFbX479/rwS3MSYcqAVEGWPmAXcAjwCZwOLj72tZ1jhgHEBsbOyfbsdzssA91uMudLoJOoMed05ODhERERQWFrJ27VoAhg8fzpw5c/jmm2+Iiori4MGDxMbGFodkvXr1CAsLY/Xq1Xz++edkZmYyYMAAunTpwuDBg3nuued49NFH2bZtG+Hh4VSvXh2AyMhI6tatS//+/Xn++edZunQpgwcPJjk5mYKCApo1a0Zqaiq//vor7du3p06dOlx99dV8++23ALz99tvUqFGDPn368OCDD7JhwwYA+vbty7hx44iPj1dw+5FjveYukeE0qVWV/dn57M/KJy07r+j60a+07KJjmbmFhIcG4zBQ6Cp627jcLnbt3kHrJtCGLMJNBmHBmVR3ZFIlOJsqrmyCC7MJKjxEQP4h3LkZWEcO4caQRTVqhtWmSmhNCK52wld1CKpdfH17tuF/S3Zy2BWEyxHM6N6dadKkHgQEQ2AVCDx6WXw7BOMIIBhYf+J7ddCpv1d3pWTwxPGPvSGGyFN87IGUDMaNX0Gh6+hj+8bQ7BQem5mSwT+PPWeAg1u7x4AP/DXvraGSTKDz0ZufHb0c5Y3nAohuWosv7on5W38uHs/hcLBo0SKaNGlCVlZW8fGgoKASH9Ydu33sWFRUFNWqVcPtdmOMweVyUVBQgGVZ9OnTh6CgICZPnsxdd931lzW43W7y8vJISkqiatWqWJZV5geFJz5XWe3iuw4dKSRhTxYJe7JYsjmNBZvTiv/CrlcjmMhaodSvUYV6NUKoXyOE8yOrERFQQENyqOPcT7W8vQRk7+LQvu3sLthObQ5Rm2xMZnUC3Q2gWv2inl71+hBaF6o2hyrhUDW86LJKGAFVw/l1v8UvqYdP633TDBjU5o/hmban8X47k/eqHY/1VLZ4ms+fx32qopvWOuMf6urVq9m0aRPt2rUrDsvPP/+crKwsbrzxRt555x127NjBq6++yoMPPojD4SAtLY0dO3YAcOutt7Jo0SImT57Myy+/TFpaGm+88QbDhg3j2WefpU+fPuTk5JCamkpOTg6ZmZkkJycDsGrVKgD27NnD/PnzqV+/Pk2aNCE5OZlevXpx5MgR4uLi2LVrF9nZ2bzyyis888wzLFy4kJdeeomdO4uWNdi8eTM5OTnFt6X8lDXWbFkWOzOO8PvuLDYeDeqNu7PIzC2gbcMatGtUk8AABwaoSh4tHHt5oJXhioaH4WASHNgKW1IgLxNqNISwSAiLgJpNoFFnws66mt3ZVZmZFkCnNq3o2qLBadXctQZ0bXn6r/VM3m/+9lhPZIunGV/aLDg2NtY6cZGphIQE2rVrV+61bN++nebNm5ORkVE8CSU2Npb4+HgmT55c7vXYwa6fvT86/gPyAIehd7sGpOcUkLAni+pVAmnfqCbtGtWkfeOatKtfhaauHTjSNsLe9WSlrOPIrg3U5DCpNKR+s/bUimwHtVtCnVZQqylUbwCOALtfppQjY8wYy7Jiy2qrMD1uT0tJSQEgOTmZrl27Fh/bu3cvubm5hIZ65axG8TOWZRG/I5OXZmwkr7BoeMrtssgrdPFw79a0axBKrZytsHM17FwDS9YV9aRrNYMGHaDB2dS85EG2Fzbmu33BxLSsd1pDD1I5KbhPomfPnpz418iECRNsqkZ8iWVZJO7LZkr8bqb+tpsgh4NuzWqxYVcWwa5cYgITeb7mIRou2QB74ouGNSLOhcjz4PxhUK8dBJWcDdwJ6NTBntcj/kfBLXKKth84zNR1u5mybje5BS6u6dyI9wd0pL07EbNtOs80mU+VtA3k1etIjdq9oPMj0CQaqqoHLZ6l4BY5wfEfMjYOr8K0dXuYsm43ew7lcXXHhvz3mqZ0yluNY/NE+HIu1G4BzXtS4/InISqGoOBqdr8EqeAU3CLHiUvJ4NYPiz5kNAZCgwO4umNjnr6kPufmLiJg0//gu7XQ7AJo2xeu+DfUOL0zOUTOlF8sMuWrFi1ahDGGzMzMMtuPrS9yquuT3H333fTr1++Unjs3N5cOHTrwxhtvnPQ+iYmJuN1uWrRowd69e0/p+1ZmeYUuXp+zmXynG4ui0/NeOWsLL+e/QMy0ywjY8QucNxxGJcKtX0P0nQptsYV63GegadOmJ2374YcfWLp0Ka+99topr0/StGlTMjJObSGb0NBQ6tWrd9L2119/HWMMbdu2Zdq0aTRs2PCUvm9ltWhzGs/+tIH6NULoEpjCP8xsrnKsgpxucN5tcNNHEFLD7jJFAPW4i61cuZIXX3yRK6+8khkzZrBjxw6MMbz00ku0bNmS1NTUUvc5xul0EhkZyYsvvsjixYv59NNPef3111m3bh2//vorPXr0YNq0aSQkJPDYY4/Rs2dPtm3bxnPPPceECRPo3r072dnZpWpKTU1l6NCh9OvXD6fTyQMPPMC4ceMYMWJEiTNeXnjhBc455xy++OILwsPDyczM5P3332f58uUsWrSIDh06EB8fz8SJE4mNjeX2228nIyOD559/nsjISEaPHs2AAQPK5efsa/YcOsKIL+J4/sc1vNchgW8cT/NV2Ns0a9me7QMWEDZ0CnT+h0JbfIp/9bhjwzz0fQ6VOnT8OiVz5syhb9++QMm1P6pVq1biPu3btwcgMDCQkSNHMm/ePM4//3xuv/12FixYQHh4ePGCVFAUsBdddBF9+vQhNDQUYwxRUVGsWLGieAbl8TIzM/nqq68YO3YsiYmJfPzxxxw+fJioqCgGDx5cfL+IiAgAmjRpAkB4eDiNGjWiR48e9OzZs/h+o0aNYurUqbz11lt89NFHREZGUqdOHfr27Uv//v0987P1E4UuNxOWbWPSgrW8GLGS7uYHTEY0XPIEIa16E6PJLuLD/Cy4Sweupxy/TklZs0ndbvef3mfw4MHExsYyaNAgHA5HmeuL5OXlsWbNGu69917cbjfp6enFa6KU9ZyNGzfmgw8+YPjw4Sxbtqz4e564PonDUfoPp9Nd36QyrW2yMjmdt3+cxx3WNOYGLyKg3nVw/Syo29ru0kROiYZKjjq2TklQUBApKSnF648cv/bHifc5NrsyNTWVevXq0a9fPzp3LlpbKyYmhoSEBFJSUorXJxk9ejSzZ8/mySefJD09nXnz5rF8+XLatGlDUlISKSkp7Nu3j/z8fAA2bNjAunXruOWWW2jatCm33XYbH3/8Mb1796Zz586kpaWRmppKdHQ0ubm5zJkzh/z8fA4cOMAFF1zAqlWrStQ4duzY4iGeIUOGsHPnTjIzM9m2bVvx4yqytOx8YifOZefE+/gofxS9OzcnYORKuO4thbb4Fa1V4iF79uxh/fr19OnTx+5SPMZffvZ/xeW2+G5JPPkLX+Nmx0ICut1B8MWPQrU6dpcmclJaq8TLcnJyuPbaa3nnnXfsLkWOOjaJpm6oIWvRu/wj71vc7a6j6pUroWZju8sTOSMKbg+oXr06a9assbsMOerYJJrz3PH8K/BzgmpHUOOu2Zj6Z9ldmohH+EVwn2wzAfEeXxpCO13L1m7gDfMq7QJTeck1iM6dBnJ/fY1hS8Xh8x9OBgUFkZeXZ3cZlU5eXh5BQUF2l3F6LIv0JR8xKP5WkojgqsKXWew4l5iWnt/zT8ROPt/jrlu3Ltu3b7e7jEqpUaNGdpdw6jJSOPL9CPbv2kNSj4/o3vZcjI9tNyXiKT4f3OHh4cU70IiU6bdvcc58nI8LrqL+lW9z8/nNARTYUmH5fHCLnFR+NkwfRX7qaoYUPkH/q/vSv2uE3VWJeJ3Pj3GLlGl3PLx/ERkF0DvnOW659hqFtlQaCm7xP/GTYGJ/UrqM4vKtN/HU9dFc11nnZkvloaES8R+uQvj5adg6h41XTOL2KVmMvbETl7fXmthSuajHLf7h8AH47HrI2MavV/zI7VOyeOWWzgptqZQU3OL70pPgo8sh8jx+Of9dhn6zhf8N6EKvtvXtrkzEFhoqEd+2YzV8fRtc8iRLwq7h4UnxvH1rV7q31AJRUnkpuMV3JUyFqQ9Bv/dZ4D6HUV/F8/7t0ZzbrLbdlYnYSsEtvmntRJj3PAz6ntkZjXjyh3V8eGc3ukZpUo2Iglt8z6oPYenrbOjzBeMWwaLNvzFxyPl0jPDQ1nUifk7BLb5l+Vuwahzr+0yi/6RdFLosggMdFLgqz9ZqIn9FZ5WI71j8X4j7BO6eyfQdwRS6ipaWdbncrEhOt7c2ER+i4BbfsOxNWPc13DUDV40mLN1ygECHIcBAUKCDmBY6i0TkGA2ViP1Wjy/6unsm1GjABwu3EhoSyKRhMazadlBLs4qcQMEt9oqfBEteg7umQ1gT4ndk8vHSbUwZeSGNw6vq1D+RMii4xT4bp8DcZ+HOqVC7OTn5Th76ai3PXX82jcOr2l2diM9ScIs9UpbDtIfh9h+hXlsAnv3pd7q3qEPfjn60846IDfThpJS/tM3wzR3Q/0No1BmAn+J3sTY1g39d297m4kR8n3rcUr5y9sMXN0HvWGh1GQA7DuYyZupGPht8HqHB+icp8lc8/i4xxvQC2gCRlmX90xhzCXApEAU8b1lWkqefU/xEwWH48hboPBC6DALA6XLz0Fdrua9nS85uopmRIqfCG0Mlw4BEYKAxJhSwKAruLOCgF55P/IHbDT8Mg3rt4JInig+/OX8r1UICGXJhcxuLE/Ev3gju4OOuG6AqMBG4GDjrxDsbY4YZY9YYY9bExcV5oRzxCYv/UzRMcu0bYAwAq7Yd5MuVqbx6c2ccDmNzgSL+wxvB/T7QBVgMTAGqAQ2AhcD+E+9sWdY4y7K6WZbVLTo62gvliO0SpsGvn8E/PofAEAAO5RbyyNfxvHxjR+rXrGJzgSL+xeNj3JZlzQHmnHD4e08/j/iJ/Qkw9UG47Vuo0RAAy7J4avJ6ererz2XttPWYyOnS6YDiPUcyYNJA6PMiNPnjr6lv43aydV8OT/ZtZ2NxIv5L516Jd1gW/HgftLkCzhlYfDg5LYexMzcxaWgMVYICbCxQxH8puMU7lr8Fh/fDLZ8VHypwunnoq3ge7t2atg1r2FiciH9TcIvnpa6E5W/C0PkQ+MdJRq/OTqR+jRBuj2lqY3Ei/k/BLZ51OB2+GwzXvQXhUcWHl245wOT4Xcx48CKM0al/ImdCwS2e43bD5Hvh7Bug7VUAxKVkMH/TPr5cmcqbA7tQp3qIzUWK+D8Ft3jOyveLziS57FmgKLRvG7+CvEI3AQ6jdUhEPESnA4pn7PsdlrwC/cdBQBAAK5LTyS88usmvZWnfSBEPUXDLmSvMg++HwuXPQe0WxYe7RIZjAQ7tGyniUfrbVc7c/OehTgs457YSh1dsO8gFLevQo1Vd7Rsp4kEKbjkzyQthw/dw77LixaMAdmbk8tkv25n+4EU00TZkIh6loRL5+45kwOT74fq3oVrJYZB/z9jEXT2aKbRFvEDBLX/frCeLTvtr1bvE4RXJ6cTvyGT4xS1tKkykYtNQifw9m2cXbfh73/ISh11uizFTN/LEVWdRNVhrkYh4g4JbTl/eoaId2vu9CyHVSzR9tTqVGiGBXNNJO7WLeIuGSuT0zX4GWl8OLS4pcfhQbiGvz9nMv65tr2ntIl6kHrecnqQFsHUejPilVNMb8zZzefuG2vRXxMvU45ZTl59TtJvNtf+DKjVLNG3Zl83ktbsY1aeNTcWJVB4Kbjl188ZAs4ugdcmzSCzL4rlpG7m/VystIiVSDjRUIqdmZxxs/AlGrCjVNC9hP7szj3Bnj2blX5dIJaQet/w1l7PoLJLLn4fQ2iWa8p0uXpi+kWeuaU9QgP45iZQHvdPkr60aB1XCoNMtpZomLNtOy3rVuaRtfRsKE6mcNFQif+7QLlj8Xxgyu8RaJAD7s/P4YFESP4y4wKbiRCon9bjlz816As4bCnVbl2r676xEbukWSfO61WwoTKTyUnDLyW3+GfZtgAsfLdW0bkcmCzenMfLSVjYUJlK5aahEylaQCzNGFZ2zHVSlRJNlWcRO/Z3H+rSlRpUgmwoUqbzU45ayLXkFIs6FlpeWapocvwuny+Km6AgbChMR9biltIPJsGYC3LesVNPhfCcvz0zkndu64HBoPRIRO6jHLaX9/DT0eABqNi7V9N7CJGJa1Ca6ae0yHigi5UHBLSVtmQv7E6D7/aWaUtNzmbgyhcevOsuGwkTkGAW3/MFZUHT635VjIbD0miMvzUhgyAXNaRSm7chE7KTglj+s+gBqNYM2V5Rq+mTZNpYnHaBbM+3ULmI3BbcUyd4HS16DK/9daobk6u0HGTN1I9l5Tu7+ZDVxKRk2FSkioOCWY+aNgS6Dypwh+fHSbViABRQ63axITi/38kTkDzodUGDnmqJdbUauLtWUV+hi9faDBAc4cLndBAU6iGlRx4YiReQYBXdlZ1kw83Ho/WypXW0AJq5I4ZzIWtx3SUtWJKcT06IO0U01zi1iJwV3ZbdxMrgKoNOAUk1ZeYW8tzCJScNiaNOghgJbxEcouCszZwHMew6ufg0cpT/u+GBREpeeVZ82DWrYUJyInIw+nKzM4j6BWs2hZa9STfuy8pi4IpVHLtfmvyK+xuM9bmNML6ANEGlZ1j+PHhsIHARWWJZ1yNPPKX9DXlbRBgm3/1Bm8//mbeEf50bSOFyTbUR8jTd63MOARGCgMSbUGNMTuAPoBWgLcF+x/E1odRk07FiqKSkth5nr9zDikpY2FCYif8UbwR183HUDRAO/AgHAoBPvbIwZZoxZY4xZExcX54VypJSsPbB6PPR6uszmV2cnMvTiFoSHBpfZLiL28kZwvw90ARYDU4DvgfpAITDtxDtbljXOsqxulmV1i46O9kI5UsrCf0OX2yE8slRT/I5M4lIyuLtHcxsKE5FT4fExbsuy5gBzTjg81NPPI3/T/k2waTo8sKZUk2VZjJ2ZwMO921A1OMCG4kTkVOiskspmbixc+DBULX1O9qLNaezPzudm7Wwj4tMU3JXJ9mWw73c4t/QfQG63xcuzEhl9RVsCA/TPQsSX6R1aWVgWzPkXXPZMqc1/Aaas201IoIMrOjS0oTgROR0K7spi409FU9vPvqlUU77TxatzEnniqrMwRvtIivg6TXmvDJwFRcu2nmRq+5crU2lZr7pW/RPxE+pxVwZxnxTtbFPG1PbsvELeWZDE6Cu0j6SIv/jL4DbGdDvhtk7t8yfHprb3HlNm84dLtnFx67q0b1x6SVcR8U2n0uMeYYzpZYy5yhjzC0UTbMRfHJva3qhTqaa07Hw++2W7FpIS8TOnMsa9DXgLaA8sBx7yakXiOcemtg9fUmbzW/O30L9LBJG1Q8u5MBE5E6fS4x4DJAHnHb1+jlcrEs/5k6ntKemHmbpuNyMvbWVDYSJyJk6lx32rZVlfHbthjCn0Yj3iKfs3waZp8EDZC3e9Mnszgy9oTu1qWkhKxN/8ZY/7+NA+enuh16oRz5k3Bi58pMyp7Rt2HWJlcjpDLtJCUiL+SKcDVkTbl8HeDWVObQd4edYmHrisNaHBOo1fxB8puCuaY1PbL/1nmVPbl245wI6DuQw4t/S4t4j4BwV3RbPxJ3DlQ8ebSzWt2X6Q//s2nv5dIwjSQlIifkvv3orEVVg0tn35c6WmtselZDDwwxXsy8rn3YVbiUvJsKlIETlTCu6KpHhq+6WlmpZtTaPQZQFQ6HSzIjm9fGsTEY/Rp1MVRV4WLPoPDPq+zOYDOQU4TNEmoEGBDi0oJeLHFNwVxfI3i3raZUxtT8/JZ+q63bwx4Bx2HDxCTIs6RDctfZqgiPgHBXdFUDy1fXGZzf+bt4Xrz2nCdZ2blHNhIuINGuOuCBb+G7oMgvCoUk1b92cz/bc9PHRZaxsKExFvUI/b36Ul/unU9henJ3DfJS2ppantIhWGetz+bm7sSae2L96cxrYDh7mje7NyL0tEvEfB7c9Slp90arvLbfHi9ASeuKodwYH6NYtUJHpH+yvLgtnPnHRq+zdrdhAWGsQVHRrYUJyIeJOC21/9ydT2nHwnr83ZzDNXt9eu7SIVkD6c9EfHprZf/WqZu7a/t3ArF7euR8eIMBuKExFvU4/bH/3J1PadGbl8sTKVx65oW+5liUj5UI/b3/zF1Pb/zErkzu7NaBhWetxbRCoG9bj9zZ9Mbf81NYNV2w4yvGcLGwoTkfKiHrc/ydp90qntlmXxwrSN/F+fNtrZRqSCU4/bnyx4EbreUebU9unr91DgcnNj1wgbChOR8qSumb/Y9zts/hlGrinVlFfoYuzMTfz3ps44HDr9T6SiU4/bX8z5F1z0f1A1vFTThGXbad+oJt1bao1tkcpAwe0PkhZAehJ0G1Kq6UBOPuMWJ/Fk33Y2FCYidlBw+zq3G+Y8A72fhcDSK/y9PmczN3SJoHndajYUJyJ2UHD7ut++hsAq0L5fqabN+7KZtWGv1toWqWT04aQvKzwC81+Amz6CMtYceXF6AiMvbUVYaJANxYmIXdTj9mUr3oPG50BUTKmmhYn72XEwl0ExTW0oTETs5PHgNsb0MsYMN8a8cNyxMGPMUmNMM08/X4V1+AAsfwt6jynV5HS5eXF6Ak/2bUdQgP7vFalsvPGuHwYkAgONMaHGmHCgK6CZIadj0X+g401Qt1WJw3EpGdz3xa+EBDro3a6+TcWJiJ28EdzHn/pggIHAnUBdoNeJdzbGDDPGrDHGrImLK3vfxEonPQnWfws9Hy9xOC4lg9s+XMGcjfvYvD+HX1MzbSpQROzkjeB+H+gCLAamWJb1nmVZdwEHgAUn3tmyrHGWZXWzLKtbdHS0F8rxQz8/DRc8BNXqlji8IjmdPKcbAJfLzYrkdDuqExGbefysEsuy5gBzyjjezNPPVSElzYe0TXDLp6WawqoWnT0SYCAo0EFMC82UFKmMdDqgL3E5YdZT0OcFCAwp0ZRX6OLjpdsYfUVbLCCmRR2im5be2V1EKj4Fty+JmwDV68FZV5dqemPuFto1qsmIXq3KeKCIVCYKbl+RexAWvQx3/FRqss1vOzP5Lm4HMx+62KbiRMSX6CRgX7HoZWh3LTToUOJwgdPN6O9+459Xt6dejZCTPFhEKhP1uH1BWmLR6X/3ryrV9N7CJBqHV+X6cxrbUJiI+CIFty/4+amitbZPOP0vcW82n/6ynekPXogpY60SEamcNFRit82zIWM7nDu0xGGny83o79bx2BVtaRRW1Z7aRMQnqcdtp8I8mDka+v631FrbHy/bRrWQQAacG2lTcSLiqxTcdlr+ZtGHka0vL3F424HDvLcwiZ/u1xCJiJSm4LZLxnZY8S4MX1zisNtt8fj3vzHy0tZE1Qm1pzYR8Wka47bJZUIAAAnVSURBVLbLzCeg+0gIjypx+IuVKThdbu7q0cyeukTE5ym47ZA4E9K3QI8HShzemZHL63O38J+bOhHg0BCJiJRNQyXlrfAIzHwcrnm9xHoklmXx1I8bGHJhc1rVr2FjgSLi69TjLm9LXy/ajqzVZSUOf//rLg5k5zPs4hY2FSYi/kI97vKUngSrPoR7l5Q4vD8rj3/PSODTwedpKzIR+UtKifJiWTD1Ibh4FIRFHHfY4pmfNjDwvCjObhJmY4Ei4i8U3OVl7UQoyIHz7y1xeMb6vSSlHeaBy7Rcq4icGg2VlIfsfTA3Fu6YDI6A4sMZhwuInfo77w+KJiQw4OSPFxE5jnrc5WHW49D1dmjYscTh56Zt5NpOjbWTjYicFvW4vS1xFuxZB/3eKz4Ul5LBlytTWL71AAseu8S+2kTELym4vSk/G6b/H/R7F4KKVviLS8ng1g9XkO90ExzgIGFPtnrcInJaNFTiTXPHQItLoEXP4kPLkw6Q73QD4HK7WZGcbk9tIuK31OP2luSFkDgD7ltW4nDCniwcBgwQFOggpkUdW8oTEf+l4PaGvEPw00i47k2o+scwyDerd7BpTzaf3H0e63cdIqZFHQ2TiMhpU3B7w89PQaveRV9Hrdl+kJdnbeKbe7vTsl51Lm5Tz8YCRcSfKbg9LXEWbFtSYohkd+YRRnzxK6/c0pmW9arbWJyIVAQKbk/KPQjTHoYbx0NI0Qp/RwpcDPt8DUMubE6vtvVtLlBEKgKdVeIplgXTH4UON0CzC48eshj9/W+0rl9Dq/6JiMeox+0p8V/A/k0lJtq8uzCJ1PTDfD28u/aOFBGPUXB7woEtMOdfcOe04ok2czfu4/NfUph8/wVUCdI6JCLiORoqOVPOfPhuMPR6Chq0B2DLvmwe//433h3UlYZhVWwuUEQqGgX3mZo7pmjD325DAMjMLeCez9bwZN92dI3SOdoi4nkaKjkTW+bAxp+KdrQxBqfLzcgv13J5uwbcFB3x148XEfkbFNx/V2YqTB4BN38CobUBeGnGJhwOwxNXnWVvbSJSoSm4/w5nPnxzB1zwIDS7AIBv1uxgQeJ+Jo+4gEDtGykiXqSE+TtmPg5hkdB9JFC0VOvLMzfx4R3dCAsNsrk4Eano1OM+XfFfwvalMHQ+GMOeQ0cY8UUcr9zcmVb1NZ1dRLxPwX069vwGs/8Jd02HKjX5JekAD38dT5/2Del1lqazi0j50FDJqcreB1/dCn1fgfrtWLUtndvGr2RfVj7fxu0gLiXD7gpFpJLweHAbY3oZY4YbY144erupMeZHY0yCMcY/1zItzCsK7S6D4Oz+5BW6eOrHDbito81O7WQjIuXHGz3uYUAiMNAYEwrssSzrBuAQUOqTO2PMMGPMGmPMmri4OC+Uc4YsC6aMhFpNoefjZOcVcteEVTSoEUKVIAcBRjvZiEj58sYYd/Bx141lWQXGmHOAZ4FS3VLLssYB4wBiY2MtL9RzZha/AulJcPcM0g8XcNeE1XSKCOO5688mfkcmK5LTtZONiJQrbwT3+0AXYDEwxRjzMPA1kAc8CczwwnN6x/rvIO4TGDqPXYfh9vG/cHWnRjx6eRuMMUQ3raXAFpFy5/HgtixrDjDnhMNtPf08Xrd1Hsx6Au6YwtYj1bjjo+UMuagFQy5sbndlIlLJ6XTAsuyKgx+GwT8msq6gMUM+XcmTV53FjVp/RER8gIL7RAe2wqSBcN1bLCtszYOTVvPyjZ3o3b6B3ZWJiAAK7pIO7YSJN8Cl/2SWswtPf7uWd2/ryvk6Y0REfIiC+5is3fDptXDeML5yXsJrM3/n08HncXaTMLsrExEpQcENkL23KLS73sn7hX2ZuGQrXw/vTvO61eyuTESkFAV3zn749FqszgMZm3UF8zft5Lt7e2jLMRHxWZU7uA/thM9vwNXhRp7a34fEfQf5Znh3alUL/uvHiojYpPIGd3oSfNaP7a1uY+ja7oQGZ/Hl0BiqhVTeH4mI+IfKuTrg3vUwoS/xzYfSa3lHtuzPIXFfNpv2ZttdmYjIX6p8wb1tCdbnNzCn2aMMiGsNWuFPRPxM5Qru+C9xfXMnsUGPMD69E6/dfA4hWuFPRPxM5RjQtSzc818gZ/Uk7sr/J9dfdCnPxjTF4TA0CKuiFf5ExK9U/ODOz+Hwt/eyK2UL/wl/ldeHXkzTOn+cn60V/kTE31To4Hbv38yhT//BosNNOdTrE8ZddBYOh7G7LBGRM1Jhgztt1XcEz3qUr6rdyZX3P07zetqBXUQqhgoX3O6CPBK+fIza22ew7Nx3GHbVNQSoly0iFUiFCO64lAxWJKfTKXg3kQseIjugIaGDF9AvKsru0kREPM7vgzsuJYPbxv/CAPcs2gf+wOKo+7nurscJCKhcZzqKSOXh9+m2IjmdGPda+gUs5ebCWHa3uFmhLSIVmt/3uGNa1OEtRxduKeyIIzBIk2hEpMLz++COblqLL+7prkk0IlJp+H1wgybRiEjlosFgERE/o+AWEfEzCm4RET+j4BYR8TMKbhERP6PgFhHxM8ayLLtrKGaMGQ/s/JsPjwbiPFiOL9NrrZj0Wiumv/taIyzLuqesBp8K7jNhjFljWVY3u+soD3qtFZNea8XkjdeqoRIRET+j4BYR8TMVKbjH2V1AOdJrrZj0Wismj7/WCjPGLSJSWVSkHreISKWg4BYR8TMKbhERP1MhgtsY08sYM9wY84LdtXjb0de6yhjzo921lAdjTDtjzEq76/A2Y0ygMeYBY8x5xphgu+vxFmPMKGPMncaYn4wxxu56vMUY09AYM8UY09MY819jzGPGmHae+v4VIriBYUAiMNAYE2p3MV62DLgSaG53Id5mjGkNtAQa2F1LOXgaaANcDThtrsWbAoFeQLJVgc+MsCxrL1ATOAs4F0gD7vDU968owX18D6XC/i8OYFlWAXAZ0N8YE2J3PV42ErgLqGuM6WhzLd4WDUwB+gLn2FyLNwUBk4GbjDFV7S6mHHglmyrE1mXA+0AX4AfLsg7bXYw3GWMGAc8B+UD3o5cVkmVZDxljmgHdLMtab3M53vYeRb/PNcBGm2vxpppACDAXKLS5Fq8xxoQDtYAsin6fEcAXHvv+FfivFRGRCqmiDJWIiFQaCm4RET+j4BYR8TMKbhERP6PgFjmBMaaBMaYynDsufkrBLXIcY0xt4Aegkd21iJyMglukpBFAD+B+Y0xLu4sRKYuCW6SkpUcv37EsK8nWSkROQsEtIuJnFNwiJbmPXTHGBNhZiMjJKLhFSvoN2ALcCVSGRZDED2mtEhERP6Met4iIn1Fwi4j4GQW3iIifUXCLiPgZBbeIiJ9RcIuI+BkFt4iIn1Fwi4j4mf8HEztpKRhl7nMAAAAASUVORK5CYII=\n", "text/plain": [ "