Numerical Integration Tutorial - R.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numerically solving differential equations with R\n",
    "\n",
    "*This is a brief description of what numerical integration is and a practical tutorial on how to do it in R.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Software required\n",
    "### R\n",
    "*In order to run the R codes in this notebook in your own computer, you need to install the following software:*\n",
    "\n",
    "* [R](http://www.r-project.org/), along with the packages from [CRAN](http://cran.r-project.org/):\n",
    "  * [deSolve](http://www.vps.fmvz.usp.br/CRAN/web/packages/deSolve/index.html), a library for solving differential equations\n",
    "  * [ggplot2](http://www.vps.fmvz.usp.br/CRAN/web/packages/ggplot2/index.html), a library for plotting\n",
    "  * [reshape2](http://cran.r-project.org/web/packages/reshape2/index.html), for manipulating data.frames\n",
    "\n",
    "To install R, download it from its homepage (Windows or Mac): http://www.r-project.org/. On Linux, you can install it using your distribution's prefered way, e.g.:\n",
    "\n",
    "* Debian/Ubuntu: `sudo apt-get install r-base`\n",
    "* Fedora: `sudo yum install R`\n",
    "* Arch: `sudo pacman -S r`\n",
    "\n",
    "To install the packages, all you have to do is run the following in the `R` prompt\n",
    "\n",
    "    install.packages(c(\"deSolve\", \"ggplot2\", \"reshape2\"))\n",
    " \n",
    "The R code presented here and some additional examples are available at https://github.com/diogro/ode_examples (thanks, [Diogro](https://github.com/diogro)!).\n",
    "\n",
    "### Running commands in  R\n",
    "You can also follow the instructions in this notebook and run the commands in R in two ways:\n",
    "* Copy the R codes snippets in this notebook and paste them in the R promt. Alternatively you can save them in a pure text file with .r or .R extension and run directly in an R shell.\n",
    "* Download the [R script](https://github.com/diogro/ode_examples) and run the commands in a R shell.\n",
    "\n",
    "### Running R commands from Jupyter\n",
    "*To run the R commands in this notebook from IP[y] you need also:*\n",
    "\n",
    "* the [Jupyter notebook](http://jupyter.readthedocs.org/en/latest/install.html) and\n",
    "* the [R Kernel for Jupyter](http://irkernel.github.io/).\n",
    "\n",
    "To install the notebook on Windows and Mac, we recommend installing the [Anaconda distribution](https://store.continuum.io/cshop/anaconda/), available at http://continuum.io/downloads. On Linux it should be available from your distro's repositories.\n",
    "\n",
    "The R kernel can be installed running, on the R shell:\n",
    "\n",
    "    install.packages(c('rzmq','repr','IRkernel','IRdisplay'),\n",
    "        repos = c('http://irkernel.github.io/', getOption('repos')))\n",
    "    IRkernel::installspec()\n",
    "\n",
    "### Running from the web\n",
    "If you for some reason don' want to install anything on your computer, 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 on http://nbviewer.ipython.org, but no computation can be performed (it just shows saved pre-calculated results).\n",
    "\n"
   ]
  },
  {
   "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 does 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 the derivative value is a good approximation to the function 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": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAADAFBMVEUAAAABAQECAgIDAwME\nBAQFBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQVFRUW\nFhYXFxcYGBgZGRkaGhobGxscHBwdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYnJyco\nKCgpKSkqKiorKyssLCwtLS0uLi4vLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6\nOjo7Ozs8PDw9PT0+Pj4/Pz9AQEBBQUFCQkJDQ0NERERFRUVGRkZHR0dISEhJSUlKSkpLS0tM\nTExNTU1OTk5PT09QUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1e\nXl5fX19gYGBhYWFiYmJjY2NkZGRlZWVmZmZnZ2doaGhpaWlqampra2tsbGxtbW1ubm5vb29w\ncHBxcXFycnJzc3N0dHR1dXV2dnZ3d3d4eHh5eXl6enp7e3t8fHx9fX1+fn5/f3+AgICBgYGC\ngoKDg4OEhISFhYWGhoaHh4eIiIiJiYmKioqLi4uMjIyNjY2Ojo6Pj4+QkJCRkZGSkpKTk5OU\nlJSVlZWWlpaXl5eYmJiZmZmampqbm5ucnJydnZ2enp6fn5+goKChoaGioqKjo6OkpKSlpaWm\npqanp6eoqKipqamqqqqrq6usrKytra2urq6vr6+wsLCxsbGysrKzs7O0tLS1tbW2tra3t7e4\nuLi5ubm6urq7u7u8vLy9vb2+vr6/v7/AwMDBwcHCwsLDw8PExMTFxcXGxsbHx8fIyMjJycnK\nysrLy8vMzMzNzc3Ozs7Pz8/Q0NDR0dHS0tLT09PU1NTV1dXW1tbX19fY2NjZ2dna2trb29vc\n3Nzd3d3e3t7f39/g4ODh4eHi4uLj4+Pk5OTl5eXm5ubn5+fo6Ojp6enq6urr6+vs7Ozt7e3u\n7u7v7+/w8PDx8fHy8vLz8/P09PT19fX29vb39/f4+Pj5+fn6+vr7+/v8/Pz9/f3+/v7////i\nsF19AAAACXBIWXMAABJ0AAASdAHeZh94AAAgAElEQVR4nO3deWATZf748acFyiUIiAcqnrhe\n+1UUUfdU1uOnK0VROWQBL1QExRMvUFRUvPDCAw9gV1lUdIVVxANFQRBRAQXkvqGU9rOriHK3\nnV+epE2TNsnM5JnZNpn3+49kZvL042S3L9qmaaIsIjJO1fQJEGVDQCLyICAReRCQiDwISEQe\nBCQiDwISkQcBiciDvIH0o12/lG61XZNGv2z3Y+rmUn/G7vRj6o+lu3wZu9uXqbtK/Bn7kx9T\nd5Rutluy2WNIYtcW61fbNWn083Y/pv5o+TN2lx9Txdrty9hSX6aWlPkydvd//Zi6w/rRbsmP\nQEoRkARIAiTjgCRAEiAZByQBkgDJOCAJkARIxgFJgCRAMg5IAiQBknFAEiAJkIwDkgBJgGQc\nkARIAiTjgCRAEiAZByQBkgDJOCAJkARIxgFJgCRAMg5IAiQBknFAEiAJkIwDkgBJgGQckARI\nAiTjgCRAEiAZByQBkgDJOCAJkARIxgFJgCRAMg5IAiQBknFAEiAJkIwDkgBJgGQckARIAiTj\ngCRAEiAZByQBkvgMqWzN5rhrIDkfC6SgQ9p4b8XW8ssvzH9yd+U1kFyMBVLAIe145LqKr0d9\nn965tsf70WsguRkLJG8grZu7qcqRqpCWvfPvValnFE8e8XZB6iWLx77wsd2puIL0TJf8Ckjz\nOm2xrLH9o9dAcjMWSLJhjc2Cmd1POPO54lQr5p+Xo/YYVBR3rAqkBxsp1WxEqiE/tFdKHTY1\n1ZKXm4aWnLU+9em6grRh+UsVkCZcH7qY26mk4hpIbsYGHtK7v1HqwH+kWvF+XuizV12aYkVB\nW71CDYo7GA9pbHiFmpRiSofwikPWJl/xVcPwkj6pztb1z0gTKiCNGhy6WJn/34rr0OVHzzzz\nzMvb7Npp7bJdk0Y7dvsy1SrxY+z2Uj+mbrP8GVvm+iNeO+Oov36YasH3e+hPzQafplhyRMTA\nlBT/lciKRv+NPVi6PXbvtMiSzsmHLIqsUK9sjLTyh2gLZ0a6IrKi8a+p7lHof/30ID13X+hi\nXf76iuvQ5aB27dqd6WQIZXd3hz/xXkux4srI5+ZZyVcUlH+C35d8yZDyJf/6duaUKVPeGT/+\nzRdDPRxq0O23337b1aGaRlY0PfPMM09rpzvxsHCHNI9UX7no55R3ujS65Q7S2NtCF4vzt1Rc\nhy5XzJ49e85mu7Za223XpNHWnX5M/cXyZ+xuP6Zutkp8GVtWZX/1BzP/k2r9nMjnXZONyZf8\nLrLkwOQrviv/7L3w788+Nezeu2684bLLulyQf/rp7du2PeSQ1s2aNXJDIFqTZrrmh0T6bdu2\nbY8tv+W400OddUG0bpdFax9ZsU/q/5W2pAnpg96hi88uLqu4rliQ+htJ4Wek8NgM/hmp+IbQ\nDy+Hp/qpY0T55+b7yZecF1lxfGhzxYLZn0x8fdQTw+6+sW/vzh1PO+X4w1qnZLJHs2YHHXzo\n8cefcGrdiMbefQYMGHDX3XffM3z48CdfHz1q9Nuh3v0k1IzHIx8ybvnypI9adA+vOGFj8pNd\ndkB4ySMp7rKk/zPSj+cvtqxhw6PXQHIzNoMh3R/+rGqxIPn658o/4z9IcFvB4lkfjh81PPLZ\nq/ber0l1KHnNDjj4+BNPOyk3/IXi7vuHPzVq1FtvT/5k+rffLl0e95jAq/o7t2MXxf0Hqjxq\nd2/D0Dd2T6a6d2v0ufx5Xqol004I/SQ2OOUDiGlCunOSZY3sM3lkt1WV10ByMTaDIe0T+XS/\nI/n6uXkRaxvCewULp08c/fig/pece+rR+8d/pclp1vro9qd16n7FgEEPDn/ltbc/nv7tkpXR\nMTO6H39G6oe/ZfEzg16r8oukqr9HWjr+7RU292/xpLk2K2Tl/EK7JS4hffG0vhw82bLKJg8b\nudKqvAaSi7GZC2lDOYJLUnxA5ItWz5t6/L+TDo39ktN4/6NPPafb1QOHPj367VFDhk7242Qz\n5JkN9tmeE5BqN6S5M6v+uBD/Fal5BMUtCT50ybRxTwzs/f9OaFUn+jVnryNOOfdvN94/4rVJ\nM36If/pAsJ8iBCRvxtZaSB8crdSeVX6sjod0S+SLyzeVR1ZM++cjN1x0ykF55Xrq7HvcWZfc\n9ODId6Yvin/WQVxAApL52NoKaX6LMIWX4w7GQyrUP5y3HBvaWjt97EPXnHvMHuV+9v7t2b1v\ne+r1z39IoScmIAHJfGxthXRjxMRv4g5WfYrQVyMeeXZQj9/tV/4gW5sOPW8fMeFrm+d9VgtI\nQDIfW1shnR/BUS/uYAyk799++MrT9o8savX77rc/N2mBzaPCSQMSkMzH1lZIfSJG9o87GIa0\n5J1hvdtHHoLb+w+9h/x9us2ToW0DEpDMx9ZWSFMiTz67LfZY4fwX+v+5pT6ce+i5A575YLnp\nfyQSkIBkPra2QpLh+k8GLqz47ePGqY/3PrFB+GtUh+tGfGL6RSguIAHJfGythSTzRz72qb4u\n/uqFq9prQ3WO7HHv20vNB1cNSEAyH1t7Iek2vDf4bP04eO6R3R/5aH1t/lPz6gHJMCA5gLRx\n3LCxNg9Tr3itX3v9m9X9Og6ZuDpyCEhAMi27IH2l/+r00OlJb189rt/xdUJfiI65/PnYp3EC\nCUimZRWkTf8X+W1rwq9Jhe/d1K6uUnXb3fBm1VfdARKQTMsqSB+VP5fnnWq3zH38vCZK1Wnb\nf9zqBB8HJCCZllWQxpVDejHuaOHEa9uEDh7Ya0yyv9sBEpBMyypIX5ZD+qTy0OrRXZor1eAv\nQ2em+DggAcm0rIJU/ky6syqeHbf0yQ55SrX826vrUn8YkIBkWnZBWtk9V+VctCy8vfzpDvWU\n+s2AyfZ/9AAkIJmWXZBEVn0e/kFo1YgzQ1+Ljr7zK0djgQQk07INkm7T+IsbKnXU7al+LIoL\nSEAyLfsgfXH9fkq1vvkLF2OBBCTTsgzS6sePV2qP7hOd/TF4RUACkmlZBemzy/ZQuac/n+J9\nFxIHJCCZlj2QCkadptQ+A+akMRZIQDItWyB9P6CZyjltdIoXuE4RkIBkWnZA+uSieqp5/9np\njgUSkEzLAkhFY0Pf0x3yoOufjCoDEpBMy3hIBU8dqtQfx7p7mK5KQAKSaRkFaeXMqq/SWPBI\na1Wva8o3GnYQkIBkWiZBureRUvuPizlQMLyVyuv6tfFgIAHJtAyC9GL4qd0NZ1Tsr7unpWpw\n1fceTAYSkEzLIEiRvyNXl0f2Nj6+n2p83aLUH+IwIAHJtAyC1CwC6S96u/iVw1T9/l69GB2Q\ngGRaBkE6KgKpZ2jz7bYqt1M6z2FIHJCAZFoGQXo07Kj+xzL9NJXT0fHfSDgISEAyLYMgFV+b\np1TTEcv61FV/+NDTyUACkmkZBElk7vh3Fz/bUh3wrMdzgQQk0zIKklifHqsaDtzg9VggAcm0\njIL0Q3eV09mLXxxVCUhAMi2DIBU/3Uy1fc/7uUASIBmXOZC+PU01eHiH52N1QAKSaZkCqfDB\nRuoPs716f6QqAQlIpmUIpI+PUXs9790bjVUJSEAyLSMgbby1ruqqnw8EJCAZF1xIs05QrcaH\nt4AEJOOCCql4eCOVX/70VCABybiAQlp0tmoyvGIHSEAyLpiQxrRQf/ouugckIBkXREgFfVT9\noTGvawIkIBkXQEhzT1BtpsUeABKQjAsepH82V+fFv+0rkIBkXNAgFQ7MzXuwyjEgAcm4gEGa\nf4o6aErVg0ACknHBgjR5H3Xu8mpHgQQk4wIF6dm83LuLqx8GEpCMCxCk4oGq8auJbgASkIwL\nDqTV56qDE78FLJCAZFxgIH1zpDplceKbgAQk44ICaUJzdVmyN94DEpCMCwikl/PqPpL0RiAB\nybhsh7T+3Ve/FXkwt/H45GuABCTjshzShAOUUt1vVi0mp1gEJCAZl92QFrSIvET+QV+lWgUk\nIBmX3ZCGRBzlLki5CkhAMi67IV0VgaRSv+URkIBkXHZDuj/iqMmmlKuABCTjshvS4r3CkG5N\nvQpIQDIuuyHN3U87urIw9SogAcm4rIY0e3/V/9PxC+2WAQlIxmUzpFmt1G1O1gEJSMZlMaSZ\n+6o7HS0EEpCMy15In7XIqfriDEkCEpCMy1pIM/fKedThUiABybhshTTvQHWf07VAApJxWQpp\n/sFqkOOxQAKScdkJaemR6hrnY4EEJOOyEtLK49UVLsYCCUjGZSOkdaeqrkX2y6IBCUjGZSGk\ngr+ov9o8KSg+IAHJuOyDVNRJnVHgaiyQgGRc9kG6VrVf724skIBkXNZBekwdkuTl65IGJCAZ\nl22QxtZpkfL1GRIFJCAZl2WQPm7UINXrBSUOSEAyLrsgfdsyd4z7sUACknFZBWnp4eqhNMYC\nCUjGZROk9Sep/umMBRKQjMsiSMUXqPPdPKEhGpCAZFwWQRqkTtyQ1lggAcm47IE0vs4+36c3\nFkhAMi5rIM3as957aY4FUhZB+tGuX61ttmvS6Jcdfkz92fJl7ObdSW9ae6R6Kt2xVvKxJpX6\nMrWkzJ+xm/2YutOyHbvZY0i77CqxSmzXpFFJqR9Td1v+jC1LdsuOv6rr0h5rJR1rlE9T7T9V\n0hrry9RSa7fdkp0eQ7L9Ksm3dqm+tbteneruGd+x8a1dFn1rZ3tOQEoBaXRO6yXpjwUSkIzL\nCkifNWg0zWAskIBkXDZAWnVYzmiTsUACknHZAOliNy8ZlCAgAcm4LID0sDoh/QcadEACknGZ\nD+nzBs3mmI0FEpCMy3hIKw7O+YfhWCABybhMh1TcMb0/nYgNSEAyLtMhPaDamf2AJEDSAcmw\nDIf0cV6zucZjgQQk4zIb0vLWOa+bjwUSkIzLbEid1QAPxgIJSMZlNKTn1fEbPRgLJCAZl8mQ\nvm/WcJYXY4EEJOMyGFLRH9VjnowFEpCMy2BIg9Vfij0ZCyQgGZe5kKbm7fWDN2OBBCTjMhbS\nhqPVWBEvHmsAkgDJuIyF1EddWjz8kJyW168xHgskIBmXqZDeyjlk9VClO894LJCAZFyGQlq2\nb92P1jYIQ1ITTccCCUjGZSikruo2mRpxpB4wHQskIBmXmZBeV8dslNnlkJ4yHQskIBmXkZBW\nHVD3k9DVsWFHjeebjgUSkIzLSEi91E366ouWIUf1XzQeCyQgGZeJkCbkHBF5/5aVD11+1zfm\nY4EEJOMyENK6Q3Lf93QskIBkXAZCukpd6+1YIAHJuMyDNLnOQeZPZogLSEAyLuMgFRyZ87bH\nY4EEJOMyDtL16gqvxwIJSMZlGqSpdQ9c7fVYIAHJuAyDVNROjfN8LJCAZFyGQXrcgyd7VwtI\nQDIusyAVtmho/nqQ1QISkIzLLEg91b0+jAUSkIzLKEjv5xzlyd+WVwlIQDIukyAVHptj/Ed8\niQISkIzLJEiD1WU+TAWSDkiGZRCkeY2aG7+DS8KABCTjMgjSOeqZxG/GbBqQgGRc5kB6VZ34\nHyABCUhmrWtdd1ridzU3DkhAMi5jIN2m+iZ+V3PzgAQk4zIF0veN9loBJAESkMy6UA1P+GbM\nXgQkIBmXIZAm5xy7CUg6IAEp/YpOVPrPYoEEJCCZNEJ1DI8FEpCAlH5r9sv7OjwWSEACUvrd\noG6MjAUSkICUdnPq770qMhZIQAJS2nVUI8rHAglIQEq3Ceq4ovKxQAISkNKs6NicyRVjgQQk\nIKXZs6pzdCyQgASk9NpwYN430bFAAhKQ0utedXXlWCABCUhptaLFHosrxwIJSEBKq+vVXTFj\ngQQkIKXT/Ib7ro0ZCyQgASmduqsnYscCCUhASqOZddsUxo4FEpCAlEZnqX/EjQUSkIDkvvdU\nu+K4sUACEpDc115Nih8LJCAByXWj1blVxgIJSEByW2GbOjOqjAUSkIDkthGqe9WxQAISkFy2\n8eB631YdCyQgAcllj6ve1cYCCUhAcldB67xq77sMJCAByWXDVJ/qY4EEJCC5akOr+vOrjwUS\nkIDkqvvUtQnGAglIQHLT2paNFiUYCyQgAclNg9QNicYCCUhActGq5o2XJBoLJCAByUW3qVsT\njgUSkIDkvBXN9lyecCyQgAQk590U+4onsWOBBCQgOW55kxarE48FEpCA5LiBalCSsUACEpCc\ntqZF0xVJxgIJSEBy2uDED9kJkHRAApKz1rVstDTZWCABCUgOe0Bdn3QskIAEJGcVtKq/MOlY\nIAEJSM56TF2VfCyQgAQkR208qF61P4ytHAukgEMqW7O5fGvjgnA/Wev01WogVWlE9VdqiBkL\npGBDWn75hflP7g5vjswPN9UapK8GAym+TW3qfF2+uaDn4W0u/SF+LJACDams79M71/Z4P7xd\nsivUnMu2Wn0+C22UACm+F6OvZbdkXxXqgGVxY4EUaEjzOm2xrLH9o/u7+82ySs5fFrfG9pwC\nAan46NwvyjcvU+HiHnkAUrAhTbg+dDG3U/TLz1v3Wdam/PUz524FUnx/VxdUbB4bgdQ2biyQ\nAg1plP5RaGX+f8t3f+0e+mL0Xf7fbujZa6Hef/3222+/f4ddu6zdtmvSaFeJH1N3WumNPSnn\nm4rNdhFIp8SNLTU8r8RZ/owt82eq5c/YnX5MLbXsx7qA9FzoK5C1Ln99+e4/9e53Q9dbux++\nUn+VGtSuXbszbYcEoU9Vx+j24Aik+2vwdOh/UGl0yx7S2NtCF4vzt0T2Sv42q+KGtfnrQpf/\n3bBhw8Yf7frV2ma7Jo1+2eHH1J+ttMaeoSZFtzcepx2dWBh7++bdhueVOMufsaW+TC0p82fs\nT35M3WlttluyOcrEHtIHvUMXn11cFtmbdYl+IPwnzerH/NUVa2y/3QzAz0jTc06I2St4sGPH\nYQXxY/kZKdA/I/14/mLLGja8fO+xx/XlyCt3WtaU7ruBVFkXNcZmLJACDcka2WfyyG6rLOvO\nSaGdXuFfKBVeestbL100JbrE9pyyH9L39Q4vshkLpGBDKps8bOTK0PXgyaGfiO5aFz62eezQ\nZxdWLrE9p+yHdI0abjcWSMGG5CDbc8p6SMsbt1xvNxZIQAKSTXcme8mTmLFAAhKQUlewT9K/\nMK8cCyQgASl1jyV6H5eqY4EEJCClrOiwevPsxwIJSEBK2WjVzcFYIAEJSCk7KeczB2OBBCQg\npWqiOsvJWCABCUipOltNdDIWSEACUopm5x7vaCyQgASkFF2hXnA0FkhAAlLyljfer8B+FZB0\nQAJS0u62f3ZQZCyQgASkpBUe0MD22UGRsUACEpCS9pK61OFYIAEJSElrlzPT4VggAQlIyZqs\nznA6FkhAAlKyOqnxTscCCUhAStLcukcVOx0LJCABKUn91JOOxwIJSEBK3Oo997J7qYbKsUAC\nEpAS95C6xflYIAEJSAkralNvvvOxQAISkBL2qpO/jI2OBRKQgJSwP6tPXYwFEpCAlKiZOae4\nGQskIAEpUVeol9yMBRKQgJSg1U32dfSHSBVjgQQkICXoITXQ1VggAQlI1Ss+wsVj3wIkHZCA\nVK3xqrO7sUACEpCqd46a5G4skIAEpGrNrXOsy7FAAhKQqnW98+d9l48FEpCAVLX1LfZc63Is\nkIAEpKo9rfq7HQskIAGpam1zZrsdCyQgAalKk9XZrscCCUhAqtJF6k3XY4EEJCDFtzjv0CLX\nY4EEJCDFd6ca6n4skIAEpLg2HdhwhfuxQAISkOJ6TfVIYyyQgASkuM5UH6UxFkhAAlJs89w+\nzS4yFkhAAlJsN6rh6YwFEpCAFNPG/fZYnc5YIAEJSDGNUlekNRZIQAJSTH9WU9MaCyQgAamy\nb3LapzcWSEACUmX91LPpjQUSkIAUrWCvPdelNxZIQAJStOdU3zTHAglIQIp2spqR5lggAQlI\nFc1Qf0x3LJCABKSKrlQvpzsWSEACUnlrm+7l5oXz48YCCUhAKu9JdUPaY4EEJCCV1y7327TH\nAglIQIo0XZ2e/lggAQlIkfqoV9IfCyQgASnchuYt0n2oAUg6IAFJ94K61mAskIAEpHB/SPdZ\nDeGxQAISkHTf5JxsMhZIQAKSboAaYTIWSEACUqjCfZu6fEuk+LFAAhKQQo1RlxuNBRKQgBTq\nDPWp0VggAQlIIt/XOcZsLJCABCSRgepRs7FAAhKQpKh1A/fvQBE3FkhAApK8obobjgUSkIAk\nHdUkw7FAAhKQFtVrU2w4FkhAAtIQdZ/pWCABCUi/qfuD6VggASnwkCar84zHAglIgYfUS401\nHgskIAUd0vo9994Yu794mvunrwIJSIGH9Ky6PmZv7ulK1bt2Y9VFdmOBBKSgQ/qDmlm5U/Bb\npbvO7VggASngkObEvbfY6LAjVW+Vy7FAAlLAId2inozZuzsCSU13ORZIQAo2pKID456v+kzE\nUc5il2OBBKRgQxof/3zVpS3DkM51OxZIQAo2pM7q33H7E/YNOWq/xO1YIAEp0JBWNDioyvNV\nV495eKLrp7ACCUjBhvSwutOTsUACUqAhHZ87z5OxQAJSkCFNVx28GQskINm02a6t1nbbNWm0\ndacfU3+x4sb2VaM9GbtltydjqmaV+DK21JepJWX+jP3Zj6m7rC12S7Z4DGm7Xbus3bZr0mhn\niR9Td1ixY39u2fwnb8aWejKmapY/Y8v8mWr/qZJOpTv8mFpi2Y/1GJLtV8kM/tZujLrSo7F8\na8e3dgGGdI7Z66vGjAUSkIILaXG9o7waCyQgBRfSA+per8YCCUjBhfTbugu9GgskIAUW0nR1\npmdjgQSkwELqq17xbCyQgBRUSIX7NF3v2VggASmokMaqy7wbCyQgBRVSR/WBd2OBBKSAQlqW\nd5jhK+fHjgUSkAIK6RE1yMOxQAJSQCGdmPudh2OBBKRgQvpKneblWCABKZiQBqjnvRwLJCAF\nElLRAY3XeDkWSEAKJKTxqoenY4EEpEBCulC96+lYIAEpiJBWVXs5O8OxQAJSECE9pW71diyQ\ngBRESL9TX3k7FkhACiCkubntbVe6GwskIAUQ0p3qMY/HAglIAYR0RN5Sj8cCCUjBg/Sh6uj1\nWCABKXiQrlT/8HoskIAUOEgbWzQv8HoskIAUOEivqis8HwskIAUOkpd/Y14xFkhAChqk5XmH\nevr0oPBYIAEpaJAeV3d4PxZIQAoapJNzvvV+LJCAFDBIS3JO9WEskIAUMEiD1RM+jAUSkAIG\n6Yj6K3wYCyQgBQvSTNXJj7FAAlKwIPVVY/0YCyQgBQpSYYuWG30YCyQgBQvSa6qfH2OBBKRg\nQcpXX/gxFkhAChSklQ3abPNhLJAESIGC9IS6z5/HMIAEpCBB+r1aBiQgAcmwebknW0ACEpAM\nG6QeAxKQBEiGHV1vOZCAJEAya6o650cgAUmAZFY/9QqQBEgCJKOKWjVZDyQBkgDJqLdUDwGS\nAEmAZFR39Q6QdEACkkHrm+63CUg6IAHJoBfV9QIkHZC8gbQ7mJDOUtOAFA5IJpCu+iJyvfGi\nzwMJaUm9owRI4YBkAqlb7o3bLKvs5WYNFwQS0jB1jwApHJBMIJU81uiIL5adrs5Z5cpR1kBq\nl/udACkckMx+Rlp5Zm79/d5wxyhrIH2t/qivgCRAEkNI80/J2aP1pIBCuk09ra+AJEASI0jb\n76rX5vOC81S3TYGE1CbyupBAEiCJEaROdW7ZFrr6e7PmSwII6ePyt40FkgBJjCD1nx25Ljgv\niA9/X63+Hr4GkgBJPPqF7I7gQdq0T7PI28YCSYAkPEUo3carXpENIAmQBEjp1lX9O7IBJAGS\nACnN1jdpVRTZApIASYCUZi/pJ36HA5IASYCUZmfrJ36HA5IASYCUXsvyjqzYBJIASYCUXo+q\nwRWbQBIgCZDS6+ScuRWbQBIgCZDSam7O76LbQBIgCZDS6i41PLoNJAGSACmtjspbGt0GkgBJ\ngJROU9W5lTtAEiAJkNKpnxpVuQMkAZIAKY02tWqyvnIPSAIkAVIa/Ut1j9kDkgBJgJRGPdRb\nMXtAEiAJkNy3Yc+WhTG7QBIgCZDcN0ZdHbsLJAGSAMl9HdVHsbtAEiAJkFy3qsEhxbH7QBIg\nCZBc95QaGLcPJAGSAMl1f1Kz4va/nzSjwIu5VQISkLIa0oI6bWN3Cy5RSh3+gQeDqwQkIGU1\npPvU0Njd/kq379Jky9MOSEDKakjH5c6P2StoGIakhnkwOT4gASmbIX2pTovdXRBxpPqbT64S\nkICUzZBuUc/E7hY04CuSDkhActdhkfdyiXZt2NE+S8wnVwlIQYe0btKM7eWbX4wLNTX+WGZD\n+kDlxx8o6BZydOhk48HVAlLAIX1y4R1X9v8psj3wqltvvfXF+GOZDalP+Xu5xPTdv6fxeyQg\neQ1pZ49J1o7+r0Z2en9V/VhGQ9q0d9MNVY/xzAYBkngOafrFOy1rYu8IoPw11Y5lNqTxqme1\nY0ASIInnkN66OXSxIH+n3l6X/+XIfyyOP5bZkLqpCdWOAUmAJJ5DGjkkdLEmv0hvf5Pf7+UH\nOk2OPfZcz549r91tV4lVarsmjUpMp25pcsDO6lN9OtkyP6butnwa68vUMp/G+jPV/mR3uYD0\nQgRNgd5eNWGHZb3ZZUfMsUHt2rU703ZIbe0NdUtNnwJlcKXRLXtIr+tPtYX5lQ92F+Uvr3rM\n9qtkbf3W7hz1afWDfGsnfGsnnn9rN61riWVN6hXe/m5h6GJTfmHssUyGtDzviARHgSRAEs8h\nbe861Sq5eXToy1iZ9UaX/1jWqL5l0WMZDulxdWeCo0ASIIn3v5Cd2uXRG67bbFkXjLa23XPR\nkAG9llQey3BIv1OzExwFkgBJfHiK0JpJM7aGrt6YZ1ll30+ctiXmWGZD+i73pESHgSRAEp60\n6ry7Ez/HG0gCJAGS846puyjRYSAJkARIjpup/pLwOJAESAIkx92onkt4HEgCJAGS04oParAq\n4Q1AEiAJkJz2vuqc+AYgCZAESE67XI1NfAOQBEgCJIcVtmyW5O9ggSRAEiA5bJy6NMktQBIg\nCZAcdpF6L8ktQBIgCZCctXFw84MAABd9SURBVLbxAUVJbgKSAEmA5KwX1IBkNwFJgCRActZZ\nalqym4AkQBIgOWppvaOT3gYkAZIAyVGPqLuT3gYkAZIAyVHtc+YmvQ1IAiQBkpPm5vwu+Y1A\nEiAJkJx0lxqe/EYgCZAESE46Ki/FW1sCSYAkQHLQVHVuiluBJEASIDmonxqV4lYgCZAESPYV\n7d9kfYqbgSRAEiDZ947qnupmIAmQBEj29VBvp7oZSAIkAZJtG/bcd1Oq24EkQBIg2TZaXZvy\ndiAJkARItv010Xu5xAQkAZIAya7E7+USE5AESAIku55I+F4uMQFJgCRAsuv3Cd/LJSYgCZAE\nSDbNz21vswJIAiQBkk33qIdtVgBJgCRAsunYuottVgBJgCRASt0X6ky7JUASIAmQUne9Gmm3\nBEgCJAFSyooOaLTGbg2QBEgCpJS9o7rZrgGSAEmAlLIearztGiAJkARIqbJ74nc4IAmQBEip\nesXmid/hgCRAEiCl6hybJ36HA5IASYCUomV5v3GwCkgCJAFSih5RgxysApIASYCUovY5cxys\nApIASYCUvJSv+F0ZkARIAqTk3Z7qFb8rA5IASYCUvDapXvG7MiAJkARISftIdXS0DkgCJAFS\n0q5SYxytA5IASYCUrMK9mxU4WggkAZIAKVmvq97OFgJJgCRAStZF6l1nC4EkQBIgJWl1w9bF\nzlYCSYAkQErSCHWTw5VAEiAJkJL0J/Wlw5VAEiAJkBI3v86JTpcCSYAkQErcYDXM6VIgCZAE\nSIk7uu4ip0uBJEASICXsU3WO46lAEiAJkBJ2jRrteCqQBEgCpEQV7tN0veOpQBIgCZAS9brq\n5XwqkARIAqREdVbvOZ8KJAGSAClBqxw/PUgHJAGSAClBT6tbXEwFkgBJgJSgP6hZLqYCSYAk\nQKrevNyT3EwFkgBJgFS9QeoRN1OBJEASIFXvqHqOXj2oIiAJkARI1Zqi/upqKpAESAKkajl9\n9aCKgCRAEiBVbeNezZ29elBFQBIgCZCq9qq6wt1UIAmQBEhVO0995G4qkARIAqQqLcs73OVU\nIAmQBEhVGqbujtlb0H3/Fud+kfojgCRAEiBV6fjc7yp3Vh6kQjWenfIjgCRAEiDFN0N1iNm7\nSYU7N+WHAEmAJECK7zo1MmavQwTSfik/BEgCJAFSXEX7N1kXs3tuBNIhKT8GSAIkAVJcb6qe\nsbtPRSBdk/JjgCRAEiDF1VlNit0tPk87OnZtyo8BkgBJgBTbyoYHxf+NefGoXt0et3nGEJAE\nSAKk2Iar291PBZIASYAUW/ucb91PBZIASYAU09c5v09jKpAESAKkmG5Wz6QxFUgCJAFSZcUH\nNVydxlQgCZAESJVNUF3TmQokAZIAqbJu6u10pgJJgCRAiraq4f6b0pkKJAGS1AikX+zabu2w\nXZNG23elunWEuj2tqVutlGPTbWuJH1N/sfwZW+bL1FKfxv7qx9Td1la7Jb96DOlXu3ZYO23X\npNH2XaluPSVnQVpTt1m70/o4u7Elfkz91Sr1ZWyZL1NLfRq71Y+pIUh2S7Z6DMn2q2RNfGv3\nVc4f0pvKt3bCt3bCz0gVXa+eTW8qkARIAqTyCvdrvCa9qUASIAmQyvunm3e7jAtIAiQBUnkd\n1eQ0pwJJgCRAirQ073AX73YZF5AESAKkSEPVPelOBZIASYAU6di6C9KdCiQBkgAp3BR1dtpT\ngSRAEiCFu8LleyLFBiQBkgBJt6FZC3fviRQbkARIAiTdS6pv+lOBJEASIOk6qM/TnwokAZIA\nKdT3ddoaTAWSAEmAFOoO9bDBVCAJkARIIkUHNVhuMBVIAiQBksj49F70pCIgCZAESCKd1L9N\npgJJgCRAMnm+ajggCZAESHJf+s9XDQckAZIASY6su9BoKpAESAKkSaqj2VQgCZAESN3VG2ZT\ngSRAksBDWpnm66tWBiQBkgQe0qNqoOFUIAmQJPCQjsudazgVSAIkCTqkaeoM06lAEiBJ0CFd\nbvCnseUBSYAkAYe0fs+9N5pOBZIASQIOaYS63ngqkARIEnBIp6iZxlOBJECSYEOalfN786lA\nEiBJsCFdq54znwokAZIEGlLBXnuuM58KJAGSBBrSs+paD6YCSYAkgYbUXs3wYCqQBEgSZEgz\n1B+9mAokAZIEGdKV6mUvpgJJgCQBhrRuz73Sf8HvmIAkQJIAQ3pS3ejJVCAJkCTAkNrmfO3J\nVCAJkCS4kD5Tf/FmKpAESBJcSD3VP7yZCiQBkgQW0uo99jX+A4pIQBIgSWAhPWz8Wg0VAUmA\nJIGF9Ns68zyaCiQBkgQV0vvqXK+mAkmAJEGF1NX0ZSErA5IASQIKaUXDAw1fFrIyIAmQJKCQ\n7lODPJsKJAGSBBNS8WF5izybCiQBkgQT0jizN7uMD0gCJAkmpDPUh95NBZIASQIJ6Zvc4zyc\nCiQBkgQSUl81wsOpQBIgSRAhrWveYoOHU4EkQJIgQhqubvByKpAESBJESMfUmePlVCAJkCSA\nkN5Vf/V0KpAESBJASOertz2dCiQBkgQP0sJ6hxV7OhVIAiQJHqSBapi3U4EkQJLAQdrYqvHK\niv11Yx4Ya/735kASIEngIL2orqjYnXKgUuo3xq/JBSQBkgQO0slqevne+oOVrq3pT0xAEiBJ\n0CB9VvnC+W+oSJ8ZTgWSAEmCBukS9feKvWfLIY03nAokAZIEDNLi+q2jf2I+uRzSXMOpQBIg\nScAgDVQPRPeKO4QddTedCiQBkgQL0uZ99lhZubv0whxV5/K1plOBJECSYEEaqfrGHVj1BW/G\n7E1AChSk4+p84/1UIAmQJFCQ3lUdfZgKJAGSBArSOeo9H6YCSYAkQYI0O7etD1OBpANSgCD1\nUaN9mAokHZCCA2lF4/1+9n4qkMIBKTiQhqi7/fmMBxKQJDiQNh2UtwJIQBIgmfWK6vUzkIAk\nQDKrvZoGJCDpgGTQh+p0ARKQdEAyqKN6E0gCJB2Q0u/bOkcXA0mApANS+l2qnhcgCZB0QEq7\nxQ323wgkHZCAZNBNaqgASQckIKXf2hbNVguQdEACUvo9oG7UV0ACkg5IaVbYOm+hvgYSkHRA\nSrMXVO/wNZCApMsQSOsmzdhesb1y8gcFoasvxoWaWmOQ/i93ZvgaSEDSZQakTy6848r+P0W2\nx15w180XfGZZA6+69dZbX6wpSOPVeZENIAFJlxGQdvaYZO3o/2p4e1n+V5b1arddVu+v4tbY\nnpO3kP6sJkc2gAQkXUZAmn7xTsua2Du8PfHK0EVh/sqd+WtqENLnOaeWbwEJSLqMgPTWzaGL\nBfk79faGRaGL2flF6/K/HPmPxTUFqbMaW74FJCDpMgLSyCGhizX5RRX7hX2GWt/k93v5gU6T\n9e4DHTp06Fxmm2W/xGkr6x5d4t20RFkenm3sWH+mcrI1d7IlLiC9EIFUENkrm9Rt8FZr1YQd\nlvVml9CFNbxTp069SuwqtUpt1zjtajUqOrbMs6kxlVr+jPVlaok/J1ti+TK1zKex3n12xU61\nP9ndLiC9fkvoYmF+5AHwH+/qOaWs/Iai/OUVa2y/Snr4rd3C+gcWVGzzrR3f2uky4lu7aV1D\nX78m9Qpv77xh8BZ9/d3C0MWm/MKagHS1ejS6DSQg6TIC0vauU62Sm0dbVmmZ9UmXreFjb3T5\nj2WN6lvxtel/CWlp433WR3eABCRdRkCypnZ59IbrNlvWBaOtJ7rfpSvads9FQwb0WhJdYntO\n3kG6Qd1fuQMkIOkyA5K1ZtIM/YXojXnW1HHhNltl30+ctqVyhe05eQZpZdMWayr3gAQkXYZA\nss/2nDyDdJsaFLMHJCDpgOS2tS2arojZBRKQdEBy293qlthdIAFJBySXbdi30ZLYfSABSQck\nlz2k+sftAwlIOiC5a2Pr+gviDgAJSDoguWu46hN/AEhA0gHJVYWH1JsXfwRIQNIByVVPq79V\nOQIkIOmA5KaNB9ebU+UQkICkA5KbHleXVT0EJCDpgOSigtZ531U9BiQg6YDkomFVH7ITIOmA\nBCQ3bWhVf361g0ACkg5IzrtfXVv9IJCApAOS49bu3WhR9aNAApIOSI67W92Q4CiQgKQDktPW\ntGy8JMFhIAFJBySn3Rn/d0gVAQlIOiA5bGXzpssSHQcSkHRActhAdXvC40ACkg5IzlratPmq\nhDcACUg6IDmrrxqS+AYgAUkHJEfNzWu1Tl8Xz11e5RYgAUkHJEddrEboqyf3UuqU6XG3AAlI\nOiA5aVrubwpDV68oXaulsTcBCUg6IDmpg/qnvjo8DEndGXsTkICkA5KD/q3C7xhblBuB1CX2\nNiABSQck+4pPKn8L8xYRSNfE3ggkIOmAZN9LqmNk49qwo/qfxN4IJCDpgGTbxkPrzIxsbeig\nHT0adyuQgKQDkm3DYl7xZOK9T8yNvxVIQNIBya41ezf4PsXNQAKSDkh2DVQ3p7oZSEDSAcmm\nBY1brEx1O5CApAOSTd3UIylvBxKQdEBK3ae5RxamXAAkIOmAlLLiU9X41CuABCQdkFL2ojrH\nZgWQgKQDUqrWt877ymYJkICkA1KqblP97JYACUg6IKVofqO9VtitARKQdEBK0cXqcds1QAKS\nDkjJ+yDnqNQPfeuABCQdkJJWfLJ6y34VkICkA1LSnlfnOVgFJCDpgJSs1a3yvnGwDEhA0gEp\nWX3VjU6WAQlIOiAlaXq9A9c6WQckIOmAlLii9mqso4VAApIOSIl7Qp3rbCGQgKQDUsKWtWg4\nx9lKIAFJB6SE9VCDHa4EEpB0QErU+zmHFzhcCiQg6YCUoMLfqnecrgUSkHRAStB9qqvTpUAS\nIOmAVL2FTZssdLgUSDogASlRHW1eOCguIAFJB6RqjVEnFjlbqQMSkHRAqtqyffKm26+KBiQg\n6YBUta7qDkfrygMSkHRAqtIb6hinv0IKByQg6YAU36oD605xcg+jAQlIOiDFd5ka4OQOVgYk\nIOmAFNe/cw5f7+geRgMSkHRAim39YbnvObuH0YAEJB2QYuun+ji7g5UBCUg6IMX0UZ3Wqx3e\nw2hAApIOSJWtPzLH5j1cEgQkIOmAVFkfdanTO1gZkICkA1K0t3IOWeP4HkYDEpB0QKpo+QF1\nP3R+D6MBCUg6IFV0vuo9sO8L9q+aXyUgAUkHpPJGqP1UqP9b5eJu6oAEJB2QIs1p0kCF+5ub\n+ylA0gEJSOUV/V79MQJpD1d3FEg6IAGpvEHqnG4RSLmbXN1TIAmQdEDSTc1rufi+CKRj3N1T\nIAmQdEAKtfpwNVbWHh6G9LbLuwokIOmAFKqLuip0Oa9TozrHjHN3R4GkAxKQdI+ptpG/Li/a\n4OpehgMSkHRAkukN9vzW9T2MBiQg6YC05gg1xv09jAYkIOmAdLG62v0drAxIQNIFHtKj6gRX\nL79VNSABSRd0SNMa7OnwrfmSBCQg6QIOaU2bHJMfkARIOiAFHVLxeeqa9O5hNCABSRdsSLer\nkzemdw+jAQlIukBDejV33wVp3sNoQAKSLsiQvmxa/6N072E0IAFJlzWQttm109oVf2BjG/WC\n7UfZj91tPqN6O6wSP8ZuL/Vj6jbLn7Flvkwttf9USWvsdj+mllj2Yz2GtMWubdaOuP3NZ6vr\nbD/Ivm27PBhSrV8tf8aW+DF1i+XP2DJfppb6M7bkFz+m7rZ+tVvyi8eQbL9KVv3W7jr1Z9ev\ndJIgvrXjWztd1nxrZ3tOVSCNyWm91OAeRgMSkHRBhfRJw8Zu3ik2eUACki6gkObua/qMhoqA\nBCRdMCGtOkYNNruH0YAEJF0gIW08TfU0vIfRgAQkXRAhFXdTZ3jxgF04IAFJF0RIN6vj0njb\niSQBCUi6AEJ6WrX63vgeRgMSkHTBgzQhr8k083sYDUhA0gUO0mdN8t7x4B5GAxKQdEGDNKtl\nznNe3MNoQAKSLmCQvmut7vPkHkYDEpB0wYK0uI2605M7WBmQgKQLDqRZE2cs+63Za9glCkhA\n0gUF0g9/Uko1UpcUe3QPowEJSLqgQOoQfs+W5mm8Sr5NQAKSLiCQZkTeRUx5+sh3OCABSRcQ\nSOPLIY3w6h5GAxKQdAGB9GU5pIle3cNoQAKSLiCQ5IywoxM9e9J3NCABSRcUSIv/EnJ08lyP\n7l9MQAKSLiiQROa8N9vzx74FSDogBQlSijdjNglIQNIByTAgAUkHJMOABCQdkAwDEpB0QDIM\nSEDSAckwIAFJByTDgAQkHZAMAxKQdEAyDEhA0gHJMCABSQckw4AEJB2QDAMSkHRAMgxIQNIB\nyTAgAUkHJMOABCQdkAwDEpB0QDIMSEDSAckwIAFJByTDgAQkHZAMAxKQdEAyDEhA0gHJMCAB\nSQckw4AEJB2QDAMSkHRAMgxIQNIByTAgAUkHJMOABCQdkAwDEpB0QDIMSEDSAckwIAFJlzWQ\nNtu1/stVtmvS6Jdtfkwt/nKRH2O37PBj6uYvv/Nl7C5fps6Z5cvYnT/7MXXhl8V2S7Z4DMm2\nj9r983/zH/Kile3ur+lTcNFJl9b0Gbio659q+gxcdGe7AueLgVQ9IPkWkAwDkm8Bya+AZBiQ\nfAtIhm2asvZ/8x/yol+n/FDTp+CiKbNr+gxcNGtqTZ+Bi+ZP2eZ88f8IElF2ByQiDwISkQf5\nDWnNnd0GLk2wXSsrHtqz91M/R7bvyw81uGbPJ2X/0id4QWS7bOyVl48qqdnzSdXM/HBPhXdi\nT7w29slyfTljwCUP/VR+JHY7WT5D2tb9yUXPX7S52natbOulgxd+c809kZ1+Y+bMmbO8Zk8o\nZc8PCZ3g3Mj2Gz2+/Lr3qJo9n1T9FDrVOV/3+Cy8E3vitbAtvWaGLhd2emfBHTdFjsRuJ81n\nSO9dXmqV9X2j2natbOrF2yxrfv5/9HbZRQtr+nRsGvJ6dLOk94eWNa3b9ho8Gwe9U/5PVMyJ\n174Kn+6ZryE99Jhl/adT5NHb2O2k+QzpoadDFy8NrrZdK/v48dDF2vwVevvHfNmyo4bPJ3V9\nP9te8UyvNfliWb/k1+4H7Yt7bIpsxJx47at4woTOGlIP/Tj99ZF/9WO3k+YzpIH697D/6l9t\nu/b2yt926qtF+bfkdxryn5o+m+SVdb6lU36/ReHtuZ1KQ5cXz6jZM7LpqZGR69gTr5V1CUHa\nlT8/tHVf+JRjt5PnM6S+74QuJvestl1b2z6yU+TTcVrXKVtXDrijhk8nRf+56JWf5NEe4R85\np3XRl70m1+wZpa6gS/m/SrEnXivTkH7M1z8eP/qI3o/dTp7PkG4ZF7r4V99q27W076+4ekHM\n7sL8oho7FUft6PKpvvq2U1no8uJpNXw2KXvu0di98hOvlWlIOyNfhZ7V+7HbyfMZ0lD9nx91\nZ7Xt2tn7nd/YHbu/Ob9WP1ofqt+/9OXK/B8ta1v4/+/a2s7u38TtR068VqYhWd30v0o3jgsf\niN1Oms+Q3r0q9G/lgDeqbdfK1pz/VXR71NDQxbzzt9bc2dg0o9/PlrX14vCnZ0nP0L/vX3at\nzY/azexW8Wuu2BOvlYUh6UfGtnSKPHIbu500nyFt6TGmcFxXsYrHbohu19rGXD5/Qagd1j/n\nWPPPf2nJrD4v1/QpJW9LzyHzFg6+qdT6JPSj0bjLlyy76pWaPqVUPRt57Dt0stETr62FIX3X\n+fMNQ2+0rAVjd0e3U+b3MxtW3xF+NsOi/G+i27W2+yO/f19rXTA69L/gHd36vlmLnyxgFT/U\n89Knt1jW4Jssq+zVKy9/pfZ+aoa6OvKdkT7ZihOvrYUhWdMHXPLgZst6M397dDtlPNeOyIOA\nRORBQCLyICAReRCQiDwISEQeBCQiDwISkQcBiciDgJTpHXNPTZ8BWUDK/Pa4rqbPgCwgZXqb\nP2/Y+fPa/LTvoASkzO5zFWpFTZ8FASnj41u7WhGQMj0g1YqAlOkBqVYEpEwPSLUiIGV6QKoV\nASnTA1KtCEiZ3h5X1/QZkAWkzO/YlpdvqulzICBlfDN6nV/LXw82EAGJyIOARORBQCLyICAR\neRCQiDwISEQeBCQiDwISkQcBiciDgETkQUAi8iAgEXkQkIg8CEhEHvT/Ab7Ci8a05mPcAAAA\nAElFTkSuQmCC",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 420,
       "width": 420
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# time intervals: a sequence from zero to ten at 0.5 steps\n",
    "time <- seq(0, 10, by = 0.5)\n",
    "# initial condition\n",
    "x0 <- 0.1\n",
    "## The function to be integrated (right-hand expression of the derivative above)\n",
    "f <- function(x){x * (1.-x)}\n",
    "\n",
    "## An empty R vector to store the results\n",
    "x <- c()\n",
    "## Store the initial condition in the first position of the vector\n",
    "x[1] <- x0\n",
    "\n",
    "# loop over time: approximate the function at each time step\n",
    "for (i in 1:(length(time)-1)){\n",
    "    x[i+1] = x[i] + 0.5 * f(x[i])\n",
    "}\n",
    "\n",
    "## plotting with ggplot2\n",
    "library(ggplot2)#load each library once per R session\n",
    "p <- ggplot(data = data.frame(x = x, t = time), aes(t, x)) + geom_point()\n",
    "analytic <- stat_function(fun=function(t){0.1 * exp(t)/(1+0.1*(exp(t)-1.))})\n",
    "print(p+analytic)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Why use scientific libraries?\n",
    "\n",
    "The method we just used above is called the *Euler method*, and is the simplest one available. The problem is that, although it works reasonably well for the differential equation above, in many cases it doesn't perform very well. There are many ways to improve it: in fact, there are many books entirely dedicated to this. Although many math or physics students do learn how to implement more sophisticated methods, the topic is really deep. Luckily, we can rely on the expertise of lots of people to come up with good algorithms that work well in most situations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Then, how... ?\n",
    "\n",
    "We are going to demonstrate how to use scientific libraries to integrate differential equations. Although the specific commands depend on the software, the general procedure is usually the same:\n",
    "\n",
    "* define the derivative function (the right hand side of the differential equation)\n",
    "* choose a time step or a sequence of times where you want the solution\n",
    "* provide the parameters and the initial condition\n",
    "* pass the function, time sequence, parameters and initial conditions to a computer routine that runs the integration.\n",
    "\n",
    "### A single equation\n",
    "\n",
    "So, let's start with the same equation as above, the logistic equation, now with any parameters for growth rate and carrying capacity:\n",
    "\n",
    "$$ \\frac{dx}{dt} = f(x) = r x \\left(1 - \\frac{x}{K} \\right) $$\n",
    "\n",
    "with $r=2$, $K=10$ and $x(0) = 0.1$. We show how to integrate it using the desolve package below, introducing key language syntax as necessary. We'll use the [deSolve](http://desolve.r-forge.r-project.org/), the most popular R library to run  numerical integration. The basic workflow recommended by the authors of the package is\n",
    "\n",
    "#### 1. Declare initial conditions, parameters, time interval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "library(deSolve)# loads the library\n",
    "\n",
    "## time sequence\n",
    "time <- seq(from=0, to=10, by = 0.01)\n",
    "# parameters: a named vector\n",
    "parameters <- c(r = 1.5, K = 10)\n",
    "\n",
    "# initial conditions: also a named vector\n",
    "state <- c(x = 0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2. Define a function in R for the ODE to be integrated\n",
    "\n",
    "Let's define the right-hand side of the differential equation. \n",
    "To be recognized by the integration routines of *deSolve*\n",
    "it must be an R function that computes the values\n",
    "of the derivative on a time $t$. \n",
    "There are many ways to do this, but the recommended format is:\n",
    "* Make a function with three arguments: time sequence, state variables and parameters, in this order.\n",
    "* The function should return a list with results of the function to be integrated. \n",
    "To do this use `with(as.list(c(state, parameters){ ... }` inside the R function.\n",
    "Include between brackets the function(s) to be integrated\n",
    "and then close returning the list of the calculated values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "## The logistic ODE to be integrated\n",
    "logistic <- function(t, state, parameters){\n",
    "    with(\n",
    "        as.list(c(state, parameters)),{\n",
    "            dx <- r*x*(1-x/K)\n",
    "            return(list(dx))\n",
    "        }\n",
    "        )\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3. Integrate the function\n",
    "\n",
    "Now call the R function `ode`, to perform the integration the basic arguments of 'ode' are\n",
    "* `y`: the vector of initial conditions\n",
    "* `times`: the vector with the time sequence\n",
    "* `func`: the R function as described above\n",
    "* `parms`: vector of parameter values (named)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "out <- ode(y = state, times = time, func = logistic, parms = parameters)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The resulting object has the values of the integration\n",
    "at each time point in the vector of times"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "head(out) # first 6 lines"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which you can plot with"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "#### For jupyter notebook only\n",
    "options(jupyter.plot_mimetypes = 'image/png', repr.plot.height=5)\n",
    "####\n",
    "\n",
    "plot(out, lwd=6, col=\"lightblue\", main=\"\", ylab=\"N(t)\")\n",
    "curve(0.1*10*exp(1.5*x)/(10+0.1*(exp(1.5*x)-1)), add=TRUE)\n",
    "legend(\"topleft\", c(\"Numerical\", \"Analytical\"), lty=1, col=c(\"lightblue\", \"black\"), lwd=c(6,1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We get a much better approximation now, the two curves superimpose each other!\n",
    "The get the same plot with *ggplot2* package use:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "## Ploting with ggplot2\n",
    "p <- ggplot(data = as.data.frame(out), aes(time, x)) + geom_point()\n",
    "analytic <- stat_function(fun=function(t){0.1*10*exp(1.5*t)/(10+0.1*(exp(1.5*t)-1))})\n",
    "print(p+analytic)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### A system of equations\n",
    "Now, what if we wanted to integrate a system of differential equations? Let's take the Lotka-Volterra equations:\n",
    "\n",
    "$$ \\begin{aligned}\n",
    "\\frac{dV}{dt} &= r V - c V P\\\\\n",
    "\\frac{dP}{dt} &= ec V P - dP\n",
    "\\end{aligned}$$\n",
    "\n",
    "All that you have to do is to write an R function that returns the values of both derivatives at each time, and to define the values of the parameters and of the initial conditions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# time sequence \n",
    "time <- seq(0, 50, by = 0.01)\n",
    "\n",
    "# parameters: a named vector\n",
    "parameters <- c(r = 2, k = 0.5, e = 0.1, d = 1)\n",
    "\n",
    "# initial condition: a named vector\n",
    "state <- c(V = 1, P = 3)\n",
    "\n",
    "# R function to calculate the value of the derivatives at each time value\n",
    "# Use the names of the variables as defined in the vectors above\n",
    "lotkaVolterra <- function(t, state, parameters){\n",
    "  with(as.list(c(state, parameters)), {\n",
    "    dV = r * V - k * V * P\n",
    "    dP = e * k * V * P - d * P\n",
    "    return(list(c(dV, dP)))\n",
    "  })\n",
    "}\n",
    "## Integration with 'ode'\n",
    "out <- ode(y = state, times = time, func = lotkaVolterra, parms = parameters)\n",
    "\n",
    "## Ploting\n",
    "out.df = as.data.frame(out) # required by ggplot: data object must be a data frame\n",
    "library(reshape2)\n",
    "out.m = melt(out.df, id.vars='time') # this makes plotting easier by puting all variables in a single column\n",
    "\n",
    "p <- ggplot(out.m, aes(time, value, color = variable)) + geom_point()\n",
    "print(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An interesting thing to do here is take a look at the *phase space*, that is, plot only the dependent variables, without respect to time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "p2 <- ggplot(data = out.df[1:567,], aes(x = P, V, color = time)) + geom_point()\n",
    "print(p2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Congratulations**: you are now ready to integrate any system of differential equations!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### For more info:\n",
    "\n",
    "* [Introduction to R](cran.r-project.org/doc/manuals/R-intro.html)\n",
    "* [Crash course in R](http://www.r-bloggers.com/a-crash-course-in-r/)\n",
    "* [ode package](http://desolve.r-forge.r-project.org/)\n",
    "* [Some additional example codes](https://github.com/diogro/ode_examples)\n",
    "* [ggplot2 package](http://ggplot2.org/)\n",
    "* [A R graph gallery](http://www.sr.bham.ac.uk/~ajrs/R/r-gallery.html)\n",
    "* [Another tutorial on numerical integration in R](http://www.r-bloggers.com/learning-r-parameter-fitting-for-models-involving-differential-equations/)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.2.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}