{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cruise Control Example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, import necessary libraries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import *\n",
    "init_printing()\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import linalg\n",
    "import numpy as np\n",
    "from numpy import linalg as LA\n",
    "from matplotlib.backends.backend_pdf import PdfPages\n",
    "from collections import deque"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function allows us to ensure the correctness of simbolic manipulations. This is helpful to obtain the analytical solution of the cruise control example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def equals(a,b):\n",
    "    assert np.isclose(float(simplify(a - b)), 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following defines the equations for the example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAM4AAAAqBAMAAAAJ5HljAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAVO8Qq4lmdpnNRCIyu925KXARAAAACXBIWXMAAA7EAAAOxAGVKw4bAAADBElEQVRYCe2Vz2sTQRTHv9vpJummTUNBL4rtRa2KuIgiHiQ5WIJKcQ9i9SI5CV5swENLKZiDtYhQc2xUJPYP0BQ8iIKtZ5FqD/bgwQitZ0PQQz3om53Z7dt1oaGwOfXBy877zmfmza/MAB20AbsTycTI5Y7kAcZ38+xoQ3fXbUfLtnvedrZsuNuZ/+m9E/NFPcK0V9h+xG2iRh2QHrRzPOxdLfNQlcVqQxUC6P+cUj4O4xKVpAfsSSDqzwZCFfSU1TeIRoCu1MRZ+krnlijzCC/ZvhklXaWTh1Bd6VM6tlowqaF0bpk6j3CDRX4POnkI1aRP6dhsQAzBddYX3vEAmGeh34NOHkI16VM6lov8mMrSmU2ycnpmuokzs8gUXdHvYR57H5LCUDFZwg/VUlNaeT22kgWuUZV0Zl9YeaHe14KzgEzeFb08opl885sUhprWENZVS00pxTqGbzZwi6qkM/u+VRaHkKqks3MwHFf08hitO6JICkMXUw3scSFoSikni3hA8gHX/3r2k8KjtGVV15ZTLWQcgcOwaOpmtfrsSLWaJyQ1TD9kDLW7SnhO0hallEc2Nkhe0U4fz6ixZ12fkKuBzmVSKd58MoU1V2CoPP+lIOUqxyE2SY7IwxYj52CQhtiAGewhV1uqSYWh8vyrxfXWzVXELxhykV5ol620sc3NZXHBqHfnsV/VefMZRy4rO2UozsNYDlJSEU10V/ZFngN2/PpKiXUTZkWERjqK/mUpMRTXcVql8efjKmu4nSd0huqkM5NT1CYK06+mgJGbOvbm8xnJWZs0hqL7ajFEucqpiQ+FGvCU6qQzi75MXMDL49HRaJgiWjSUew3lN1HmUaAs5CSYRaNhihpE3aP4ynraptguGvUuoK3HSw2gXVSO578xtfkYy0xtook6bQe5deXtRKF37KIaZHy/Zt8aclMYiC+D6vl9poIlB+Nx57Hpthy0MRp3HnlbzoHu/7iNnqINWH9E3HkOylciVVmMOU96U74SPY4Tcx6rgWQZxn36K8Vp/wD5v7zjjMuUXwAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$\\frac{1}{m} \\left(- d v + k \\left(- v + vd\\right)\\right)$$"
      ],
      "text/plain": [
       "-d⋅v + k⋅(-v + vd)\n",
       "──────────────────\n",
       "        m         "
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t, v, m, k, vd, d, v0 = symbols('t, v m k vd d v0')\n",
    "dv_dt = (1/m)*(k*(vd-v) - d*v)\n",
    "dv_dt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rewrite the previous equation into the form $av + b$, so that we can define the analytical solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJYAAAArBAMAAABoe2x6AAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAELvv3c2ZVESJZiJ2qzKF9OxMAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAC0klEQVRIDa1WTWgTQRh9kybZTdLUteJBqaTgTQsWEfHvEPAgFQK5RARRFloPFcFAvYhiChU8dg/Wg1TIQT14MQexBxGjh1CwlCAIihGKFw8emp5SEdRvdmdjJjOJGfA7zM73vve93Z3deQwwegQd8fVoR2I+PSe1tKTMNJG67S3T9k6+3J3IdtZM53J3zDXt7+THXOtJhdWqOAu8WN4sdtb8+dhVRMsKqgUyxXeL5YQ9iRzskzjkKCRvAdG6gmqB0pizFyvxLG5gVxmXFU6qOAPLU2AtMO0S7AxVcRfXHEwpHIZTsNU3V3gcyL2t0rjooIrDYNsckoM+9LCM9MrYdor33wM89gNWU+XRh06oqA5JNlmLqF9gbbAW0vmbCildx20F1ALpcfwirSXsASbwqK6uciLPVFCrFSmi9gZIPygD+9a+N+ZV1qenbSxFLH1YFT3eGx3lpZH9Lr+0I3Caz+18wMkFn7ejKNN9p9ktY//Mkq5PoR+Ixy1/pMF3mkQ1TAe7RoNFoR+IR6gVOA2bDNBBx+cBUbhnqCWc5vygKgGvRpfU7Lpwz1DLdxrgvpnWAaIvVCLCe0Mt32mAx2ZaZwB2AvF80BVq+U4DXNdqsYIfG0rxGBDfQtSjwlKhcLpQuMgp0y4fsYnf3eHjf4d2uUkYaQ2NIzMflMPnCpyGtIyC3jHjoSR6hJZwGtwxkgKtfaaID2LvCS3hNKZrT589Uk3mhJ0JLeE0mDV7LnoN1lh/tRp0CS3hNLhkpiX2kGgSWiJjWTMtsbdF006p2XRv46DULiWmngPfCyWJdtLnNm2ONOnt0cmKRPwvif3w5VpjZPmjTqxfTcdPRCaQWcWcaU3H/xbN47WHZ6Y1Hd+hnV9y8N60puPznT8DOs/ool9Nx6dT1xTsn8y0puMf54eleH7FtKbh06mJDksxj3tvd/SrdXN5bmcx7MK6ovt1+9W6tP4ASAXBENmDg2oAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\frac{k vd}{m} - \\frac{v}{m} \\left(d + k\\right)$$"
      ],
      "text/plain": [
       "k⋅vd   v⋅(d + k)\n",
       "──── - ─────────\n",
       " m         m    "
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = -(1/m)*(k+d)\n",
    "b = (1/m)*(k*vd)\n",
    "dv_dt_norm = (a*v + b)\n",
    "equals(dv_dt, dv_dt_norm) # Proves correct rewrite.\n",
    "dv_dt_norm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define analytical solution, and prove that it is indeed a solution by differentiating it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATEAAAAyBAMAAADGqPFBAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAELvv3c2ZVESJZiJ2qzKF9OxMAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAGXElEQVRYCa1ZW2gcVRj+Z2+zl9mLRoVKyy5G7EvBULyhFRbFYpXIgm4FQRlQUdvYLBTB1mqWVhAE6T6YolLoKipYCllQvDwUt0qD2DRdlD6IUZfqkz4keTG9Rs9/LjPnzJzZ3dnkwJ75/st3zj/nzJn58wcgRHsvhK/qukkV11tKNYYeMVkZjmqMliVibtSWJAmelTBClXbhLo9ZET9RpMGFuC375muy5GBzzIEcqLQVr1mWI0OutxrL4YI8poNH6g7kQKEll71mWTbukaXBsRrLu3ribT61QkuVfXZZsVcWBsdqLDu1RNO/JgpN3VrfELHhtnMn7N8NGycg1oTMS/P6ByZm+2ZDmjHbhkcAvjy6WPPZZYV5RZYGxcaK9fVlaByCWAcOtSL+xcGBZlre4SgtlRyDcUjeB7frn06HtN1BIYC5fMpoZmrPgdkw7oVERUt9wKeltK8SZdgHI0143mdXFTPayE8yJ+9d5zg3ga8iA7ZBspZYBv0TkbuoTkQkSitE2/AWTBZgh8+uKoptVWbSCXpJYT/ZpRjg+zaYdYZjC1sIIMfegmgJilzLbKJPjAnkXBmNHNA2bAXDH7njSEFctxWp7U00/opdETtssTrAbnIlrVj/tg5Ajn0Kig2YojpvF+14NZxGDmjDuATmks+uKiztESijU66B/Svkl0VQLADwD+0MFGsNyHbgDQLgvNlCu6ctNj0KciYo7Tcwu8YKZCuv+RwUReayIjLBqOA128X+CPlFELxJfnEWw0+Q7zYgVTEaEGmnx+muo4vczvnDZbRp2ACwBT7uNGR3Db5fozPtJtHGyO+muSdJTyM7kjsLFuoBNoM1UQD45Tg5BwvzJ3+gSk836pEdWvaDJln8uX8W6n4PRXOrIjEhM4E3/AWZ97EcvqxoZE+9WiyYHY23XrVVrw6hPRboe5Ccczs3Ruw0smtdchxLgd4eg/GvRxFenOk6nIO3kFaA/2hbgncAojWrA+erj1erzdy2zeRQVBznAHCA63NXAxy0arOtUS/qlNTvAjmPrciHBbZmVunPDCT77qaIzNR/TDUBoEobWdEO8IYbyJq1Fv8mZtzNSPsvC9KNIGehF5GlV4VmkKs2snzgMsQKkFw4sIuMjJHFW9l5SNT7zSMi078ng9giMpa+cK98iYINc/z9DpB+4fiefV0Aq8Zd2AlAIS40gVcRWUqfgATwRGQ0fbGOsp2JVtA7+Sh8I1jmabjbRGNJaHDNaNsvQOBVREa+9CEaj4ylL5/B65QaLeNlpGbsoSJ2vxtlum1vC02GA8MGfnLZ+ZV7dJmuVh+sVp9GmF3CHm6sYnsG8rKrhNErVa0+cWe12hHpy8OQLyA3Nob9zyd+pBJiKJt2zCRX3/fQcn2oo6YTa8Yj03joVGI3MX0xVsnLCp1YZNckf6Nj1bNnJEUYKCIbajdp+pK7gkk8abEy9pewW48mIktpc5igGcSaYfpCvoexNjqyE7CDHMUgWii9iMwa6n2G6QtZM5L2ksYyvNNgOG+NUJF4nUVk6X4pq0IUa4bpC3nO8nW05m3srV0koxm0eQsZ19/hMkVkpj8ybyFDojmR4UDibJJ0OXTL07Pj0na68DoOc5qM1PspdWmGfOhn+fss+IvuTudF5I8MbGJ9wDslGv1ZkChk9KQRpvW+jQPAlOf+qbJPR/7IwCamEFMyLe8fUiQUUmWm6kljLrQ/1pKEASHfBjGFmFJhTyoSCqKQ0ZMmsf6Q8EDQLWSIKeK2+anv/mbYjjtDuoWMnjTHH8jfpCGbW8gQUxRrZw43vaMsqk+JVMjoSXOHyVx18UBIKmSIKaY2Fm72cdl70lGPuIWMnjSHAOEyPMIThYxpN6941nbHc1BiyYEIJnkhox/NJUXKLh4IRUtOIUPc/Ph3bT81qb5qpUJGT5o7UF4zqGvVIKmQwacwLmbUKBhLqX/JhYzeNGfKc54j5BiCgFTI4FOkl4yVlN9dOZxyIaM3zRlIVzxwjDogFTL4FNkSrGoiizZkulTI6E0TJDNUGoUsqZDBp4jUYPaUGNC9qhVkqZDRmyYGEO9lIYe68imCOO4HW/XoQ+POe8M+ZvIkIq+QdRL+vCsJEuxDY57D/qdCmicYmuVgW19L1u7rsgaHl9fA/WgN3P7UbLu/T4BHshJgWCf1i0OPsyn88/8/VfGjP7zNOTcAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\frac{k vd}{d + k} + \\left(- \\frac{k vd}{d + k} + v_{0}\\right) e^{- \\frac{t}{m} \\left(d + k\\right)}$$"
      ],
      "text/plain": [
       "                        -t⋅(d + k) \n",
       "                        ───────────\n",
       " k⋅vd   ⎛   k⋅vd     ⎞       m     \n",
       "───── + ⎜- ───── + v₀⎟⋅ℯ           \n",
       "d + k   ⎝  d + k     ⎠             "
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v_t = exp(a*t)*(v0 + b/a) - b/a\n",
    "equals(dv_dt.subs(v, v_t), diff(v_t, t))\n",
    "equals(v_t.subs(t, 0), v0)\n",
    "v_t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Explicit_Euler:\n",
    "    def simulate(self, f, x0, t0, h, max_steps):\n",
    "        x = x0\n",
    "        xs = [x]\n",
    "        for i in range(1, max_steps):\n",
    "            x = x + f(x, t0 + h*(i-1))*h\n",
    "            xs.append(x)\n",
    "        return xs\n",
    "\n",
    "class Implicit_Euler:\n",
    "    def __init__(self, max_its):\n",
    "        self.xp = None\n",
    "        self.max_its = max_its\n",
    "    \n",
    "    def simulate(self, f, x0, t0, h, max_steps):\n",
    "        x = x0\n",
    "        xs = [x]\n",
    "        for i in range(1, max_steps):\n",
    "            converged = False\n",
    "            it = self.max_its\n",
    "            xc = xn = x + f(x, t0 + h*(i-1))*h # Use Euler to get better initial guess.\n",
    "            while (not converged) and (it>0):\n",
    "                xn = x + f(xc, t0 + h*(i-1))*h\n",
    "                en = xn - xc\n",
    "                converged = np.isclose(LA.norm(en),0.0) if len(en) > 1 else float(sqrt(np.dot(en,en)))\n",
    "                xc = xn\n",
    "                it = it - 1\n",
    "                if (not converged) and (it == 0):\n",
    "                    print(\"Warning: implicit euler did not converge.\")\n",
    "            x = xc\n",
    "            xs.append(x)\n",
    "        return xs\n",
    "\n",
    "def convert_autonomous_ode(A, B, x0):\n",
    "    n = np.shape(A)[0]\n",
    "    x0 = np.concatenate((x0, np.matrix([[1]])), axis=0) # vertical\n",
    "    A_expanded = np.concatenate((np.concatenate((A,             B)            , axis=1),\n",
    "                                 np.concatenate((np.zeros((1,n)), np.zeros((1,1))), axis=1)), axis=0)\n",
    "    return (A_expanded, x0)\n",
    "\n",
    "def analytical_solution_matrix(A, B, x0, h, max_steps):\n",
    "    # Computes solution to dx/dt = A*x + B\n",
    "    n = np.shape(A)[0]\n",
    "    (A_expanded, x) =  convert_autonomous_ode(A, B, x0)\n",
    "    exphA = linalg.expm(h*A_expanded)\n",
    "    xs = [x]\n",
    "    for i in range(1, max_steps):\n",
    "        x = exphA * x\n",
    "        xs.append(x)\n",
    "    return [x[0:n] for x in xs]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4lHXW8PHvyWRSCCmUJASSkFCkQ8CAIGDFroCriBVd\nFbDsqrurK+76WFb3edm17Pq4lhV0xS6iLooF7KAgECAQioJIIKGFTihJJpPz/jEDhpAyCZlMkjmf\n68qVzMxdTqLc575/7YiqYowxJniFBDoAY4wxgWWJwBhjgpwlAmOMCXKWCIwxJshZIjDGmCBnicAY\nY4KcJQJjjAlylgiMMSbIWSIwxpggFxroAHzRtm1bTUtLC3QYxhjTpCxZsmSnqsbXtF2TSARpaWlk\nZWUFOgxjjGlSRGSjL9tZ05AxxgQ5SwTGGBPkLBEYY0yQ83sfgYg4gCxgs6peLCKtgbeBNCAXuEJV\n99T2uC6Xi/z8fIqKiuoz3KAVERFBcnIyTqcz0KEYYxpYQ3QW3wmsAWK8rycBX6jqZBGZ5H19b20P\nmp+fT3R0NGlpaYhI/UUbhFSVXbt2kZ+fT3p6eqDDMcY0ML82DYlIMnARMLXc26OAad6fpwGj63Ls\noqIi2rRpY0mgHogIbdq0sacrY4KUv/sI/gn8ESgr916iqm71/rwNSKxsRxGZICJZIpK1Y8eOSg9u\nSaD+2N/SmEYobxHMe8Lz3Y/8lghE5GKgQFWXVLWNeupkVlorU1VfUNVMVc2Mj69xPoQxxjQMXy7O\n9bFN3iKYNhK+/Kvnux+TgT/7CIYCI0XkQiACiBGR14DtIpKkqltFJAko8GMMxhjjkbcIcudB2nBI\nGVS3bY5cnN0l4AiD6z84bjvdtBBeGQXuEtThZONFb7G7dT8OFrs5VFLKwWI3LQqWEr/8VpaFO+j/\n3ZMsS32Ctc4eHHK5KSpxc9jl5oK9b9DLoSyNiiKz2EVG7ryq4z5BfksEqnofcB+AiJwB3K2q14rI\nY8D1wGTv95n+isEY0wwE6AK+9rzX2RbTl8KiUgqLStlf5KLn+vcYWlqMgzLcpcXMeOcN3gh3caDI\nxcFiNweLSxnnfpezWijLoqLof7iEL2a8ybPuQ8eca0zUm3yVHIdLBKcq5236kG81nsgwB5FOz9fi\n6ESeCo/HJRCmMKVVEhn19TetIBBLTEwGpovITcBG4IoAxHDCJk2aREpKCrfffjsADz30EC1btuTu\nu+8OcGTGNCH1cAGvuI2Om0lhwgD2HCxhzyEXew6WEL/8I3qWFhPivYB/8sF0ZrSAfYdd7DvsYv9h\nF1cWv8OZkeUu4P+dzrPu4mNOdXJIIuHRkWRHhNK/qJRljt7ERITSIS6CqLBQosJD0aLuTDw4H5eA\nMw5+174/09IHERXmoEVYKFHhDmYun49rUw5lIpQC6UO7sGD42ceca2rOar5ZGkIZiktCyNJDTTsR\nqOrXwNfen3cBZ1e3fW09/OEqVm/ZX5+HpGf7GB68pFeVn48dO5a77rrraCKYPn06s2fPrtcYjGnS\n6uMinzsPdZcg6kbdJWxeNoeV+1LYeaCEXQdK2HmgmMy8t7nYe5deWlrMP194kX+V7jzmMAMkjnuj\nI8mODCWjqJSP9ndmV1kJsZFO2sdGEhPpxFXcnYlFv1zA70k5mXe7DiE6wklMhJPoiFDW7evA+Nmv\nUlLmIizEyZTzMslIOPbyPDVnAS7vBbxUQjjc3snpJx3bz3la7zG8uvkTXGWlOB1OMrtdetyfJzMx\nkzBHOK4yF84QJ5mJmXX77+CDJrHoXGPUv39/CgoK2LJlCzt27KBVq1akpKQEOixjGkZ9NMVsmAfe\ni3yZu4SV387im4RWFBQWU1BYREFhMfF7I3mqzIETxaUO7lgQxdL5S48eI66Fk70RXWkf3oJlEQ4G\nFLtp1/kc7m/Xg7gWYbSOchLXIoztxW24/ft2uI5ewIdUcgFfeMwFvDDRwckdWx+zTdb2LErK3JQB\nLi0ja3vWccfx5QKekZDBlPNeImt7FpmJxyeTo9ucO6XabepLs0gE1d25+9OYMWOYMWMG27ZtY+zY\nsQGJwZh6V98X+dISsud+yJzWMWzbd5gt+4rYtq+IdvtDmeb45SL/0IpWLNW1xEY6SYgOJyEmnKjO\np/LPsvsodn1Px7anc1eXK2jTMoy2LcNpHRWG0xFCdkE242cn/nKXfnpXMhI6HRPP1JyVuOrhAu7Z\nJqzmi7wPF/CMhIwaL+6+bFMfmkUiCJSxY8cyfvx4du7cyTfffBPocIypmR+aa9Yt/JhvNrQlf88h\n8vYcJn/PIdrscfKSeC/yOHh0ZWtyQn4mMSaC9rGRZKTEkRR3Lv8pFvYfmk+PDiN4qs/VxEeHE+F0\nHD1VdkE24+e8SYmWELZ7A1PihtKrwoXR97v0E7+A1+dFvjGxRHACevXqRWFhIR06dCApKSnQ4RhT\nPR8v8rhLwHuRL1r7NWu0K7k7D3q+dh0ifFsMf1EHTvXcyU9aEstSXUN0eCjJrVvQsU0UyV3O5qVS\nKDw8n97JI3i+/7W0jQonJOSXiYuei/yrlLhLCMtbz5QemaQ4K7nIu0soowxXmcuvF/kj2zWWu/SG\nZIngBOXk5AQ6BGM8arrbL3eRx13ieZ0yiFJ3GZt2H2JdwQEKd6UzklBCvM0113zuZOln8wEIEWgf\nF0l62778Pf5eytyLOKndmfyl99WktGpBTGTo0Rnqx1zkN61nSvdMEqIb/0U+WFkiMKYpqIcmndLU\noYQ4nFAKbgnlqbUJfJY1lw07D1LiPrIKTDizox9iRORa9iaewoWpp3BbmyjS2kaR0jqS8FCH9yL/\nNiWUELZjA1OcQ4htYRf5pswSgTGNXS2bdHCXUPTTN6wo7cKarftZs3U/q7fu58dthfRyT2JwyBq+\n1x7s2p1Ml/hIzugeT5f4lnRNjKZzfBTr968ia3sYpyT2OK7TFewi3xxZIjAm0OrYpHPE3kMlbAjp\nTW8JJcTbbn/1nFCW6gIAWrVw0rN9DNcN7kjL2FAOiJNJnU5lUPsBx53Kc7c/3tOk4whjyrlT7CIf\nBCwRGBNIvtztpw33fOYuQR1h/BjRj+++3cCK/L0sz9tL7i7P8gUD5D4uiF7PvsRTOLvTEH6bFEOP\npBgSY8IREe9F/k+UuEuYmTut0ou8L3f7dpFvfiwRGBNINdztb9tXxKLdyWxL/ydh+d/x0f7OLH63\nGFhNu5gI+qXEcsXAFDKS4yAikdV7ohmR2IeMhC7Hnaq+mnTALvLNjSUCY/yppmafCnf7W+JO5tvF\nm1i4YTeLc3eTt/swAC3C4uifeg2DMuIYnxxHv5Q4EmMijh7Gc7d/a4M16ZjmxRKBMf5SQ7OPqpIb\n2Ys1J/+bknVzmbkvna9ePwjk0DoqjIFprbh+SBqD0lvjCt3Ash1LvBfndsedypp0zImwRNDI5Obm\n0qNHD7p163b0vd///veMGzeuyn1s5dMAqUMn757WGXy3fiffrtvJvHU72bz3MBBKh7iLOOWk1vy/\n9NYMTGtN5/ioY8bk3zJnwgnf7YNd5E3lgisR+LKueSPQuXNnsrOz/Xb80tJSQkOD6z99vfOxk1cd\nYai7BLeE8qesGGZ8/BmqEB0Ryqmd23DLGZ0Z3qUtHdu0qLJcaH3e7RtTmeC5GvjyD7cWalOP4LHH\nHmP69OkUFxdz6aWX8vDDD7N48WJuuukmFi1ahNvtZtCgQbz99tu0bNmyynO2bNmSAwcOADBjxgxm\nzZrFyy+/fMw269ev5/bbb2fHjh20aNGCKVOm0L17d2644QYiIiJYtmwZQ4cO5cknn6zz726otpN3\n32EXX/9YwGerQ9nr+jN9S3NYTE9I6MVdZ8czrGtb+iXHEurwVIrNLsjm85VVX8Dtbt/4m98SgYhE\nAHOBcO95ZqjqgyLyEDAeOFKR/k+q+rG/4jiqhtEZteVrPYI5c+awbt06Fi1ahKoycuRI5s6dy2mn\nncbIkSO5//77OXz4MNdeey29e/cmNzeX9evXk5Hxyz/op59+muHDh/sU14QJE3j++efp2rUrCxcu\n5LbbbuPLL78EID8/n/nz5+NwOGo4iqlNJy+OMAraDOTj7zbw2ZrtLPx5N6VlStuWYZzd53Qyeozh\n1s5tiI5wHncYX8bt292+8Td/PhEUA2ep6gERcQLfisgn3s/+oaqP+/Hcx6vwD5c03y6sVfG1HsGc\nOXOYM2cO/fv3B+DAgQOsW7eO0047jQceeICBAwcSERHB//3f/x3dp65NQwcOHGD+/PmMGTPm6HvF\nxb9UWBozZowlAV/48vSYMoiNF7/FpmWzeW93Ou+/UgispnN8FDcP78Q5PRPJSIkjZ+dysrbPImF/\nJhkRx1/AfWn2AbvbN/7lz5rFChzwvnR6v9Rf56tRyiDPP+h67CPwpR6BqnLfffcxceLE4z7btWsX\nBw4cwOVyUVRURFRUVLXnK9+GXFRUdNznZWVlxMXFVZlEajq+8arm6TFv9yE+WL6FD5dv4YdtBwmR\nYWR2bM2fBiYwokcineJ/adqrr1m6xvibX/sIRMQBLAG6AM+o6kIRuQD4rYiMA7KAP6jqHn/GcVTK\noHrtJPalHsF5553H//zP/3DNNdfQsmVLNm/ejNPpJCEhgYkTJ/LII4+wYcMG7r33Xv71r39Ve77E\nxETWrFlDt27deP/994mOjj7m85iYGNLT03nnnXcYM2YMqsqKFSvo169fvf3OzUItm332JJzCzO82\nMHP5FpZt2gvAyR1b8ZdRvbiwTxJtW4ZXehrr5DVNhV8Tgaq6gQwRiQPeF5HewHPAI3ieDh4BngBu\nrLiviEwAJgCkpqb6M8w686UewbnnnsuaNWsYMmQI4Onwfe211/j0009xOp1cffXVuN1uTj31VL78\n8ks6dep0XB/BjTfeyB133MHkyZO5+OKLiY+PJzMz82jHcXmvv/46t956K48++igul4srr7zSEkF5\nPjb7HLrqfX5a/Anv7UrnlZf3Uab76N4umnvP784l/ZJIbtWC7IJs/rvhc+vkNU2eeFpwGuBEIg8A\nh8r3DYhIGjBLVXtXt29mZqZmZWUd896aNWvo0aOHHyINXkHxN533BHz5V0+zjzjgrD/D8D8Anma8\n5fn7eGvRJj5YvoVDJW5SW7dgVEZ7RvZrT9fEX57AfGn2ObKd3e2bQBGRJapaY3ujP0cNxQMuVd0r\nIpHAOcDfRCRJVbd6N7sUWOmvGIw5TiWDBvYddjEzezNvLNzED9sKiXQ6uKRfEmMHpjIgNa7S8f3W\nyWuaE382DSUB07z9BCHAdFWdJSKvikgGnqahXOD4XtQmKCcnh+uuu+6Y98LDw1m4cGGAIgpSNbX/\newcN6IZ5/BDRjykLwvg453OKXGX07hDDo6N7MyqjfaVDPcuzTl7TnPhz1NAKoH8l719XyeZNXp8+\nffw6G9j4wIf2/8Mlbt7dnMi0rEzWFRygZfh2fjUgmasGptInOfbodjU16Vgnr2lOgmdmsWn+qhn2\nuaOwmFcX5PLq9xvZc8hFnw6x/P2yvlzUN4mo8GP/Gfja/m/NPqa5sERgmo9K2v/XbS/kxW838N6y\nzZSUljGiRwLjh3diUHrrE1rbx5jmxBKBaVqq6wMo1/6f4+zDPz6Hr36cS3hoCJefnMxNw9LpHF/1\nWk5HWPu/CTaWCEzT4cP6/rP3pfJ09iBWbdlPm6h9/G7ESVw7OJU25SZ9Wfu/MceyROBnW7Zs4Y47\n7mDGjBk+73PDDTfwzTffEBvr6bxs0aIF8+fPr3af8iuTNltV9AGoKl/+UMCTn61l1Zb9pLeNYvKv\n+jC6fwcinMeurWTt/8YcL6gSQSAm97Rv375WSeCIxx57jMsvv9wPEXnunFWVkJAQvxzfbyr0AWjH\nYcxdu4MnP1vL8ry9pLZuweNj+jE6o/3RJZ4rsvZ/Y47XxK4EdXfkTvDppU8zfs54sgtObKjnpEmT\neOaZZ46+fuihh3j88eMXVM3NzaV3b8/E6ZdffpnRo0dzzjnnkJaWxr/+9S+efPJJ+vfvz+DBg9m9\ne3e156x4jiPLVlf02GOPMXDgQPr27cuDDz54NI5u3boxbtw4evfuTV5eXl1+7cA6snDgWX9mxdmv\nMOYjN9e/tIidhcVM/lUfvvjD6Vx+cnKVSQB+af93iMPa/43xCppEUNmd4IkYO3Ys06dPP/p6+vTp\nVa5AWt7KlSt57733WLx4MX/+859p0aIFy5YtY8iQIbzyyitHt7vnnnvIyMggIyODa665xue4ytc/\nyM7OZsmSJcydOxeAdevWcdttt7Fq1So6duxYi9+2geQt8iwBkbeoyk0Wu7tw5ZpTGTnTRd6eQzwy\nqhdf3n06Vw5KxekIIbsgm6k5U6tM9Efa/3/T/zdVNgsZE2yCpmmovkeC+FqPoKIzzzyT6OhooqOj\niY2N5ZJLLgE8E9JWrFhxdLu6Ng1VVf8gNTWVjh07Mnjw4Fofs0HU0BGct/sQj360mtmrttO2ZTgP\nXNyTq09JPaYPwNr/jamboEkE/hgJ4ks9gorCw38ZvRISEnL0dUhICKWlpdXuGxoaSllZ2dHXldUk\nqKr+QW5ubuOuR1BFR/DB4lKe/fonpszbgEOEP5xzEjcP70Rk2PEFdqz935i6CZpEAPV/J+hLPYL6\nlJaWxqxZswBYunQpGzZsOG6bquofNHqVdAT/d1k+kz/5ge37ixmd0Z5JF/SgXWxElYew8f/G1E1Q\nJYL65ks9grq65557ePTRR4++XrRoEZdddhmvvPIKvXr14pRTTuGkk046br+q6h80+hKV5SrIrYvM\n4N5ZpSzdtJy+ybE8e80ATu7YusZD2Ph/Y+qmweoRnAirR9AwAv03LSgs4u+f/siMJfm0bRnOH8/v\nxuUDkgkJ+WUpCFvf3xjfBbwegTHHqGZpiLIy5Y1Fm5j8yQ8Ul7qZeHonfnNml+OWgva1M9gYUzuW\nCOqJ1SOoRjUjgjbtOsS9765gwc+7GNqlDY+O7kN628o7ta0z2Bj/aNKJQFWrXEGyoTX1egR+bSKs\nZERQWYeBvDw/l8dm/0hoiDD5V30YOzCl2v+e1hlsjH/4s1RlBDAXCPeeZ4aqPigirYG3gTQ8Fcqu\nUNU9tT1+REQEu3btok2bNo0mGTRVqsquXbuIiKh6RM4JqTAiKD/uZO769wKyNu7hjG7x/O+lfWgf\nF1njYawz2Bj/8FtnsXiuzlGqekBEnMC3wJ3Ar4DdqjpZRCYBrVT13uqOVVlnscvlIj8/v9Kx9Kb2\nIiIiSE5O9t9Q07xFuH+ey4f7OnHvokginA4evKQnl/bvcDSRW0ewMfUr4J3F6skwR5bDdHq/FBgF\nnOF9fxrwNVBtIqiM0+kkPT39hOM0DWNtWA/uySlhef4+zu0Zz6Oje5MQ88sTiHUEGxM4fu0j8Bau\nXwJ0AZ5R1YUikqiqW72bbAMS/RmDCSxV5ZUFG/nrR2toGRHK01f15+K+Scc151lHsDGB49dEoKpu\nIENE4oD3RaR3hc9VRCptmxKRCcAEgNTUVH+GaU5UFUND9x12MendFXyychtndU/gscv7HlMgpjzr\nCDYmcBpsQpmIPAAcAsYDZ6jqVhFJAr5W1W7V7VtZH4FpJKoYGro8by+/eXMpW/cW8cfzu3HzsE7H\nTAyrjPURGFO/At5HICLxgEtV94pIJHAO8DfgA+B6YLL3+0x/xWAaQIWhobphHv/ZGM//+2QNCdER\nvD1xCCd3bOXToWxVUGMCw59NQ0nANG8/QQgwXVVnicgCYLqI3ARsBK7wYwzG38oNDVWHk7/90Jbn\nf17NiB6JPD6mL3EtwgC72zemMfPnqKEVQP9K3t8FnO2v85oG5l0sbkv2HB7Kac2XufHcf1F3bhqW\nfsywUBsRZEzjFTQVyox/qCpTc9ty2oIBrArpzju3DOHm4Z2OGRVU39XhjDH1q0kvMWECq6S0jPve\ny+Hdpfmc2zORxy7vR2yL4yek2YggYxo3SwSmelUMDd1zsISJry1h0Ybd3DWiK3ee3bXKpT5saQhj\nGjdLBKZqVQwNXb/jADe9vJgt+4p46soMRmV0qPFQNiLImMbLEoGpWiWrhs4v6cQtry7B6QjhzfGn\n+FQ5zBjTuFlnsanakaGh4gBHGJ8f7sq4FxeRGBPBf28fejQJZBdkMzVnKtkFTXcZbmOCmT0RmKp5\nh4aWbZjHa1tTeOBLB8O7tuGZawYQ460eZkNDjWn67InAVOtQ4gBuyT2dB5ZFce3gVP5zw8CjSQBs\naKgxzYE9EZgq7TpQzA3/WcyqLft44OKe/Hpo2nEjg2xoqDFNnyUCU6mC/UVcM3Uhm3YfYsq4TM7u\nUflq4TY01JimzxKBOc6WvYe5ZupCtu8v4j+/HsipndtWu70NDTWmabNEEMwqmSy2adchrp76PfsO\nuXj1pkE2PNSYIGCJIFhVMllsfURPrpmykKJSN2+MH0yf5NhAR2mMaQCWCIJVhcliO3K+YOzS/agq\nb44fTI+kGMCWjzYmGFgiCFbl6giUhTi5e3FLHE54/eYhdEloCdgcAWOChc0jCFbeyWKbB/yece77\n+Sm8F9Mn/pIEwOYIGBMs/JYIRCRFRL4SkdUiskpE7vS+/5CIbBaRbO/Xhf6KwVRvoasz5y7OJC+q\nN9NvGULHNlHHfH5kjoBDHDZHwJhmzJ9NQ6XAH1R1qYhEA0tE5DPvZ/9Q1cf9eG5Tg6Wb9nDDfxbT\noVUkr998CokxEcdtY3MEjAkONSYCEUkAhgLtgcPASiBLVcuq209VtwJbvT8XisgaoOb1io3frdte\nyI0vLyYhJpw3xp9CQvTxSeAImyNgTPNXZdOQiJwpIrOBj4AL8BSj7wncD+SIyMMiEuPLSUQkDU/9\n4oXet34rIitE5CURaXUC8Zta2rz3MONeWoTTEcKrN1afBIwxwaG6J4ILgfGquqniByISClwMnAO8\nW90JRKSld5u7VHW/iDwHPAKo9/sTwI2V7DcBmACQmprq0y9jyqlkstjugyWMe3EhB4pLmT5xCKlt\nWgQ4SGNMYyCq6r+DiziBWcBsVX2yks/TgFmq2ru642RmZmpWlo1Y8Vklk8UOJgzg6qkL+WHrfl69\n6RQGpbe2OQLGNHMiskRVaxzl4UsfwZ3Af4BCYCqeJp5Jqjqnhv0EeBFYUz4JiEiSt/8A4FI8fQ6m\nPlWYLFb681xumaOs3LyP5689+WgSsDkCxhjwbfjojaq6HzgXaAVcB0z2Yb+h3m3PqjBU9O8ikiMi\nK4Azgd/VMXZTlXKVxdQRxlM/JTJv3U7+36/6cE5PzyqiNkfAGHOEL8NHjyxAfyHwqqqukoqL0ldC\nVb8tt295H9ciPlMX3sliumEeU/M78PSKaCZd0J0rMlOObmJ1BIwxR/iSCJaIyBwgHbjPOyeg2qGj\nphFIGcTTa1vx5Iq1TDitE7ec3vmYj22OgDHmiCoTgYg4VdUF3ARkAD+r6iERaQP8uqECNHXz2vcb\nefKztVw2IJlJ53evdBubI2CMgeqfCBaISD7wKfCpqu4FUNVdwK6GCM7UzbfrdvLAzJWc1T2ByZf1\nISSkxpY8Y0wQqzIRqGqmd3jn+cA/RaQD8C3wCfCNqhY3SISmVjbtOsRv3lxK14Ronr6qP06HrSto\njKletVcJVc1V1edVdTRwKvAhMAKYJyIfNUSAxncHi0uZ8GoWZWXKC+NOJircVhk3xtTM5yuFt7/g\nS+8X3icE00ioKvfMWM7a7YW8/OtBx60kaowxVamx3UBELhaRZSKyR0T2i0ihiOxX1c0NEaDxzbNf\nr+fjnG1MuqA7p50UT3ZBNlNzppJdkB3o0IwxjZwvTwT/BH4F5Kg/16MwvquwjtBXPxTw+JwfGdmv\nPeOHd7JZw8aYWvElEeQBKy0JNBIV1hHaPOot7ni3mB7tYvjbZX0RkUpnDVsiMMZUxZdE8EfgYxH5\nBjg6UqiyReRMAyi3jpC6S/h01gycjlG8MO5kIsMcgM0aNsbUji+J4K/AASACCPNvOKZG3nWE1F2C\ni1A+OdCFZ24cQHKrX5aUtlnDxpja8CURtK9pmWjTgLzrCH33+fs8uTaBkRePYkjnNsdtZrOGjTG+\n8mW20ccicq7fIzE+m70/lWt/HEanAWdx/alpgQ7HGNPE+ZIIbgU+FZHD5YeP+jswU7m83Ye4e/py\n+iXH8ujo3viwEKwxxlSrxqYhVY1uiEBMzUrdZfzubc+8gH9dPYAIpyPAERljmoPqitenVbejeCTX\nd0Cmas99vZ6sjXt4ZHRvUlpbvWFjTP2ormnoMRF5V0TGiUgvEUkQkVQROUtEHgG+A3pUtbOIpIjI\nVyKyWkRWeUteIiKtReQzEVnn/d6qnn+nZik7by///GIdozLaM7q/re5hjKk/1RavF5GewDV4yk4m\nAYeANXiqjM1Q1aJq9k0CklR1qbeYzRJgNHADsFtVJ4vIJKCVqt5bXZDBXrz+YHEpF/7fPErdysd3\nDmdD4SobGmqMqVG9FK9X1dXAn+sSgLdA/Vbvz4UisgboAIwCzvBuNg34Gqg2EQS7hz9cRd7uQ7w1\nYQgbClfZ8hHGmHrVIIvVe/sb+gMLgURvkgDYBiQ2RAxN1Sc5W5melc+tZ3RmUHprKzpvjKl3fl+w\nXkRaAu8Cd6nq/vLDHVVVRaTStikRmQBMAEhNTfV3mI1HuQXltsX0ZdJ7OfRNjuWuEScBtnyEMab+\n+TURiIgTTxJ4XVXf8769XUSSVHWrtx+hoLJ9VfUF4AXw9BH4M85Go9yCcuoI47lW/0tJaXueuvKX\nSmO2fIQxpr5VV7x+QHU7qurS6j4Xz63/i8CaCgvUfQBcD0z2fp/pc7TNXfkF5UpLiNryPQ+Oup/0\ntscWmbHlI4wx9am6J4InvN8jgExgOSBAXyALGFLDsYcC1wE5InKkOsqf8CSA6SJyE7ARuKJuoTdD\n5RaUKy5zoGlDGTswJdBRGWOaueqK158JICLvAQNUNcf7ujfwUE0HVtVv8SSOypxd60iDQcogiq9+\nn9fefp35pT14/OqrbAkJY4zf+dJH0O1IEgBQ1ZUiUuVEMnNi/jcnmmn7LuDVmwbRKspW/TbG+J8v\niWCFiEzBvKTpAAARsklEQVQFXvO+vgZY4b+QgteSjXuYtmAjN5yaxvCu8YEOxxgTJHxJBL/GswLp\nnd7Xc4Hn/BZRkHK5y/jTezm0j43gnvO6BTocY0wQ8WX10SIReR74WFV/bICYgtLUeRv4cXshU8dl\nEhXu9+kdxhhzVI0zi0VkJJANfOp9nSEiH/g7sGCyadchnvpiLef3akfbtluZmjOV7ILsmnc0xph6\n4Mut54PAIDxrAqGq2SKS7s+ggomqcv/MlYSGhHDFMDfj59xi6wgZYxqUL2sNuVR1X4X3gmOmbwOY\ntWIrc9fu4O5zT+LnAytsHSFjTIPz5YlglYhcDThEpCtwBzDfv2EFh32HXTz84Wr6Jsdy3ZA0cnba\nOkLGmIbnSyL4LZ6lqIuBN4DZwKP+DCpY/P3TH9h9sJiXfz0QR4jYOkLGmIDwZdTQIeDPIvJX78+m\nrsqtLLqkrCuvL9zEzcPS6d0h9ugmto6QMaah1ZgIRORUYCrQEkgVkX7ARFW9zd/BNSsVVhZ9Nfwv\ntI/tzO/OOSnQkRljgpwvncX/AM4DdgGo6nLgNH8G1SxVWFm0/d4s/jKqt80ZMMYEnE8VylQ1r8Jb\nbj/E0rwdWVlUHBSrA00bxoieVpzNGBN4vtyO5nmbh9RbaOZOPAXsTW2kDELHzeSdd99i5t5OPHGF\nrb5tjGkcfHkiuAW4HU/h+S1Ahve1qaVZe1L44/YRnHPuxbSLjQh0OMYYA/g2amgnnhVHzQk4XOLm\n0Y9W06eDZ86AMcY0Fr6sNdRJRD4UkR0iUiAiM0Wkkw/7veTdfmW59x4Skc0iku39uvBEf4GmYuq8\nn9m+v5gHL+mJI8SKzRhjGg9fmobeAKYDSUB74B3gTR/2exk4v5L3/6GqGd6vj30NtCnbUVjM89+s\n57xeiWSmtQ50OMYYcwxfEkELVX1VVUu9X6/hqWNcLVWdC+w+4Qibgae+WEtRaRkjT3HZyqLGmEbH\nl1FDn4jIJOAtPIvNjQU+FpHWAKpa24v9b0VkHJAF/EFV99Ry/yblp4IDvLkoj/MHFPHQojtsZVFj\nTKPjyxPBFcBE4Cs8S1HfClwJLMFzMa+N54BOeEYebQWeqGpDEZkgIlkikrVjx45anqbx+NunPxDp\ndHBSWoGtLGqMaZRqTASqml7NV42dxhWOtV1V3apaBkzBU+egqm1fUNVMVc2Mj2+a9XsXbdjNZ6u3\nc8vpnTg9dTBhjjAc4rCVRY0xjYovaw2NAT5V1UIRuR8YADyiqstqezIRSVLVrd6XlwIrq9u+KVNV\n/vfjNbSLieCmYZ2IDHPYyqLGmEbJlz6C/1HVd0RkGDACeAx4Hjilup1E5E3gDKCtiOTjqXR2hohk\n4OlryMXT5NQsfZSzley8vfz98r5EhjkAW1nUGNM4+ZIIjqwrdBHwgqp+JCI11iNQ1asqefvF2gTX\nVBWXuvn7pz/SvV00lw1IDnQ4xhhTLV86izeLyL/5ZbRQuI/7Ba3Xvt/Ept2HmHRBd5s8Zoxp9Hx5\nIrgCz8Swx1V1r4gkAff4N6wmyFt05kDSEJ7+8gDDurTl9JOaZie3MSa4+Fqh7L1yr7fiGfppjihX\ndCZcQulUdB/3XXgjIvY0YIxp/KwqSn0oV3RGypSbkjfTq31szfsZY0wjYG399cFbdMZNCC5CGXTG\nqEBHZIwxPrNEUB9SBrH+gjd4wjWGGb2fJb7n8EBHZIwxPrOmoXryUHYUK8Mv55tLzgx0KMYYUyv2\nRFAPlm3aw7x1O5l4emdiIpyBDscYY2rFEkE9eOarn4hr4eTawR0DHYoxxtSaJYITtGrLPj5fU8BF\nmSW8tfZlqzVgjGlyrI/gBD3z1U9Ex+Qze/cUXDtcVmvAGNPk2BPBCVi3vZBPVm6jf7fduNwuqzVg\njGmSLBGcgGe/Xk9EqIObM0dYrQFjTJNlTUN1tHHXQWZmb+bGoemc1rEnUyKt1oAxpmmyRFBHz329\nnlBHCBNO8xRps1oDxpimypqG6mDz3sO8uzSfKwemkBATEehwjDHmhPgtEYjISyJSICIry73XWkQ+\nE5F13u+t/HV+f3rhm/WowsTTOwc6FGOMOWH+fCJ4GU8dg/ImAV+oalfgC+/rJqWgsIg3F+dx2YBk\nOsRFBjocY4w5YX5LBKo6F9hd4e1RwDTvz9OA0f46v79MnbeBUncZt55hTwPGmOahofsIEr2FbQC2\nAYkNfP4TsvtgCa99v5GR/dqT1jYq0OEYY0y9CNioIVVVEdGqPheRCcAEgNTU1AaLq1LeMpSfFqRy\nqKQFt5/ZJbDxGGNMPWroRLBdRJJUdau39nFBVRuq6gvACwCZmZlVJgy/85ahVHcJl5Y5yO/8BF0T\nowMWjjHG1LeGbhr6ALje+/P1wMwGPn/tectQirpxUsq49nmBjsgYY+qVP4ePvgksALqJSL6I3ARM\nBs4RkXXACO/rxi1tOOoIo5QQ3CFO2vU9J9ARGWNMvfJb05CqXlXFR2f765x+kTKID/o9x4/ff8Ko\nS6+gW8qgQEdkjDH1ymYW16DUXcbfVsawrOONdMscEehwjDGm3lkiqMFnq7ezreRHUjstsKIzxphm\nyRJBDZ5d8AVRHafy6eaXGT9nvCUDY0yzY4mgGmu27ufHfdmIuK3ojDGm2bJEUI1p83NxlHQhzOG0\nojPGmGbL6hFUYc/BEv6bvZlL+5/K2GGDreiMMabZskRQhbez8ihylXH9qWl0T4ixBGCMabasaagS\n7jLl1QUbGdypNd3bxQQ6HGOM8StLBJX4fM12Nu89zA2npgU6FGOM8TtLBJWYNj+XDnGRjOjRpFbJ\nNsaYOrFEUMGP2wqZv34X1w7uSKjD/jzGmObPrnQVTFuQS3hoCFcOTAl0KMYY0yAsEZSz75CL95du\nZlRGe1pFhQU6HGOMaRCWCMqZnpXHYZeb662T2BgTRCwReLnLlFe+z2VQWmt6tY8NdDjGGNNgLBF4\nffVDAXm7D9vTgDEm6ARkZrGI5AKFgBsoVdWAL+AzbUEuSbERnNvLhowaY4JLIJeYOFNVdwbw/Ef9\nVFDIvHU7uee8bjhtyKgxJsjYWkN5i9jw8QwGhXbgyoFWgcwYE3wClQgU+FxE3MC/VfWFgESRtwid\ndglnuko43ekkbM8QaGk1iY0xwSVQ7SDDVDUDuAC4XUROq7iBiEwQkSwRydqxY4d/osidh5aWECpl\nOCmF3Hn+OY8xxjRiAUkEqrrZ+70AeB847jZcVV9Q1UxVzYyPj/dPHB2HkRUeyb9jY1ke2QLShvvl\nPMYY05g1eCIQkSgRiT7yM3AusLKh4wB4r9DFze0SeLZ1LOOTEsgOt9nExpjgE4gngkTgWxFZDiwC\nPlLVTwMQB9NXfk2ZlFEGuLTM6hEbY4JSg3cWq+rPQL+GPm9FRS43azYk4EhyIuK2esTGmKAVtMNH\nv1hTwP69HXj43CdwOX+yesTGmKAVtIlgxpI8kmIjuDbjdBwhZwQ6HGOMCZignEa7fX8R36zdwWUD\nknGESKDDMcaYgArKRPD+ss2UKVx2cnKgQzHGmIALukSgqsxYkk9mx1akt40KdDjGGBNwQZcIlufv\n46eCA1xuTwPGGAMEYSJ4JyuPCGcIF/VNCnQoxhjTKARVIihyuflg+RYu6J1EdIQz0OEYY0yjEFSJ\n4LPV2yksKrVmIWOMKSeoEsE7S/LpEBfJkE5tAh2KMcY0GkGTCLbtK+LbdTu4bEAHQmzugDHGHBU0\nieC9Zfk2d8AYYyoRFIngyNyBQWmt6djG5g4YY0x5QZEIlm7ay887DnJ5pj0NGGNMRUGRCGYsySfS\n6eDCPjZ3wBhjKmr2iaDI5WbW8i1c0KcdLcODdrFVY4ypUkASgYicLyI/ishPIjLJn+eavWobhcU2\nd8AYY6oSiJrFDuAZ4AKgJ3CViPT01/lmLMknuVUkg9Nt7oAxxlQmEE8Eg4CfVPVnVS0B3gJG+eNE\nXy+aBnv/wuj05TZ3wBhjqhCIRNAByCv3Ot/7Xr3KXvkGd696jBVtNvLWwRfJXvlGfZ/CGGOahUbb\nWSwiE0QkS0SyduzYUev9s36ejUugTIRS8bw2xhhzvEAkgs1ASrnXyd73jqGqL6hqpqpmxsfH1/ok\nmZ3OI0zBoYpTPa+NMcYcLxDjKRcDXUUkHU8CuBK4ur5PktH7aqbgeRLI7HQeGb3r/RTGGNMsNHgi\nUNVSEfkNMBtwAC+p6ip/nCuj99WWAIwxpgYBmWGlqh8DHwfi3MYYY47VaDuLjTHGNAxLBMYYE+Qs\nERhjTJCzRGCMMUHOEoExxgQ5UdVAx1AjEdkBbKzj7m2BnfUYTkOwmP2vqcULFnNDaWoxVxdvR1Wt\ncUZuk0gEJ0JEslQ1M9Bx1IbF7H9NLV6wmBtKU4u5PuK1piFjjAlylgiMMSbIBUMieCHQAdSBxex/\nTS1esJgbSlOL+YTjbfZ9BMYYY6oXDE8ExhhjqtGsE4GInC8iP4rITyIyKdDx1EREXhKRAhFZGehY\nfCEiKSLylYisFpFVInJnoGOqiYhEiMgiEVnujfnhQMfkCxFxiMgyEZkV6Fh8JSK5IpIjItkikhXo\neGoiInEiMkNEfhCRNSIyJNAxVUdEunn/tke+9ovIXXU6VnNtGhIRB7AWOAdPOczFwFWqujqggVVD\nRE4DDgCvqGrvQMdTExFJApJUdamIRANLgNGN/G8sQJSqHhARJ/AtcKeqfh/g0KolIr8HMoEYVb04\n0PH4QkRygUxVbRJj8kVkGjBPVaeKSBjQQlX3BjouX3ivd5uBU1S11nOumvMTwSDgJ1X9WVVLgLeA\nUQGOqVqqOhfYHeg4fKWqW1V1qffnQmANfqg/XZ/U44D3pdP71ajvhkQkGbgImBroWJorEYkFTgNe\nBFDVkqaSBLzOBtbXJQlA804EHYC8cq/zaeQXqaZMRNKA/sDCwEZSM28zSzZQAHymqo095n8CfwTK\nAh1ILSnwuYgsEZEJgQ6mBunADuA/3ia4qSISFeigauFK4M267tycE4FpICLSEngXuEtV9wc6npqo\nqltVM/DUyx4kIo22GU5ELgYKVHVJoGOpg2Hev/MFwO3eps/GKhQYADynqv2Bg0Cj71cE8DZjjQTe\nqesxmnMi2AyklHud7H3P1CNvO/u7wOuq+l6g46kN76P/V8D5gY6lGkOBkd729reAs0TktcCG5BtV\n3ez9XgC8j6e5trHKB/LLPR3OwJMYmoILgKWqur2uB2jOiWAx0FVE0r0Z80rggwDH1Kx4O15fBNao\n6pOBjscXIhIvInHenyPxDCb4IbBRVU1V71PVZFVNw/P/8Jeqem2Aw6qRiER5BxDgbWI5F2i0o+FU\ndRuQJyLdvG+dDTTaQQ8VXMUJNAtBgGoWNwRVLRWR3wCzAQfwkqquCnBY1RKRN4EzgLYikg88qKov\nBjaqag0FrgNyvG3uAH/y1qRurJKAad5RFiHAdFVtMkMym5BE4H3PvQKhwBuq+mlgQ6rRb4HXvTeO\nPwO/DnA8NfIm2XOAiSd0nOY6fNQYY4xvmnPTkDHGGB9YIjDGmCBnicAYY4KcJQJjjAlylgiMMSbI\nWSIwxpggZ4nAmEp4lyS+rZrPI0XkG+98hKq2+VxEWvknQmPqjyUCYyoXB1SZCIAbgfdU1V3NNq/W\ncAxjGgVLBMZUbjLQ2Vvw47FKPr8GmAmeugwiMte77UoRGe7d5gM80/+NadRsZrExlfAuqz2rsgJB\n3iUINqlqO+/rPwARqvpXb1NRC299BkRkHTBYVXc1WPDG1FKzXWvIGD9qC5QvWrIYeMm7Eut/VTW7\n3GcFQHvAEoFptKxpyJjaOwxEHHnhrSx3Gp5lzl8WkXHlto3wbm9Mo2WJwJjKFQLRlX2gqnsAh4hE\nAIhIR2C7qk7BU05ygPd9AdoBuQ0RsDF1ZYnAmEp42/S/83b+VtZZPAcY5v35DGC5iCwDxgJPed8/\nGfheVUv9Ha8xJ8I6i42pAxEZAPxOVa+rZpungA9U9YuGi8yY2rMnAmPqQFWXAl9VN6EMWGlJwDQF\n9kRgjDFBzp4IjDEmyFkiMMaYIGeJwBhjgpwlAmOMCXKWCIwxJsj9f44jvfeunbY3AAAAAElFTkSu\nQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x73ae790>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "v0_num = 0.0\n",
    "v_t_example = v_t.subs(m,1576.0).subs(vd,40).subs(v0,v0_num).subs(k,1000).subs(d,0.5)\n",
    "\n",
    "#expr.subs(vd, 40.0)\n",
    "#expr_num = expr_num.subs(mc, 1576.0).subs(dc, 0.5).subs(kc, 1000)\n",
    "\n",
    "dv_dt_example = dv_dt.subs(m,1576.0).subs(vd,40).subs(v0,v0_num).subs(k,1000).subs(d,0.5)\n",
    "\n",
    "h=0.2\n",
    "t_range = np.arange(0,7,h)\n",
    "num_steps = len(t_range)\n",
    "\n",
    "xs_c = [v_t_example.subs(t, x) for x in t_range]\n",
    "xs_ee = Explicit_Euler().simulate(lambda v_val,t: dv_dt_example.subs(v,v_val), \n",
    "                                      v0_num, 0.0, h, num_steps)\n",
    "xs_ie = Implicit_Euler(10).simulate(lambda v_val,t: np.matrix([[dv_dt_example.subs(v,v_val)]]), \n",
    "                                        np.matrix([[v0_num]]), 0.0, h, num_steps)\n",
    "\n",
    "plt.figure(1)\n",
    "plt.plot(t_range,xs_c, label='v')\n",
    "plt.plot(t_range,xs_ee, 'o', markersize=3, label='v_exEuler')\n",
    "plt.plot(t_range,[x[0,0] for x in xs_ie], 'o', markersize=3, label='v_imEuler')\n",
    "plt.xlabel('t (s)')\n",
    "plt.ylabel('speed (m/s)')\n",
    "\n",
    "plt.legend()\n",
    "\n",
    "\n",
    "pp = PdfPages('car_example.pdf')\n",
    "pp.savefig()\n",
    "pp.close()\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
