stochastic_processes_data_analysis/2_Euler_Runge_Kutta.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# The Euler method for numerical integration\n",
    "\n",
    "- [Introduction](#intro)\n",
    "- [Ordinary differential equations (ODE)](#ode)\n",
    "    - [1st order ODE](#1ode)\n",
    "    - [Numerical calculation](#calculation)\n",
    "- [A very simple problem](#problem)\n",
    "    - [Euler Method](#euler)\n",
    "    - [Leap-frog method](#leapfrog)\n",
    "    - [Runge-Kutta 2nd order](#runge2)\n",
    "    - [Runge-Kutta 4th order](#runge4)\n",
    "- [Extra](#extra)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div id=\"intro\" />\n",
    "\n",
    "## Introduction\n",
    "\n",
    "The Euler method for numerical integration is the most basic explicit method for numerically solving ordinary differential equations. Although many other methods with higher accuracy have been proposed, the simplicity of the Euler method makes it ideal when studying the subject for the first time. We will use it to simulate the motions of Brownian particles."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<div id=\"ode\" />\n",
    "\n",
    "##  Ordinary differential equations (ODE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div id=\"1ode\" />\n",
    "\n",
    "### 1st order ODE\n",
    "\n",
    "Consider the following 1st order differential equation.\n",
    "\n",
    "$$\n",
    "\\frac{dy(t)}{dt}=f(y(t),t) \\tag{A1}\n",
    "$$\n",
    "\n",
    "Assume that the initial conditions are $y=y_0$ at time $t=t_0$. We need to determine $y(t)$, for any $t \\ge t_0$.\n",
    "\n",
    "#### Formal solution\n",
    "\n",
    "Integrate Eq.(A1) over time, from $0\\rightarrow t$, to obtain the formal solution for $y(t)$\n",
    "\n",
    "$$\n",
    "y(t)=y_0+\\int_{{{t}_{0}}}^{t}{}\\,d{t}'f(y({t}'),{t}') \\tag{A2}\n",
    "$$\n",
    "\n",
    "Let us consider an ordinary differential equation, shown here in Eq.(A1). We want to know the temporal values of $y$ for any time $t$, under the initial condition $y=y_0$ at $t=t_0$. \n",
    "\n",
    "First, let us write down the formal solutions of the differential equation by integrating Eq.(A1) in time, from $t_0$ to $t$. Eq. (A2) gives the exact value of $y$ at any given time $t$, starting from the initial condition of $y_0$ at $t=0$. However, an analytical evaluation of the time integrals is only possible for very limited cases. In most situations, we must solve these differential equations numerically. The Euler method is the simplest method for performing this numerical integration."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<div id=\"calculation\" />\n",
    "\n",
    "### Numerical calculation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Discretization\n",
    "\n",
    "Divide the total time span $t_0\\rightarrow t$ into $N$ equally spaced segments, each describing a time increment $\\Delta t$.\n",
    "\n",
    "![](images/discretize.png) \n",
    "$$\n",
    "t_i\\equiv t_0+i\\Delta t,\\ \\ \\ \\ \n",
    "y_i\\equiv y(t_i),\\ \\ \\ \\ \n",
    "f_i\\equiv f(t_i,y(t_i))\\tag{A3}\n",
    "$$\n",
    "\n",
    "To solve differential equations numerically using computers, we must first discretize the time axis. We divide the total time span, from $t_0$ to $t$, into $N$ equally spaced segments, each describing a short time increment $\\Delta t$, and define the following discrete variables."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Advancing the solution forward a small step $\\Delta t$\n",
    "\n",
    "Integrate Eq.(A1) over a small time interval, from $t_i \\rightarrow t_{i+1}(=t_i+\\Delta t)$,\n",
    "\n",
    "$$\n",
    "y_{i+1} \n",
    "=y_i+\\int_{{{t}_{i}}}^{t_{i+1}}{}\\,d{t}'\\,f(y({t}'),{t}') \\hspace{50mm}\\tag{A3}\n",
    "$$\n",
    "$$\n",
    "= y_i+\\int_{0}^{\\Delta t}\\,d{\\tau}\\,f(y({t_i+\\tau}),t_i+\\tau)\\ \\ \\ \\ \\ \\ \\ \\ (\\tau\\equiv t'-t_i)\\tag{A4}\n",
    "$$\n",
    "$$\n",
    "= y_i+\\int_{0}^{\\Delta t}\\,d{\\tau}\\,\\left[f_i+\\mathcal{O}(\\tau)+\\mathcal{O}(\\tau^2)+\\cdots)\\right]\\hspace{14mm}\\tag{A5}\n",
    "$$\n",
    "$$\n",
    "= y_i+\\left[\\tau\\, f_i+\\mathcal{O}(\\tau^2)+\\mathcal{O}(\\tau^3)+\\cdots)\\right]_{0}^{\\Delta t} \\hspace{23mm}\\tag{A6}\n",
    "$$\n",
    "\n",
    "$$\n",
    "= y_i+\\Delta t \\,f_i+\\mathcal{O}(\\Delta t^2)+\\mathcal{O}(\\Delta t^3)+\\cdots \\hspace{23mm}\\tag{A7}\n",
    "$$\n",
    "\n",
    "Now consider how to perform the integration over a small time increment $\\Delta t$,  from $t_i$ to $t_{i+1}$. Using the Taylor expansion, we expand the integrand around $t_i$, in increasing powers of $\\tau$, where $\\tau = t^\\prime - t_i$. For simplicity, we only explicitly write the zero-th order term, proportional to $f_i$. The integrals can be performed analytically, with each term one order in $\\Delta t$ higher than the previous one."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div id=\"problem\" />\n",
    "\n",
    "## A very simple problem\n",
    "\n",
    "Numerically solve the following differential equation and determine $y(t)$ for $0\\le t\\le 10$ with the initial condition $y=1$ at $t=0$. Then compare it with the analytical solution $y=\\exp(-t)$.\n",
    "\n",
    "$$\n",
    "\\frac{dy(t)}{dt}=-y(t) \\tag{A22}\n",
    "$$\n",
    "\n",
    "Let us consider a very simple problem here and try to apply the Euler method to numerically solve it. We want to numerically solve the following differential equation and determine $y(t)$ for $0\\le t\\le 10$, with initial condition $y=1$ at $t=0$. Then we compare it with the analytical solution $y=\\exp(-t)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "plt.style.use(\"ggplot\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<div id=\"euler\" />\n",
    "\n",
    "### Euler method\n",
    "\n",
    "- Difference equation $\\rightarrow$ 1st order in $\\Delta t$\n",
    "\n",
    "\n",
    "$$\n",
    "y_{i+1} = y_i+\\Delta t \\,f_i \\tag{A8}\n",
    "$$\n",
    "\n",
    "- Simulation procedure $\\rightarrow$ Explicit method\n",
    "\n",
    "$$\n",
    "y_0, f_0 \\xrightarrow[\\Delta t]{\\rm Eq.(A8)}\n",
    "y_1, f_1 \\xrightarrow[\\Delta t]{\\rm Eq.(A8)}\n",
    "\\cdots\n",
    "y_i, f_i \n",
    "\\cdots\n",
    "\\xrightarrow[\\Delta t]{\\rm Eq.(A8)}\\tag{A9}\n",
    "y_N, f_N \n",
    "$$\n",
    "\n",
    "- Forward difference approximation (1st order)\n",
    "\n",
    "$$\n",
    "\\left.\\frac{dy(t)}{dt}\\right|_{t=t_i}\\simeq\\frac{y_{i+1}-y_i}{\\Delta t}\\tag{A10}\n",
    "$$\n",
    "\n",
    "In the limit when the time increment $\\Delta t\\rightarrow 0$, we can neglect all terms in Eq.(A7) of second order or higher. This first order approximation is known as the Euler method, shown here as Eq.(A8). By repeatedly applying Eq.(A8) $N$ times, we can obtain an approximation of $y$ at a given time $t=N \\Delta t$. This process is schematically shown in Eq.(A9).\n",
    "\n",
    "The Euler method is a first-order method, which means that the local error (the error per step) is proportional to the square of the step size, and the global error (the error at a given time $t$) is proportional to the step size. Although simple and not very accurate itself, the Euler method serves as the basis to construct more accurate methods. \n",
    "\n",
    "From Eqs.(A1) and (A8), it is straightforward to solve for $f$, and to recognize that the Euler method is a forward difference approximation to the derivative of $y$, as defined in Eq. (A10)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABB3klEQVR4nO3deXhcZdn48e85M5OZyZ6ZyZ60adOU0p2SUihaWhoWEZBLWWRRuPjJq61Kiy8/rBWQrb+3KrWIFkGpRV99sfoqyCIIBVkD0oXuW9qmS9K0WSZ7MsnMnPP740xmMk1K0jTJZCb357pyZc4690nTO8885zn3o+i6riOEECLqqZEOQAghxOCQhC6EEDFCEroQQsQISehCCBEjJKELIUSMkIQuhBAxwhzJNz9+/PiAjnO5XNTW1g5yNCObXPPoINc8OpzNNefk5Jx2m7TQhRAiRkhCF0KIGCEJXQghYkRE+9CFENFP13U8Hg+apqEoyhkff/LkSTo6OoYgspGrr2vWdR1VVbHZbGf0M5WELoQ4Kx6PB4vFgtk8sHRiNpsxmUyDHNXI1p9r9vl8eDwe7HZ7/8/b1w5PPfUUW7ZsISUlhVWrVvXYrus669at49NPP8VqtbJ48WLGjx/f7wDOhO7zQnUVnrIdaPt2oRRNQTln2pC8lxCifzRNG3AyF6dnNpvP+JNLn33o8+fPZ/ny5afd/umnn3LixAmefPJJ/uM//oNnn332jALoL+2NF9C+cyOHnv07B3/zB3wv/Rl9+6YheS8hRP8NpJtF9M+Z/mz7/LM6efJkqqurT7t906ZNzJs3D0VRmDhxIq2trdTX15OWlnZGgfQpKZVO1c7uibeComLyd+BsP05GWQfpWWYSElX5xRJCjGpn/TnJ7XbjcrmCy06nE7fb3WtC37BhAxs2bABg5cqVYcf1xTt5Gp+kzwLF+FDhN1mpNo2jeks7AInJZvLGxJOTH092XjxxcbE1gMdsNp/RzysWyDVHh5MnT551l8vZHp+dnc25554bXL7uuuu4++67z+qcXXbu3MmJEycoKSkZlPN16c81W63WM/p9GNaOr5KSkrAfypk8KaXbEkltqWSX1kqeYiVFCQ+9pcnH3p1N7N3ZhKJAmstEepaFjEwzKQ5T1Lfe5Wm60SEar7mjo+OsbmqazWZ8Pt9ZxWCz2XjjjTfC1p3tObts27aN7du3M3/+/EE5H/T/mjs6Onr8PnzWk6JnndAdDkfYG9bV1eFwOM72tD0oVhsucxMH2o7wkc1BMiaWHn0f//lXU9sch7/bz0bXwV3jx13jZ98OsMQppGeaSc8yk55lwR4fW613IUYC/13XDuy4M9jX9JuX+r1vU1MTX/ziF1m3bh0TJkxg8eLFXHzxxdx6660sW7aMbdu24fF4+OIXv8i9994LwNatW3nwwQdpa2vDarXy/PPP8/jjj+PxePjkk0/4zne+w5e+9KUzvMLhc9YJvbi4mNdff52LL76YsrIy4uPjB7//vEtOPrltNdTYHDThx9NUxkUpe2DhPNx1PmpO+Kiu8tHUEP4r4u3UOX7My/FjXqCdxGSV9CwL6VlmnOlmzObobr0LMdp5PB4uu+yy4HJX4l2xYgX33HMP3/jGN2hsbOTWW28F4Pvf/z5paWn4/X5uuukmdu/ezYQJE1i0aBG/+tWvmDlzJs3Nzdjtdu699162b9/OihUrInV5/dZnQn/iiSfYvXs3zc3NfOtb3+LGG28MflS4/PLLOe+889iyZQt33303cXFxLF68eMiCVbLyyT1cw1bHOQAcj0+H48dQTQquDAuuDAvnTocOj2Yk9xNeak746OwInza1pUmjpamD8v0dqCo4XGZcWWbSM82kpEV/94wQo43NZuPNN9/ssX7evHm88sorLF++PGz7yy+/zB//+Ef8fj8nT56krKwMRVHIyMhg5syZACQlJQ1X+IOmz4S+dOnSz9yuKArf+MY3Biuez5aTT+7u0FDFyvgM9BPHeuxmtankFcSRVxCHrus0NfipOWm04N01PjQttK+mQW21j9pqH3uBOKuCK8PonnFlWohPkO4ZIaKVpmmUlZVht9tpbGwkJyeHo0eP8swzz/Dqq6+SmprK0qVL8Xg8kQ51UETV0wBKdh65ba8Flyvi0+F46WcfoyikpJlJSTMzYRL4fDruGiO515zw0tykhe3f2RHePZOQpAb63y04081Y4qT1LkS0+PWvf01RURHLli3je9/7Hi+99FKwKyU5OZmamhr+9a9/cdFFF1FYWEh1dTVbt25l5syZtLS0YLPZSExMpKWlJdKX0i9RldDJzie3LTQmvjI+A72mCt3nRTFb+nUKs1khI9tCRrYFsNPeplF70kjuNSd7ds+0Nmu0Nndy+EAnigKpDlOw9Z7mNKGqkuCFgDO7YdndYIxyObUPfcGCBdx44408//zzvPrqqyQmJjJnzhx+/vOfc++99zJ16lTmzZtHTk4Os2fPBiAuLo5f/epX3H///Xg8Hmw2G+vXr2fu3LmsWbOGyy67bMTfFFV0Xdf73m1oDGSCC99/fp3bpn+PdrMNgLWlj+D84UqUnDFnHY/RPaNRe9JI7nU1PrTPuAVvMoMz3eh7d2VaSEoZuoebonE429mSa44ObW1txMfHD/j4wUjo0aa/19zbz3ZIhy0ONyXQSj+QbCTwyvgMnFUVMAgJ3eieMZGSZqJwEvj9Ou5aH7UnfNSc9NFYH57d/T6orjJG1oAHq03BlRlK8DI8UggxnKIzoZ+oCUvo06qOMRTtYpNJIT3TQnqmhXOBjg6NupNGcq856aO9Nbz/vcOjU3nES+WRQP97ooor02x8ZZiJs0qCF0IMnahL6GTncdXm1/h89VZy2mpI99RDVdywvLXVqpIzJo6cMcb7tbb4qTnho/akMUrG23lK/3uLRmtLJ0cOdgKQkmYKJniHS8a/CyEGV9QldCU7n6Lm8KGKelXPoYvDISHRRMIEEwUTrOiaTmOD37jBetKHu7Zn/3tjvZ/Gej8H9xrj39OcJlyZFpwZZrnBKoQ4a1GX0MnO77nuRCW65kdRI1ckX1EVUh1mUh1mJpxr9L/X1xrJvfakj4Z6P3RrwGsa1NX4qasxsr7JbDzglJ5pxplhJiXVhCIJXghxBqIvoaekoSQkobc2h9Z5O6GuBtKzIhfXKUwmBVemBVemMZzS26lRV+On9qSX2pO+HuPf/T4CY+ONO9+WOAVneqj/3emM2GAkIUSUiLq7dIqiYM4bC4BXMXE0PpMWsx0i1O3SX5Y4laxcC1NnxTP/C8lcdm0y510YT/64OOzxPVvi3k6dE5Vedm5p553Xm1m/7jCbP2rlyMEOWlv8RHC0qRAjTm5uLg8//HBw+emnn+51hrWhtG3bNh544IEBHXv99dezbdu2s44h+lrogCmvgJ/r5/B29mw0ReV7u/7I56uOoUyfHenQ+s1mV8kbG0feWOMGa1uL3yhBELjB2uEJT9jt7X7aj/o5ftRrHB9vlChwZZhxZkiJAjG6Wa1WXnvtNb773e8OSbXXvvh8PmbMmMGMGTOG/b27i8qEbs4vIKG8HC0w2cXx+PQR30LvS3yiiTGJJsaMt6LrOi1NWrDGTF0vI2g8bToVh71UHDaGSMYnqDiDCd4sY+DFqGIymbj11lv59a9/zbJly8K2LV26lJKSEq6++moAioqKKCsro7S0lFWrVpGcnMzevXu55pprmDRpEmvXrsXj8bB27VoKCgqoq6tj2bJlVFZWAvDwww8ze/ZsVq1axeHDhzl69Ci5ubncdtttPP300/z+97+ntbWV+++/n+3bt6MoCvfccw9f/OIXg2V7Ozo6uOqqq4JlewdLdCb0vAJy2/4dXK6MT0c/vjmCEQ0uRVFISjGRlGJiXJGR4FWSObC/htqTRoGxUx8ya2vVaCvv5Fi5MUQyIdFI8F1J3maXBC+Gx/Pba/jTjrp+7Xv5hBSWXBw+0GHNv6t440BjcPmr05zcPD29z3PdcccdlJSUnFHF1927d/POO++QmprK3Llzufnmm3n11Vd59tln+e1vf8sjjzzCgw8+yF133cUFF1xAZWUlt9xyC++++y4AZWVlvPDCC9jtdkpLQ3WlnnjiCZKSknjrrbcAaGhoAEJlexVF4Stf+Qq7d+9m8uTJ/Y63L9GZ0PPHkdteE1yuiM+AQ8fQNQ1Fjb3EpSgKTpcVXbFReA5omk5jvT/YenfX+PCfMkSyawz80UOBBJ+kGjdZA0leEryINUlJSVx//fWsXbsWu93er2NmzJhBZmYmAGPHjuWSSy4BYNKkScEE/f7777N///7gMS0tLbS2tgJGCfHe3uv999/nqaeeCi6npqYCvZftHfUJXXVlkudvCi5Xxmfg7/Cg1lWPqJEuQ0VVFdKcZtKcZorOBc2v0+DuluDreo6B7yoyJglexLJvfOMbXHnlldx0003BdWazGS1QM1vTNLxeb3BbXFzooURVVYPLqqoGa61omsbLL7+MzWbr8X5nUsOme9lel8vFd77znUEv2xuVCV1RVZIyM0jraKLemkynycJJu5O8isOjIqGfSjUpONLNONLNMCUwBr7Oj7vG6IOvrw2vAQ+9JPiuLpp06YMXZ+fm6en96iI5nW/Pyebbc7IHdGxaWhrXXHMNzz//PF/96lcByMvLY8eOHVx77bW88cYbYQm9Py655BLWrVvHokWLAGPS6KlTp37mMfPmzeO5557jkUceAYwul+5le6urq4NlewdT1P6vVfIKGNN6Irh8NCETvfJw5AIaQUwmYwTMxCk25i5I5Movp3DRgkQmTrHhTDfRW69Ua4vG0UOdfPrvNja83MRbrzax9ZM2jpV30tZ6JrM+ChFZ3/zmN3G73cHlW2+9lY8++oiSkhI2b958xpUhH330UbZt20ZJSQnz58/nv//7v/s8ZsmSJTQ2NnLppZdSUlJCaWkpU6ZMCZbtXbRoUbBs72CKuvK5YJQYrf7f3/Pbjyt4OX8eAF8t/yc3utoxfWtZH0dHp8Esq9rVgq+rNkoE99aCP5U93njQyRFowSckDl2p4C7RWEr2bEXjNUv53DMn5XNPoeQWMKY1NB3d0YQsqHgncgFFka4WvCvD+Of3+3Ua6vzU1Zy+D769TafiiJeKI8bHVatNCXbPONPNJCYPfYIXQny2qE3o5I49pcslC/YcR+/oQLFaIxhY9DGZlOAQR6aEbrLW1RgteHetD/8pjYkOT/ep+kKlChzpJpzpZpJTpdiYEMMtahO6Ep9AvtXoLUrubMHV0WA8Dn/8KIwrinB00a37TdYiQsMkgy34Wh++U+4rdZUqOFFpbOgqNuZIN1rwqQ4TJpMkeCGGUtQmdABbbi7PffgQyd624Dq9ohxFEvqg6j5McsIk0DWdpsZAH3ytMZrm1LlYTy02pqqQ6jThcBkJPs1lxmKRBC/EYIrqhK7kFpC87ZPwlZVHIhPMKKKoCilpZlLSzIw/h2Cpgroa4yGnuhofnvbwBK9p4K7x467xc2BPByiQnGLCmW4KtuKttqgddCXEiBDVCZ28gh6r9IrDwx7GaNe9VEHBBKNUQXurUS7YXeOjrtZHa/Mpw2h0aGrw09Tgp7wsNBbekW7G4QokeSkZLMQZieqEruQV0OO/fOVhdF2XERcRpCgK8Ykm4hNN5I8znrzztGu4a7ta8EYiP1VXuYJj5cay3d5GqlMNJni50SpOJz8/n0mTJuH3+8nPz+fJJ58kJSXltPvv3LmTkydPsnDhQgDeeOMN9u/fz3e+853hCnlIRHVCJyMbvzmOg/ZMjiRmU2dN4auH34RGN6Q6Ix2d6MZmV8nJjyMn30jw3k4dd52R4N21Phrq/D3Gwre3+2mv8FNVEbrRmuYMtOBdRp++WfrhBWCz2XjzzTcB46Ge5557jiVLlpx2/127drF9+/ZgQr/88su5/PLLhyXWoRTVCV0xmdByxvLDgjvxB6afu/bYeyRWHJaEPsJZ4hQysy1kZhszOvn9Oo1uP3WBVnx9rR+vt+eN1trAlH4Q6ofvasE7XFKyQMD555/Pnj17APj000958MEH6ejowGaz8bOf/YwxY8bw+OOP4/F4+OSTT4I1VbZv386KFSs4duwY3/ve96ivr8fhcLB69Wpyc3MjfFX9E9UJHSAubwy5bdUcTTRqPxxLyGRSxWGUqedHODJxJkzd69Gca9xoNSnJHDpQG2zFt7ed0sHWrR/+8AGjH94erxjDJV3GSJrkFFXmZh1GL69vGLJzX3NTap/7+P1+PvjgA26++WYAJkyYwAsvvIDZbOa9997jxz/+Mb/5zW+49957gwkcYP369cFz3H///dxwww3ceOON/OlPf+KBBx7gt7/97ZBc02CL+oROXgFjdpwMJvSjCVlMkhujUU9RFBwuKxpWCiYYD4q1t4X64d01Ppoae9YraG/TqTzqpTIws5PZQnDIpcNlkm6aGOXxeLjssss4ceIERUVFzJtnlARpampi6dKllJeXoyhKvwpzbd68mWeffRaAr3zlKzz22GNDGvtgivqEruQVMObjfwWXjyZkoVec/dx8YuSxx6vkjokjd0yoH76+zmi9u2v91PdSssDnDR8P372bJs1lJHl7vJQtiHZdfejt7e3ccsstPPfcc/yf//N/+OlPf8rcuXNZu3Ytx44d4/rrr490qEMq6hM6vVRdpLwC3duJYon7jANFtLPEKWRkW8gI9MNrmk5TvT+Y4N21Pedm7a2bxmpTAl00JhxOM8lp8lTrQPWnW+RUg1mcy2638+ijj3LnnXdy++2309zcTFaWUVL7z3/+c3C/xMREWlpaej1HcXExf//737n++uv529/+xpw5cwYltuHQr4S+detW1q1bh6ZpLFy4kOuuuy5se21tLWvWrKG1tRVN07jllluYNWvWUMTbg5KUwhhTqEj8kYRsdL8fKo5ICYBRRlUVUp1mUp2hB57a2zTjgaZaH/V1PpoaenbTdHh0qiq8wdE0qgopDiO5pwW6aWQCkOgxdepUzj33XF588UUWLVrE0qVL+fnPfx4c0QIwd+5c1qxZw2WXXdZjqOJjjz3GPffcw9NPPx28KRot+iyfq2kaS5Ys4f7778fpdPKDH/yAJUuWkJeXF9znmWeeYdy4cVx++eVUVFTwX//1X6xZs6bPNz+b8rndS4x6f/kYtyVdg8ds9LU+W/oYzhtuRZ3/hQGdfySKxrKqZ2sortnrNbpp6gMt+Ia6nvOz9saeoOJwGt00aU7TkI2Jj8Z/Zymfe+YiVj73wIEDZGVlBefdmzt3Lhs3bgxL6Iqi0NbWFgwgLS2tz0AHk2lMIeOOHWdP6jgADiXl4jx6cFhjENHBYlHIyLKQkWV00+iaTnOTcbO1vtZHfZ2f1pZebra2alS2asGbraoJUoOteCPJS+kCEWl9JnS3243TGRrT7XQ6KSsrC9vnhhtu4LHHHuP111+no6ODBx54oNdzbdiwgQ0bNgCwcuVKXC7XwII2m8OO7Zh2HuP2fBBK6Ik5XFh5BOcAzz8SnXrNo8FwXXN6BoyfEFpub/NRfcJD9QkPNSc81FZ34PefUpvGH6pNAx0AJCWbSc+0kZ5lIyPLhsNpRT3Dvvho/Hc+efIkZvPZ3Y472+OjUX+u2Wq1ntHvw6D8FD/88EPmz5/PNddcw/79+/nFL37BqlWrUE+Z66ykpISSkpLg8kA/Wp76sVRPS2dcS6j7pjwxB9/ed6mpqkKxWAb0HiNNNH4UP1uRvOaEZBiXDOMm2tD8Vpoa/Ljr/IFWfC9j4oHmJh/NTS0cKjNutqkmSEkzWvGpTqMvvq8Hn6Lx39nj8WAymQZ8vHS5nJ7H4+nx+3BWXS4Oh4O6urrgcl1dHQ6HI2yft99+m+XLlwMwceJEvF4vzc3Nn1lLYTApqU6K9EZmuvcxvrmSyY3lxmOFx4/A2Al9n0CIz6CaQjdbmWjcp+mqTVMfSPKN9T1LF2h+qK/1U18bGktps3eVIjYSfEqaCZM5ukfUqKqKz+cbla3soeTz+Xo0ivvS579AYWEhVVVVVFdX43A4KC0t5e677w7bx+VysXPnTubPn09FRQVer5fk5OQzi/4sjclI4cHta8PW6UcOokhCF0MgVJvGWNb8Oo0Nfurr/DTU+XDX+Wlv7dkX72kPH1GjKJCcaiLNaSLVYcZi6oy64nI2mw2Px0NHR8eA4rZarXR0dAxBZCNXX9es6zqqqmKz2c7ovH0mdJPJxJ133smKFSvQNI0FCxaQn5/P+vXrKSwspLi4mK9//es888wzvPrqqwAsXrx42H8hlbGF6Ds2ha88IjdGxfBQTaFJQCDUim9wB7pp3Eai95/y4JOuQ2O9n8Z6P9DJ1k+OYolTSHWEknyq04TVOnJvuCqKgt1uH/Dx0djNdLaG6pr79Rlp1qxZPcaV33TTTcHXeXl5PProo4Mb2RlSxhb2KKWrHzkQkViEAKMVn5WrkpUbevCpubGrFW882dpyap14jCdgQ0+3Gq24+ESVNIeJVKeZNIdJHn4SvYqdTq8x4V0rOqBUHkb3eVHMsXFjVEQ3tdtMTwR+Xb2dWqD1biT4xnqNDk/PJN/WotHWEho2qahGCYPurfjEJClhMNrFTkJPc9KYls2fXBdyKCkHHYWfbPklHD8GY8ZHOjohemWJU8nIUoPj4p1OJ0ePVAf74uvrjDIFp95w1bXwrhowCpGlOowJuY0uG3nCdbSJmYSuKApxufm84ZyDrqioukaHasF25ACKJHQRJRRFISHRREKiibyxRi0iv1+nqSHQinf7aHD7e07ph1GILFQv3mCzK6Ek7zSRmmbGEiet+FgVMwkdIH5MATm1tVTGZ6ApKkcSsjhHnhgVUc7U7YbruMAN184OjYb6UFdNg9tPZ0fPsfGedp0TlV5OVIbKxiYkqYFWvJHoU1Kjf+ikMMRUQlfGFjKu/AiV8RmAUQJgoox0ETEozhreVWMUItNpCCT3erePRre/x6gagNZmjdZmjcojoaGTSSmmYFdNqsOY8Fvmb40+MZXQGTuBcS2lfJA5E4DyxFw4tFlujIqYpygK8QkK8Qlx5Iwx1mmaTkuTRkOgm6bBbfTHn1qOT+9WUvjoIWOdaoKUVCO5p6QFbromyuxPI11sJXSHi/H+huDiocQco2Ox4jAUSCldMbqoqkJyqlEZsus2kt9n9MfXu/3BRN9bf7zmx3gKti5009VkNkoZpKYFumocJhISZWTNSBJTCV1RFManJwaXjyZm4VNULIf2oUhCFwKTWTGqQ7pCD0B5O3Ua60Ot+AZ377Vq/L7uBckMZgukpplJ6equSTNhT5AkHykxldABkseNI/2EmxqbA69q4UhCNhMO7YNLr450aEKMSJY4BVemBVdmqFuyw2M85dpYH2rJ95j9icDImmoftdW+sPOlpHV11xjfZZq/4RFzCV0Zfw5F+/dSYzMKiJUl51N4aF+EoxIiulhtKpk5Kpk5oZuunnY9LME3uP14O3smeW+n3mP4ZFc5g5Q0SfJDKeYSOuOKmNj8JqUZMwAoSx7DlXs/Rm9uREkanuqPQsQaRVGwxyvY40OlDLqm+Au15P00uv14vb0n+bDJugm15LNzIc7WSUqaiXjprjkrMZfQFVs8RXGdweVKe7rx4tB+mDE7QlEJEXuMkTUm4hNMwaqTuq7T1qrR6PbTUO8PfPfh8/Y8PtSSbwius1gUktOMvviu1nyClDTot5hL6ACFOaks3f0/FDUfI6vdqOWuH9qHIgldiCHV/UnXruGTXUm+qyXfGPjea0veq1NX7aOuW5981+ialFRToBaOicRkVcbJ9yImE7p1fBHz3n89bJ1eLv3oQkRC9ySfe0qSb6z30+mJo6qyxUjyvfTJh4+uMT59qyajOFlKt5Z8UopUoIzJhK6MP6dHKV3K96NrfhR14FNlCSEGR/ckb9QGDz3t2jWEsqnBf9qSBpqf4I3Z0DkhMVkNJHgzKalGmWGLZfQk+ZhM6GTlgT0e2ttC6zztRuXFvIKIhSWEOL3uT7tm5xnrukbXNDWEhlA2Nvjx9DJOXtehuVGjuVGj4nCo0z4+UQ101xhfyammmK1CGZMJXVFVGDcR/+5tVMZnUJ6Uw7yTnxr96JLQhYga3UfXdA2hBGOcfGODn6ZACeHGej+tLT2feIVQLfmuaf8ArDbjKdquvvnktNh46jUmEzoA487hmylX4bamAjC5oZyMQ/tg3hWRjUsIcdastvDiZGDcUO1qyTfW+2iq99PcpPWoXQPQ4ek5jNJkNuZ3TQmUS4jGfvmYTehq4SRyKqqCCb0sOZ90ecBIiJhlsSg4080400NlDfx+Y9q/xnp/MNk3NfRehdLvg/paP/W1p/TLJ6kkp4USfXKqCattZHbZxGxCZ1wRE1/+NzvTjLm+ypLGMPfQDvTWFpSExD4OFkLEApOpa4KPUKrTNZ2WFs2oMFnvpzGQ6Hu7+arr0Nyk0dwUKjcMRpdNV39819dIqEYZswldSUymSG0NLpclB558OLBHHjASYhRTVIWkZBNJyeHDKDs8eqglH0j2p+uX7/DoVFf5qK4Kddl0DaXsnuSTU4d3lE3MJnSAiTmhR/0PJuXhV1SUsp3ygJEQIoyiKNjsCjZ7+M1Xn1enqVuXTVODn6ZGP1ovXTa9DaUEiE9QA8ldDSZ5p7OXjv1BENMJ3VFUhGt3PbW2NDpMcRxOyKZw/65IhyWEiBJmi4LDZcbhOk2XTbd++d6qUQK0tWq0tWqcqOx+3hbOvyiejOzBnXgnphO6MnEKkz56mw9saQDsSSmg8OjH6J52FJs9wtEJIaJRb102YAylDCb5wPeW04yy8Xl17PGDf2M1thO6M4NzfTV8EFjenTqOqys/hEP7YPLMSIYmhIgxVptKepZKerehlH6/TkuTn6YGLdhd09Tgx+8zJusebDGd0AGmuGzB13tSxqEDetkuFEnoQoghZjIpgYJioXW6rpOU6KCltX7Q329kDqYcRPmFY0jwGiUAGuOSqLK70KUfXQgRIcYN2KGpKRXzLXTTOVOY9ulHNFoSObexHIvmg0P70L1eFMvg3pAQQohIivmETmYu//foSyjNjeHrD5dB0eTIxCSEEEMg5rtcFEVBmTilx3p9/84IRCOEEEOnXy30rVu3sm7dOjRNY+HChVx33XU99iktLeUvf/kLiqIwduxYlixZMtixDphSNBV9c2nYOr1M+tGFELGlz4SuaRpr167l/vvvx+l08oMf/IDi4mLy8vKC+1RVVfHiiy/y6KOPkpiYSGNj42eccfgpE6f0nPDiwF50vx/FJBNeCCFiQ58J/cCBA2RlZZGZmQnA3Llz2bhxY1hCf+utt7jiiitITDSKXqWkpPR6rojJHcMxx1jeTp3KnpQCzm08zO2HXjX60QsnRTo6IYQYFH0mdLfbjdPpDC47nU7KysrC9jl+/DgADzzwAJqmccMNNzBz5swe59qwYQMbNmwAYOXKlbhcroEFbTaf8bHbJs7i77aLAPCqxmXbj5aROOdzA4phuA3kmqOdXPPoINc8iOcdjJNomkZVVRU/+tGPcLvd/OhHP+Lxxx8nISEhbL+SkhJKSkqCy7W1tQN6P2MOwjM7dsKYDNSTGpqiciQxm1aTDTaV4llwzYBiGG4DueZoJ9c8Osg1n5mcnJzTbutzlIvD4aCuri64XFdXh8Ph6LFPcXExZrOZjIwMsrOzqaqqGlCwQyVh6gzGthgxaYrKntRxcHAfuqc9wpEJIcTg6DOhFxYWUlVVRXV1NT6fj9LSUoqLi8P2ueCCC9i1yxg10tTURFVVVbDPfcTIymV6e6jc2fbUCcYUJTLaRQgRI/rscjGZTNx5552sWLECTdNYsGAB+fn5rF+/nsLCQoqLi5kxYwbbtm3jnnvuQVVVbrvtNpKSkoYj/n5TFIXpaSb+HljeHpjJSN+9FWVa8ekPFEKIKNGvPvRZs2Yxa9assHU33XRT8LWiKNx+++3cfvvtgxvdIJt8Th7mAz58qpmjidk0WBJJ3b010mEJIcSgiPknRbuzT5nOxKajweUdaRPg+FH0BncEoxJCiMExqhK6kpzGdH/oznKw22XPtkiFJIQQg2ZUJXSAGZmhmYq2p00wniCVbhchRAwYdQm9aPJ4bL4OAHyKmWZLPPqebei9zRMlhBBRJPbL557CPHEqS//0MDkt1eS2VaMANLZB5RHIK4hwdEIIMXCjroWuWK1ckKaQ15XMA/QdmyMWkxBCDIZRl9CBXsed69s3RiASIYQYPKMzoU+f3XPlwb3oLU3DH4wQQgySUZnQycqF9CxO2tL4R+5cPsiYAbqGvnNLpCMTQogBG3U3RcF4svWTaVewUjOmpitqOsLnqrfB9o1w4fzIBieEEAM0OlvowJQphai6BsCBpHwaLQnou7ag+3wRjkwIIQZm1Cb0pMlTmNh8DABdUdniOAfaWuHg3ghHJoQQAzNqE7pitnB+XEtwebPzXEBGuwghoteoTegAxeNDU0BtdUzEr6iS0IUQUWtUJ/SCWdNwehoAaDPb2Zs8Fk5UoFcfj2xgQggxAKM6oaspDs73hqbK2+ycBID+6ceRCkkIIQZsVCd0gOIsW/B1sB9904eRCkcIIQZs1Cf06bOnYdG8ABxLyKLalgaHy9DrqiMcmRBCnJlRn9DtuXlMbQ/1mW92BLpdNpdGKiQhhBiQUfmk6KkWuHQyD5RyUc0OJjeWA6BvKYXLr4tsYEIIcQYkoQPzLpzM5157Knzlwb3o7loUh6v3g4QQYoQZ9V0uAErOGMjO77Fe//SjCEQjhBADIwk9QDl/bo91+mYZ7SKEiB6S0AO6Enqbycr7GTNpNtvhwB70BneEIxNCiP6RhN4lt4A/TL2BOy7+Easn38Inrimg6+ibPoh0ZEII0S+S0AMURSE5OxOfatwn/ih9GgD6x+9EMCohhOg/SejdXHT+xODr7WlFtJptcOQAetWxCEYlhBD9Iwm9m6zCAsZ31ADgU8184jRmNNI/+lckwxJCiH6RhH6Ki52h1+9lngeA/u930TUtQhEJIUT/SEI/xbwLzw2+3pE2AXdcMrhroGxXBKMSQoi+SUI/RUami6mdJwHQFJX3M2YC0u0ihBj5+pXQt27dypIlS/jud7/Liy++eNr9Pv74Y2688UYOHjw4WPFFxCV5oZK6wW6XLaXonR2RCkkIIfrUZ0LXNI21a9eyfPlyVq9ezYcffkhFRUWP/drb23nttdcoKioakkCH09yLphMXKKlbnpTLkYRMaG+TiS+EECNanwn9wIEDZGVlkZmZidlsZu7cuWzc2HPezfXr1/OlL30Ji8UyJIEOp8REO7P12uDy+xmBVvr7b0QqJCGE6FOf1RbdbjdOZ2joh9PppKysLGyfQ4cOUVtby6xZs3jppZdOe64NGzawYcMGAFauXInLNbBKhmazecDH9tfV5xfgeWsrl5zYwuy63cbKfTtI7WjDnDtmSN+7N8NxzSONXPPoINc8iOc92xNomsbvf/97Fi9e3Oe+JSUllJSUBJdra2s/Y+/Tc7lcAz62v86dXMAP//RTqAnvXnK/9CfUG+4c0vfuzXBc80gj1zw6yDWfmZycnNNu67PLxeFwUFdXF1yuq6vD4XAElz0eD8eOHePhhx/m29/+NmVlZfzkJz+J+hujiqKgzLuix3q99C10b2cEIhJCiM/WZ0IvLCykqqqK6upqfD4fpaWlFBcXB7fHx8ezdu1a1qxZw5o1aygqKuK+++6jsLBwSAMfDsrcS8F8yj2Blmb0LVInXQgx8vSZ0E0mE3feeScrVqzgnnvu4aKLLiI/P5/169ezadOm4YgxYpSEJJTii/ErKp84J/OJczIA+nuvRzgyIYToqV996LNmzWLWrFlh62666aZe933ooYfOOqiRpPz8K3lMmYvbmkpuWzWz63aj7N+FXnUMpZdZjoQQIlLkSdE+5Ew5B4/ZDkBlfAa7UscDoP/r1UiGJYQQPUhC74PdYmJeYntw+c3sOQDoH76F3tocqbCEEKIHSej9cMVF5wRff5Q+jSZLPHR2oL/3zwhGJYQQ4SSh98P47DQm0gQYddLfypoNgP72K+g+byRDE0KIIEno/XTFlMzg63/kXoxPUaHBjb7pwwhGJYQQIZLQ+2netHxSNA8AdbZUPkqfDoD+5t/RdT2SoQkhBCAJvd/iTCpX5VuDyy/lfR4d4OhB2L8zYnEJIUQXSehn4Mo5E4jTfAAcTM5nT0oBANo//jeCUQkhhEES+hlItVuYn2bcBI3zezkWH+hX3/0p+sG9EYxMCCEkoZ+xay+exM3Vpfz64xVcUfXv4Hrt1T9HMCohhJCEfsby0+zcOGcsyd628A07NqEfLuv9ICGEGAaS0AdAueASSM/qsV57ZX0EohFCCIMk9AFQzGaUL1wfXG41BSaV3vYJ+tHorgMvhIhektAHSLloASeyJvDUxK9w10U/pNaaAoD24h8jHJkQYrSShD5AitnCL6d/nQ05c/CYrbyQP9/YsGMT+r4dEY1NCDE6SUI/C9fPGRd8vSHnAtxxyQBo//ucPD0qhBh2ktDPwqz8ZIqsxvyiXtXCX8cuMDYcLpMaL0KIYScJ/SwoisJNFxYEl9/IvpDjdhcA+gu/l0qMQohhJQn9LBXnJjEl0Q+AXzXxh/FXGhtqTqC/81oEIxNCjDaS0M+SoijccfH44PLH6dPZmzwWAP2l59Gb6iMVmhBilJGEPggmuux8zhX6Uf6u8ItGJcb2VvT//V3E4hJCjC6S0AfJ1+YWYNY1APalFPBR+jQA9I/eRi/bHcnQhBCjhCT0QZKVFMdVhUnB5R2phcHX2v88je73RyIsIcQoIgl9EN04K4fxlg6W7XiOb5a9GNpQcRj97VciFpcQYnSQhD6IkqwmVn1lChckdvbYpr/43+jVxyMQlRBitJCEPshUkxn11kWgKOEbOjvRfvcLdE2LTGBCiJgnCX0IKOOKUOZ/AYAO1cyhxBxjw/5d6P96NXKBCSFimjnSAcQq5cu3s6m8jmcz59NmtvHzT1aR6m1B/9vv0Kedj5KRE+kQhRAxRlroQ8RrsfLriV/mpN1JsyWBX53zFWNsemcn2trV6D5fpEMUQsQYSehDJM6k8p3PjQkub3RN4Z3M842FQ/vQ/y5104UQg6tfXS5bt25l3bp1aJrGwoULue6668K2v/LKK7z11luYTCaSk5NZtGgR6enpQxFvVJmZncBV45P4x6FmANYWXcu0hgO4OhrRX/8r+qTpKFPOi3CUQohY0WcLXdM01q5dy/Lly1m9ejUffvghFRUVYfsUFBSwcuVKHn/8cS688EL+8Ic/DFnA0eb22dlkWY3a6G1mO09Ougk/xggYbe3P0Bul1osQYnD0mdAPHDhAVlYWmZmZmM1m5s6dy8aNG8P2mTp1KlarFYCioiLcbvfQRBuFbGaVpfPGohg96OxMm8BfCkqMjc2NaL/+ifSnCyEGRZ9dLm63G6fTGVx2Op2UlZWddv+3336bmTNn9rptw4YNbNiwAYCVK1ficrnOMFyD2Wwe8LGR8HkXfL1B43cbjU82fxm7kHMby5lRfwD278L69/8m+Zv/9zPPEW3XPBjkmkcHueZBPO9gnuy9997j0KFDPPTQQ71uLykpoaSkJLhcW1s7oPdxuVwDPjZSvjQhgS2HLOyo86IrKk+cezOrNv0cR2cT7a+/gMeVjXrJlac9Phqv+WzJNY8Ocs1nJifn9EOe++xycTgc1NXVBZfr6upwOBw99tu+fTsvvPAC9913HxaLZUCBxjKTqvCfl4wl1Wx0vTTGJQUrMgLozz+Dvn9npMITQsSAPhN6YWEhVVVVVFdX4/P5KC0tpbi4OGyf8vJyfvOb33DfffeRkpIyZMFGuzS7mf+8ZAxmdO7a/wJfrOw276jfj7bm/6EfPxq5AIUQUa3PLheTycSdd97JihUr0DSNBQsWkJ+fz/r16yksLKS4uJg//OEPeDwefvaznwHGx4nvf//7Qx58NJqelcCvvlSI629J6KfW6mprQXviIdRlP0FxjK4+RSHE2VN0Xdcj9ebHjw+s+mAs9LnpXi/aqh/Cwb09N+aMQb1vJUpCYnBVLFzzmZJrHh3kms/MWfWhi6GhWCyoi5dDRjbtpjhWTv06ewJzkXL8KNrPH0Jvb4tskEKIqCIJPYKU5FRav/0jHj7vW3zimsqK6XdSnpBtbCzfj/bEjySpCyH6TRJ6hLUmuTiZlgcYT5I+MuMbHLcH+s8P7ZOWuhCi3yShR1h2UhwPlRQQbwoNZ3xoxl1U2QMPcx3ci/azB9CkRIAQog+S0EeA8Q4bD146ljjFSOq1tjTun7mIivhAgbPDZbiXL0KvPRnBKIUQI50k9BHi3Ix4frhgTDCp11uTeWDmtziSkAWA//hRtJXfR68oj2SYQogRTBL6CDIzO4EHF47Bpoa6Xx6Y+S12pxQYOzS60VYuQ9/6ceSCFEKMWJLQR5hpmQk8VDI22KfeYonnoRn/wc7U8cYOHe1oa/4f2qt/JoKPEAghRiBJ6CPQuenxPHrZOFIDJXHGtlQxoelY2D76i39Ae3oleltLBCIUQoxEktBHqAlOGz+5ajznp1v4YfU/sWnenjtt+QjtkaXo5fuHP0AhxIgjCX0Ey0yM48HLC3Hc9whx588NrteBGmuqsVBXjfbj76O99ld0zR+ROIUQI4Mk9Cig2ONJ/cGPUb7wFQD+nn8JSy74TzZkzTbmQfL70f/2O7QfL0OvqvjMcwkhYpck9CihmEyoX76dnbf9kD+M/wIek5WnJt3Aj6d8nSZLvLHToX1ojywxWuu+XrpohBAxTRJ6lEmZPoOcxFDV40/Sp7K0+HuUpk8zWus+r9Faf2Qp+t7tEYtTCDH8JKFHmYI0Gz+7egJXTUgOrmuwJvP4lK+xcurt1FoDE4xUHUNbdT/aMz9BrzkRoWiFEMNJEnoUsppVvjknhwfm55FmCY1F3+iawt2z7+WvYxbQqRqteH3TB2gPLEb702/Qm5siFbIQYhhIQo9ixbmJ/PK6iVw5Pim4zmO28sfxX+Du2fdS3TUSxu9Df+tltB/chfbX36E3N0YmYCHEkJKEHuUS40wsuiiX/7psDGPiQ+uTvK24Ok5J3B3t6K//FW3ZN9D+vBa9rmZ4gxVCDClJ6DFickY8T3zpHL5ZnEGy6ufOqn+hEl4aoNqaioYCnR3ob/4dbfldaL/+KfrBvVJGQIgY0Ock0SJ6mFSFq85xsLAwlTjfA+hvvID+xovQ4cGvqDxw3rew+Tu59th7fL56K3GaD33j++gb34e8ApR5V6LMuQQlPiHSlyKEGACZJDpKDPSa9QY3+ut/5b3dVaw+58bg+uTOFhae2Mjlxz8m09Nt8gxLHMrMOSgXLYDJ56GYTIMR/oDIv/PoINd8Zj5rkmhpocc4JdWB8tW7aNlagX1XA+2Bf/KmuEReGLOAF/MvYXr9ARac2MSc2l1YvZ2hVntiMsp5F6IUXwwTp6GY5ddFiJFM/oeOElfPzGP+5Gz+uesEr+2po0aPA0BXVLY5JrLNMRG7z8Oc2p1cWfkRE5uPQUsT+vtvoL//BsQnoEw9H6bPRpk6CyUhqY93FEIMN0noo0hinImvnJfLdTNy2HS0kde2HGZrmxVdUQBoN9t4J6uYKQ2HjITeXVsr+ifvwSfvoSsqjC1EmTwT5dwZMO4cFKs1AlckhOhOEvooZFIV5hSkMqdgJjWtXv619QhvH26iinhMmp8LaneF7e9TVFZOvYNJjeXMrC9jXHMlpsNl6IfL0P/xFzCZYVwRyoTJKIXnwPhJKMmpkbk4IUYxSeijXHqChRsvnsANc3UO1LRQvquMpOkzYcdm8HYCUJaUzxbnJLY4J/E/fIF4XzuTG8qZ3FjOuY3ljGs5TtyBPegH9oQGSjozYOwElLGFKGMKIX8cJKeiBD4NCCEGnyR0AYCiKBRlJFGUMQsWzEL3tKPv2Iy++QO2NaSG7dtmtrPJNZlNrskAmDUf41oqKWo6xrSGg8yp3QV11VBXjb6lNJTkk1IgdyxKdj7k5BvfM3MgxSGJXohBIAld9Eqx2VFmfw5mf44vNLWTva2MbUfdbPMmUm9JDNvXp5opSx5LWfJY3NZkI6F3szd5LC1mO3lt1WTs3YEaqAIZTPRWG2RkQ3oWiisL0jNRnBn4Cieiq2YUWzxCiL5JQhd9ciTbufTz07kU0HWd40er2LmrnN3VbezzJ1BlcwT3LWyu7HH8K3mfozRjBgAWzUtWex3ZbbXktNeS4XGT4XGTWVuPq3IT1sBUezpQ13UCewKkOSHViZLqgJQ0SHEYXTjJqZCcCknJEJ+IosrDz2L0koQuzoiiKOSOzSF3bA5XYCT4psrj7NtTzsETTcy0NIOigq4Fj6mMzwi+9qoWjiVkcSwhq9fzL9n9PJdUfxq2rjRxPHEdPlKP1ZBysJxkbwtWzQcQXtxAUSExCRKSICHRGEcfnwDxiYGveLAnotjjwR4Ptniw28FmB6sdrFYUNXIPUglxtiShi7OiKAopeblckJfLBQBche5ph2Pl6EcOwNGDTGuvJrGxnWP2DJriEj/zfI7OniV+n5n4ZZot4eUIbL4OUrwtJHlbSfK2kehrJ9HbxvVH3yKtOTQNnw+F8qRc4n0ebP4O7P4ObP5OTttjHxcHcTajGyjOanxZbYH1VhSL8R1LHFgsxnezxXht7v5lRjFbjBFA5sCXyQQmS+B713K316oJzW41fn4mk/EHymSS+wui3/qV0Ldu3cq6devQNI2FCxdy3XXXhW33er388pe/5NChQyQlJbF06VIyMjJ6P5mIeYrNDkWTUYqMm6Z3AbrPByeP01JxjOPHa6lyt3KizcdJn5mTcSnUWlOptabi6mgIO5dXMfVI5mCUCfaYrZy0O8PWf+nYu2HLLZZ47jv/7vD4dA2bvxOrvxObFvgeWH5w+1rUzk5oMf6w1FpT+GfOhcR1eLFofiyaF4vmwaL5sOg+zJoPs+bHrPux+zs4p+lo8H10wKNaaLIkYNY1VN2PSdcCX8ZrVdfDiqj1Wv9SUUFVwaSCYjJeqyooSrfXavh6pfvrrmUFUELHBbd1++LU5cD7K4Rvh1P26dp2ymtAoZf9g8dAo82G5uno2rnbPt337basnLq+l4XT/RHs7x/H0+53uvOe2TmbbDY0jwfl6q+ipKT1L6Z+6DOha5rG2rVruf/++3E6nfzgBz+guLiYvLy84D5vv/02CQkJ/OIXv+DDDz/kj3/8I/fcc8+gBSmin2I2Q+4YknLHcA5wTmC9ruvQ6Iaak/iqq1CSLkNxn0Svq0FtqMPb0MiCqo3UW5NpiEuiyZJAoyUR/2m6RpJ87WHL7aaeDzzpikq72Ua72Ra23qz5elSorLOm8NexC/t1jekeN898vDJs3RbnJB6f8rXPPE4JJPlJjYd5ZNuvw7a9mX0B6wsuQ9U1FHTjD4Cuo6Kh6HpwnYLOTPc+vn7otbDjX8u5iA/SZ6Jg7NN1jALB1wRez6ndxeVV/w47/qW8z7M7dRzoXfsDuo4S+DEp3X5el5zcwuy6PWHH/++YS6lIyAi+R9h1d7134EPZ5cc/DvuDCPDHcVfQEHf6p5K7n/PaivfIawv/k7h2wrXByV56HHvK8k2H3yCtsyW47FVM/HbCtad971M6/Ljj4CvBrkCABksC6wsu63FUYXMlJSc2GjFceo1xT2iQ9JnQDxw4QFZWFpmZmQDMnTuXjRs3hiX0TZs2ccMNNwBw4YUX8tvf/hZd1+WjouiToiiQatzwtARa9F1cLhdUV7OkuREa6qC+1ig21nCQ1uZmmlraaWr30uzx0eLTadFM2PwdYefQFJVxzZW0m620m4yvTlNcr7HEaT0n1u5ULf2+FrOm9VjnV/ruk9cVFZ+ioik9b+i2m6y4u6YV7ENuW3WPdSftDvakjuvX8Xm9HH8wKY9PXFP7dfyE5ooeCX17WhE70wr7dfx0d1mPhP5R+nSOx6f36/jPV2/tkdDfziru8Yf7dK6ueD8softUE//MvahfxwJ87dBrQCiht5nt/DN3bo/9Lq7eGkzog63PhO52u3E6Qx9rnU4nZWVlp93HZDIRHx9Pc3MzycnJYftt2LCBDRs2ALBy5UrjP+xAgjabB3xstBqt15yekQEZGUBRn/vr3k605kb05ia05ia05kaSW1t4prUJveUkWlsrelsL3tZW2j2dtHs6aevw0uHz4+n04/P5epwz3VPPzeWv06la8KpmvKqZTtWCTzHhVc34VBM+xfie1tGz/9+ieXF56vErJvyKGvgyoSkqPlVF65bw1V4Kn/aW5E+nt+P1fvUFDJazK9yqnOXx0SgtLRXzIP6/HtaboiUlJZSUlASXB1o+Usptjg4Du2YFElKMr6z80+5lD3w5uq3TNQ06O6DTAx4PdHaQ09nBTR0e6OwEbwd6Z6fxBK3Pa3zv7ARfW2DZC67LwOcDnxfd5+NCv48LG/8Bfh/4/ca2rtd+H7rfj+bX0DQNXddQ4hPRfV7Q/KBpXHH8Iz5XvRUdBU1RAt+N1ryGgq4YaVBTVBJO6W4C+EJlKXNqd6KjogU+MXedp+s1gU6YTI+7x/HXVLzPRTXb6Uq3erdP3aHOG+PHXtBS1eP4Lx99m4UnPgn7w6J3HQDo3f7eTDyldQ5wS/nrtJlCLWz9NJ/6dRRy23regbjzwEv4evuU1Mt5Uru1zgEsmo//2P+3077fqSxaeIMgxdvCXftf6LFfdnvod7q+vgHFdma/42dVPtfhcFBXFxwRTF1dHQ6Ho9d9nE4nfr+ftrY2kpKkGp+ILoqqGkMYbXZIPs0+QxzDqX/EEnSdBF0DvxZM8sEvveu7HljnD73WddA1cjSNHHTQjGVjPca+EDiHDuiBbYGvwLaJOkwMrtNDx9N9HeGvCZ3jvGCju/t5u7XEdZ2kpCSam5uBeaFGfmCfi4PHdjsP4fuEXs/u8fMs6b4Q9gmmt08Ds8I2WYAvnPb4sA2B7zPC1iYCV53m7RITr6WlpdV4hmIQ9ZnQCwsLqaqqorq6GofDQWlpKXffHT5q4Pzzz+edd95h4sSJfPzxx0yZMkX6z4UYBIqiBEa2mDBSTHTpTxawu1y0jrJPn/EuF21DcM19JnSTycSdd97JihUr0DSNBQsWkJ+fz/r16yksLKS4uJhLL72UX/7yl3z3u98lMTGRpUuXDnqgQgghPptMQRcl5JpHB7nm0WGopqCTwhdCCBEjJKELIUSMkIQuhBAxQhK6EELEiIjeFBVCCDF4orKFvmzZskiHMOzkmkcHuebRYaiuOSoTuhBCiJ4koQshRIyIyoTevcDXaCHXPDrINY8OQ3XNclNUCCFiRFS20IUQQvQkCV0IIWLEsE5wMRj6mrA61tTW1rJmzRoaGhpQFIWSkhKuuuqqvg+McpqmsWzZMhwOx6gY1tba2srTTz/NsWPHUBSFRYsWMXHixEiHNaReeeUV3n77bRRFIT8/n8WLFxMX1/v0gNHqqaeeYsuWLaSkpLBq1SoAWlpaWL16NTU1NaSnp3PPPfeQmJg4KO8XVS30rgmrly9fzurVq/nwww+pqKiIdFhDymQy8bWvfY3Vq1ezYsUK/vnPf8b8NQP84x//IDc3N9JhDJt169Yxc+ZMnnjiCX7605/G/LW73W5ee+01Vq5cyapVq9A0jdLS0kiHNejmz5/P8uXLw9a9+OKLTJs2jSeffJJp06bx4osvDtr7RVVC7z5htdlsDk5YHcvS0tIYP348AHa7ndzcXNzunlOFxZK6ujq2bNnCwoULIx3KsGhra2PPnj1ceumlgDGXakJCQoSjGnqaptHZ2Ynf76ezs5O0tLRIhzToJk+e3KP1vXHjRi655BIALrnkkkHNYVHV5dKfCatjWXV1NeXl5UyYMCHSoQyp5557jttuu4329p5zZMai6upqkpOTeeqppzhy5Ajjx4/njjvuwGbr32z10cjhcHDNNdewaNEi4uLimDFjBjNmzOj7wBjQ2NgY/OOVmppKY2PjoJ07qlroo5nH42HVqlXccccdxMfHRzqcIbN582ZSUlKCn0pGA7/fT3l5OZdffjk/+clPsFqtg/oxfCRqaWlh48aNrFmzhmeeeQaPx8N7770X6bCGnaIogzpdZ1Ql9P5MWB2LfD4fq1at4vOf/zxz5syJdDhDat++fWzatIlvf/vbPPHEE+zcuZMnn3wy0mENKafTidPppKioCIALL7yQ8vLyCEc1tHbs2EFGRgbJycmYzWbmzJnD/v37Ix3WsEhJSaG+vh6A+vp6kpNPMyP5AERVQu8+YbXP56O0tJTi4uJIhzWkdF3n6aefJjc3l6uvvjrS4Qy5W265haeffpo1a9awdOlSpk6d2mNS8liTmpqK0+kMTsm4Y8cO8vLyIhzV0HK5XJSVldHR0YGu6+zYsSPmbwR3KS4u5t133wXg3XffZfbs2YN27qh7UnTLli387ne/C05Y/eUvfznSIQ2pvXv38uCDDzJmzJjgR7Obb76ZWbNmRTiyobdr1y5efvnlUTFs8fDhwzz99NP4fD4yMjJYvHjxoA1lG6n+/Oc/U1paislkoqCggG9961tYLJZIhzWonnjiCXbv3k1zczMpKSnceOONzJ49m9WrV1NbWzvowxajLqELIYToXVR1uQghhDg9SehCCBEjJKELIUSMkIQuhBAxQhK6EELECEnoQggRIyShCyFEjPj/KeWMKynkS/AAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Euler method\n",
    "dt, tmin, tmax = 0.1, 0.0, 10.0 # set \\Delta t,t0,tmax\n",
    "step = int((tmax-tmin)/dt)\n",
    "# create array t from tmin to tmax with equal interval dt \n",
    "t  = np.linspace(tmin, tmax, step)\n",
    "y  = np.zeros(step) # initialize array y as all 0\n",
    "ya = np.exp(-t) # analytical solution y = exp(-t)\n",
    "plt.plot(t, ya, label=\"Exact\", lw=5)\n",
    "\n",
    "y[0] = 1.0 # initial condition\n",
    "for i in range(step-1):\n",
    "    y[i+1] = y[i] - dt*y[i] # Euler method Eq.(A8)\n",
    "    \n",
    "plt.plot(t, y, ls=\"--\", lw=3, label=\"Numerical\") # plot y vs t (numerical)\n",
    "plt.plot(t, y/ya, lw=3, label=\"Ratio\") # plot y/ya vs. t\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<div id=\"leapfrog\" />\n",
    "\n",
    "### Leap-frog method\n",
    "\n",
    "- Central difference approximation\n",
    "\n",
    "$$\n",
    "\\left.\\frac{dy(t)}{dt}\\right|_{t=t_i}\\simeq\\frac{y_{i+1}-y_{i-1}}{2\\Delta t}\\tag{A11}\n",
    "$$\n",
    "\n",
    "- Difference equation, Substitute Eq.(A11) in Eq.(A1)\n",
    "\n",
    "$$\n",
    "y_{i+1} = y_{i-1}+2\\Delta t \\,f_i \\tag{A12}\n",
    "$$\n",
    "\n",
    "- Simulation procedure $\\rightarrow$ Explicit method\n",
    "\n",
    "$$\n",
    "y_{-1}, f_0 \\xrightarrow[\\ \\ \\ 2\\Delta t\\ \\ \\ ]{\\rm Eq.(A11)}\n",
    "y_1, f_2 \\xrightarrow[2\\Delta t]{\\rm Eq.(A11)}\n",
    "\\cdots\n",
    "y_i, f_{i+1} \n",
    "\\cdots\n",
    "\\xrightarrow[2\\Delta t]{\\rm Eq.(A11)}\n",
    "y_N  \\tag{A13}\n",
    "$$\n",
    "\n",
    "$$\n",
    "y_{0}, f_1 \\xrightarrow[2\\Delta t]{\\rm Eq.(A11)}\n",
    "y_2, f_3 \\cdots\n",
    "\\xrightarrow[2\\Delta t]{\\rm Eq.(A11)}\n",
    "y_{N-1}, f_{N} \\tag{A14}\n",
    "$$\n",
    "\n",
    "Instead of using the forward difference approximation of the Euler method, we can use a centered difference, to design a higher-order (more accurate) integration method.\n",
    "\n",
    "The central difference approximation given by Eq. (11) estimates the derivative of $y$ at step $i$, in terms of the $y$ values at steps $i-1$ and $i+1$. Substituting Eq.(A11) into Eq.(A1), we obtain the difference equation referred to as the Leapfrog method, shown in Eq. (12).\n",
    "\n",
    "This integration scheme gives a method that is accurate to second order in $\\Delta t$ (this means the error is $\\mathcal{O}(\\Delta t^3)$). The simulation procedure corresponding to the Leapfrog method is schematically represented in Eqs. (A13) and (A14).\n",
    "\n",
    "In contrast to the Euler method, the difference equation for the Leapfrog method (Eq. (A12)) involves values of $y$ and $f$ at distinct times. In this case $i$ and $i-1$, respectively. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<div id=\"runge2\" />\n",
    "\n",
    "### Runge-Kutta 2nd order\n",
    "\n",
    "- Difference equation\n",
    "\n",
    "$$\n",
    "y\\,'_{i+\\frac{1}{2}} = y_i+\\frac{1}{2}\\Delta t \\,f_i  \\hspace{50mm}{\\rm (Euler)} \\tag{A15}\n",
    "$$\n",
    "\n",
    "$$\n",
    "y_{i+1} = y_i+\\Delta t \\,f(t_{i+\\frac{1}{2}},y\\,'_{i+\\frac{1}{2}})\n",
    " = y_i+\\Delta t \\,f\\,'_{i+\\frac{1}{2}} \\hspace{15mm}{\\rm (Leapfrog)}\\tag{A16}\n",
    "$$\n",
    "\n",
    "- Simulation procedure $\\rightarrow$ Explicit method\n",
    "\n",
    "$$\n",
    "y_{0}, f_{0} \\xrightarrow[\\ \\ \\ \\frac{1}{2}\\Delta t\\ \\ \\ ]{\\rm Eq.(A15)}\n",
    "y\\,'_{\\frac{1}{2}}, f\\,'_{\\frac{1}{2}}\n",
    "\\hspace{15mm}\n",
    "y_i, f_{i} \n",
    "\\xrightarrow[\\frac{1}{2}\\Delta t]{\\rm Eq.(A15)}\n",
    "y\\,'_{i+\\frac{1}{2}}, f\\,'_{i+\\frac{1}{2}} \n",
    "\\hspace{60mm} \n",
    "$$\n",
    "$$\n",
    "y_{1}, f_{1} \\xrightarrow[\\ \\ \\ \\frac{1}{2}\\Delta t\\ \\ \\ ]{\\rm Eq.(A15)}\n",
    "y\\,'_{1+\\frac{1}{2}}, f\\,'_{1+\\frac{1}{2}} \n",
    "\\hspace{8mm} \\cdots\n",
    "\\xrightarrow[\\frac{1}{2}\\Delta t]{\\rm Eq.(A15)}\n",
    "y\\,'_{N-\\frac{1}{2}}, f\\,'_{N-\\frac{1}{2}} \\hspace{20mm}\\ \n",
    "$$\n",
    "\n",
    "$$\n",
    "y_{0}, f\\,'_{\\frac{1}{2}} \\xrightarrow[\\ \\ \\ \\Delta t\\ \\ \\ ]{\\rm Eq.(A16)}\n",
    "y_1, f\\,'_{1+\\frac{1}{2}} \\xrightarrow[\\Delta t]{\\rm Eq.(A16)}\n",
    "\\cdots\n",
    "y_i, f\\,'_{i+\\frac{1}{2}} \n",
    "\\cdots\n",
    "\\xrightarrow[\\Delta t]{\\rm Eq.(A16)}\n",
    "y_N  \\tag{A17}\n",
    "$$\n",
    "\n",
    "Here we have written down the difference equations for a 2nd order Runge-Kutta method. The first Eq.(A15) represents an intermediate operation, to estimate $y_{i+\\frac{1}{2}}$ and $f_{i+\\frac{1}{2}}$ at the midpoint, using the Euler method. The second equation represents the main operation, where $y_{i+1}$ is calculated using a Leapfrog like method, with $f$ evaluated at the midpoint. This particular method is second order in $\\Delta t$, like the Leapfrog method, but higher order Runge-Kutte methods can be easily derived."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAA3IElEQVR4nO3deXhU5dn48e8zW/ZtJpAACQRCABEFMYDgwhbQqihVwbq0+PKzfQt1gb7WqsWtSktVikVRXBDbahWtdcU1oqCkKovsCAHCEggEEsg+ycyc8/tjkkkmCWZhksnM3J/ryjVznrPM/US888x9znmO0nVdRwghRMAz+DsAIYQQviEJXQghgoQkdCGECBKS0IUQIkhIQhdCiCAhCV0IIYKEyZ8ffuTIkXbtl5iYyIkTJ3wcTdcmfQ4N0ufQcCZ97tmz52nXyQhdCCGChCR0IYQIEpLQhRAiSEhCF0KIICEJXQghgoQkdCGECBJ+vWxRCCFCgebSqanRqanWqanRKD9VTmSsjsGgfPo5ktCFEKINXK7axFyt46jRapN0XZvW4L2Oo8bd5nQ2PkoFWVNiiYiUhC6EEGdM13VcTmoTsDsROzzJuXGido+sHdU6LpdvPt9RoxMR6Ztj1ZGELoQIeLqu43A0TMgNkrDXsjtZ17VpWufEpxSYLQpLmPsnOjoc1QFnMCWhCyG6FE3T65NwowTcsA7taJCkHTU6nfXsNWUAi0VhsSjMYQqLxeBJ1BaLwhLmXjZbFGG17SazQqn68kpHTXcgCV0I0SF03V2eaK7W7GhQ5kCvoaK82pOonY7Oi9FgrE3OtYnZ7EnKtYnZUp+ozbVtJhNeybkrkYQuhGiRrunUNCppNJuk21XS8E0GN5nA3CgJu0fKjUfQ9W0mU9dMzO0lCV2IENLcqLk+EdeWLxpcqVG3zlHTuc+SNzdJyo1KG3VJucF2BmNwJef2kIQuRIBqXGt2NBwdn669E08EAhgM1CfkBqNnc4MyR2K3OOz2Mk+5w2xRXbak0dVJQhfCz3Rdx+lwJ1/dZed4oaNBMm46ina/79xaM4DZrLxqzGZLg5OAlsbr3KNpo7HlenNiYhQnTlR1Ui+CmyR0IXykSTmj4XXNnvdak1G19xUaZR0eZ+NRc31irh8hN07SZovy+V2NwvckoQvRjLq7AetGw47GNWVPovZe15nlDPCuNZstTWvLDZN1W0bNIjBJQhdBzeWqHwXXNKoxN03Smmc0rfnobsDWMprcl89FRllQBlejWnODUXSjUbMkZtGQJHQREDSX+07A050APN3I2dVkDo2O9WPljLr3XiUNi/cVGqH4fE3hO5LQRadqPGI+bWKu0dFclVRVOqjxQ2JWimaTsPdy3fv6OwOD7bpmEVgkoYt2cTkbj5i1JqPk+sSteZb9kZjNlsa15AZ3BZq9Sxh1l9IZu/DdgEKcjiT0ENZwtjmvhNzg1dHw8rkG6zq7xozCq0Thea27trm2tuxVa7YYMJklMYvQIQk9COha7UxzDu/Rsferd8Kue6938lUZnlnnvK7KaHBlRm27OUzRvXsCFZUlWCQxC9EqktC7EM+Jv5rGiVnj4N4iSk5VNj+K7uTbsqF+xjnvEbOhvrzRpP5cO2JuQykjMTEcTpR3cE+ECB6S0H1M13WcTmoT7WnKGF4lDc096VGL9WV7h8RbN9ucuUFSbly+aPgq1zIL0XVJQj8Nr9Fyk4TcfPmi6V1/ncdkBrPFgNncoIxhbnr3X8N2s1lhlKsyhAgaLSb0Z555ho0bNxIXF8fChQubrNd1neXLl/P9998TFhbG7Nmz6devX4cEqztqoOAQVVu+Rdu5DTVsJGrgOaffXtdxOvA6oeeVfOtqzo62jpY7RuMrMhq+j4uPwuG0NzNiltuyhRBuLSb0cePGcdlll7FkyZJm13///fccPXqUxYsXk5uby4svvsif/vQnnweqffgmNSvfpjSyF0fMUThMUTgNpThrqpqMkD3LDh38MFo2GmlapmiUgBvf8ddSfTkx0SY3nAghflSLCX3w4MEUFhaedv369eu55JJLUEoxYMAAKioqOHnyJAkJCT4NlLgESqJ6893599S3acCOat9+Ti3PaNmsTjtqNjc8CdhgnVHmZRZC+MEZ19CLi4tJTEz0LNtsNoqLi5tN6NnZ2WRnZwOwYMECr/1a4hhyHoff/aTN8ZnMirAwI2FhBizh7tewcCOWMIO7PdyAJaxxuwGzxdClTvqZTKY2/b6CgfQ5NEiffXhcnx/xR2RlZZGVleVZbksJQY+MIcZezDGtmmp07OiMP7KWqAmTsESFNztSbrm27Kr9cU8srQPVNe6friYU5/iQPocG6XPb9OzZ87TrzjihW61Wr8CKioqwWq1netgmlNlCTHw4W0q3cyDa3aFJBV8xIHIgasDpT4wKIUSoMJzpATIzM1mzZg26rrN7924iIyN9Xz+vpVL70rf8iGd5f3QP9EN5HfJZQggRaFocoT/55JPs2LGDsrIyfv3rXzN9+nScTvc1fZMnT+a8885j48aN3HHHHVgsFmbPnt1x0aakkfZdfQLfH90T8iWhCyEEtCKhz5kz50fXK6W49dZbfRXPj39WSl/SVuV4lvdH9UDP39Qpny2EEF3dGZdcOlVKGmkV9SWXQ1HJuA4fQnd19tR/QgjR9QRWQo+3EhtmxFpdAkCN0UyBJQ4Kj7SwoxBCBL+ASuhKKUjpS1p5gadNTowKIYRbQCV0AJWSRlrDK12iesiJUSGEIAATeuMRemG4FT3/gB8DEkKIriHgps9VqWmcc2oPD2x+gbTyAuId5RBv83dYQgjhdwGX0OmRSpzLzrCTufVtp4rQy0tR0bH+i0sIIfws4EouymzB2Kt30xX5+zs9FiGE6EoCLqEDmNP6N2nT5cSoECLEBWRCN9UmdLvBzK7Y3myNT4dD+/0blBBC+Fng1dABU1oGuTEp3DP8NnRloE95AYvy3/V3WEII4VcBO0JPqTzuWT4U1Z3qggJ0h8OPUQkhhH8FZEI3WhOJiI6iZ21S15SR/RHd4PB+/wYmhBB+FJAJHYA+6fQvy/cs7olJQT+w148BCSGEfwVsQld9+pPeIKHvjUmBA3v8GJEQQvhXUCV0XRK6ECKEBWxCp086fcuPYNA1APIju1NVUIDu6IJPeBZCiE4QsAldxVsJj40hteIYALoykBeZBDJRlxAiRAVsQgegSdklFf1A7o/sIIQQwSugE3rjOvqemBSQK12EECEqsBN6Wn/6l+UTX1NG5okdnFWSh75fTowKIUJTQN7679Ennf5lh1iW8wiqrs1oRK+pRlnC/BmZEEJ0usAeoccmoBIS65M5gMslU+kKIUJSQCd0APo0M5WuXI8uhAhBAZ/QVZ/0po2S0IUQISiwa+i4T4zaDSa+TDqfvTEpnAyL5Q/7V/k7LCGE6HStSuibNm1i+fLlaJrGxIkTmTp1qtf6EydOsGTJEioqKtA0jRtvvJHhw4d3RLxN9emPQddZ3v8qaoxmAE7uegubvRIVHtk5MQghRBfQYslF0zSWLVvGfffdx6JFi1i7di35+fle27z11luMHj2axx57jDlz5rBs2bIOC7gxFROH2ZbodT36rthUyJMbjIQQoaXFhL5nzx6Sk5NJSkrCZDIxZswY1q1b57WNUorKykoAKisrSUhI6JhoT0P1G8TA0v2e5d2xfdD37erUGIQQwt9aLLkUFxdjs9k8yzabjdxc79HvtGnTePTRR/n444+prq7m/vvvb/ZY2dnZZGdnA7BgwQISExPbF7TJ5LVv5dDzGfj2J57lH+L6YD60mYR2Hr8ratznUCB9Dg3SZx8e1xcHWbt2LePGjWPKlCns3r2bp556ioULF2IweH8ByMrKIisry7N84sSJdn1eYmKi1756UgoDS+on5dobk0LF9//Cefw4SqnmDhFwGvc5FEifQ4P0uW169ux52nUtllysVitFRUWe5aKiIqxWq9c2q1atYvTo0QAMGDAAh8NBWVlZu4Jtl15pxKsakqrccToMZvKIhWNHOi8GIYTwsxYTenp6OgUFBRQWFuJ0OsnJySEzM9Nrm8TERLZt2wZAfn4+DoeD2NjYjom4GcpkgrQMBpXs97TtiuuDvu+HTotBCCH8rcWSi9FoZObMmcyfPx9N0xg/fjypqamsWLGC9PR0MjMz+cUvfsFzzz3HypUrAZg9e3anlzpUv0EM2HKQ1cnnA7Artg/s3QVjJnZqHEII4S+tqqEPHz68yXXl119/ved9SkoKjzzyiG8jayOVPpBBX+d4lnfF9kHft9aPEQkhROcK+Fv/PfoNonfFMcKd1QAUhcdTdPwkelWlnwMTQojOEfC3/tdRsfEYE7szY+8HxDnKGVh6gISacsjbDYOH+Ts8IYTocEGT0AFU+iAu/eZLrzZ93w8oSehCiBAQPCUXgH6DmjTpe+WOUSFEaAiqhK7SBzZt3LcLXdM6PxghhOhkQZXQ6ZUGYeEAOJSRnbF9OKpZ4MhB/8YlhBCdILhq6EYjpA9iZUkk/+x3OTVGC9ftz+amXdtQKWn+Dk8IITpUcI3QATVgCAk1ZdQYLQBsj09H373Nz1EJIUTHC8qEPvjUPs9ybmwq1Xt2oeu6H6MSQoiOF3QJnbQM4nGQWnEUAKfBxG5DPBQc8m9cQgjRwYIuoSuzGfoN9Bqlb5OyixAiBARdQgd32WVIg4S+Pb4f7N7ux4iEEKLjBdVVLnXUwCEM/uR9z/Lu2N7YN/+bCF0PmgdeCNFV6LqO3W5H07R2/f917NgxqqurOyCyrqulPuu6jsFgIDw8vE2/06BM6PQbSIJWTa+KYxyOSnLX0VUcQ48dhuQUf0cnRFCx2+2YzWZMpvalE5PJhNFo9HFUXVtr+ux0OrHb7URERLT6uMFZcjFboN+AJmUXqaML4XuaprU7mYvTM5lMaG28yz0oEzq46+hnlzSqo++SOroQviZlzI7T1t9tUCf0waf2EeG0M7zoB0Yd34a+e5tcjy6ECFrB+z2p3yCsrir+sfYhjHqDry1SRxci6KSmpjJoUP1sq1dffTW33XabT469bds2jh07xsSJXf9xlkGb0FVYGPQdgHHPDq92fccmlCR0IYJKeHg4n332WYcce/v27WzZskUSur+ps4ehN5PQmXClfwISIoi5fnlV+/Zrw7bGF95r9balpaVcccUVLF++nP79+zN79mwuvPBCbrrpJu655x42b96M3W7niiuu4K677gJg06ZNPPDAA1RWVhIWFsZrr73GE088gd1u57vvvuO2227j6quvbmMPO09wJ/TB56G/+y+vNu2HbRicDpTJ7KeohBC+ZrfbmTRpkme5LvHOnz+fuXPncuutt1JSUsJNN90EwO9//3sSEhJwuVxcf/317Nixg/79+zNr1iyeffZZhg0bRllZGREREdx1111s2bKF+fPn+6t7rRbUCZ20/hAZxU5zN9Z0H8Ym60Cm7/+Mift2wYAh/o5OCOEjpyu5XHLJJXzwwQfcd999Xuvff/99Xn31VVwuF8eOHSM3NxelFN27d2fYsGEAxMTEdFb4PhO0V7kAKIMRzhrKjri+fNJrDMcibGy2DnCXXYQQQU/TNHJzc4mIiKCkpASAgwcP8txzz7FixQqys7OZOHEidrvdz5H6RlAndAA1eBjDind7ljcnZOCShC5ESHj++efJyMhgyZIl/Pa3v8XhcHhKKbGxsRw/fpwvvvgCgPT0dAoLC9m0aRMA5eXlOJ1OoqOjKS8v92MvWi+4Sy6AOmsYff/5LLE15ZRaoimxxHDgeDnpFWWoqMD7SiVEV9WWE5YNmUwmnE7nGX124xr6+PHjmT59Oq+99horV64kOjqaUaNG8be//Y277rqLIUOGcMkll9CzZ09GjBgBgMVi4dlnn2XevHnY7XbCw8NZsWIFY8aMYcmSJUyaNKnLnxRVuh/vtDly5Ei79ktMTOTEiROt3t71h1+zyDaWr5LOA+DmvR9y3eUjUZkXtevz/aGtfQ4G0ufAUFlZSWRkZLv390VCDzSt7XNzv9uePXue/rit+fBNmzaxfPlyNE1j4sSJTJ06tck2OTk5vPnmmyil6NOnD3feeWdrDt0p1NnDGLpztyehb7ZmcO2OTQGV0IUQoiUtJnRN01i2bBnz5s3DZrNx7733kpmZSUpK/c05BQUFvPPOOzzyyCNER0d7Tj50FWrweQxd+5RneWdcX6p2fEikTKcrhAgiLZ4U3bNnD8nJySQlJWEymRgzZgzr1q3z2ubzzz/n0ksvJTo6GoC4uLiOiba9Bp6DzVlO7/ICwP1Yui16HBzN93NgQgjhOy2O0IuLi7HZbJ5lm81Gbm6u1zZ1tfD7778fTdOYNm2a51rOhrKzs8nOzgZgwYIFJCYmti9ok6nN+xYPOofMop0cjO4BwLrEwUzYs52oc85rVwydrT19DnTS58Bw7NixM54+NxSn321Nn8PCwtr078Env0VN0ygoKODBBx+kuLiYBx98kCeeeIKoqCiv7bKyssjKyvIst/fkT3tOHGlnDWPEZ1/wnz4TANhgO4vSnLeouviydsXQ2QLxZNmZkj4Hhurq6jN6QIWcFD296urqJv8efuykaIslF6vVSlFRkWe5qKgIq9XaZJvMzExMJhPdu3enR48eFBQUtBhsZ1JDR9K/9BCxNe7rSROqSzl16DB6eamfIxNCCN9oMaGnp6dTUFBAYWEhTqeTnJwcMjMzvbYZOXIk27e7Hx5RWlpKQUEBSUlJHRNxeyWnYOyWxP/teJXn/zufhRv+hrW6BH3rBn9HJoQ4Q7169eLhhx/2LC9dupSFCxd2agybN2/m/vvvb9e+1113HZs3bz7jGFosuRiNRmbOnMn8+fPRNI3x48eTmprKihUrSE9PJzMzk6FDh7J582bmzp2LwWDg5ptv7nLzICilUENHck62980P+uZvYfR4P0UlhPCFsLAwPvroI26//fYmFYTO4HQ6GTp0KEOHDu30z26oVTX04cOHM3z4cK+266+/3vNeKcWMGTOYMWOGb6PzMTV0JHqjhM7279Fl9kUhAprRaOSmm27i+eef55577vFaN2fOHLKysrjySve02RkZGeTm5pKTk8PChQuJjY3lhx9+YMqUKQwaNIhly5Zht9tZtmwZaWlpFBUVcc8993D48GEAHn74YUaMGMHChQvZv38/Bw8epFevXtx8880sXbqUf/zjH1RUVDBv3jy2bNmCUoq5c+dyxRVXeKbtra6u5vLLL/dM2+sroXVquf9giIiCqor6NnsV7N4GgwPjahchurrXthzn9a1FLW8ITO4fx50Xpnq1Lfm2gE/31N/L8rNzbNxwbrcWj3XLLbeQlZXF7NmzWx3rjh07+PLLL4mPj2fMmDHccMMNrFy5khdffJGXXnqJP/7xjzzwwAP88pe/ZOTIkRw+fJgbb7yR1atXA5Cbm8vbb79NREQEOTk5nuM++eSTxMTE8PnnnwNw6tQpoH7aXqUU1157LTt27GDw4MGtjrclIZXQlcmEGjIcbd1XHIhKZr1tMCWWaG7d9B1KEroQAS0mJobrrruOZcuWERER0ap9hg4d6jnf16dPH8aOHQvAoEGDPAn6q6++Yvfu+gn+ysvLqahwDwonT57c7Gd99dVXPPPMM57l+Ph4oPlpeyWhn4mhIzm5eTO/HfFbAIyai59te4YYuWtUiIB36623ctlll3mVhE0mE5rmfq6wpmk4HA7POovF4nlvMBg8ywaDwXNZoaZpvP/++4SHhzf5vLbMYVM3be/KlStJTEzktttu8/m0vSGX0NWQ87E6F5Fels/emBRcBiPrDYmMz98PqX39HZ4QAe+Gc7u1qkRyOr8Z1YPfjOrRrn0TEhKYMmUKr732Gj/72c8ASElJYevWrVx11VV8+umnXgm9NcaOHcvy5cuZNWsW4H5o9JAhP/6AnEsuuYSXX36ZP/7xj4C75NJw2t7CwkK++OILRo8e3Y5enl7Qz4femIqKhoyzueD4Vk9bTrdz0dev9WNUQghf+d///V+Ki4s9yzfddBP//e9/ycrKYsOGDW2eGfKRRx5h8+bNZGVlMW7cOP75z3+2uM+dd95JSUkJEyZMICsri5ycHM4++2zPtL2zZs3yTNvrSyExfW5j2qoPOPLOW/xm1O8BMGlOlue+QMzDT3bZsksg3kF4pqTPgUGmz227jpo+N+RG6ABq+Bh62IvpW+a+DMlpMLFOS4D8/f4NTAghzkBoJvR4K2QM5sLj9XdmSdlFCBHoQjKhA6jzL2RM4RbP8ibrAMo3rsOPFSghhDgjoZvQh48hufok6WXuOdGdBhPf6VYpuwghAlboJvTassuYwvqyy9ru56Kv/9qPUQkhRPuFbEKH2rLL8fqyS0FEIo71/5WyixAiIIXcjUUNqeFjSHr9BW7Z8z5nleynf9khFMDBfdAn3d/hCSFaKTU1lUGDBuFyuUhNTWXx4sU/+ijMbdu2cezYMSZOnAjAp59+yu7du7nttts6K+QOEdoj9Nqyy1X5X5FRl8wB/Zsv/BqXEKJtwsPD+eyzz1i1ahXx8fG8/PLLP7r99u3bWbVqlWd58uTJAZ/MIcQTOoAacUmTNv3b1eghdqODEMHi/PPP5+jRowB8//33TJkyhcmTJ3PVVVexZ88eampqeOKJJ3jvvfeYNGkS7777LitWrOAPf/gDAIcOHWLatGlkZWUxffp0z7S5gSCkSy4AasTF6CtegIYJvKwEtn8PQ31/a64Qwez9Fac67NhTro9vcRuXy8XXX3/NDTfcAED//v15++23MZlMrFmzhr/85S+88MIL3HXXXWzZsoX58+cDsGLFCs8x5s2bx7Rp05g+fTqvv/46999/Py+99FKH9MnXJKFHRcPQkbAhhzJTBGu7D2VLQn/u+u/nmCWhCxEQ7HY7kyZN4ujRo2RkZHDJJe5v3qWlpcyZM4e8vDyUUq2amGvDhg28+OKLAFx77bU8+uijHRq7L4V8yQXAMHoiLhR3jvg/nh9wDd90O5edeYXoFeX+Dk0I0Qp1NfTvvvsOXdc9NfTHH3+cMWPGsGrVKl5++WWqq6v9G2gHC/kROgBnn4cxJpZRJ7bxSa8xAHzZbRhD1n2FGvcT/8YmRABpTVmkMV9OzhUREcEjjzzCzJkzmTFjBmVlZSQnJwPwxhtveLaLjo6mvLz5AVtmZibvvvsu1113Hf/5z38YNWqUT2LrDDJCp/ZJRqPGMv7oBk9bTrdzsX+7xo9RCSHaY8iQIZx11lm88847zJo1iz//+c9MnjzZ64/GmDFjyM3N9ZwUbejRRx9lxYoVZGVl8dZbb3nmNA8EITl9bnP0g/twPTKH20f+jiOR7sn5b9/5OhNv+yUquZdPP6s9AnFa1TMlfQ4MMn1u28n0uR0ttS+qVx8mHF3nafq05wXoX33qx6CEEKL1JKHXUkqhxkxkQsE6TJr7L+euuDTyvt+C7qjxc3RCCNEySegNqDETiNdruOD4Nk/bJ/EyT7oQP0bmPuo4bf3dSkJvQEXHojIv4tIj//W0rU4aTuWaz/wYlRBdm8FgCLkaeGdwOp0YDG1L0XLZYiNq3E8Y/M3dpFYc5VBUMnZTGKsro7k8Pw+V0tff4QnR5YSHh2O326murm7XM3nDwsKC/vrwxlrqs67rGAwGwsPD23TcViX0TZs2sXz5cjRNY+LEiUydOrXZ7b755hv++te/8uc//5n09ACdrbDfQFRKGpOPfMOyjKnE1pSjo9BXf4y6aZa/oxOiy1FKERER0e79A/HKnjPVUX1uMaFrmsayZcuYN28eNpuNe++9l8zMTFJSUry2q6qq4qOPPiIjI8PnQXYmpRRq7E8Y9/pyYhyVjD6+FbPuQi+KQL92Biq8/ZdnCSFER2qxQLNnzx6Sk5NJSkrCZDIxZswY1q1b12S7FStWcPXVV2M2mzsk0M6kLhhLlElxSeEmzLrL3VhdhZ6z6sd3FEIIP2oxoRcXF2Oz2TzLNpuN4uJir2327dvHiRMnGD58uO8j9AMVHom6YGyTdj37PXTN5YeIhBCiZWd8UlTTNP7xj38we/bsFrfNzs4mOzsbgAULFpCYmNiuzzSZTO3et7Wc19xM0eqPvdpKT5WSuncH4aPHd+hnN6cz+tzVSJ9Dg/TZh8dtaQOr1UpRUZFnuaioCKvV6lm22+0cOnSIhx9+GIBTp07x2GOPcffddzc5MZqVlUVWVpZnub0nBTrlJEpkLAwZjr5tI9vj+/Fu6li2x/XjubdXEJ9xTsd+djPkxFFokD6HhjPp84/d+t9iQk9PT6egoIDCwkKsVis5OTnccccdnvWRkZEsW7bMs/zQQw/x85//PHCvcmnAMPmnuLZtZFn/qzgQ7f4lflJtY/qenaj+Z/k5OiGE8NZiDd1oNDJz5kzmz5/P3LlzGT16NKmpqaxYsYL169d3Roz+M+hcVGpfrj5UP+vih70upPqz9/wYlBBCNK9VNfThw4c3OeF5/fXXN7vtQw89dMZBdRVKKdTkqVz40mJe6XcZxWHxnAqLZfV+B5MLj6C6n/6rjxBCdDa59b8FKvNizHHxXJlfP5/LW73H4/jwLT9GJYQQTUlCb4EymVBZV3HpkW+IdlQAcCzCxpf7TqEfP+rn6IQQop4k9FZQl1xKRLjFq5b+794TcKx8049RCSGEN0noraAiIlGTrubywzleo/TVMkoXQnQhktBbSU24splR+nicMkoXQnQRktBbqdlReriVnTv2oRcW+Dk6IYSQhN4mDUfp557M5YkNf+Psk3vR333V36EJIYQk9LaoG6VPPbSaBze/QN9y98hc/24Nel6un6MTQoQ6SehtpCZeiTE6hsbPZdH+/ZI8W1EI4VeS0NtIhUeirrqhSbsj9wfY/K0fIhJCCDdJ6O2gLpoMye4nNjmUkfdSLuaXo//AofffRZeH5Qoh/EQSejsokwnDtTMAeGbgdbzcfwqllmiWx41C+/LjFvYWQoiOIQm9vYaOhAFDuDL/a5SuAbDRNohvvlqPXnLSz8EJIUKRJPR2UkphmPY/pJcfZvKR+tr5S70vpfKNv/sxMiFEqJKEfgZUWgbqoknclPcxsTXlAJwIT+CNE+HoOzf7OTohRKiRhH6G1LUziA43MWPvSk/b+ykXk/fvN9EdDj9GJoQINZLQz5CKjkVdM4NxxzYw+NQ+AFwGI0ttF+H88N9+jk4IEUokofuAujALlT6IX+1+G6PmAmBXXBrvbSlAP7jXz9EJIUKFJHQfUAYDhptn0dt+nGkHPve0v5Y2iQOv/F1KL0KITiEJ3UdUSl/UZddyzcFVpJflA9DdfhLH8UL091/zc3RCiFDQqodEi9ZRV/4M05Z13L5zBat6ZHJD3ieEaU70j/+DPnQkKn2Qv0MUQgQxGaH7kDKbMfzPHHpXn+CWvSsJ02qnAdA1tBeeQK8s92+AQoigJgndx1Tvfqgrr2+6oqgQ18uLZUZGIUSHkYTeAdRl10FahmdZB7KTR7CwOh3X5x/4LzAhRFCThN4BlMmE4Ve/g4hINBSLB13PM4OmkdN9KO/+dw963m5/hyiECEKS0DuI6paMYcYdGNCJc9TXzl/peykb//kv9FNFfoxOCBGMJKF3IHX+GNT4K7hp38cMLNkPgKYMLOx9FfnPL0GvqfZvgEKIoNKqyxY3bdrE8uXL0TSNiRMnMnXqVK/1H3zwAZ9//jlGo5HY2FhmzZpFt27dOiLegKOmzcS8bxe/2/5P7j7/dorD4qkwR/KnhPE89o+lRP+/O1Cq8QPthBCi7VocoWuaxrJly7jvvvtYtGgRa9euJT8/32ubtLQ0FixYwBNPPMEFF1zAK6+80mEBBxplNmOYdS/WMAP3bv07FlcNAEciu/F4VV9qPnjDzxEKIYJFiwl9z549JCcnk5SUhMlkYsyYMaxbt85rmyFDhhAWFgZARkYGxcXFHRNtgFK2bhhm30e6/Ri3/fCmp32zdQBP5Wo45SlHQggfaLHkUlxcjM1m8yzbbDZyc3NPu/2qVasYNmxYs+uys7PJzs4GYMGCBSQmJrYxXDeTydTuff0m8WKqZt3DRU89yuG8bqzoOxmAr5LOI/bbr5mbvIWIiyacdveA7PMZkj6HBumzD4/ry4OtWbOGffv28dBDDzW7Pisri6ysLM/yiRMn2vU5iYmJ7d7Xr84dibr0GqZ/8h9OWWL4pNdoALJ7jOSKpYvp6XKhzj6v2V0Dts9nQPocGqTPbdOzZ8/Trmux5GK1Wikqqr/ErqioCKvV2mS7LVu28Pbbb3P33XdjNpvbFWgoUNf8AsMF47k19x1GF24hxlHBw5ueI7miEG3JfPTt3/s7RCFEgGoxoaenp1NQUEBhYSFOp5OcnBwyMzO9tsnLy+OFF17g7rvvJi4ursOCDQbKYEDNuB3jOeczZ+dr/GXDUwwoO+Re6ahBe/pR9K0b/BukECIgtVhyMRqNzJw5k/nz56NpGuPHjyc1NZUVK1aQnp5OZmYmr7zyCna7nb/+9a+A++vE73//+w4PPlApkwnD//4e85MPkrxnh/dKp4PjLz5F4i2zMJw3yj8BCiECktL9OFvUkSNH2rVfsNTc9MpytCcfggZTARyMTGLeebMYf2wDt1zQG9PYS4Hg6XNbSJ9Dg/S5bc6ohi46joqMxjDnYeg3EIBiSyyPnPv/KDdH8n7KxTy5oZjqd1+TGRqFEK0iCd3PVGSUO6n3P4toZyUZdfV03Jc0PnoomrLlT8s0AUKIFklC7wJURCSGOx/CMugc/m/7K1x6OMezbktCBne7zmPTg/eiF4fW11IhRNtIQu8iVHgEhtvnYbpgHL/KfYcb99XfPVoQmcjvEi/n26efQd+9zY9RCiG6MknoXYgymVEz52D4yXVcd3AVd23/J2G1c79UmcL5c/p0Xv7PWmreex1dc/k5WiFEVyMJvYtRSmG45heoGbcz5uRO/rRxCd3s9XPjvJM6lgcOxeFYeD968XE/RiqE6GokoXdRhosmYbjrT/Q1VfP4hsUML/rBsy6zaAfG3dvQHrodbc0nchWMEAKQhN6lqfRBGOYtJLZPH+7bupxf7P2AzBM7uPrQGvcGVZXo/1yCtugB9ONH/RusEMLvJKF3cSrehuGuPxEzbQZT87/i3m0vY8B7RJ6bX8T7z76C491/yeWNQoQwSegBQBmNRN/4Kwz/9ygq3ua1zoXi+QHXsCz9Sn5XkMSOBfPRN6yVMowQIUgSegBRA8/B8NBTqAsnetrWJJ3H3pgUAPJienFf+o08vuYQRx7/I/qurX6KVAjhD5LQA4yKisZwy50Y7nwIrIlcVLiZm/Z95Hm0HcDa7sO4o8e1vPDOt5z425/Qd2/3X8BCiE4jCT1AqSHDMTy8BMvkq7n28Fcs/u4JRhdu8ax3GkysTLmIWQlTeOHtHI4v/CP6tg1SihEiiElCD2AqPALDdbdgeHAx3dP78rsdr/DnjU9z1qk8zzY1RjMrUy5mVvefkvPqW2gP3oa2+mP0ajl5KkSwkYQeBFSPFAxzHsJw50MMjDXy6KZnuXfrctLL8j3b6CgGlB6EgkPorzyD9rtb0P61FP1Q3o8cWQgRSHz6TFHhP0opGDIcw+Bh6N+tYcT7r5O5YTEbrIN4My2LXpWFWGtK63eoquDQt+v57z47441FdBuRiRp5CSo2wX+dEEKcEUnoQUYZDKgLxqGPuBh9/ddkrnyD8zc+TY2h6X/qj3tewIcpF/G6rnHOD3sZu+YpRiXoRJ03EjX8AknuQgQYSehBShmNqFFj0UdcDFvWEZb9HjS4jNFuMPNl8vkA6MrAloQMtiRk8KzmZNjm3YzJ/huZsS5ihpyLOjcTeqW5vwUIIbosSehBThkMMGwUxmGj0A/uQ1/1Afq6rzA6nPxq99t8kZzJloT+6Mp9OsVpMLE+cTDrEwdj0F2clbefzHVvkVWVS/SAQXDWMNSgc1EJthY+WQjR2eSZogHCl33WKyvQv1uD/tWncHAvRZZY1iSdx9fdh5EX06vJ9krXWJ7zR2IdlfXHSEzCkDEY+g9G9RsIPVNRBqNP4qsj/51Dg/S5bX7smaIyQg9BKjIKNe4nMO4n6EcOkvjtan767Wp+umE1BRE2crqdy3eJZ5Mb2xuAfuVHvJI5wDZnNM9oIxm0bj8DPt/IwKqj9OkWjalPP+idjurdD5J6oYy+TfJCiNOThB7iVM/eqJ/+HH3qzbB/Dz035nDtxhyu3fgFJy3RbLCehUVzNNlve3w/jkXYOBZhY3VtLd7iqqHviSP0y9tLv/I1pFUVkhprIaxXCiSnoHqkQo9U6JaEMpk7u6tCBD1J6AKoveyxbwaqbwb6Nb+Ao4exbl1H1pb1sGdnk+13147eG6oxWtgVl8auuDRPm0F3MW3/51z/zSv1c0QqA05rd0xJyajEZEhMgsQklK0bWLtBbLy79i+EaBNJ6KIJpRT0SEH1SIHJP0WvtkPudvSdm9F3bYOD+7h368vsi+nF7tg+7Irtze7Y3pwIb3qZo6aMxDoqvBt1jUdTpnAwKpmksiKSC4tIsn9PN/tJuttP0s1RhjXCSHFiN7ToWIhLgDgrxMWjYuMhJh5i4iAmDmWWkb4QdSShixapsHAYcj5qSO1ljlWVhO39gUF7djAwbzdT8t6BqgpOmaPYF5PCnpgU9kf34EBUT45GWEmtPNbkmIcjunHKEsMpS4zXiL6h2JpybNUl/ObbN+lX7j6BXjfK3xXbm2hHFXHKSVSkBRUVA1ExqKhoiIqGiKja10iIiEJFREJ4JIRHuNvCIyAsXEo/IqhIQhdtpiIiYchw1JDhAOiaBseOkHBwL+cf3Mfwg3vhwDooL6XKaMHU6IHWNQYTpZaoFj+n1BJNqSUao+69v1MZuHf4bZ5lo+Yi2llJjMP9E3WiihhnJVGOIqKc+Uw99CVhmtOzvUMZORZhJcJZTQQuwkwGjGFhYAmDutfaH2W2gMUC5jAwm8FsAZPZ/d5k9nqvTCYwmmrb696bwGh0v697NRhr3xtr3xvcrwaj+3cpRDu1KqFv2rSJ5cuXo2kaEydOZOrUqV7rHQ4HTz/9NPv27SMmJoY5c+bQvXv3johXdEHKYKgv0YwaC+Ce1bH0FFH5+9GP5rvnkCnIh2NHsJQU86818ygKj+dYuJWjETaOhVs5Hp7A8fB4joclcDIsBk25r5CxVpd6fV6JOdpr2WUwUmKJocQS02x81xz8wmv5SGQic0f8n1ebxeUgTKshzFVDuKuGMM1BmKuG+Ioy7lr3qte2+ZHdWJM0HJPmxKI5MWtOTJoTs+7ErLkwaU5MugtT7R+ajAZz6gBUmMIpNUdh1FyYdA2D7n416hpVuoYBHaMCg8EAdT+qmfdKebcpVd9W914ZmmlTgAJD7atXO023a9xe9x4Fyvu9opltPf9Qarej/gWgJDwcrbqmSbvXQsNjtWqbVtwE57VNK2+aO91mbbzprjQ8HM1uR135M1Sc7+7IbjGha5rGsmXLmDdvHjabjXvvvZfMzExSUlI826xatYqoqCieeuop1q5dy6uvvsrcuXN9FqQIPEqp2tp3Aurs87zW6fYqDIUFJJ84RtKJo5x7/Bh68XEo2gn7j0NVJS5l4JQ5mqKwOKKdVV77OwwmMkoPUGqOptQcRZUp/LRxWFwOzI1G+FXGptvXGM3UGM2Umb2/Odjsp5psmx/ZnX/3mdikvTkZpQf4y8YlXm2rk4bzYsbUVu0/unALv9vxilfbm30m8GGvCzHoOgoNg667f9BQuo4B3f2qa1xS+D3XHPjSa/83+mTxvXUAQO0xan90wPNeRwETC9Zx0fHNXvv/K20y+6N7evarq4O5l2uPobuXf3I4h7NLvCeA+3u/KygOi/Usq9L6W2FUg9tiFDpTD66md6OS3fMZU6k2mGu3wWv7xm7I+9RrDiOHMp7md++9b91x/2fP+4Q1uMqrxBzF62mTm9m/6bFMusb/2/Oe15qCCBvb4/qRdXSd+3MmTHH/f+IjLSb0PXv2kJycTFJSEgBjxoxh3bp1Xgl9/fr1TJs2DYALLriAl156CV3X5VZx0SwVHgG9+0Hvfs0OeHR7JYaTxaRrDhIP5EHJSSgphtJT6KWnSC4r4S+5r0B5GegaDmWkzBxJmTmKMnMk5aYIKkwRlJsicZ3mapkelcexG8OoMoZRbTR77pRtLKyZSzYdhtbX3U3NlFBcqvXX5jeXpKqM4af9NtLYyVN7mrQVRNhOe96isXNONt1/V1wftiZktGr/zKKd0Cihr7cN4nBUUqv2H3d0Y5OEvibpPCpNEa3a/6pDa7wSuqYUn/Uc1ap9AX6+70OvfwNVxjA+6TW6VftaXI4mCf2UJYbN1gxPQve1FhN6cXExNlv9bd42m43c3NzTbmM0GomMjKSsrIzY2Fiv7bKzs8nOzgZgwYIFJCYmti9ok6nd+waqkOtzSm9MJhNhQ0ecdhNd09Ary9FKTtG9vBS9vBStrAS9ohytoty9XFmB3muc++7Yqkq0ygoGV1XwzLZn0e1VoLnQgRqDmSqjhWqjhWpD7avRjFFvmpB7VxRwQ97HOAxmagwmHMqE02DCYTDiVCacBiMOgwmXMpBa0fSEcLirmqSqIpzKiKYMuJQBp8GIhsGz7Kq969bQzI3cWhsGSqqZ+8D1tuzfzB8UvbXlidPsL+olJMRj8uH/1516UjQrK4usrCzPcntvfZVbhUNDq/scFun+sSW36rh1lVxd18HpgGo7EdXVRNTYoaYaqquhxg41NeiOGnebowYcDnDUkOZwkOZ0v8fpAGcVOB3oTge4XOB0QrUDXE73cmpf96vLBS4nk6r3MemHXPd6TQOX5n6va+5tammnSYfX7/+Mqw+tRseAphQaCk0Z0JW76KEpg3tfpYhpdIcvwPT92Vx6+L+g3MfXMaArT7HF633PquNN9r8x7xNK87+u/bz636pO/R+LuqSfUXaoyf4z9n1IRV3Zq8HfBh3v4wGkNHOF1C93v4Oz0TQTDf/INPyDlVBT5rWdSdf49a63Gu3bSIP9LS7vb2gxjkp+tfs/TWJqLg5DM4OBpKoiJhbUj85PnjyFCm/b/9dndOu/1WqlqKjIs1xUVITVam12G5vNhsvlorKykpiY1n0lFMJflFLuq1bMFog+zTadGxI2m40TxwtB0zC4NNBd7qRf96PrRGkaUZrm/gNQ24amATpounsfnQbrG7xHJ0WHFL1Be923AE+bO82j1/1Qv4zOYK9t6tfVf5nQGxyzQbqsfT+iYQrVISY6mrKyMmjU7nZRo3adcY2O14RXe6bXKhNwqdfxf+Q4ADQ8/6MTDfykyef9yO54nz+y1f5ER19FeXkFxMb/2M5t1mJCT09Pp6CggMLCQqxWKzk5Odxxxx1e25x//vl8+eWXDBgwgG+++Yazzz5b6udCtINSyj3JmcEYcBcVt/f/+IjERCpC7NtnZGIilR3Q5xb/yRiNRmbOnMn8+fPRNI3x48eTmprKihUrSE9PJzMzkwkTJvD0009z++23Ex0dzZw5c3weqBBCiB8n0+cGCOlzaJA+h4aOmj5XZkASQoggIQldCCGChCR0IYQIEpLQhRAiSPj1pKgQQgjfCcgR+j333OPvEDqd9Dk0SJ9DQ0f1OSATuhBCiKYkoQshRJAIyITecIKvUCF9Dg3S59DQUX2Wk6JCCBEkAnKELoQQoilJ6EIIESQCbILOlh9YHWxOnDjBkiVLOHXqFEopsrKyuPzyy/0dVofTNI177rkHq9UaEpe1VVRUsHTpUg4dOoRSilmzZjFgwAB/h9WhPvjgA1atWoVSitTUVGbPno3FYvF3WD71zDPPsHHjRuLi4li4cCEA5eXlLFq0iOPHj9OtWzfmzp1LdPRpJuRvo4Aaodc9sPq+++5j0aJFrF27lvz8/JZ3DGBGo5Gf//znLFq0iPnz5/PJJ58EfZ8BPvzwQ3r16uXvMDrN8uXLGTZsGE8++SSPP/540Pe9uLiYjz76iAULFrBw4UI0TSMnJ8ffYfncuHHjuO+++7za3nnnHc455xwWL17MOeecwzvvvOOzzwuohN7wgdUmk8nzwOpglpCQQL9+/QCIiIigV69eFBcX+zmqjlVUVMTGjRuZOHGiv0PpFJWVlezcuZMJEyYA7ufHRkVF+TmqjqdpGjU1NbhcLmpqakhISPB3SD43ePDgJqPvdevWMXbsWADGjh3r0xwWUCWX1jywOpgVFhaSl5dH//79/R1Kh3r55Ze5+eabqaqq8nconaKwsJDY2FieeeYZDhw4QL9+/bjlllsIDw/3d2gdxmq1MmXKFGbNmoXFYmHo0KEMHTrU32F1ipKSEs8fr/j4eEpKSnx27IAaoYcyu93OwoULueWWW4iMjPR3OB1mw4YNxMXFeb6VhAKXy0VeXh6TJ0/mscceIywszKdfw7ui8vJy1q1bx5IlS3juueew2+2sWbPG32F1OqWUTx/XGVAJvTUPrA5GTqeThQsXcvHFFzNq1Ch/h9Ohdu3axfr16/nNb37Dk08+ybZt21i8eLG/w+pQNpsNm81GRkYGABdccAF5eXl+jqpjbd26le7duxMbG4vJZGLUqFHs3r3b32F1iri4OE6ePAnAyZMniY2N9dmxAyqhN3xgtdPpJCcnh8zMzJZ3DGC6rrN06VJ69erFlVde6e9wOtyNN97I0qVLWbJkCXPmzGHIkCFNHkoebOLj47HZbJ5HMm7dupWUlBQ/R9WxEhMTyc3Npbq6Gl3X2bp1a9CfCK6TmZnJ6tWrAVi9ejUjRozw2bED7k7RjRs38ve//93zwOprrrnG3yF1qB9++IEHHniA3r17e76a3XDDDQwfPtzPkXW87du38/7774fEZYv79+9n6dKlOJ1OunfvzuzZs312KVtX9cYbb5CTk4PRaCQtLY1f//rXmM1mf4flU08++SQ7duygrKyMuLg4pk+fzogRI1i0aBEnTpzw+WWLAZfQhRBCNC+gSi5CCCFOTxK6EEIECUnoQggRJCShCyFEkJCELoQQQUISuhBCBAlJ6EIIEST+PxB7/+FNOEAyAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Runge-Kutta 2nd\n",
    "dt, tmin, tmax = 0.1, 0.0, 10.0 # set \\Delta t, t0, tmax\n",
    "step = int((tmax-tmin)/dt)\n",
    "# create array t from tmin to tmax with equal interval dt \n",
    "t = np.linspace(tmin, tmax, step)\n",
    "y = np.zeros(step)\n",
    "y1 = np.zeros(step)\n",
    "ya = np.exp(-t)\n",
    "plt.plot(t, ya, lw=5, label=\"Exact\")\n",
    "y[0] = 1.0\n",
    "\n",
    "for i in range(step-1):\n",
    "    y1[i] = y[i] - 0.5*dt*y[i] # Runge-Kutta Eq.(A15)\n",
    "    y[i+1] = y[i] - dt*y1[i] # Runge-Kutta Eq.(A16)\n",
    "    \n",
    "plt.plot(t, y, lw=3, ls=\"--\", label=\"Numerical\")\n",
    "plt.plot(t, y/ya, lw=3, label=\"Ratio\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<div id=\"runge4\" />\n",
    "\n",
    "### Runge-Kutta 4th order\n",
    "\n",
    "- Difference equation\n",
    "\n",
    "$$\n",
    "y\\,'_{i+\\frac{1}{2}}=y_i+\\frac{\\Delta t}{2}f_i,\\hspace{20mm}\n",
    "f\\,'_{i+\\frac{1}{2}}= f(y\\,'_{i+\\frac{1}{2}},t_{i+\\frac{1}{2}}) \\tag{A18}\n",
    "$$\n",
    "\n",
    "$$\n",
    "y\\,''_{i+\\frac{1}{2}}=y_i+\\frac{\\Delta t}{2}f\\,'_i,\\hspace{20mm}\n",
    "f\\,''_{i+\\frac{1}{2}}= f(y\\,''_{i+\\frac{1}{2}},t_{i+\\frac{1}{2}}) \\tag{A19}\n",
    "$$\n",
    "\n",
    "$$\n",
    "y\\,'''_{i+1}=y_i+{\\Delta t}f\\,''_{i+\\frac{1}{2}},\\hspace{20mm}\n",
    "f\\,'''_{i+1} = f(y\\,'''_{i+1},t_{i+1}) \\tag{A20}\n",
    "$$\n",
    "\n",
    "$$\n",
    "y_{i+1} = y_i+\\frac{1}{6}\\Delta t \\left[\n",
    "f_i+2f\\,'_{i+\\frac{1}{2}}+2f\\,''_{i+\\frac{1}{2}}+f\\,'''_{n+1}\n",
    "\\right]\\tag{A21}\n",
    "$$\n",
    "\n",
    "In fact, the most popular Runge-Kutta method is the 4th order method we present now. The appropriate difference equations for this scheme are given in Eqs. (A18)-(A21).\n",
    "\n",
    "The first three Eqs.(A18)-(A20) represent intermediate operations to estimate tentative values for $y$ and $f$ using the Euler method. The 4th equation is the main operation, where $y_{i+1}$ is calculated using the intermediate values of $f$ obtained previously. This method is fourth order in $\\Delta t$, which means that the error is $\\mathcal{O}(\\Delta t^5)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAA28ElEQVR4nO3dd3hUZfr/8fczJb3OJCSBhJYA0iRiEIhKkcDay9p111VW1wUb7Pq1rW11WVkVYVWsi7iu/hArFqwRBSSiFOkCoRMIhCSk15lzfn9MMmRSTGGSyczcr+uaK3POec6Z+4nyyZnnNKXruo4QQgivZ/B0AUIIIdxDAl0IIXyEBLoQQvgICXQhhPAREuhCCOEjJNCFEMJHmDz54YcPH+7QejExMeTn57u5mu5N+uwfpM/+4WT63LNnzxaXyR66EEL4CAl0IYTwERLoQgjhIyTQhRDCR0igCyGEj5BAF0IIH+HR0xaFEMLX6bqOrRZqazVqa3Rqa3TKissICdMxGJVbP0sCXQghWqHZdWprHWHs/FnTyrz697U6NLlJeTkZF0UQHCKBLoQQ7XJiL7k+aLUmAWyr1alp+LPBMrvd/TXV1ugEh7h3mxLoQgivYLc3vwdsa2GvuPW95K5jNIHZrDAHOF6hoUEYOuEIpgS6EKJLaJpj77fxUEXB0WKOH69yBrCthUDWOmEvua2UApNZERCgMDUIZnNz75uZZzC4Dq101u0OJNCFEG3S5OBe4z3lVvaS7baWtlzRJfU33ktuEsB1Px2BbXCGd0CAwmgCpdw73t0ZJNCF8BO67gjVJgfyGowhN3ugr8GwhicpRRvD2PEKMCtMDdo33kv2RRLoQniJ5gK5uRC2tRTKtTqefiS8uT5kGwRxeHgQdq2mxeGK+pA2Gr1jL9mTJNCF6CJtCeQ9pnxKiytchywa7EF7OpBNJpqOITfeO25hLNlkVs0Gsj/ePrezSKAL0Ua6rmOz0ba94mbmtS2Qqzq1D0YjTcO3QfCazL8eyP4wbOHNJNCF39A1nVpbg6BtMnaMaxA3DmybZ099AzAYaXb8uNkzL8xNg9vdVyaK7kUCXXgN59V6zQxV2BoHcDOBbGvxLIuu01ogR0WFUV1T4TjDovHwhQSyaIUEuugSDcePbY2HJloYU0avpLKitluch1zPeepbMwf3mh3GaNiuDYEcExNNfn436KjwShLook2cF4W0GMBNw7rx+/Yf0HP/LrXJXHeBSKOgbTx+7DKW3GBaxpBFdyaB7gea7B03PmjXIHxbGrpo+aKQrlN/tV793nBAsyEM5gBD3TSuy0wKJYEsfJgEuhfQ7DpVlXbKS+2tDlPYWthT9vTpbtBg/LjR8ETjoYr66ZjYKMrLS5zLvOVqPSE8RQK9kznPrGgwNGFrtIfsstfczN6xY+y42NNdqQtbXIckWgvnBvON7TygFxMTTH5+eSf1RgjfI4H+K5o779hlHLmFoYqGy7rDmRXQ+t5xk4BuNN9klr1jIbq7VgP9hRdeYP369URGRjJnzpwmy3VdZ+HChfz8888EBgYyffp0+vfv7/ZCdbsdjh5Cz9lHaWEe9p3bMJyVgRqZ3nz7VsaNWwpm13lu70aHKAUBAQaMDa/SM58YI24S0M2EtJzuJoTvazXQJ0yYwLnnnsv8+fObXf7zzz9z5MgRnn32WbKzs/nPf/7DP//5T7cXqn/0FpXLviEvJpVaUzA2U19s2YHYasobhTPd5iKQevVnVrgEbgtB3Nw8owliY2Pl8mghxK9qNdCHDBlCXl5ei8vXrl3LuHHjUEoxcOBAysvLOX78ONHR0W4tlF59KA/uwZbBN7nO39+5u9H15x03DuQWw9hlGZjkzAohRBc56TH0wsJCYmJinNNWq5XCwsJmAz0zM5PMzEwAZs+e7bJea2zDUsl5b2m7ajOaFAEBBscr0FB3j2Njg/eNlgWeWB4Q4JjXXc47NplM7fp9+QLps3+QPrtxu27f4q/IyMggIyPDOd2eIQQ9MJToqmNk28uoVIoaXeOqvZ8RftlVmCNCmz3A13oYa3WvBp8DVNc4Xt2JP96RTvrsH6TP7dOzZ88Wl510oFssFpfCCgoKsFgsJ7vZJpTJREhMBDnFv7A7IgmAy4q2k8R+VGKq2z9PCCG8zUk/pjQtLY0VK1ag6zo7d+4kJCTE/ePndVRiX/qW5zqn94UloOfs65TPEkIIb9PqHvq8efPYtm0bpaWl/PnPf+aqq67CVndy9ZQpUzjttNNYv349d955JwEBAUyfPr3zqk3sS5+DO52T+0MTQAJdCCGANgT6jBkzfnW5Uoqbb77ZXfX8+mf16kPf8uXO6X1hCeiH1nXJZwshRHd30kMuXSqxL33KjjgnD4bGYT+c47joSAgh/Jx3BXpkNOFBJizVRQDUGszkBkTB0UMeLUsIIboDrwp0pRQk9qVvg730faEJ6If2e7AqIYToHrwq0KHuTJeyw87p/WHxcmBUCCHwwkCnVx/6lDv20ENsldiVQU5dFEIIvPD2uSqxL6cXvMLLP/yTmOoiFIAl1tNlCSGEx3nfHnpCb0L0WmLrwxyg8Bh6RZknqxJCCI/zukBXgYEYExKbLjh0oOuLEUKIbsTrAh3A1CelyTwZRxdC+DuvG0MHMPVJpuyH5eSE9mBfaAKBWi1nSqALIfyc1wb6L5H9eDT1TwCklBwgPedrD1clhBCe5ZVDLua+KfRpcNfFA6EJ2HMOoGtyCwAhhP/yykA39Egg0gwxVccBqDGaOWiOhCNyCwAhhP/yykBXSkHvZPqXngjw3eGJ6Pt3e7AqIYTwLK8MdADVJ5mU0hzn9J7wXrB/lwcrEkIIz/LaQKdPCv3LTgT67rBE9AOyhy6E8F9eG+iqdzLJDYZc9oX1xH5gP7qm/cpaQgjhu7w20ImNJ9Kkux4YNYXD0cOtrCiEEL7JawNdGQzNHxiVYRchhJ/y2kAHUL37k1wmB0aFEAK89EpRpz4pnLJmCcOO76Z/WQ6phTvRVaCnqxJCCI/w6kBXfZIZXrSb4UUNhlmqQtA1zTEkI4QQfsS7U69HTwgMdp1XWQHHjjTfXgghfJhXB7rjwGi/JvN1GUcXQvghrw50ANXMvdGRM12EEH7Iq8fQAeiTTGFABB/2Hs+esEQCtFoe2f+Tp6sSQogu5/WBrnonY9A1liaeDUCAvRbbusUY5MCoEMLPtCnQN2zYwMKFC9E0jUmTJnHppZe6LM/Pz2f+/PmUl5ejaRrXXXcdI0eO7Ix6m4rvRZTBTmxVIceCLNQYzewzRjLw6CFISOqaGoQQohtodRdW0zQWLFjAAw88wNy5c1m1ahU5OTkubd5//33Gjh3Lk08+yYwZM1iwYEGnFdyYMhih3wAGFZ94SPTOiN7oe3Z0WQ1CCNEdtBrou3btIj4+nri4OEwmE+np6axZs8aljVKKiooKACoqKoiOju6calug+g9iYMl+5/SOiD4ggS6E8DOtDrkUFhZitVqd01arlezsbJc2V155Jf/4xz/44osvqK6u5qGHHmp2W5mZmWRmZgIwe/ZsYmJiOla0yeSybnXqKAatfM45vTOiN8YDa7B2cPvdUeM++wPps3+QPrtxu+7YyKpVq5gwYQIXXXQRO3fu5LnnnmPOnDkYGh2UzMjIICMjwzmdn5/foc+LiYlxWVe3xtO37DBmrZZag5mjwVbyDx9FyzmACgrpWKe6mcZ99gfSZ/8gfW6fnj17tris1SEXi8VCQUGBc7qgoACLxeLSZtmyZYwdOxaAgQMHUltbS2lpaYeK7QgVEYU5Jtblzos7w5Ngn1xgJITwH60GenJyMrm5ueTl5WGz2cjKyiItLc2lTUxMDFu2bAEgJyeH2tpaIiIiOqfiFqh+gxhYIgdGhRD+q9UhF6PRyNSpU5k1axaapjFx4kSSkpJYvHgxycnJpKWlccMNN/Dyyy+zdOlSAKZPn+54kHNX6j+IQXt/4JO6yZ0RfdD3/Ny1NQghhAe1aQx95MiRTc4rv/rqq53vExMTefzxx91bWTup/oMY9MHbzunDIVb0LTvQdb3r/7gIIYQHeP2Vok5JfbFqldy1bRHJZYfoWXEMhQ75RyE23tPVCSFEp/OZQFcmM/RJZvxu12EWfc8OlAS6EMIP+NTNTlT/QU1n7t3Z9YUIIYQH+Hygy5kuQgh/4VOBTl2g60BusJXM+FGUHs5Fr672bF1CCNEFfGYMHYDoGIiyMjvxQtbEDAUgdMsbnLl3B5xyqoeLE0KIzuVTe+hKKdTAoSSVH3XO2xbVH33nVg9WJYQQXcOnAh2AAUMZWrTHObktsh96tgS6EML3+Vygq4FDOaVkPwZdA2BfWAJl+/eh22o9XJkQQnQunwt0EpIIDg6kX92NunRlYHtIL9gvD44WQvg2nwt0pRQMGMLQ4hPDLlsjZRxdCOH7fC7QwTHsMqRor3N6W5SMowshfJ9vBvqAYQwp3ouqG0ffHd6Lyj270TW7hysTQojO45OBTlJfwkzQu/wIAJoysiMgBnL2t7KiEEJ4L9+6sKiOMhghZQhDi/awP8zxuKZtUf05LXsrqnd/D1cnhG/RdZ2qqio0TevQraqPHj1KtZ9dzd1an3Vdx2AwEBQU1K7fqU8GOjjG0YetWMeOiD4MKd7LaYU70HeWwqSLPF2aED6lqqoKs9mMydSxODGZTBiNRjdX1b21pc82m42qqiqCg4Pbvt2TLay7UgOGMuaDNxiT3+BgaPZxeeCFEG6maVqHw1y0zGQytfubi2+OoQP0TQFzgOu80mI4JOPoQriT7CB1nvb+bn020JXJDCmDm8zXf9nogWqEEKLz+WygA6jBI1ymqw1m9G0bPFOMEKLTJCUlMXnyZOfr+eefd9u2t2zZwjfffOO27XUmnx74UkNS0T94gw+SJrDeego7I3rzyrqnsdhqHXvwQgifEBQUxNdff90p2966dSubNm1i0qRJnbJ9d/LpQCepH4SGs846mF+i+gGwKbQ3E3bvgEHDPFycEL7FfsvFHVuvHW2Nr37c5rYlJSVccMEFLFy4kJSUFKZPn86ZZ57J9ddfz3333cfGjRupqqriggsu4O677wZgw4YNPPzww1RUVBAYGMiiRYt4+umnqaqq4qeffuL222/nkksuaWcPu45PB7oyGFGnnMqI/GxnoG+MHsD4bRtQEuhC+IyqqiomT57snK4P3lmzZjFz5kxuvvlmiouLuf766wG49957iY6Oxm63c/XVV7Nt2zZSUlKYNm0aL774IqmpqZSWlhIcHMzdd9/Npk2bmDVrlqe612Y+HegADEllxEef8Ha/KQBstAxA++U9DJf9zsOFCSHcpaUhl3HjxvHpp5/ywAMPuCz/5JNPeOutt7Db7Rw9epTs7GyUUvTo0YPU1FQAwsPDu6p8t/Hpg6LgODCaUppDiK0SgMLAKA7lFaOXl3m4MiFEZ9M0jezsbIKDgykuLgbgwIEDvPzyyyxevJjMzEwmTZpEVVWVhyt1D98P9Nh4jDE9GHb8xP3QN0alwI7NHqxKCNEVXnnlFQYMGMD8+fP5y1/+Qm1trXMoJSIigmPHjvHtt98CkJycTF5eHhs2bACgrKwMm81GWFgYZWXesQPYpiGXDRs2sHDhQjRNY9KkSVx66aVN2mRlZfHuu++ilKJPnz7cdddd7q61w9TgVEbsyuanWMe4+cboAVz4ywbUyLEerkwI39GeA5YNmUwmbDbbSX124zH0iRMnctVVV7Fo0SKWLl1KWFgYo0eP5t///jd33303w4YNY9y4cfTs2ZNRo0YBEBAQwIsvvsiDDz5IVVUVQUFBLF68mPT0dObPn8/kyZO9/6CopmksWLCABx98EKvVyv33309aWhqJiYnONrm5uSxZsoTHH3+csLAw51eb7kINSWXEmv84p7dEJVO77WsCPViTEMJ9Dh482Oz85cuXO98/+uijzvfz5s1rtn1qaiqffvqpy7zQ0FA+++yzk66xK7Q65LJr1y7i4+OJi4vDZDKRnp7OmjVrXNp88803/OY3vyEsLAyAyMjIzqm2o045lYSqAnpUFgJQZQpke00g+rEjHi5MCCHcp9VALywsxGq1OqetViuFhYUubQ4fPkxubi4PPfQQf/vb35xjUN2FCg1D9R3AaYXbnfPWWQajb17rwaqEEMK93HLaoqZp5Obm8sgjj1BYWMgjjzzC008/TWhoqEu7zMxMMjMzAZg9ezYxMTEd+jyTydTudctGn03aV6v4slc64bXlmHUb5u0bib7qxg7V0NU60mdvJ332DkePHj3puy36490a29LnwMDAdv3/0OoWLRYLBQUFzumCggIsFkuTNgMGDMBkMtGjRw8SEhLIzc0lJSXFpV1GRgYZGRnO6fz8/DYX2lBMTEy719VThjLsnf8ya/18BpYcwIhOjcnEsZyDqKC232/YUzrSZ28nffYO1dXVJ3U/c3ccFPU2be1zdXV1k/8fevbs2WL7VodckpOTyc3NJS8vD5vNRlZWFmlpaS5tzjjjDLZuddx3vKSkhNzcXOLi4lottksl9ScwPILBJfsxojvm2Wwgd18UQviIVvfQjUYjU6dOZdasWWiaxsSJE0lKSmLx4sUkJyeTlpbGiBEj2LhxIzNnzsRgMPC73/2u211lpQwG1Klp6Cu/cpmvb1qDOm2Mh6oSQgj3Ubqu65768MOHD3dovY5+LdU3rEab/0/XmZHRGJ5ciDJ072usvPGr+MmSPnuHiooKQkJCOry+O4ZcevXqxZ/+9CceeeQRAF566SXKy8v561//elLbbY+NGzfy3nvv8fjjj7fatnGfr7jiCh566CFGjHC95Xdzv9uTGnLxKaeMAJOZUlMIy+NO4+kh17PGFA8Hdre+rhCi2woMDOTzzz9vcgZeV7HZbIwYMaJNYd6Z/CrQVVAwDBrG0sQz+ffga8nqMYIfYoejb5LTF4XwZkajkeuvv55XXnmlybIZM2a4XCw0YMAAwHF1++WXX85NN93E2LFj+ec//8kHH3zABRdcwKRJk9i3bx/gOBHklltu4fzzz+f88893XoczZ84c7rjjDi655BLuvPNOsrKyuOGGGwAoLy9n5syZTJo0iYyMDJYuXQrAfffdx3nnnce4ceN4+umn3f578LtzhdSpo0j79DPe6eu4THi95RRsmxYRcPG1Hq5MCN+waNMx3t5c0HpDYEpKJHedmeQyb/6PuXy168TV5tcMt3LtqbGtbuvGG28kIyOD6dOnt7nWbdu28d133xEVFUV6ejrXXnstS5cu5T//+Q+vvfYajz32GA8//DC33HILZ5xxBocOHeK6665zXoGanZ3Nhx9+SHBwMFlZWc7tzps3j/DwcOeTjoqKioATt+1VSnH55Zezbds2hgwZ0uZ6W+N/gT48jf6LXiWquoSiwAhKAsLYWVTL0KICVJS19Q0IIbql8PBwrrjiChYsWEBwcNtORR4xYoTzjLw+ffowfvx4AE455RRnQK9cuZKdO3c61ykrK6O8vByAKVOmNPtZK1eu5IUXXnBOR0VFAc3ftlcC/SSo2HgMPZMYVbCNr3s6zm7Jih3OkJ9/RE0838PVCSFOxs0338y5557L1Vdf7ZxnMpnQNA1wXARZW1vrXBYQEOB8bzAYnNMGg8F50FLTND755BOCgoKafF57DgbX37Z36dKlxMTEcPvtt7v9tr1+F+gAKnU0Y3/Y4Az01THDuWndJxgk0IU4adeeGtumIZKW3DY6gdtGJ3Ro3ejoaC666CIWLVrENddcA0BiYiKbN2/m4osv5quvvnIJ9LYYP348CxcuZNq0aYDjodHDhv36E8/GjRvH66+/zmOPPQY4hlwa3rY3Ly+Pb7/9lrFj3XvHV786KFpPnZ7OsKLdhNU6vjYVBEWRnVuCXlLk2cKEECft1ltvdTnb5frrr+eHH34gIyODdevWtfsUy8cff5yNGzeSkZHBhAkT+N///tfqOnfddRfFxcWcc845ZGRkkJWVxdChQ5237Z02bZrztr3u5FfnodfTdR3tb7cy33IW3yScAcDFB5dzU3ofDOPO7fB2O5M3np98sqTP3qE7nIfubdraZzkPvQ2UUqiR6Yw9tsk574fY4Wjrsn5lLSGE6N78MtAB1OlnMvz4bkJrKwA4FmRhV04helmJhysTQoiO8dtAp28K5mgLZxRsA8Cs1XIwJBZ9w48eLkwIITrGL89ygbphl9PTOXd1FqcV7uD0gu0E26vR1wFnTW51fSGE6G78dw8dx7DLgNIczsrbSLC92jHzl43o5d7xhG8hhGjIrwOdfgOh8dWhdhv6ejk4KoTwPn4d6MpgQJ2e3mS+vvpbD1QjhOiopKQkJk+ezDnnnMMf/vAHiouLf7X9li1bnPdZAfjqq694/vnnO7vMTufXgQ6gRk8AQEOxOSqZ5wddyaa8avSCPM8WJoRos6CgIL7++muWLVtGVFQUr7/++q+237p1K8uWLXNOT5kyhdtvv72Tq+x8fntQ1KlvCsT34p2gIc47MNYaTIxY/R3qgqs8XJwQor1OP/10fvnlFwB+/vlnHn74YaqrqwkKCuKZZ56hd+/ePP3001RVVfHTTz8576myadMmZs2axcGDB/nLX/7C8ePHsVgszJ07l169enm4V23j94GulEKNmcjozG+cgf5jzFAqflpA6PlXopTycIVCeI9PFhd12rYvujqq1TZ2u53vv/+ea6913A47JSWFDz/8EJPJxIoVK/jXv/7Fq6++yt133+0McIDFixc7t/Hggw9y5ZVXctVVV/H222/z0EMP8dprr3VKn9zN74dcANTo8fQty6VPWS4ANcYAftAssH+XhysTQrRFVVUVkydPJjU1lfz8fMaNGwc4Hlp/6623cs455/D3v/+dHTt2tLqtdevWcdlllwFw+eWX89NPP3Vq7e4kgQ6omDjUgCGMP7reOW953Ej01d95righRJvVj6H/9NNP6LruHEN/6qmnSE9PZ9myZbz++utUV1d7ttBO5vdDLvXUmImcvfh//K//eejKwJaoZI79/Dk9rrChTPJrEqIt2jIs0pg7b84VHBzM448/ztSpU/nDH/5AaWkp8fHxALzzzjvOdmFhYZSVNX+9SVpaGh999BFXXHEFH3zwAaNHj3ZLbV1B9tDrqLQzsWqVDD/uGGbRlYFvwwfClnUerkwI0R7Dhg1j8ODBLFmyhGnTpvHEE08wZcoUlz8a6enpZGdnM3nyZD766COX9f/xj3+wePFiMjIyeP/99533NPcGfnn73JbYX5rNioMVzBtyHQCxVYW8UPENAXc85PbPai9vvK3qyZI+ewe5fW77ye1zu4DhzAzGHttMeN2DL44FWdh4qBS98JiHKxNCiNZJoDc09DTM0RYmHlnrnPVVwhnoK7/2YFFCCNE2EugNKIMRdfYUphx23EJX6RoaBmzff41ut3u4OiG6Jw+O2vq89v5u5fSNRtRZGfT8ZBHTdrzHiMKd9KgucizYvAZSx3i0NiG6I4PBgM1mwyRng7mVzWbDYGjfPneb/gts2LCBhQsXomkakyZN4tJLL2223erVq3nmmWd44oknSE5Oblch3YWKssKpZzB5w2qX+dqKrzBKoAvRRFBQEFVVVVRXV3foyurAwECfPz+8sdb6rOs6BoOBoKCgdm231UDXNI0FCxbw4IMPYrVauf/++0lLSyMxMdGlXWVlJZ9//jkDBgxoVwHdkWH8b9AaBTpb1qEX5KGsPTxTlBDdlFKK4ODgDq/vjWf2nKzO6nOr+/O7du0iPj6euLg4TCYT6enprFmzpkm7xYsXc8kll2A2m91eZJcbcho0Cu5aDOjffuahgoQQonWt7qEXFhZitZ54CITVaiU7O9ulzZ49e8jPz2fkyJF8/PHHLW4rMzOTzMxMAGbPnk1MTEzHijaZOrxuW5Wfexllb71MYUAEn/VK56ueo3l87f9I+0MIhuCOn3PbUV3R5+5G+uwfpM9u3O7JbkDTNN544w2mT5/eatuMjAwyMjKc0x39ytEVX9H008+Cxa/xevKFfB+XCsDHMafT9+N3MEy6sFM/uznytdQ/SJ/9w8n0+aQuLLJYLBQUFDinCwoKsFgszumqqioOHjzI3//+d2677Tays7N58skn2b17d4eK7S5UeCRq7EQuzFnpnLcy7jQKln+DrskpjEKI7qfVQE9OTiY3N5e8vDxsNhtZWVmkpaU5l4eEhLBgwQLmz5/P/PnzGTBgAPfcc4/XnuXSkMq4mIGlBxlUvA8Am8HEF4HJsLHpMQQhhPC0VgPdaDQydepUZs2axcyZMxk7dixJSUksXryYtWvXtra6V1M9e8OwkVx8cIVz3pe9xlKZ+akHqxJCiOa1aQx95MiRjBw50mXe1Vdf3WzbRx999KSL6k4Mky/hjLmP0qOykLxgC6XmUL4rC+G8fdmovt5/iqYQwnfIpf+tGZyKsVdvLjj0vXPWB70nUPPZex4sSgghmpJAb4VSCpVxMZMP/0hEjeOG+MeCLCzPrUXP2evh6oQQ4gQJ9DZQoycQFBXJxTknxtLf6zMJ26fverAqIYRwJYHeBspsRp17Becd+oGw2goAjgZbWXWwDP3wAQ9XJ4QQDhLobaTOyiA4LJQLc1YSV1nAbdvfYeyxTehLZS9dCNE9SKC3kTIHoM79LZccXM7zPz3FpCNrMeka+pqV6EdyPF2eEEJIoLeHOvs3BIaFYdS1EzN1Df2j/+e5ooQQoo4EejuowEDUlMuazNfXfo++N7uZNYQQoutIoLeTmnAeREY7pzdHJfP48KmUvf+mPIpLCOFREujtpAKDUBddC8CLAy/nkdRb+dl6Ch9Vx8Jm374VghCie5NA7wB11mSI78UpxScuLPo4aRz5H70nd2IUQniMBHoHKKMRw2U3MO7oz/QpOwxAtTGARcFD0Vd94+HqhBD+SgK9o04bgzF5EDfsPvFYumXxaez8MhO9vMyDhQkh/JUEegcppTBccROnHd9JWv42AHRl4JVek7EtedPD1Qkh/JEE+klQKYNRZ4xj6q6PMWu1AOyOSOLrnYXo+3d5uDohhL+RQD9J6sqbiNcr+e3+b53z3ur3G4oWLUTXtF9ZUwgh3EsC/SSpKCvq4mu59OB3xFU6nr1aZg7lDdMp6N9/5dnihBB+RQLdDdQ5FxIY35M/7vrYOW917DAKl7yLXnDMg5UJIfyJBLobKJMJw3V/Jq3gF8bmbWJkwXbmrXmG6NJjaP97Xq4gFUJ0iTY9U1S0Tg0ahkqfxJ2rFxOg1aLqF2z9GX1VpuNiJCGE6ESyh+5G6qo/EhgRfiLM6+jvLEAvlKEXIUTnkkB3IxUahuH3tzWZvzUgjqKFL8ptAYQQnUoC3c3UqaNQYycCYFcG3u47mYdTb+V50zC0pe95uDohhC+TQO8E6upbINLC9oi+vNN3MroysDZmCEvX7UffudXT5QkhfJQEeidQoWEYps5gaMleLshZ6Zz/3+Tz2fXW/9DLSjxYnRDCV0mgdxI1JBV13hXcsPsz+pUeAsBmMPF00oWULHhWxtOFEG4ngd6J1MXXYU4eyF+3vUWQvRqAI8ExPGMYju39NzxcnRDC17Qp0Dds2MBdd93FHXfcwZIlS5os//TTT5k5cyZ33303jz32GMeOySl6UHff9Jvvpqeq5Pbt7zjnb7AM4s1d1Wirv/NccUIIn9NqoGuaxoIFC3jggQeYO3cuq1atIicnx6VN3759mT17Nk8//TRjxozhzTfl9rH1lDUWwx//Qnr+Fq7Yl+mcv6T3BL77fKU8XFoI4TatBvquXbuIj48nLi4Ok8lEeno6a9ascWkzbNgwAgMDARgwYACFhYWdU62XUqeOQl32e67Z97Xz3ukA81MuY9PC19Hzcj1YnRDCV7R66X9hYSFWq9U5bbVayc5uea9y2bJlpKamNrssMzOTzEzHXurs2bOJiYlpZ7kOJpOpw+t6iv67WynOO8yMHxZx38jbyQmNw6zZ0CsrUM89RvQTL2OMsrS4vjf2+WRJn/2D9NmN23XnxlasWMGePXt49NFHm12ekZFBRkaGczo/P79DnxMTE9PhdT1Jv+ZWQg7s4cHNr/HU0N8zffu79CvPxQ7kP3IXhv+bhQoKaXZdb+3zyZA++wfpc/v07NmzxWWtDrlYLBYKCgqc0wUFBVgsTfckN23axIcffsg999yD2WzuUKG+TgUGYrj9IXqEmHhy3bP0K28w1HJgN9pz/0CvrvJcgUIIr9ZqoCcnJ5Obm0teXh42m42srCzS0tJc2uzdu5dXX32Ve+65h8jIyE4r1heoaCuGGY+iQsObLMvJyaNw/lMS6kKIDml1yMVoNDJ16lRmzZqFpmlMnDiRpKQkFi9eTHJyMmlpabz55ptUVVXxzDPPAI6vE/fee2+nF++tVEIShjseQnvmQaipAWB/aByPjLiVqJpS/j7/SaJvuxdVd6BZCCHaQukefPrC4cOHO7Ser4y56ZvWoL3wBJUYuO2MeygKjACgd1kuD5esIGb6X1GhYYDv9Lk9pM/+QfrcPic1hi46jzp1FIY/30Mwdm7c/SkG3fFQ6QNhCTwQPZmcebPRiwpa2YoQQjhIoHuYSh2D4dZ7GVewmTt/eRuD7rjHS16whb/1vITseXPQc3Na2YoQQkigdwvqtDEY/nQP4wq2cN/m/xJgd4yrFweE81C/q1n94qtU/7zaw1UKIbo7CfRuQo0ci+HOh0kr388jG18ltLYCgCpTIP8aeDWvvf4J9i+XyAOnhRAtkkDvRtSQVAz/N4vBehGzfn6RHpWOWyjoysCb/c/juc2laAueQa+q8HClQojuSAK9m1F9UjDc+y96hyr+tf45hhTtcS5LrMiDH5ej/eOv6Af2/MpWhBD+SAK9G1JxPTH8bQ6RKSk8svFVMg7/yMiC7VxycIWjwdFDaE/8H9o3n6BrmmeLFUJ0G269l4twHxUajuGuRwl49zWmffM+tQYTBhqMn9tqyV3yPkEbfsZyw59QsfGeK1YI0S3IHno3poxGDNfcgmHqTAIDXO+PY1MG5gy5nhmW81n13EtomR+j2+WxdkL4Mwl0L2AYOxHLUwsgsa9z3ge9J7I7IomSgDCeGnQNs7baODL7YfRdv3iuUCGER0mgewlTYl8M9z+FmngBAANLDmCpLnIuXxszhDuTrua9RZ9TveDf6AXyGEAh/I0EuhdRAYEYrrsVw18eJ9VQzL/XPMPkwz86l9cYA3iz//ncro9ixbz52N99Hb28zIMVCyG6kgS6F1KDR2B49FnCzprItOwPeGL98/QpO3Gjs7xgC8+ccg335vdk/RP/xL7kLfTyUs8VLIToEhLoXkoFhWC4fhqG+59iUJSZp9c9yx+zlxBeW+5skx3Rm6cGXEX5Fx+h3Xcz2vv/RT8uN/sSwldJoHs51W8ghr89jenaW7igeDMv/PgvLj3wLWatFoALc74n1F4FVZXoX7yPdv/NaAvmou/f5eHKhRDuJueh+wBlMKImXoA+ejxhX7zPDZmfcP6hLD5KGs9FOStdG9vt/Cc/HMt/P2KSKZ/IsyeiRp2NCgr2TPFCCLeRQPchKiQM9ds/oE84n9jP3uWP338GdptLm9xgK5/1SkdXBhZpNtLXbGLS548yLDke49gJcMqpKIPRMx0QQpwUCXQfpCyxqN9NRz//SvTP30f//muwOYZgvug5Fl05RtpsBhMr4kayIm4ksVWFjP9wFeMq3iBpyEDU6WfCwKES7kJ4EXkEnZc4mT7rJcfRv/0c/bvPqK6oICt2BF/0Gkt2RO9m2/cpy2XssU2cWbabxAH94NRRqKGnoUJCT6YL7Sb/nf2D9Ll9fu0RdLKH7gdURDTqkuvQz7ucoB+XM3HFl0xc/zy7w3qxLCGNlT1SKTOfCOv9YQnsD0ug6sB33LD6M1j9LbrBAP0Goganogaf6nhvDvBgr4QQjUmg+xEVEIg6ewqcPQV9/25SVn5J8ppMbtz1Keusg1kRl8p6y2BqjI77xowq2HZiZU2D3dv5jxqAdeOXDC19hf49wjClDEElD4J+g1DhER7qmRACJND9luqTjOozHf3qmwnctJYxq79lzJa3qdQNrLUOYYNlIAOL97usU2oK4bNeZzrH4INtVQzK2c+grT8wqORtBgRUE5bUG/oko3onQ2JfVESUB3onhH+SQPdzyhwAp6djPD0dvbKCkE1rOHvtKs7etgRwPbyyKTrFGeYAlaYgNlgGscEyyDmvZ8UxkrflMHD1l1xwaBVERkOvPqiEJOiZhIpPgrieEBGFUqqLeimEf5BAF04qOAQ1ejyMHo9eUw07NqNv/Al9y3ooyGNI8R6mb3+XrVH92RKVTEFQVJNtHA6J5XBILPvCEhyBXnwcio+jb9vA3tAEdkb2oWfFMXrZS4mODkfFxqFiEyAmDmWNBWsPsMSiAoO6/hcghJeTQBfNUgGBMDwNNTzN8WDqY0ewbNvA5O2byMj+En17EUeDLOyI6MOOyD7siOjDgdB47HWnOfYty22yzbUxg1nU71zndJC9mrjKAuIPFhC/cz+x1RvoUXWc2Krj9FA1KGs09vBIVGQ0RFogMtoxhFP/Co+E4FCUQS54FgIk0EUbKKWgRwKqRwJMOM8R8EcPkZC9jYQ9Oxi/Zz3s+ogaZWRfaAK7wxOJr2x6z5ickDiX6SpjIPvDerI/rOlpWBcfXM6Nu5cCJwZ+vu8xguMBEUTWlBJVU0ZkbRkRtgoiAgwYw8IhNBxCw1AhYRAaBiFhEBIKwSGo4FAICobgEAgKcbwPCoKAIPmDIHyGBLpoN6UUxCei4hPh7CkA6JUVBB3cw6D9uxl0YDf6wWooMUKDpyidejwbg65xOCSWQyGxVJhavt2Atbq4ybyvEkazJTql2fahtRWE2yoIr60grKyC3277lqHFe53LdWC9ZRBK1wmxVxFkrybYVk2wvZogo8IcYIbAIMcrIPDEyxyACggAc/3L7PhpMjvemxq8zGaU0eSYbzQ5XiYTGI1103U/DYa6ZUYw1L2Mjp+63Y6u63J8QXRImwJ9w4YNLFy4EE3TmDRpEpdeeqnL8traWp5//nn27NlDeHg4M2bMoEePHp1Rr+imVHAIDByGGjjMOU+31cKRHPTDByE3h0m5B5h09EfYcxi9poZScwhHgqwcCbaSF2QhLyiaY3WvuMrCJp9xPCC8xc8vN4dQbg7hSN3fiN8cXt2kzQuDrqAwMLLZ9U2ajUB7DYFaDUH2Wu7Z8ga9K446+gFoKOYPuoIAzUaAVo5Zs5146SfemzQ7Zt3GqPxtmPUTf8xqlZEDofGYdDtGXcOo2zFq9e/rpnWNirppMzoYDY7wV0bHz4YvZQCDqvtpAKVOzFeq7tWgjXNeg2WKpsug0XQz8+vfo+q2Uf/Hp24aUI3Xq9dM2+KgILTqGufsBo2bWa+tbdrwB9GlTRv/gLbUrJ1/gEuCgtCqqlAXXuMYUnSTVgNd0zQWLFjAgw8+iNVq5f777yctLY3ExERnm2XLlhEaGspzzz3HqlWreOutt5g5c6bbihTeSZnMkNgPldjPZb6uaVBUQNSxo0TmH2FQ3hEozEMv+AVyjkFRgeO890am5P7E0aBoigLCKQoIp9gcSklAGGWmYJezbwDCbRVN1q8wtnyg1WYwYTOYKCfEUXujM3xqDSa+TRjV5r6/8f0jmG2Vzun8oCj+L+2uNq0baK9h0coHwXbid/BLZF/+MXwqBl3DgO74qesoHD8Nuo4BDaXrxFUV8vCmBS7bXG8ZyKJ+56J0DYWjPVC3jbqXDqDTr+wwN+3+1GX9H2KG8U3CGShdR9WtW79e/a/KMQ1DivZw4aFVLut/FzeSddZT6j7D9ferShq813VGFu7g7LwNLut/0XPMiSubG1zc7hrJjvln5m0k9Xi2y/ofJo0nNzjmV2L7xDYzcn9iQGmOy9L/13cKJQFhLvNUCxfZX5Szkp6VrleBLki5CJtyHF9KLj1ExpE1jm2cc5HjTDA3aTXQd+3aRXx8PHFxjvHP9PR01qxZ4xLoa9eu5corrwRgzJgxvPbaa/K1UbRIGQxgiXWczTJoWJPlumaHkmKidBtF+/agFxdCUSEXlxShlxyBku2QVwylJVBdiR1FuTmYUlMoJeYQyk3BJJUfdd0mcFrhDipMQVSYgqg0BlJlDKDSGEilKRBNud6zJtBe4zJdY2jf6KRZc70pml21/Z44Br3pHzObMlJpatuZP839qys1hbI7PLGZJW3bwtFgK+utp7RpbbNmg0aBviesF6t6pLZp/aiasiaBvjWqf5vX711+tEmgr7MOZltU/zatP7xod5NAX9UjldyQmDatf1behiaBnpkwmmqj48rqM/M2OAPd3Vr9v7SwsBCr1eqctlqtZGdnt9jGaDQSEhJCaWkpERGuVw5mZmaSmZkJwOzZs4mJadsvqEnRJlOH1/VWftfnHnGYTCbMA4b8ajO9phqtpBitrAS9rASttBitrBS9vMwxXVGOXl6KVlHOA5Vb0cvKHcuqKtGqKqG6Ch1HYFbVhXy10YylpsTlcwK0Wm7b/g41BjM1BjO1BiO1BjO1BpPLy6aM2JQRU4PhFgCjbqdf6SHsyoDd4GhjVwY0ZXB5b1cGAuruZd+QXbX9wG1zfxD0du1cNd3zbM8Nnxp/uxEti46OwuTGf9ddelA0IyODjIwM53RHb04jN/PxD23vswHCohyv+LZtWwFG6r4N1FRjrKoisKaKyJpqqK6GmmqoqYHaavSaGoJqa8iorYH6l60WamvrflaArdZxzMBuB5sNThnuuHWxzfFKsNuYc+h9xzxNc/y02x3vNfuJ+fWvRoYX7eZ/Kx9GUwpNGdCUQsfxU8MxT1eOQY/mAv30gu38a92zgHKug3JEr44BXVE/8EJog6GiemOPbXZ866lfF+o+78RvVK+bF1NV1GT98UfXM6D0oGM9qBucweXLQP28xt+uAM49tJrTCnecWK+RhvMHlB5osvySg8sZd/TnFtZ1ldxo7xzgmn1fUdHgG1JLf7J0FAnNnOF1065PnH+UExrsvR8/XoQKat+/65O6OZfFYqGg4ESBBQUFWCyWZttYrVbsdjsVFRWEh7d8AEuI7kIZjHWnMYa03KYL64mJieFYXp5jnFizOwPeoGmYNQ30E/McbRr91Ov/IOig1U3rOpG6TmR9G+rbc+I91LXlRJv65bpOz7qXc75LmxPtqF9Uv6BuXoqu4zw/qX7dumbh4WGUlpY1mX+CzjBgWIPt0fi9S/MzXNYFOKNJm4bvG29nZJNNjmtaVPPbAmBEkya/aTQdFnYxZWXljusp3KjVQE9OTiY3N5e8vDwsFgtZWVnceeedLm1OP/10vvvuOwYOHMjq1asZOnSojJ8L0UHO8+KN3nkv+vb+yw+OiaHcz759hsTEUNEJfW410I1GI1OnTmXWrFlomsbEiRNJSkpi8eLFJCcnk5aWxjnnnMPzzz/PHXfcQVhYGDNmzHB7oUIIIX6dPODCS0if/YP02T901gMu5JpnIYTwERLoQgjhIyTQhRDCR0igCyGEj/DoQVEhhBDu45V76Pfdd5+nS+hy0mf/IH32D53VZ68MdCGEEE1JoAshhI/wykBveIMvfyF99g/SZ//QWX2Wg6JCCOEjvHIPXQghRFMS6EII4SO69AEX7tDaA6t9TX5+PvPnz6eoqAilFBkZGZx//vmeLqvTaZrGfffdh8Vi8YvT2srLy3nppZc4ePAgSimmTZvGwIEDPV1Wp/r0009ZtmwZSimSkpKYPn06AQEBni7LrV544QXWr19PZGQkc+bMAaCsrIy5c+dy7NgxYmNjmTlzJmFhYa1sqW28ag+9/oHVDzzwAHPnzmXVqlXk5DR9uogvMRqN/P73v2fu3LnMmjWLL7/80uf7DPDZZ5/Rq1cvT5fRZRYuXEhqairz5s3jqaee8vm+FxYW8vnnnzN79mzmzJmDpmlkZWV5uiy3mzBhAg888IDLvCVLljB8+HCeffZZhg8fzpIlS9z2eV4V6A0fWG0ymZwPrPZl0dHR9O/veLhtcHAwvXr1orCw0MNVda6CggLWr1/PpEmTPF1Kl6ioqOCXX37hnHPOARzPjw0NDfVwVZ1P0zRqamqw2+3U1NQQHR3t6ZLcbsiQIU32vtesWcP48eMBGD9+vFszzKuGXNrywGpflpeXx969e0lJSWm9sRd7/fXX+d3vfkdlZdNnW/qivLw8IiIieOGFF9i/fz/9+/fnxhtvJCgoqPWVvZTFYuGiiy5i2rRpBAQEMGLECEaMaProNl9UXFzs/OMVFRVFcXGx27btVXvo/qyqqoo5c+Zw4403EhLS8vMvvd26deuIjIx0fivxB3a7nb179zJlyhSefPJJAgMD3fo1vDsqKytjzZo1zJ8/n5dffpmqqipWrFjh6bK6nFLKrY/r9KpAb8sDq32RzWZjzpw5nH322YwePdrT5XSqHTt2sHbtWm677TbmzZvHli1bePbZZz1dVqeyWq1YrVYGDBgAwJgxY9i7d6+Hq+pcmzdvpkePHkRERGAymRg9ejQ7d+70dFldIjIykuPHjwNw/PhxIiIi3LZtrwr0hg+sttlsZGVlkZaW5umyOpWu67z00kv06tWLCy+80NPldLrrrruOl156ifnz5zNjxgyGDRvW5KHkviYqKgqr1ep8JOPmzZtJTEz0cFWdKyYmhuzsbKqrq9F1nc2bN/v8geB6aWlpLF++HIDly5czatQot23b664UXb9+Pf/973+dD6z+7W9/6+mSOtX27dt5+OGH6d27t/Or2bXXXsvIkSM9XFnn27p1K5988olfnLa4b98+XnrpJWw2Gz169GD69OluO5Wtu3rnnXfIysrCaDTSt29f/vznP2M2mz1dllvNmzePbdu2UVpaSmRkJFdddRWjRo1i7ty55Ofnu/20Ra8LdCGEEM3zqiEXIYQQLZNAF0IIHyGBLoQQPkICXQghfIQEuhBC+AgJdCGE8BES6EII4SP+P5QoH3qXiVVbAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Runge-Kutta 4th\n",
    "dt,tmin,tmax =0.1,0.0,10.0\n",
    "step = int((tmax-tmin)/dt)\n",
    "# create array t from tmin to tmax with equal interval dt\n",
    "t = np.linspace(tmin,tmax,step)\n",
    "y = np.zeros(step)\n",
    "y1 = np.zeros(step)\n",
    "y2 = np.zeros(step)\n",
    "y3 = np.zeros(step)\n",
    "ya = np.exp(-t)\n",
    "\n",
    "plt.plot(t, ya, lw=5, label=\"Exact\")\n",
    "y[0] = 1.0 # initial condition\n",
    "\n",
    "for i in range(step-1):\n",
    "    y1[i] = y[i] - 0.5*dt*y[i] # Runge-Kutta Eq.(A18)\n",
    "    y2[i] = y[i] - 0.5*dt*y1[i] # Runge-Kutta Eq.(A19)\n",
    "    y3[i] = y[i] - dt*y2[i] # Runge-Kutta Eq.(A20)\n",
    "    # Runge-Kutta Eq.(A21)\n",
    "    y[i+1] = y[i] - dt*(y[i] + 2.0*y1[i] + 2.0*y2[i] + y3[i])/6.0\n",
    "    \n",
    "plt.plot(t, y, lw=3, ls=\"--\", label=\"Numerical\")\n",
    "plt.plot(t, y/ya, lw=3, label=\"Ratio\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "<div id=\"extra\" />\n",
    "\n",
    "## Extra\n",
    "\n",
    "1. Compare the Euler method with the Taylor expansion\n",
    "\n",
    "$$\n",
    "y_{i+1}=y_i+\\Delta t \\frac{dy}{dt}+\\frac{\\Delta t^2}{2} \\frac{d^2y}{dt^2}+\\frac{\\Delta t^3}{3} \\frac{d^3y}{dt^3}+\\cdots\n",
    "$$\n",
    "\n",
    "to find the order of the method in terms of $\\Delta t$.\n",
    "\n",
    "2. Compare the Leapfrog method with the Taylor expansion to find the order of the method in terms of $\\Delta t$.\n",
    "1. Use the code example for the Euler method and examine the amplitude of the numerical errors (i.e. deviations of y/ya from 1), by changing the time increment to dt=0.01, 0.05, 0.1, 0,5, and 0.9.\n",
    "1. Perform the same error analysis for the Runge-Kutta (2nd order) method by modifing the code example for the Euler method.\n",
    "1. Perform the same error analysis for the Runge-Kutta (4th order) method.\n",
    "\n",
    "\n",
    "1. Euler method:\n",
    "\n",
    "$$\n",
    "y_{i+1}=y_i+\\Delta t f_i=y_i+\\Delta t \\left.\\frac{dy}{dt}\\right|_{i}\\tag{AH1}\n",
    "$$\n",
    "\n",
    "Taylor expansion ($+\\Delta t$):\n",
    "\n",
    "$$\n",
    "y_{i+1}=y_i+\\Delta t \\left.\\frac{dy}{dt}\\right|_{i}+\\frac{\\Delta t^2}{2} \\left.\\frac{d^2y}{dt^2}\\right|_{i}+\\left.\\frac{\\Delta t^3}{3} \\frac{d^3y}{dt^3}\\right|_{i}+\\cdots\n",
    "\\tag{AH2}\n",
    "$$\n",
    "$\\rightarrow$ 1st order\n",
    "\n",
    "\n",
    "\n",
    "2. Leapfrog method:\n",
    "\n",
    "$$\n",
    "y_{i+1}=y_{i-1}+2\\Delta t f_i=y_i+2\\Delta t \\left.\\frac{dy}{dt}\\right|_{i}\\tag{AH3}\n",
    "$$\n",
    "\n",
    "Taylor expansion ($-\\Delta t$):\n",
    "\n",
    "$$\n",
    "y_{i-1}=y_i-\\Delta t \\left.\\frac{dy}{dt}\\right|_{i}+\\frac{\\Delta t^2}{2} \\left.\\frac{d^2y}{dt^2}\\right|_{i}-\\left.\\frac{\\Delta t^3}{3} \\frac{d^3y}{dt^3}\\right|_{i}+\\cdots\\tag{AH4}\n",
    "$$\n",
    "Eq.(AH2)-Eq.(AH4) gives\n",
    "$$\n",
    "y_{i+1}=y_{i-1}+2\\Delta t \\left.\\frac{dy}{dt}\\right|_{i}\n",
    "- 0 \n",
    "+\\left.\\frac{2\\Delta t^3}{3} \\frac{d^3y}{dt^3}\\right|_{i}+\\cdots\\tag{AH5}\n",
    "$$\n",
    "$\\rightarrow$ 2nd order"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "celltoolbar": "Raw Cell Format",
  "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.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}