{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Notebook for examples of intergation of simple ODEs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Python initialization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np #initialization\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Defining growth functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Malthus(Population):\n",
    "    \"\"\"\n",
    "    this function returns the growth rate of a simple exponential growth\n",
    "\n",
    "    Malthus(Population)\n",
    "\n",
    "        Population -- current population size\n",
    "    \"\"\"\n",
    "    GrowthRate = 7.0\n",
    "    return GrowthRate*Population"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def MalthusForODE(Population, t, GrowthRate=1.0):\n",
    "    \"\"\"\n",
    "    this function returns the growth rate of a simple exponential growth\n",
    "\n",
    "    Malthus(Population,t,GrowthRate)hel\n",
    "\n",
    "        Population -- current population size\n",
    "        t -- time, unused\n",
    "        GrowthRate -- population Growth Rate\n",
    "    \"\"\"\n",
    "    return GrowthRate*Population"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def MalthusCapacity(Population):\n",
    "    \"\"\"\n",
    "    this function returns the growth rate of a simple exponential growth\n",
    "    with carrying capacity\n",
    "\n",
    "    Malthus(Population)\n",
    "    \"\"\"\n",
    "    GrowthRate = 1e-9\n",
    "    CarryingCapacity = 10\n",
    "    return GrowthRate*Population*(CarryingCapacity-Population)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def MalthusCapacityParams(Population, GrowthRate, CarryingCapacity):\n",
    "    \"\"\"\n",
    "    this function returns the growth rate of a simple exponential growth\n",
    "    with carrying capacity\n",
    "\n",
    "    Malthus(Population)\n",
    "    \"\"\"\n",
    "    return GrowthRate*Population*(CarryingCapacity-Population)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Integrator functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Euler(xPrime, t0=0.0, x0=0.0, T=1.0, dt=0.1):\n",
    "    \"\"\"\n",
    "    Solves 1-d ODE using the Euler method.\n",
    "    Euler(xPrime,t0=0.0,x0=0.0,T=1.0,dt=0.1):\n",
    "    \"\"\"\n",
    "    t = np.arange(t0, T+dt, dt)\n",
    "    x = np.zeros(t.size)\n",
    "    x[0] = x0\n",
    "    for i in range(1, t.size):\n",
    "        x[i] = x[i-1] + dt * xPrime(x[i-1])\n",
    "\n",
    "    return (t, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def RK2(xPrime, t0=0.0, x0=0.0, T=1.0, dt=0.1):\n",
    "    \"\"\"\n",
    "    Solves 1-d ODE using the Runge-Kutta 2 method.\n",
    "\n",
    "    RK2(xPrime,t0=0.0,x0=0.0,T=1.0,dt=0.1):\n",
    "    \"\"\"\n",
    "\n",
    "    t = np.arange(t0, T+dt, dt)\n",
    "    x = np.zeros(t.size)\n",
    "    x[0] = x0\n",
    "    for i in range(1, t.size):\n",
    "        fx = xPrime(x[i-1])\n",
    "        guess = x[i-1] + dt * fx\n",
    "        fxdx = xPrime(guess)\n",
    "        x[i] = x[i-1]+0.5*(fx+fxdx)*dt\n",
    "\n",
    "    return (t, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def EulerArg(xPrime, t0=0.0, x0=0.0, T=1.0, dt=0.1, arg=()):\n",
    "    \"\"\"\n",
    "    Solves 1-d ODE using the Euler method.\n",
    "\n",
    "    EulerArg(xPrime,t0=0.0,x0=0.0,T=1.0,dt=0.1),arg=():\n",
    "    \"\"\"\n",
    "    t = np.arange(0, T+dt, dt)\n",
    "    x = np.zeros(t.size)\n",
    "    x[0] = x0\n",
    "    for i in range(1, t.size):\n",
    "        x[i] = x[i-1] + dt * xPrime(x[i-1], t[i-1], arg)\n",
    "\n",
    "    return (t, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def RK2Arg(xPrime, t0=0.0, x0=0.0, T=1.0, dt=0.1, arg=()):\n",
    "    \"\"\"\n",
    "    Solves 1-d ODE using the Runge-Kutta 2 method.\n",
    "\n",
    "    RK2Arg(xPrime,t0=0.0,x0=0.0,T=1.0,dt=0.1,arg=()):\n",
    "    \"\"\"\n",
    "\n",
    "    t = np.arange(0, T+dt, dt)\n",
    "    x = np.zeros(t.size)\n",
    "    x[0] = x0\n",
    "    for i in range(1, t.size):\n",
    "        fx = xPrime(x[i-1], t[i-1], arg)\n",
    "        guess = x[i-1] + dt * fx\n",
    "        fxdx = xPrime(guess, t[i-1], arg)\n",
    "        x[i] = x[i-1]+0.5*(fx+fxdx)*dt\n",
    "\n",
    "    return (t, x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Constant intialization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "SimulationTime = 2.0 # time to solve for\n",
    "P0 = 10.0  # initial population size\n",
    "GrowthRate = 7.0 # max growth rate of the population"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'absolute error')"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEaCAYAAAAG87ApAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd8VGX2+PHPIYQOoRN6782CYAMBRYpd197X1S1iWXVX17K2dXWL69r2Z1n9quuuruu6ChqKqBQRFVCB0BFRAimEhABJIO38/nhucIgpM8nM3EnmvF+vvGDuzNx7nrkz99z2PEdUFWOMMSZYjfwOwBhjTP1iicMYY0xILHEYY4wJiSUOY4wxIbHEYYwxJiSWOIwxxoTEEkeMEZETRGSziOwXkbODeH0fEVERaRyN+PwkImtFZGId53GfiLzqx7KrmG+d1p+I3Ckifw93XN68fyci2SKSEYn5+01EtonIKX7HUR/FVeLwviiF3ka5/O8pv+Oq4AHgKVVtpapvV3wynr/sqjpcVRfG27LLichEEUkLnKaqv1fVn0RgWT2BW4FhqpocpnmKiMwUkdUiUiAiGSKyUEQuCsf8a1j2SyLyuxBeP1FEyrxtxD4R2SgiV4fw/lrtoAQx3zsrbL8KvTg7VvH6PiLykfd5bwjXtiOuEofnDG+jXP43s7IXVbYHGOpeYS33InsDa2vxvgYrHo6mYlBvYLeqZoX6xmrW1xPAzbiE1AHoDtwNTKtiPiIifm6jdqpqK6AN8EvgeREZ7GM85TsKh7ZfwB+AhaqaXcVbXgO+xH3edwFvikincAQSN3/ANuCUKp67ClgKPAbkAL+rYloj3Jf9WyALeAVI8ubRB1DgGuA7YHEVy7oW2OLNcxbQzZv+NVAGFAL7gaYV3vePCs//OmCZV3rLzAbuCnhPI+AOb967gTeA9lXENRFIw/2ws4B04OqA5xcCP6nwmX0c8FiBXwCbgX3Ag0B/YBmw11t2k4DXnw58BewBPgFGVVhXtwOrgYNA48D1ByQAd3rt2gesBHp6zz0ObPeWuRIYHzDf+4BXq2h/R+BdL54cYAnQqOJ3x5vHf4BXvWWvAQYBv/E+t+3AqVV97wJjCFh/jb3HVwPrvfluBX7qTW/prfcyb93vB7pVbA9wJm7HY4+3voZWiOM27zPNA/4NNKvkczilwrJeCnLeh62vCvMcBJQCY2r4jS4EHsL97gqBAV47Z3nrZAtwrffaZt5rOnqP7wZKgDbe498BfwWuA4qBIq89s2v6PPB+CxViywLOD3hc6fcMlwiLvGXuB1Z505OAF3C/qx1efAl12J4J7vt/ZRXPD/LWReuAaUuAn9V5W1rXGdSnP2pOHCXADbiNVPMqpv3Y+/L2A1oBbwH/8ObRB7cReAX3Q29eyXIm4zbuRwFNgScJSDDVxVjZ8wHLfN6Lb7T3ZRnqPX8z8CnQw1ves8BrVcx7otfeB4BEYAZQALTznl9IzYljFm4PbbgXxwfeZ5UErCv/knvtzwLG4ZLAlV7bmga08yugZ/nnyOEb71/hNtiDvR/QaKCD99xluD2sxrgkmMH3G4T7qDpxPAw847U9ERgPSCXLvg84AEz1lvEK8A1ujy4Rt2PwTTXr7FAM/DBxnIZLtgKc5H3+RwWsn4obs8B5DQLygSleHL/GfVebBMTxOW5D3B6XoCrdiFRcVpDzPmx9VZjfz4BtQfxGF+J2gIZ7n20isAj4Gy5RHAHsAk72Xr8YOM/7/3zchnR6wHPneP9/CfhdJb+lSj+PwPbjdr7OxCXSIwPeH9L3DHgb9/trCXT2ll2+Y3AiLiFX9XdiJZ/VBFxialXFZ3kOsL7CtKeAJ+u6LY3HU1Vvi8iegL9rA57bqapPqmqJqhZWMe1S4C+qulVV9+P2Mi+qcHh+n6rmB8wj0KXAi6r6haoe9N5/nIj0qWO77lfVQlVdBazCbUgBfoo7Aknzlncf8KNqTicUAw+oarGqpuC+mKEcnv9BVfeq6logFZjvfVZ5wBzgSO911wLPqupnqlqqqi/jEs2xAfN6QlW3V/E5/gS4W1U3qrNKVXcDqOqrqrrbW2eP4hJmMG0oBroCvb32L1Hv11aJJao6T1VLcEcfnYBHVLUYeB3oIyJtg1jmYVT1PVX92mvTItzGcHyQb78QeE9V3/fi+DNuZ+L4gNc8oao7VTUHmI3bEIdz3lWtr464DeshIpLm/QYPiEjvgKdeUtW13mebjNuo3q6qB1T1K+DvwOXeaxcBJ3nf51G402EniUgz4BjcHnZ1qvs8uonIHtxRzf+AW1T1y/InQ/meiUgXYDpws7dtyMKdybjIm9fHqtq2mr+PK5ntlcCb3naoMq1wR1KB8oDW1X8kNYvHxHF2hRXyfMBz2yt5fcVp3XCnqcp9i9vj6FLDfCp9v7fSd+PO99ZF4I+yAPelAXeu+n/liRK3V1VaId5Au70fbGXzCkZmwP8LK3kcGNetgUkct7faLeD11X2OPXF7lz8gIreKyHoRyfPmm4TbcNXkT7i96PkislVE7qjmtRXbla2qpQGPIbTPrTz26SLyqYjkeLHPCDJ2+OF3qwz3GQZ+t6r6noRj3tWtr924pHyIqvbAta0p7girsvl0A3JUdV/AtG8DlrsId3RwFO4I9H3ckdqxwBat+tx/ueo+j52q2hZ3BP0E7mzBISF+z3rjjp7SA77vz+KOPEImIs2B84GXq3nZfi/2QG1wp0HrJB4TR3Uq27usOG0n7ktQrhfu9E7ghqS6IYcPe7+ItMQd7u6oQ4zV2Y47dA9Mls1UNdjlBcoHWgQ8rsvdNtuBhyrE1UJVXwt4TXVt3Y47pXMYERmPO9d+Ae4UW1vcXpZUfG1FqrpPVW9V1X7AGcAtInJyCG2qSlCfm4g0Bf6L25vv4sWewvex17TuK363BJdga7OuazPv6uL7EOghImOCWFbgfHYC7UUkcC+5V8ByP8Ht5Z8DLFLVdd7zp+GSSjCxVR+MO1K/HRhZfot8EN+zisvbjjui7hjwfW+jqsPL51fhbqmKfxWPOs/FXfNZWE3oa4F+FT670YTh5htLHKF7DfiliPQVkVbA74F/V9hLr86/gKtF5AhvQ/F74DNV3Rbk+zNx1wyC9QzwUPmpABHpJCJnhfD+QF8B54pICxEZgLsJoLaeB34mIuO8u2daishpFb7k1fk78KCIDPTeP0pEOuAOw0tw58Ebi8hv+eFeV6VE5HQRGeBtFPfijsxKa3hbML7Cnc5M9DacP6ridU1we9+7gBIRmQ6cGvB8JtBBRJKqeP8bwGkicrKIJOLOux/EbVzrqk7zVtWNuD3s10Vkiog0F5EEDj/VVdn7tnvLeFhEmonIKNz37p/e8wW4C9PX832i+AR3ijYwcYT6u6kYRxHwKPBbb1JN37NM3OnKRt7703GnHR8VkTYi0khE+ovISd7zS/Twuz0r/lU85XYl8Eo1p1JR1U2479693md3Du503n9r+zmUi8fEMbtCJv9fiO9/EXd302LcBdEDuIvnQVHVD4B7cCsvHbfXHMp97A8Dd3uHu7cF8frHcRes54vIPtyF8nEhLC/QY7i7RTJxh8j/rOV8UNUVuOscTwG5uFNEV4Uwi7/gNmbzcRv5F3Dn3OfhrqVswp3SOED1p1ACDQQW4A7xlwF/0/D03bgHt55zgftxOw8/4J2OuRHXrlzgEty6K39+A27HZau3/rtVeP9G3AXbJ3E3YJyBu/28qK4NCNO8r8ed8vkLbm85DXfn3YW4C+JVuRh3E8FO3LWGe1X1/YDnF+FOA30e8Lg17jda7gVgmPe5/aB/VJBeBHqJyBnU/D37j/fvbhH5wvv/Fbidg3W49fsmFU7fBUNEuuNOm71SyXPPiMgzAZMuAsZ4y3sE+JGq7gp1mT9YTjUJyxhjjPmBeDziMMYYUweWOIwxxoTEEocxxpiQWOIwxhgTEkscxhhjQtIgRx3t2LGj9unTx+8wjDGm3li5cmW2qgY1cm6DTBx9+vRhxYoVfodhjDH1hoh8W/OrHDtVZYwxJiSWOIwxxoTEEocxxpiQWOIwxhgTEkscxhhjQmKJwxhjTEga5O24xhgTDw4Ul5JbUERufjG5BUUUlZYxaXCtigqGxBKHMcb4TFXJLyolN7/IJYKC4kr/v6egmJz8IvZ40wuLD68z1r5lE764Z0rE47XEYYwxYVRWpuw94DbwuQXF7Cko8jb2xV4i+P4IYU9BMTkFLhEUl1ZeG0kEkpon0q5FE9q1SKRrUjOGdm1D+5aJtG3R5ND0di2b0L5lk6i00RKHMcaEqKS0jE2Z+1mdtodVaXlsztxHTkERuflF5BUWU1ZFfbzGjcTb2LtE0KdjC45s0Za2LZoclggC/5/UPJGERlL5DH1iicMYY6qhqmzbXeCSxPY8VqftIXVnHgeKywBo06wxQ7q2YWhyG9q1dAnhB4mgRRPatkykddPGuJL29ZslDmOMCZCRd4BVaXtYnbaH1Wl5rE7LI6+wGICmjRsxonsSF4/txegebRnVI4k+HVrSKMaOCCLNEocxJm7tKSjykoM75bQ6bQ+Zew8CkNBIGNylNTNGJjOqR1tG92jLoC6taJwQg70YSopg1WuQuRZm/DHii7PEYYyJCwVFJazduZdV2/ccShbbdhccer5fx5Yc16+DSxI92zK8WxuaJSb4GHEQigvhi3/A0sdhbxp0O9JNS2we0cVa4jDGNDjFpWVszNjnTjltz2NV2h42Ze47dNG6a1IzRvVI4vwxPRndoy0jeySR1DzR36BDcXAfrHgRPnkK8rOg57Fw5uPQ/2R3G1aEWeIwxtRrZWXK1uz8Q9ckVqXtYe3OvRSVuIvXbVskMqpHW6YM6+Kdckqic5tmPkddS4V74PPn4NO/QWEu9JsIE/4Pep8QlYRRzhKHMaZeUVW+3rWfD9ZnsXjzLlZvz2PfwRIAmicmMLJ7Elcc25tRPdtyRI+29GzfvP7fyZSfDcuehs+fh6J9MGg6TLgNeozxJRxLHMaYmHewpJTPtubw4YYsPtyQxXc57trEkOTWnHlEN0Z71yUGdG4Vc30e6mRvOnzyJKz8P3ftYvjZMP5WSB7pa1iWOIwxMSlr7wE+2ugSxZLN2RQUldK0cSNOGNCR6yb0Y9KQznRvG9mLwL7J/RaW/hW+fBXKSmHUBXDiLdBpkN+RAZY4jDExoqxMSd2Zxwfrs/hoYxar0/IAdyH7nCO7c/LQzhzXryPNm8T4nU51kb0FPv4LrHodGiXAEZfCCTdB+75+R3YYSxzGGN/kHyxhyeZsPtqQxYcbs9i17yAicGTPtvxq6mAmD+nMkOTW9f8aRU0yUmHJo7D2f9C4GYy9Do6/AZK6+x1ZpSxxGGOi6rvdBXy4IZMPNmTx2dYcikrLaN20MRMGd+LkIZ05aVAnOrRq6neY0bFjJSx+FDa+B01awYk3w7HXQ6tOfkdWLUscxpiIKiktY+W3uXy4IYsPNmSxJWs/AP06teTK43szeUgXxvRpR2Is9siOlG8/gcV/gq8/hGZtYeJv3FFGi/Z+RxYUSxzGmLDLzS9i0aZdfLAhi0Ubs9h7oITEBGFc3w5cMrYXk4d0pk/Hln6HGV2qLlEseRS+XQotO8Ep98Mx10DT1n5HFxJLHMaYOlNVNmXu54MNmXy4PosvvsulTKFjqyZMHZ7MyUM7c8KAjrRuVo96Z4dLWRlsmuuOMHZ+AW26w/Q/wpGXQ5MWfkdXKzGfOESkEfAg0AZYoaov+xySMQZXtnTZ1t18uN7dMrtjTyEAI7q3YebkgZw8pDMjuyfF3cixh5SVwrq33TWMrLXQtjec8TiMvhga1+9rOL4kDhF5ETgdyFLVEQHTpwGPAwnA31X1EeAsoDuQA6T5EK4xxpOTX8SC9ZnMX5vJ0i3ZFBaX0jwxgRMHduSGyQOYNKQzXerrcB7hUloMa/7jTknt3gIdB8E5z8GI8yAh5vfVg+JXK14CngJeKZ8gIgnA08AUXIJYLiKzgMHAMlV9VkTeBD6IfrjGxK+dewqZvzaDeWsz+XxbDqVlSrekZlwwpgeThnTm2H4dYn8U2WgoOeg67C39K+z5zvXuPv9lGHomNGpYF/59SRyqulhE+lSYPBbYoqpbAUTkddzRxnagyHtNKcaYiPt6137mpmYwf20Gq7yOeAM6t+LnJ/Vn6vBkRnRv0/D7VgSruBBW/B988gTsS4cex8CMP8PAU6M68GA0xdJxU3dckiiXBozDnbp6UkTGA4urerOIXAdcB9CrV68IhmlMw6OqpO7Yy9y16cxbm3noltnRPZL49bTBTB2eTP9OrXyOMgZtXgApt0HuN9BnPJzzLPSd0GATRrlYShyVfdKqqgXANTW9WVWfA54DGDNmTBWl4o0x5UrLlOXbcpibmsH76zLZsaeQhEbC2D7tufzY3kwZ1oVuDXUsqLrauxPm/sZd/O4wAK54xw1xHidiKXGkAT0DHvcAdvoUizEN0oHiUj75Opu5qRksWJ9FTn4RTRo3YsLAjtx8ykBOGdqFdi2b+B1m7CotgeXPw4cPQWkRTLobTrix3t8lFapYShzLgYEi0hfYAVwEXOJvSMbUf/sOFPPRxl3MW5vBwg1Z5BeV0rppYyYN6cy0EcmcNKgTLZvG0qYgRqWthHdvhozVMOAUmPEnaN/P76h84dftuK8BE4GOIpIG3KuqL4jITGAe7nbcF1V1rR/xGVPfZe8/yIJ1mcxbm8HSLbspKi2jY6smnHlEN6YOT+a4/h1o2tjuhApKYS588IC7AN46Gc5/CYad3eCvY1THr7uqLq5iegqQEuVwjGkQ0nILmLfWJYsV23IoU+jRrjmXH9ebqcOTObp3u4ZV5CjSVGH1GzD/LijYDcf+3I0p1ayN35H5zo5PjamnVJUtWe622XnrMkjdsReAwV1aM3PSAKaOSGZYV7tttlZ2bYL3boFtS6D70XDZf6HraL+jihmWOIypR8rKlNU78g71sdianQ/Akb3acsf0IUwdnkzfeBs8MJyKC12P74//6saROu0vcPRVrqiSOcQShzExLvC22XlrM0jPO0BCI+G4fh24+oQ+TBmWTHJSnA/zEQ6b3/f6ZGyDURfCqb+DVp39jiomWeIwJgYVl5bx6dbdzPGOLLL3l98224nbTh3MyUM707aF3TYbFnt3wtw7YN070GEgXDEL+p3kd1QxzRKHMTHiYEkpH2/OZk5qBgvWZ7KnoJgWTRKYNNjdNjtpSGda2W2z4VNaAp8/Bx89BGUlMPluOD7++mTUhn0LjfFRYVEpizZlkbImgw83ZLH/YAmtmzXmlKFdDvWxsAEEIyBthdcnYw0MmOL1yejrd1T1hiUOY6Js34FiPtyQxdzUDBZu3EVhcSntWiQyY2Qy00d25YT+HWnSuGGNphozCnNhwf2w8iXXJ+OCV9zotXbnWUgscRgTBXsKinh/XSZzUzNYsjmbotIyOrVuynlHd2f6iK6M69uexvFUczvaVGH1v2HeXVCY4/pkTLqz3pVsjRWWOIyJkF37DjJ/XQZzUzNY9vVuSrw6Fpcd25vpI5M5qpd1yIuKw/pkjIHL37I+GXVkicOYMMrIO8Dc1HTmpGaw3Ou93btDC64Z35fpI7oyukeSdciLluJCWPxnWPq465Nx+mNw1FUNrqiSHyxxGFNH23MKmOMliy+/2wPAwM6tmDlpANNGdGVo19aWLKJt03zXJ2PPtzDqIjj1QeuTEUaWOIyphfIKeSlr0lm70w31MaxrG247dRDTRnRlQGcreuSLvB2uT8b6Wa7W95WzXWElE1aWOIwJgqqyIWMfc1IzmJuazqZMVyHviJ5t+c30IUwf0ZVeHVr4HGUcKy2Bz5+Fj37v9cm4x+uTYZ0kI8EShzFVKC+nmpKaztzUDL7JzkcEjunTnnvPGMbU4clWIS8WbF8O7/4SMq1PRrRY4jAmQFmZ8uX2PYcucKflFh4aF+qaE/ty6vAudG5t40LFhMI9sOBeWPkytO5qfTKiyBKHiXulZcqKbTneaagMMvYeIDFBOHFAR26cPJApw6ycaszZ8B68ewvkZ8Gxv4BJv7E+GVFkicPEpZLSMj7dmsOc1HTmrc0ke/9BmjRuxEmDOnH7yMFMHtKFpOaJfodpKsrPhjm/htT/QufhcPFr0P0ov6OKO5Y4TNwoKilj6dfZzF2Twfx1GeQWFNM8MYHJQ2wQwZin6pLFnF/Dgb0w8U448Zd28dsn9isxDdqB4lKWbM5mzpp03l+fyb4DJbRq2piTh3Zm+oiunDSoE82b2CCCMW1vuuv5vTEFuh0FZz0NXYb5HVVcs8RhGpyCohIWbtzFnNQMPlyfSX5RKW2aNebUYcnMGJnMCQM62oiz9YEqfPkPmHc3lB6EKQ+66xkJttnym60B0yCUjzg7Z00GCzdlcaC4jPYtm3DmEd2YNqIrx/XrYCPO1ie522D2TbB1IfQ6Hs56Cjr09zsq47HEYeqtvIJiFqzPZE5qOos3fT/i7AVjejJtRDJj+9iIs/VOWRksf94NfS4CM/4MY66x8aVijCUOU6/k5Bcxf20Gc1IzWLol+7ARZ2d4I842shFn66fszTDrBvhuGfQ/Gc74K7Tt5XdUphKWOEzMy9p3gHlrM5mzJp3PvsmhtEzp1b4F15zYl+kjbcTZeq+0BJY9CR89DInN4Oz/B6Mvto58McwSh4lJ6XmFzE3NYM6aDJZ/m4Mq9OvUkp+f1J9pI5IZ3q2NJYuGICMV3rke0r+CIafDaY+6ynwmplniMDFje06BG3E2Nf3Q8ORDkltz08kDmTGyKwM7t7Jk0VCUFMGSP8OSR6FZWzj/JRh2th1l1BOWOIyvvsnOd7Us1mSwZkceACO6t+FXUwczbUQy/TvZ8OQNTtpKd5Sxaz2MvACmPQItO/gdlQmBJQ4TdZsz3fDkKWvS2ZCxD7DhyeNCUQEs/D0sexpaJcMlb8CgqX5HZWrBEoeJuEO1LNakk5KawZas/YjAmN7tuOf0YUwbkUx3G568Ydu2FGbNhJytcPRVMOUBaJbkd1SmlixxmIhQVdbsyGNOagZz1qSzbXcBjQTG9m3PFccNZ+rwZLq0seHJG7yD++D9e2HFC9C2N1wxC/qd5HdUpo7qReIQkZbAYuBeVX3X73hM5QJrWaSsyWDHHlfL4vj+HbhuQn9OHd6Fjq2a+h2miZYtC2D2zZCX5oYKmXw3NGnpd1QmDHxJHCLyInA6kKWqIwKmTwMeBxKAv6vqI95TtwNvRD1QU6PSMmXlt7mkrEk/rJbF+IGduOmUgUwZarUs4k5BDsy7C1b9y9X9/vE86DXO76hMGPl1xPES8BTwSvkEEUkAngamAGnAchGZBXQD1gF2XiNGlJSW8fk3OV5JVatlYQKsnw3v3erqZoy/FSb82nXqMw2KL4lDVReLSJ8Kk8cCW1R1K4CIvA6cBbQCWgLDgEIRSVHVsorzFJHrgOsAevWyYQrCrbi0jE++3s2cNenMX5dJTn4RzRMTmDSkE9NHdLVaFvFufxak/ArWvQ3JI+HS/0DX0X5HZSIkln7p3YHtAY/TgHGqOhNARK4CsitLGgCq+hzwHMCYMWM0sqHGh4MlpXy8OZuUNRm8vy6DvV4ti8lDOjNjZDInDepstSzinSqs+Y8rsFSU765jnHAzJNgRZ0MWS4mjsi6jhxKAqr4UvVDi14HiUq+WRTofrM9i/8ESWjdrzJRhXZgxoisnDrRaFsaTtwPe/SVsngc9joEzn4LOQ/yOykRBLCWONKBnwOMewE6fYokr+QdL+Gijq2Xx0cYsCopKadcikRkjk5kxsivH9+9otSzM91Thi5dh/j1QWgxTfw/jfgaNbIciXlSbOLwL1jeq6mNRiGU5MFBE+gI7gIuAS6Kw3Li090AxH67PYk5qOgs37uJgSRkdWzXhnCO7M2NkV8b1tVoWphI538DsG+GbxdBnPJz5BLTv53dUJsqqTRyqWioiZwFhTRwi8howEegoImm4/hkviMhMYB7udtwXVXVtOJcb7/IKinl/vRuefMlmV/ioS5umXDy2F9NHJDOmT3sSrJaFqUxZKXz+HHzwAEgCnP5XOOpKK7AUp4I5VbVURJ4C/g3kl09U1S9qu1BVvbiK6SlASm3na36ovPBRSmoGn3iFj7q3bc7lx7nCR0f2tMJHpga7NrlBCdM+h4GnwumPQVIPv6MyPgomcRzv/ftAwDQFJoc/HBMOVRY+Gt+XGSO6MsoKH5lglBbDJ0/Awj9AYnM451kYdaENfW5qThyqOikagZi6OVT4KDWD5dsOL3w0fWQyw7pa4SMTgvTV7igjYzUMPdPV/m7dxe+oTIyoMXGISBJwLzDBm7QIeEBV8yIZmKlZWq5X+GhNOl94hY8Gd7HCR6YOSg7C4j/Bx49B83ZwwSsw7Cy/ozIxJphTVS8CqcAF3uPLgf8Dzo1UUKZq27Lz3YizqemsTnO5e3g3K3xkwiBthVdgaQOMugimPQwt2vsdlYlBwSSO/qp6XsDj+0Xkq0gFZH5oS9b+Q7Us1qfvBWC0FT4y4VJUAB89BJ/+DVp3hUv+A4NO9TsqE8OCSRyFInKiqn4MICInAIWRDSu+qSobM/eRssbVstictR9whY/uPm0o00Yk06OdJQsTBts+hndmQu43cPTVXoGlNn5HZWJcMInjZ8Ar3rUOgFzgysiFFJ9UlbU79x6qv701Ox8RGNunPfef6QofJSfZKKMmTA7shQX3uQJL7frAlbOh74Sa3mUMUHPP8UbAYFUdLSJtAFR1b1QiiwOqylfb97gL3KnpbM9xhY+O7deea8b35dRhyXRqbYWPTJhtXgCzb4K9O+DY62HyXVZgyYSkpp7jZV5v7jcsYYRHWZmy8rtc5qzJYG5qOjvzXOGjEwZ05IZJAzllWBfaW+EjEwmHFVgaDNfMh55j/Y7K1EPBnKp6X0Ru44c9x3MiFlUDU1qmfP5NDnNSXZW8rH2u8NGEgZ24bepgTh5qhY9MhB1WYOk2mPArK7Bkai2YxPFj79/rA6YpYCObVaO4tIxPt+4mZU0G89dmsDu/iGaJjZg0uDPTRiQzeUhnWjezZGEibP8uSLnNCiyZsArmGsdlqro0SvHUa0UlZSzdkk3KmnTeX581sQCJAAAXfElEQVTJnoJiWjRJ8AofdWXi4E60aBJLI9mbButQgaXboWi/FVgyYRXMNY4/A8dFKZ5650BxKYs37WJuagbvr89k34ESWjdtzCnDujB9RDITBnWywkcmuvbudAWWNs2F7mPgrKetwJIJq2B2f+eLyHnAW6pqJVmBgqISr0peBh+uzyS/qJSk5olMG+4VPhrQgaaNLVmYKLMCSyZKgkkctwAtgVIRKcSVeFVVjateQvsPlvDB+kzmproqeQeKy+jQsglnHtGdGSOTObZfBxKt8JHxS+42mHUjfLPIFVg643Ho0N/vqEwDFczouK2jEUgsyissZsG6TOakZrB48y6KSsro1LopF4zpyfQRXTmmTzurkmf89YMCS4/BUVdZgSUTUcGMjivApUBfVX1QRHoCXVX184hH54Pc/CLeX5dJSmo6S7dkU1yqdE1qxqXjejFjZFeO7mWFj0yM2LUJZs2E7Z/BgClwxl+twJKJimBOVf0NKMMVbnoQ2A88DRwTwbiirrColGtfWcGyrbspLVN6tm/O1Sf0ZfqIZEb3aGvJwsSO0hKvwNIjVmDJ+CKYxDFOVY8SkS8BVDVXRBpc1+bmTRJo2TSBn07ox4yRXRnezQofmRiUscYNfZ6+ygosGd8EkziKRSQB1+kPEemEOwJpcJ69fIzfIRhTuZKDsPjP8PFfXIGl81+G4Wf7HZWJU8EkjieA/wGdReQh4EfA3RGNyhjzvbSVXoGl9e6U1LRHrMCS8VUwd1X9U0RWAifjbsU9W1XXRzwyY+JdcaErsLTsaWiVDJe8AYOm+h2VMUEdcaCqG4ANEY7FGFNu21J3x1TOVjj6Kq/AUlKNbzMmGmzgJGNiycF9rsDS8r9D295wxSzod5LfURlzGEscxsSKLR+4Akt5aTDu53DyPVZgycSkoBKHiPQGBqrqAhFpDjRW1X2RDc2YOFGYC/Puhq9ehQ4D4cfzoNc4v6MypkrB9By/FrgOaA/0B3oAz+Aulhtj6mLDe/DuLZC/C078JZx0hxVYMjEvmCOO64GxwGcAqrpZRDpHNCpjGrr8bJjza0j9L3QZAZe8Dt2O9DsqY4ISTOI4qKpF5b2oRaQxXmdAY0yIVF2ymPNrOLAXJt3lCiw1bnCDMZgGLJjEsUhE7gSai8gU4BfA7MiGZUwDtDcd3rsFNqZAt6NcgaUuw/yOypiQBZM47gCuAdYAPwVSVPX5iEYVQETOBk4DOgNPq+r8aC3bmLBQhS9fhXl3QelBmPIgHPsLSLCbGk39FMyg/Teo6vOqer6q/khVnxeRm+qyUBF5UUSyRCS1wvRpIrJRRLaIyB0Aqvq2ql4LXAVcWJflGhN1e76Df5zjOvN1GQ4//wROuNGShqnXgkkcV1Yy7ao6LvclYFrgBG8gxaeB6cAw4GIRCTyOv9t73pjYV1YGnz8PfzsOtn/uRrG96j2rymcahCp3e0TkYuASoK+IzAp4qjWwuy4LVdXFItKnwuSxwBZV3eot/3XgLBFZDzwCzFHVL6qJ9zrcbcP06tWrLuEZUze7v4Z3ZsJ3n0C/Sa6Ma7vefkdlTNhUd7z8CZAOdAQeDZi+D1gdgVi6A9sDHqcB44AbgFOAJBEZoKrPVPZmVX0OeA5gzJgxdteXib6yUjcg4UcPQUJTd/H7iEutwJJpcKpMHKr6LfAtcFyUYqns16Wq+gRuaHdjYlfWejf0+Y6VMHgGnPYXaNPV76iMiYhgeo7v4/t+G02ARCBfVduEOZY0oGfA4x7AzjAvw5jwKimCjx+DxX+CZm3gvBdgxHl2lGEatGDqcbQOfOzdHjs2ArEsBwaKSF9gB3AR7hqLMbFp55fuWkZmqksW0/8ILTv6HZUxERfMXVWHUdW3gcl1WaiIvAYsAwaLSJqIXKOqJcBMYB6wHnhDVdfWZTnGRETxATf0+fMnu6FDLvoX/OhFSxombgRzqurcgIeNgDHUccgRVb24iukpQEpd5m1MRH33mbuWsXszHHkZnPo7VwPcmDgSTC+kMwL+XwJsA86KSDTGxKqifPjgQfjsGUjqAZe9BQNsgGgTn4K5xnF1NAIxJmZtXQSzboA938Ix18Ip90LT1jW/z5gGqroOgE9SzSkpVb0xIhEZEysO5MH8e+CLl6F9P7gqBfqc4HdUxviuuiOOFVGLwphYs2kezL4Z9mfA8TfCpDshsbnfURkTE6rrAPhy4GMRae0m6/6IR2WMXwpyYO4dsPrf0GkoXPgq9Dja76iMiSnB3FU1AvgHrnSsiMgu4Aq7VdY0OGvfhpTbXA3wk26H8bdC46Z+R2VMzAnmrqrngFtU9SMAEZkIPA8cH8G4jImefZmQciusnw1dR8Pl/4PkkX5HZUzMCiZxtCxPGgCqulBEWkYwJmOiQxVWve5OTRUXwin3wXE3WK0MY2oQzC9kq4jcgztdBXAZ8E3kQjImCvLS3MXvLe9Dz3FuJNuOA/2Oyph6IZjE8WPgfuAt3Ai2iwDr22Hqp7Iy+OIlmP9b0FKY9gcYey00SvA7MmPqjWA6AOYCN8KhKn0tVXVvpAMzJuxytsKsG2HbEuh7kiuw1L6v31EZU+/UOMihiPxLRNp41zXWAhtF5FeRD82YMCkvsPS34yF9FZzxBFzxjiUNY2opmNFxh3lHGGfjBiDsBVwe0aiMCZesDfDiVJh3J/Q7CX7xKRx9pdXLMKYOgrnGkSgiibjE8ZSqFouIlWY1sa20GJY+Dov+AE1awbnPw8jzLWEYEwbBJI5ncSPirgIWi0hvwK5xmNiVvsoNfZ6xBoafA9P/BK06+R2VMQ1GMBfHK9b8/lZEJkUuJGNqqeQgLPqjK+XasqMbLmToGTW/zxgTkmCGHOkA3AuciBst92PgAWB3ZEMzJgTbl7ujjOyNcMSlMPUhK7BkTIQEc3H8dWAXcB7wI+///45kUMYEragA5t4JL0xxxZYu/S+c/TdLGsZEUDDXONqr6oMBj38nImdHKiBjgvbNYldgKXcbjLnGDRnSrI3PQRnT8AWTOD4SkYuAN7zHPwLei1xIxtTgwF5YcC+seBHa9YWr3oM+J/odlTFxo7oKgPtw1zQEuAV41XuqEbAfd93DmOja/L4bY2rfTjhuJky6C5q08DsqY+JKdYWcrKiyiR0FOa4T36rXoNMQuOB96DHG76iMiUtBjR8tIu2AgUCz8mmqujhSQRlzmHWz4L1boTAHJvzK/VmBJWN8E8ztuD8BbgJ6AF8BxwLLgMmRDc3Evf1ZriLfuncgeRRc9l/oOsrvqIyJe8EccdwEHAN8qqqTRGQIbph1YyJDFVa/AXNvd7fYTr4HTrgJEhL9jswYQ3CJ44CqHhARRKSpqm4QkcERj8zEp7wd8O4vYfM86DEWznoKOtnXzZhYEkziSBORtsDbwPsikgvsjGxYJu6owhcvw/x7oKwEpj4M435qBZaMiUHBjFV1jvff+0TkIyAJmBvRqEx8yfkGZt/oOvT1GQ9nPmm1MoyJYUHdVVVOVRdFKhATh8pK4fPn4IMHQBLg9L/C0VfZ0OfGxLiQEocxYbNrE8yaCds/g4GnwumPQVIPv6MyxgQh5hOHV7L2b0ARsFBV/+lzSKYuSkvgk8dh4R9cj+9znoVRF9pRhjH1SDCj44adiLwoIlkiklph+jQR2SgiW0TkDm/yucCbqnotcGbUgzXhk7EG/j7ZnZoaPA2u/xxGX2RJw5h6xpfEAbwETAucICIJwNPAdGAYcLGIDMN1PNzuvaw0ijGacCk5CB8+BM9NhL074YJX3F+rzn5HZoypBV9OVanqYhHpU2HyWGCLqm4FEJHXgbOANL7vtV5lohOR64DrAHr16hX+oE3tpK1wBZZ2bYBRF8G0h6FFe7+jMsbUgV9HHJXpzvdHFuASRnfgLeA8Efl/wOyq3qyqz6nqGFUd06mT1Zf2XVEBzL/bFVg6uA8u+Q+c+6wlDWMagFi6OF7ZiW5V1Xzg6mgHY+pg21J3x1TOVjj6apjygBVYMqYBiaXEkQb0DHjcA+uhXr8c3AcL7oPlf4d2feDK2dB3gt9RGWPCLJYSx3JgoIj0BXYAFwGX+BuSCdqWBa7AUl4aHPsLmHw3NGnpd1TGmAjwJXGIyGvARKCjiKQB96rqCyIyE5gHJAAvqupaP+IzISjMhXl3wVf/hI6D4Jr50HOs31EZYyLIr7uqLq5iegqQEuVwTG2tfxfeuwXys2H8rTDh15DYrOb3GWPqtVg6VWXqi/xsSPkVrH0LuoyES96Abkf4HZUxJkoscZjgqcKaN2HOr6Fov7uOccLNVmDJmDhjicMEZ+9OePcW2DQHuh8NZz0NnYf6HZUxxgeWOEz1VOHLf8C8u6H0IJz6EBz7cyuwZEwcs8Rhqpa7DWbfBFsXQu8T4cwnoEN/v6MyxvjMEof5obIyWP48LLjfjVx72l9cD/BGsTRCjTHGL5Y4zOGyt7jhQr5bBv1PhjMeh7Y9a36fMSZuWOIwTmkJLHsKPvq964tx9v+D0RdbrQxjzA9Y4jCQuRbe/gWkfwVDTofTHoXWyX5HZYyJUZY44llJESx51P01S4LzX4JhZ9tRhjGmWpY44tWOlfDODZC1FkaeD9P+AC07+B2VMaYesMQRb4oL3XWMZU9Bqy5w8esweLrfURlj6hFLHPHk22WujGvO13DUFTDlQWje1u+ojDH1jCWOeHBwP3xwP3z+vLu19vK3of8kv6MyxtRTljgauq8/hFk3Qd52GHsdnPxbaNrK76iMMfWYJY6GqnAPzL8LvnwVOgyAq+dA7+P8jsoY0wBY4miINqS4Akv7M92w5xPvgMTmfkdljGkgLHE0JPnZMOd2SH0TOg+Hi/4F3Y/yOypjTANjiaMhUHXV+FJ+BQf2wsTfwIm3QOMmfkdmjGmALHHUd3vT4b1bYeN70O0oV2CpyzC/ozLGNGCWOOorVfjqnzD3TldgacoDcOz1kGCr1BgTWbaVqY/2fOcKLH39IfQ6Ds58CjoO8DsqY0ycsMRRn5SVwYoXYMF97ohj+p/gmJ9YgSVjTFRZ4qgvdn8N78yE7z6BfpNcgaV2vf2OyhgThyxxxLqyUlj2NHz0ECQ0daeljrzMhj43xvjGEkcsy1znBiXc+QUMmg6nPwZtuvodlTEmzlniiEUlRfDxY7D4T9CsDZz3Aow4z44yjDExwRJHrNn5pbuWkZnqksX0P0LLjn5HZYwxh1jiiBXFB2DRI7D0CZcoLvwnDD3d76iMMeYHYj5xiMjZwGlAZ+BpVZ3vc0jh992n7ihj92Y44jKY+jto3s7vqIwxplIR7QAgIi+KSJaIpFaYPk1ENorIFhG5o7p5qOrbqnotcBVwYQTDjb6ifDco4YvToOQAXPYWnP20JQ1jTEyL9BHHS8BTwCvlE0QkAXgamAKkActFZBaQADxc4f0/VtUs7/93e+9rGLYuhFk3wp5v4Zhr4ZR7oWlrv6MyxpgaRTRxqOpiEelTYfJYYIuqbgUQkdeBs1T1YeAHJ/VFRIBHgDmq+kUk442KA3kw/x744mVo3w+uSoE+J/gdlTHGBM2Paxzdge0Bj9OAcdW8/gbgFCBJRAao6jOVvUhErgOuA+jVq1eYQg2zjXPh3V/C/gw4/gaYeCc0aeF3VMYYExI/EkdlnRG0qher6hPAEzXNVFWfA54DGDNmTJXz80VBjruWseYN6DQULnwVehztd1TGGFMrfiSONKBnwOMewE4f4oiOtW9Dym1QmAsn3Q7jb4XGTf2Oyhhjas2PxLEcGCgifYEdwEXAJT7EEVn7MiHlVlg/G7qOhsv/B8kj/Y7KGGPqLKKJQ0ReAyYCHUUkDbhXVV8QkZnAPNydVC+q6tpIxhFVqrDqdZh7BxQXwsn3wvE3WoElY0yDEem7qi6uYnoKkBLJZfsiLw1m3wxb3oee49xItp0G+R2VMcaEle0Gh0NZGXzxEsz/LWgpTPsDjL0WGiX4HZkxxoSdJY66ytnqOvJtWwJ9J8AZT0D7vn5HZYwxEWOJo7bKSuGzZ+CDByEh0VXkO+pKG/rcGNPgWeKojV0bXYGltOUwcKorsJTU3e+ojDEmKixxhKK0GJY+Dov+AE1awrnPw8jz7SjDGBNXLHEEK301vPMLyFgDw86GGX+CVp39jsoYY6LOEkdNSg7Coj/C0r9C8/ZwwT9g2Jl+R2WMMb6xxFGd7cth1kzYtQFGXwxTfw8t2vsdlTHG+MoSR2WKCuCjh2DZ09CmG1z6Jgyc4ndUxhgTEyxxVPTNEph1A+R+A2N+DKfcD83a+B2VMcbEDEsc5YryXYGlFS9Auz5w5bvQd7zfURljTMyxxFGuUSJs/wyOvR4m3+VutzXGGPMDljjKNW4C135otTKMMaYGjfwOIKZY0jDGmBpZ4jDGGBMSSxzGGGNCYonDGGNMSCxxGGOMCYklDmOMMSGxxGGMMSYkljiMMcaERFTV7xjCTkR2AXuAvIDJSdU8Dvx/RyA7DGFUXF5tX1vVc5VNr66NFR9bm+OrzeFqb1Ux1eZ14WpzVc/FS5vD9b3uraqdgnqlqjbIP+C5YB9X+P+KSCy/tq+t6rnKplubrc1VtTlc7Q2lzTW9Llxtruq5eGlzJL7XNf015FNVs0N4XPG5SCy/tq+t6rnKplubrc0VH/vZ5ppeF6421/R5hEMstzkS7a1WgzxVVRciskJVx/gdRzRZmxu+eGsvWJsjqSEfcdTWc34H4ANrc8MXb+0Fa3PE2BGHMcaYkNgRhzHGmJBY4jDGGBMSSxzGGGNCYokjSCLSSEQeEpEnReRKv+OJBhGZKCJLROQZEZnodzzRIiItRWSliJzudyzRICJDvXX8poj83O94okFEzhaR50XkHRE51e94okFE+onICyLyZl3nFReJQ0ReFJEsEUmtMH2aiGwUkS0ickcNszkL6A4UA2mRijVcwtRmBfYDzYifNgPcDrwRmSjDKxxtVtX1qvoz4AIg5m9fDVOb31bVa4GrgAsjGG5YhKnNW1X1mrDEEw93VYnIBNwG8BVVHeFNSwA2AVNwG8XlwMVAAvBwhVn82PvLVdVnReRNVf1RtOKvjTC1OVtVy0SkC/AXVb00WvHXRpjaPAo3bEMzXPvfjU70tROONqtqloicCdwBPKWq/4pW/LURrjZ773sU+KeqfhGl8GslzG2u8/arcV3eXF+o6mIR6VNh8lhgi6puBRCR14GzVPVh4AenKEQkDSjyHpZGLtrwCEebA+QCMV+QPUzreRLQEhgGFIpIiqqWRTTwOgjXelbVWcAsEXkPiOnEEab1LMAjwJxYTxoQ9t9zncVF4qhCd2B7wOM0YFw1r38LeFJExgOLIxlYBIXUZhE5F5gKtAWeimxoERNSm1X1LgARuQrviCui0UVGqOt5InAubucgJaKRRU6ov+cbgFOAJBEZoKrPRDK4CAl1PXcAHgKOFJHfeAmmVuI5cUgl06o8b6eqBUBYzg/6KNQ2v4VLmPVZSG0+9ALVl8IfStSEup4XAgsjFUyUhNrmJ4AnIhdOVITa5t3Az8Kx4Li4OF6FNKBnwOMewE6fYokWa7O1uaGyNkexzfGcOJYDA0Wkr4g0AS4CZvkcU6RZm63NDZW1OYptjovEISKvAcuAwSKSJiLXqGoJMBOYB6wH3lDVtX7GGU7WZmsz1mZrc6TiiYfbcY0xxoRPXBxxGGOMCR9LHMYYY0JiicMYY0xILHEYY4wJiSUOY4wxIbHEYYwxJiTxPOSIMVElIvfhRjjNBuarakPv2WwaKDviMCb6rgK6+R2EMbVlicOYCBKRu7xCOwuAwd7kMcA/ReQrEWnuY3jG1IqdqjImQkTkaNz4QUfifmtfACuBFcBtqrrCx/CMqTVLHMZEznjgf96Q/IhIQx90z8QJO1VlTGTZYHCmwbHEYUzkLAbOEZHmItIaOMObvg9o7V9YxtSNjY5rTASJyF3AFcC3uMI764BvgN8DhcBxqlroX4TGhM4ShzHGmJDYqSpjjDEhscRhjDEmJJY4jDHGhMQShzHGmJBY4jDGGBMSSxzGGGNCYonDGGNMSCxxGGOMCcn/B3c9uguhcbBHAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "dt = np.logspace(-1, -6, 10)\n",
    "PopulationE = np.zeros(dt.size)\n",
    "PopulationRK2 = np.zeros(dt.size)\n",
    "t_endE = np.zeros(dt.size)\n",
    "t_endRK2 = np.zeros(dt.size)\n",
    "\n",
    "\n",
    "\n",
    "for i in range(dt.size):\n",
    "    t, P = Euler(Malthus, 0.0, P0, SimulationTime, dt[i])\n",
    "    PopulationE[i] = P[-1]\n",
    "    t_endE[i] = t[-1]\n",
    "    t, P = RK2(Malthus, 0.0, P0, SimulationTime, dt[i])\n",
    "    PopulationRK2[i] = P[-1]\n",
    "    t_endRK2[i] = t[-1]\n",
    "\n",
    "PopErrorE = PopulationE - P0*np.exp(GrowthRate*t_endE)\n",
    "PopErrorRK2 = PopulationRK2 - P0*np.exp(GrowthRate*t_endRK2)\n",
    "\n",
    "plt.loglog(dt, np.abs(PopErrorE),label='Euler error')\n",
    "plt.loglog(dt, np.abs(PopErrorRK2), label='RK2 error')\n",
    "plt.title('Error of the numerical simulation for GrowthRate=' + str(GrowthRate))\n",
    "plt.xlabel('dt')\n",
    "plt.ylabel('absolute error')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}