{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cruise Control Example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, import necessary libraries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "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 symbolic manipulations. This is helpful to obtain the analytical solution of the cruise control example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "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 symbolic equations for the example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAKkAAAArCAYAAAADraj8AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAGfUlEQVR4Ae2djXXUOBSFMzkpgIUOQgewdAAdhN0KFjogJXCyHcBWsJvtIKGCQDpIqICfDsL9PJIiazyOZcsTO346x7EtPT09XV09SR7LWd3c3OxZKIfAarU6lLaXwvVjOa0PX5NwO1ItL4XbdVrb/TSi5L0KfqbjTMcPHR9K6i6tS/Y9imy96qMfHcr3wQjajp5w2uCFMPtfuU6URievhZUSaxFj3KhgCnntDBmjiGI6ZesXKfssW9/mKnV5qeeGN8jVtQT5lBe6p5N/En7P4/qP6kkpSAW/dAWexwWXvlY5RzreFdD7TDrOcvWo7DfKA7knT9CCWOXCFORlwwYvhN1PCZyn7XgQco138Uqqr50B45VSQHMTcBlqjyVLXS10Q2AbL94r+1cdf3s1o3tSFUSPGdWL+soUOG8DrlW1yF3No+bgRVsrstvERl44Z3YdOYy9op7UNRYeBbf9TQfkZPikd1RBMsw7GBqf6riSUaHHKI0Fy4ni7mtlXANO9rDifKyDORJ2bRvKkRulI04Nrz72KA+duJUXSo8DWL7WscZUwO+VOKSQBv6h49Dr0zVzOxZNj6I4GltRa3kf7+Ig6Wkc1/Va+SDKu67yqZzy0nkw7Ig0Bea3gEu9iHiT5vH3SqOeW9O9XJ+z9BbHSzp7Y5Vrj+Q78SLGRnlwYl98XJHh3vWuUyk+luLY23DNsy88K4soGv2CawV6yvfq6vYPvS3Of5sy/hVgEvzE/aOry2fFvdV1m3enXmld0DUoTA2vXHu68qIBJLDEaVThgL9O2ad1VOe/f6nhLp30ic4o/S/JTcPz/MuH78rj7/9QZJgGOAHI7Ens89TOspWyPKHiNMp/rPQ/40h3zcKNTtEWmI9SPp4rPH7SNXFtBEUnUwLkNsJAbAfhNQJWufZ05UWKG46Kjr8OaoTBw700Mcyfxbp0D2kYJvn1pVYGcS4tTAOQURgyXPcewlzZ1VQDG3QwfHe2RbLUf6Oeab373qNbBwAVwUt6hmLVyR6Vk8ULj4/ysY7R7Zo3B7oZFJyngJDeq3p9VISSmhYU21bRT3zmXZ5dHei5DOsM93j7K52ZqjTZn5pXG57SxAL3k8JL9bnTnp688FAxMoVp376PLXBOh2mG3Yq4MpgH7bfue+3Ka6RGRvL/FrCjjwqmHnGHgnQEOh/TIYattgCgADtWALsp4ZVjTw4vPH7oD9OnfR/b9yxPgzK8DYqroEZldYbL9r3hheT8NTJcx/KQAZlaQyC4o4BnCB7T1YmiY5vbTMHu2k95bcI90qaG15329OSFhwYsQ3sU+e1epIRk/+jwvYYCIC8eiLiNYVN5ThVPZZkL4sXuWpwgtjVIH56Yx1/huetW4SRBeZmD8mKIX9ThPZmbVs9ydWalT30ag2SZ2pAf+VGCyiiGl3T1xspXros9ksnmBfqVj/cnwsK8CEm94fd5LgH8EPtVPp3tlYja1fsOKW5Q3vvGqs142QaxeUYaOvxDIimVwyNv9Xht4AxNE7hMcZ6r/PD4aqjOsfI7ItwbVm31km2Mvt+EYxgRHwxJ2yq+qzQBzLSBJwST96a7wiSnHNd5dv+qXo6RD0CWHwwm/XL3xDFm3r3xo4t50sKtJm/AUwvbPpKJq3BjunTeNAoZSTPBNPHdIzD4OenuTbYSl4aAkXRpLT7D+q5kc/VmxwxtN5MXgoDNSRfS0HOupg33c269hdhuJF1IQ8+5mkbSObfeQmw3ki6koedcTSPpnFtvIbYbSRfS0HOuppF0zq23ENuNpAtp6DlX00g659ZbiO2DtzQvBKdQTb1Sxn4mNu4d6uCLKwT2DBFe6LjgrXLJsRGRHbN8E4t49kCFzWW6t9ARASNpR6AiMfYxHYuE1eY9xfNRjGqrg+LYwsJXrRFn82FFYt1DYl7o/Y0EC3kIGEkz8BLZ8I4XLgue9KcnqIvze++fKD72msRXe7CcnJ0yEDCSZoAlUb4pdek8JiRl2I8DJCak37iq7etfi9jfrgjYwqkrUpITQf1OVOal3MfekijI2PRVa+SZHljogYCRtAdoygIZm762AhnDByZQLa9LHEN9Fe+8MEkWOiJgJO0IVCLGt6NqXjSaAqQek92PLKL4xDZThMoLJ/rstgUBI2kLOE1Jjmh4xpSMFfkapgC/S9YTmj35NU/bVIbF1REwktbx6HKHN2Te6Ynn8xDf9D0rFlH8IzO27NqefI9WxvkX/4Nyhllp2tkAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$\\displaystyle \\frac{- d v + k \\left(- v + vd\\right)}{m}$"
      ],
      "text/plain": [
       "-d⋅v + k⋅(-v + vd)\n",
       "──────────────────\n",
       "        m         "
      ]
     },
     "execution_count": 3,
     "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": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAI8AAAArCAYAAABB7ttrAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAGOUlEQVR4Ae2cj3XcKBDGvX4uwC/pwOkgfzpIOnAuFSTpICnhntNBchXcSzqwr4Ik7sC+Cnx2B77vh4EFCa2EvCuxu/AeFhoGmG9mGJAsdnF/f39Q03o0sFgsTtTTa+n023p6nLcX4TmVBJfCc52S5NARxXisfK58q3zl6Ju6aoznwXhfNzXOVP0Ky7HG+jrEcebEjpyB3lfaWVh+CNOZ+JkU7SSGgzCL47cySojom7rXWIS+0031P1W/Vm8nOePNid3K22tn8TEpfqdw+cgjBpeeq3DubjZ5lUe/tv1fbHKcTfctHB80xi8pOBneU+M/Brvanip/SvWbQRtkZ2G6U58XqfGOwsEeAyjsJ6P8RrzXVsCMZsWxfpZEYMlJs2EfYec/Bexf5S8hwGbkmRoQkWfbo47ZD+REHWuAObFn2dlO7uvA6QyEKPKIEgESM7vtJ8ovlM+cgkRnHSRUP1O+Et17pOrYhMEbPXGIjpKZoYTBG2WchtCJV8+axuAJBEZHKydAgdgH2TnASBGMb+2Ve+1W7cZYRRzCb15VZk3F4AwE/UPAi3OIZOpuHd3ScJ7vDRp93Cr7DaXK7Kvo5DjknaMsGbLwhDJaHF43YR1lpbVjV5847KfmWEPu1W6wncP+1I5gEW2cj0RwCZAktzn6psZ3mjX/ifZRZRNJ7Cz6aTgfPJH6MBFdXjmCndXfdf9ZfYQbSsq8QyASzZZy8SQEZYI1dWDYCsU+yM4JnGDE8ZZJxjPRRxTetRAdVj6+qd5HCssfzQDREM4/eqtMv60II5pb3nz0c7JMeZUcWXiasqk9OuPFYAuH6I/CrvZnyrw6aWZ0x7hNOvdR1G/KpfpBdk60Y4sh8hKnB6wKs9zoynLFkhI5RdiIshJOknKKqJ14AHkette9C51JpYe8U5Ul0yA8TXksviQOW7d27Or3MctWlp0dXo3Zcp4jEQ9seCX8sjyxbPFm8UpXlpWuzWDXjv0pfZJsvzjKpSEs/2AoHLer7yXndKVePB2itMO5GEvEbmXKtbODzYNTuO04OLQ1f3ANjOnWcLPGaVDCZzMhROQU4mNG/N1k1L3bI7mqdyqYtrRRpq+5Uw6eUFYUimK7UknYx9jZ4UI/d+6G66G9Ydb5KCAnckyRp1led6HOG10OgKO9UlvvULYf+g352LUTAl3ftHFlkWdLvXg6JAMvrzKiVCj2MXZ2uMDofQTiQiAJsexx2CizXJkkGnsf8x5HV/Pk9VCz/CsenqJQOusokSt6twPNOtVfKroZiAA4J9EM2qqlUdXTpSF4mtKoDUswukNXUdoUdvVLhOe1h3+/Fg3ccaN2o+xMd2rLZvy9xvTBwThPx1iVPFADUiyT540UO0kEHes8A+G02DQeqwrveKIJUp2npap8gpTLUvxCyv2Y3zq/hTUmkd5tL/I7yWih8VglbjReFOmq82QocRWrFMySwNPqJNFnlSzrrLOO+o9wtfZ1bsO8zvH2tS/+78MLuF1L7GvB1ko18rRUMp6gWcpT5S59hspyfNEVTavzjPeVvW9Zl629d4HxCqjOM153e99yIQ2Y/3LOrQmtq8jSmdyuv5MhXRG91EqxqN8i8KdkK51W9zylW6hg+eqyVbBxShetOk/pFipYvuo8BRundNGq85RuoYLlq85TsHFKF606T+kWKli+6jwFG6d00arzlG6hguWrzlOwcUoXzRy96RNSr/D5TpePp0+UORFK4jtaEqdDf+rfC1/Ex4ftnIy4sXS+7Y0+mha9JquBrderjMvnjCuzsLqz3HwtR/YH+1Tm+1b+P8QH8/7wm8o4V3SOvW+cfauXfrZar72Rx0YTd/KByHMnI4ffsrozS09FD6MMdByrpoQGdkGvvc4j3Pz40qXA4gg4D8tXmFiqSM2fSonOCD2w1L+BBrZer70bZjmO+0KffQ/LWxhdIOEkqV/3gp8lrqaEBnZBr73OE+DGSfyBr4COk/jDgtDtRpBIZeg2alFVU1sDW6vXHOfhnHMUdYKlrBlh+Nre/H6veFjqTNRq661SpIGt1esg57EOQCRpOknXUvZSvM7ROMsURabqMg8a2Ha9DnIeQSV6sK9xDuHsD711Pl00Ns/8WDRHN3bxLJNgrSVttV7/B/9NG/o62/GKAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$\\displaystyle \\frac{k vd}{m} - \\frac{v \\left(d + k\\right)}{m}$"
      ],
      "text/plain": [
       "k⋅vd   v⋅(d + k)\n",
       "──── - ─────────\n",
       " m         m    "
      ]
     },
     "execution_count": 4,
     "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": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASUAAAAyCAYAAAAaw83gAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAN6UlEQVR4Ae2djbXVNhLH8zgpIJtUsKQDAhUEOiCkAqADOFSwh3RAqCCBDshWAKEDSAWB1wH7/ykar6xr+/pDsmVfzTl+19bHzGhmNBp92O/qy5cvX1XIK4Grq6vnovBGsv4jL6WKfQ0JSJ/fiM5L6fOnNehdGo0bl9bgtdvrHdJX1SGtLfl89KTLa2F/Id3+mY/K5WIe7ZQYHXS90fVZ14ecIhP+WwGtFzlp5cStNjwR/psy4qc56fThFv2qsz7hLEz3g0x1TAvl2FX9SsLtSu9Nk6EzOrxTvce9hRJliBbM/SRarxOhXA2NeL8vYoT4/1qNaA+hqrMewSRIlmzfCM176XmTgSdBE4pDMTpSCji/pXsUkRWk7LuewO7WYcS7W3MQ/6WsOVSd5bNWdPwosNd8lC4E8ySnFAh+DUdxTzr4qBGI+fve4JUY/l28ryGnQdlUnQ2KZ3Gmt8+HQoTOKySQwCSnJHprOgoipc079VQZywkwbYP3UsL5qrOpSpxYXo6J5YVP0j27rLsEBi9dRNQnoHTWznrXkZU32G5f/2RTQOlPdDGraMFUp9RyFEJ4XxehK0zfBLN+WVyFGGks9Dag5w+6HjUJ/kZpN3VR/rku6iKcVaYcMS8JnlHQrwVFeGd1Rpsl80l6U/kj6SyB2t0ghO26fpACYWocXmcnEZ3n+ZZs9n0PTeoMBQgnjiXCw3LPuyiNHelflHbq0JRB5tlLlSHMwvN9ygvczpJ+MXoSHvn05/6X9M8hbj3jbV9Faa6c0tmlcnzonkaA8xtL28Ov+MXhwnfTli35Rn6en0GdwaNgtN5U9jA6S6kfyQX7fpESZ0pc4g37fBLjVFqrT3bk4zhc/47zeBYMtpl8Xc4G4/qkx3lnnZEh8ZWdo9A9Dsk5DH51mUNilLAOACMfrL5nHgac8ftn6n7W1Wqwnqn7Z1h3D/fiGaPkkORoueYsK16Q96DOoC8YrTeVPZTOUspfsrFBqbjBVLzhWJzThE9rt+7R54lToYwu+jkX9XoHWuWd1Df83r5cfZVjMANvq7yeW05xdOcBkS4cSAthRLxRhi/b8spKg6nGY4JLl+s0ER4a0TivMK/Ue/HLdJO23C2FR/FyVmfwKhitN4/zEDpLrSfkqAvZtOw+NZ25+MRXK0jwuqdPtvjVcxw8xDMeHBy2ZZc5O3tu+rin4QIPlWfwQ0aNU/T5rYF8ilOCMHNLPCfTq1ZDQkEpj4Z2GW6rjsrAbJuh/yu2mM4dtq3vXu1wDqAvf4t08TRaZ/AnOKs3lTmMznLoRPJhQbe4KF884QxO+FIakUvsRGiDi4z0y2DbimRiuSl/KFDBpsCH3+gsR36I84YSzoIWwmgQXg6kLE5xcJJFaQh2Qd+Oz3dW2OMEb7y45nCKztDCmqEp6feBmCmG5xk6Q5aDejugznLYz29CyhsJ9JeSgH7lFpvFG47I4KNuYl7pl598Aeq9VZ2+vu6L9f7g1H7zfuM2NqSLSCwEo+XSRjkllaTDEVVZpzMkMM/ODeFcCDSy5Ww8IygshrdRws96dnWpoysWWFR8+0fxiOCRRfZDpRNaO1VnoB6rt93rbIIcpxa1tw/mduKp9MaWtz5FxGL9mLo4qjvcBEDQ8Ux2bc7rez3P7YfUNZkwm3im61pXCK3nr8OcgXtG0KYhck7XYpjieNkuIL1RisrSYe+onjGHgwMHOJvGeiHQwc2hteoovVSwtjYyKoDRqTqD5UG9HUxnWVQkGX2UHdPJkP+vWYjMQApfqnbyapjXaQuj0rDjJLYsXA1N3Z/IQ7IiamoFK6PefVNFIgCmbo1TURoeFy/IusXJuRzls/6EIMjHCXUxhLN6qctGXgSBQom8SOOdoiTCEa5s4OVzW7xu/p6bNXKOzqh7Tm/KP4TOTE45fr3sWZOhfxQP4pfAgHUllmYmg+oPnXEaxKe6+JXGcVF4lFMaxFoz6cgs/vKSMqNjhQuXgOyBQZVNHTed2IM4xLOL9sXzakGAaBLYnAQ0Y6dve5DrJjxKsEQOXESFFaoEkICbHcg2ZkcQKcUoPtzW6hicKjumWMoybJi18F2sU5IgmMsSYs8KWQMp2pqYM8Qgvd5mkEBCvWXgrkHJ4jFA9GHroy5h7h+1G1xj3qfkUz8sgTSwp4gNpi/WKTUaW35jTqlGSstleRQMZgvJ1pT8tGq1qdWWirixJfGD0P7Wt8MM8SDNqs2YK4EgUjHbmIvqIuvVSGm52llPAloh8z9J9W8pEtD0Bz1x7sbtGMtxNNN25TH15rWmkx3ihfybbZxF4/njDM/fvjB8wtPFDXbVKZ01l7MF7JT6p7Mla4EtJfBMHfypX5vhuErjlDxTzbkilWFXCL3iIL6jni8z5YdBalSk5HniYOE9c0JKY1mArfI5tFUtHYgXnPkPujgaZG3iGZ7hk8OXnNpOsn5WnZKkuRBGj4YL6dTqMyXgO7idhePztfEAQsd3p5p9B2wckZ55q4B3v6g3Fc7ahnBTBif5MHBIpMETnb4EQF5Ekzh2nJEdg+G91ddqA0k/60rjlIRw9HYhlEsDCaW9nxgxqPZxZuRulMwjiv9W+QgzBj7DO9YI3cih8tcxkq5n0YPuf7vyBtIw2JbCq94m6e2T5GcHfx9Izv+JZI3uzGnhDJrDfNSTrF+ht7E69rjpyBZV+KTOHw4PA9giERoAP08n0nMVc/wxGQj3Q/AjC/2EMiXKxLGmAREc/aWAI5WV9DgS0PpqwZz2CQfKUNXLlOPa7V6iN9VlcOr9eoXy6Gzkc76o6Rd65nBs6036ML/rXuWJLFqf/OgpB73Bt/C76q2dRnuMpu7pO5zEdjKyPNItbclv3X2TJBcCI6KNHgtR1eqZJcCI3vXPKFg/AlgfAYhUQhgb9YR1uHe2ESd2PBe9mO0jo/cB38gpjIze+ygvLBMUn3Zb15SmyaurtBkwobrdd5WradtLgM7U6jjqTCcvhKpMlzMhipoC2MMYZ9Nbho6uiCNekJ/CQ5Ky4gG7bpYzYp703OSlIFgjpeVStC3cqUa7nHLFMFUCOACLhiy65UsULUelMvFaUPw8lm6vwwkQsJjdWvMkMtFFuq2DBcWPf1sjpeU6tuhoruEu56BiGCUBOR+OBLBozeYHaz6sibCwbWBOJB5geI4dl9Xp+6VOV8TVKi/6v4gfHCROiC80OlB6s9huaZfyO+iUJKhbEgQKvK2Lf664qqBEH8Uyd4U+q/3Jju0LXyo4a3ipCI3Fs6XeSteZbKh3qqE8vvGFY+oaYOx9trNq8DKg3PXZwioguptP0cbwuVaZQackYTE63JOQ2SHg4NSqgJGIIPQZQUYbxUgmU71HZKNrMy0YST9bsS31lllnyCyV3vrkT8SC43J0ZHusOb327eqrE6ebLdjUPs6vzwMSOLumJKXYfHeyMaBQXXb2YoCNs1lEbEmdIkY20dA6mRQOC+vdobLOQhskztVbyTpDjKn01qcS4Sdq4Z+m8kkNZgmcsu6NrnrwmFOyAaunWE3uksBgpOQr9G2jduFLnja3cyVnZBghxscUsyTYTG870Vmvrrxj6s0fkXHHl7EBa0SVWsQkcDZSUkEipclRkhFI8LtZ55rAO/IhmisJttTbHnSWU1cuUpJzq5HSDCm3IiWNcAiT3QjWcpgPW2eLj+UrazVodS7xyByfhUimS6W8Re12TcTb1FcRkgixQL3tQWdJZN+DpNX+njI1uUcCjVOSYSNIdrp+MA+vNFvH2SRSopOLHyIQ5xT1zPoUZzdwnraFmvpzE0I9GWwRHhnC32pQmt52pLMsOlL7sU3s1vpOFjpHRuqckjckHBIvAYYhJ/f8RxEipy2ATg78IR5xSO4j47pnG/6x+CrBIbH4yjF7ZMTLvas5pUL1tgudSU+5gEgeWM0O/iF3nL/OKak57DLg3X+PmoaBnRWuOgf1zRhDFOBc8iY+axN0dqZpzRkp7ySLcEhBY5Gd/fO+IDnr7Wy9VZ1l04utp4WDezZiR0Ts/sWSDPTkXwT5UZh0tkRnTd+Eg1Fj9sf5VZ9Tt+xgvIUPXXy/pciDZuIVp0zIPlteqjsJRDO53i5JZ5OEPbKw5MeZPmYcRdrpyGZsWuxrCZFohivevnSRz1yHtLRVni/m50zTmL4RsXF+hOnkLCe5lKeh+p5HRsfm4N1Q+aV5JeptbzpbqoO4vtpvZ/JKi+JjVot+vhFwRzQSAusjzlFJ2ByCxEGsCQ8gFjgg1pEAHCjvCjF1KQ3g6ZF4czyuxFxJetujzlKq6ZmQTT39nZL+IXDdUKe/VkuIPBqno07F2gi7XjYv5k1qu1fyKsB0rYmIPJ8QXpuP0Y0Vj4yQyBPjzAqF6m13OkulJPUZlioYjMIXfFOhvyg8Fikx5bgjwT7h0v07XRgYEQnPW2xvomC2/UNA4Y89T1uenQp5iu8fKgE5wn9uKE1ve9VZCj0xELE7XOygmaKRa+BwC925CKljLlrozsVXbrxqt3uBWAba7BjmppkK/6XqbIn8vMxeCse/fQS7BN3F183tlFy0cGmKkpEyFcYx/ai2xxsIRRudeL9Inc1VipfXX6rPv8tulhvm4qv1NDuTIKscMkjAj56cryrxG1AZWnyZKKVnDh3zra/dRcWlaszWlErlb7d8yUg5wsBRBoy2wgElIN2yIcT7jtUhJdRvdUoJhRmjMmOV8bJZUOFAEpBOOcfHGTq3IXSgpm3elOqUMqtARssO2ffeiDNTq+jXkIB0ybob0dGPa9C7NBr/A6KTYmS/JtiKAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$\\displaystyle \\frac{k vd}{d + k} + \\left(- \\frac{k vd}{d + k} + v_{0}\\right) e^{- \\frac{t \\left(d + k\\right)}{m}}$"
      ],
      "text/plain": [
       "                        -t⋅(d + k) \n",
       "                        ───────────\n",
       " k⋅vd   ⎛   k⋅vd     ⎞       m     \n",
       "───── + ⎜- ───── + v₀⎟⋅ℯ           \n",
       "d + k   ⎝  d + k     ⎠             "
      ]
     },
     "execution_count": 5,
     "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": "markdown",
   "metadata": {},
   "source": [
    "The following represents the explicit euler method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "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)*h\n",
    "            xs.append(x)\n",
    "        return xs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following appplies the previous numerical solver to the cruise control example. \n",
    "The results are plotting against the analytical solution `v_t` obtained earlier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl4VPXZ//H3nclKIIR9S0ICouybIYqIdUGq1orWpbhbrWi1j9rWVm37a221fWytba1Pqw8ulbohohbriru4sYOgURBBEraw79nv3x8ZfCJkGUImJ8l8XteVK5mZc+Z8EuXcc77nu5i7IyIisSsu6AAiIhIsFQIRkRinQiAiEuNUCEREYpwKgYhIjFMhEBGJcSoEIiIxToVARCTGqRCIiMS4+KADRKJz586enZ0ddAwRkRZl/vz5m9y9S33btYhCkJ2dzbx584KOISLSopjZl5Fsp6YhEZEYp0IgIhLjVAhERGJc1O8RmFkImAescffTzSwHmAp0BBYAF7t76cG+b1lZGYWFhRQXFzdu4BiVnJxMRkYGCQkJQUcRkSbWFDeLrwfygbTw4z8Af3H3qWZ2H3AFcO/BvmlhYSHt2rUjOzsbM2u8tDHI3dm8eTOFhYXk5OQEHUdEmlhUm4bMLAP4FvBA+LEBJwLTw5tMAc5syHsXFxfTqVMnFYFGYGZ06tRJV1ciMSra9wj+CvwMqAw/7gRsc/fy8ONCoFdNO5rZJDObZ2bzNm7cWOObqwg0Hv0tRZqhgjkw666q71EUtUJgZqcDRe4+v/rTNWxa41qZ7j7Z3XPdPbdLl3rHQ4iINI1ITs6NsU3BHJhyBrzxu6rvUSwG0bxHMAY4w8xOA5KpukfwVyDdzOLDVwUZwNooZhARqVIwB1bNguyxkJnXsG32nZwrSiGUCJc+d8B2vno2/GsCVJTioQS+/NZUtnQcxu6SCvaUlrO7pII2RQsYN/f7hCrLqIhL4OHD/sayhAHsKauguLSCvWUVnLrtcS4oLyFEZdXxVs2qPfchilohcPdbgFsAzOx44EZ3v9DMngLOoarn0KXAjGhlEJFWIKAT+LJvPsb6tKHsLC5nZ3E5O4rLGLjiGcaET84V5SVMf+pxHk8qY1dxGbtLKthdUs4lFU/zo1AJ8VZJRZnz1PQn+EfFnq8d65rQ85wcX0acVVJZUcbe5W/zbmIXUhJDpCRUfX2ROoKKPVOJoxwLJVb9blESxBQTNwFTzex2YCHwYAAZDtlNN91E7969ueaaawC49dZbadeuHT/5yU8CTibSgjTCCXz/bfySGezsOpKtu0vZuqeMrbtL6bL4BQaWlxAXPoG/9Nw0preB7XvL2L63jB17y5hY8hQ3xP3fCfy5f0/jHxUlXzvUkXHdGJUYTwLllBPPwtBg0pLj6ZWeTGpiPKlJ8fQqOxn/+N9UVpZDKIHjTjyTo3qMIjUxRJvEeFKTQrTblEboqefwilLi4xO57tLvcd0Bv/8xUDCo/iLYCJqkELj7W8Bb4Z+/ABr1N/rNfz7mk7U7GvMtGdgzjV9/e1Ctr0+cOJEbbrjhq0Iwbdo0Xn755UbNINKiNcZJftUsvKIU8wq8opQ1C2eydHsmm3aVsnlXKZt2lZBb8CSnhz+ll5eX8NfJD/I/5Zu+9jYjLZ3HwifwMovnhR192VxZSvuUBHq2TyEtJYGu5Sfjn1adwC0+gZPGn81JGaNol5xAWnIC7ZLjaZN4GlY4BlbNIj57LHfUeHIeBHm9YdUs4rLHcnRN23QaW/X71neSz8yLagHYp0VMOtccjRgxgqKiItauXcvGjRvp0KEDWVlZQccSaRqN0RSzchaET/KVFaUsffd53u7agaKdJRTtLKZoZwldtqVwd2WIBJwyD3HdB6kseH/BV++R3iaBbcn9OIUEoIxKS6D7sJP5ZfcBpLdJpGNqAultEunQ5njKt+aRtO5DUnLGcm+NJ9chUND7q9/ryEM5OTfWNk2kVRSCuj65R9M555zD9OnTWb9+PRMnTgwkg0ija+yTfHkpi975DzM7prF++17Wbi9m/fZiuu+IZ0ro/07yt37UgQW+jPYpCXRtl0TXtCRS+x7DU/ydAaWL2dNzNDdkHkWntol0bptEx9REEkLhjo8FI2HVLELZY7motpNr5zHQb0zdv3szOjk3pVZRCIIyceJErrzySjZt2sTbb78ddByR+kWhuWb57Bd5e2VnCrfuoWDrXgq37qHT1gQesvBJnhC3L+3Ikrgv6JaWTM/2KQzPTKdH+njeqOxNv72LqOx9LHf3PYYu7ZJITgjtF2o4cF7dv1eMnsAbiwrBIRg0aBA7d+6kV69e9OjRI+g4InWL8CRPRSmET/LFy94i3/uxatPuqq/Ne0han8ZvPUSCV32Sv3l+exZ4Pu2S4sno2IbenVLJOOwkXrMsjihehGWP5b4jxtA5NYm4uP2HEg2ggZMLSCNSIThES5YsCTqCSJX6Pu1XO8lX75deXlHJ6i17WF60i52bcziDeOLCzTUXvpbAglffByDOoGd6Cjmdh/JI13sYXrGEyt5j+O1hY8js0Ia0lPj9RqgPQif5lkGFQKQlaIQmnfKsMcSFEqAcKiyeu5d15dV577By025KK/bNApPEK+1uZVzKMrZ1O4rTso7imk6pZHdOJbNjCknx+5ptjgJ0X6y1UCEQae4OskmHilKKP3+bj8oPI3/dDvLX7eCTdTv4bP1OBlXczNFx+XzoA9i8JYPDuqRwfP8uHNalLf26taNvl1TaJWsq8lijQiAStAY26eyzbU8pK+MGM9jiiQu3218wM54F/gEAHdokMLBnGhcf3ZsBPYZwRPd2XNe1bQ03ZSVWqRCIBCmST/vZY6teqyjFQ4l8ljyM995dyUeF21hcsI1Vm6umLxhpt3BquxVs73YUJ/UZzX/1SGNAjzS6pSVpdlmpkwqBSJDq+bS/fnsxc7ZksD7nryQWvscLO/oy9+kS4BO6pyUzLLM9543KZHhGOoMzxpOmZh1pABUCkWiqr9lnv0/7a9OP5N25q5m9cgtzV22hYMteANokpjMi60LyhqdzZUY6wzLT6ZaW3MS/jLRWKgQi0VJPs4+7syplEPlH/i+ly99hxvYc3nxsN7CEjqmJjMruwKWjs8nL6cjAHmnEh6K9jpTEKhWCZmbVqlUMGDCAI4444qvnfvzjH3PJJZfUus+tt95K27ZtufHGG5siouzTgJu8WzsO570Vm3h3+SZmLd/Emm17gXh6pX+Low7vyH/ndGRUdkf6dklVu740mdgqBJHMa94M9O3bl0WLFkXt/cvLy4mPj63/9I0uwpu8HkrEK0qpsHh+Pi+N6S++iju0S47nmL6duPr4vow9rDO9O7XRiV8CEztng0j+4R6Eg1mP4M4772TatGmUlJRw1lln8Zvf/Ia5c+dyxRVXMGfOHCoqKsjLy+PJJ5+kbdu2tR6zbdu27Nq1C4Dp06fz/PPP8/DDD39tmxUrVnDttdeyceNG2rRpw/3330///v257LLL6NixIwsXLmTkyJHcddddDf7dhTpv8m7fW8ZbnxXx6ifxbCv7BUPLlzCXgdB1EDec1IVj+3VmWEZ7NfVIsxG1QmBmycA7QFL4ONPd/ddm9jDwDWB7eNPL3D16H3/3qad3xsGKdD2CmTNnsnz5cubMmYO7c8YZZ/DOO+9w3HHHccYZZ/DLX/6SvXv3ctFFFzF48GBWrVrFihUrGD58+Ffvcc899zB2bGSrE02aNIn77ruPfv36MXv2bK655hreeOMNAJYtW8Zrr71GKKT+4/U6iJu8hBIp6jSKF99byav5G5j9xRbKK53ObRM5acg3GD7gXH7Qt5MGakmzFc0rghLgRHffZWYJwLtm9lL4tZ+6+/QoHvtA+/3DPdRl3yJdj2DmzJnMnDmTESNGALBr1y6WL1/Occcdx69+9StGjRpFcnIyf/vb377ap6FNQ7t27eL999/n3HPP/eq5kpL/W2Hp3HPPVRGIRCRXj5l5fHn6VFYvfIVntuTw7L92Ap/Qt0sq3x/bh5MHdmN4ZjqhAyZZE2l+orlmsQO7wg8Twl8erePVKzMvshWBDkIk6xG4O7fccgtXXXXVAa9t2bKFXbt2UVZWRnFxMampqXUer3obcnFx8QGvV1ZWkp6eXmsRqe/9JayOq8eCLXt4bvFa/rN4LZ+u302cHUtu7478fFRXxg3oRp8utTftiTRXUW2kNLOQmS0CioBX3X12+KXfmdlHZvYXM0uKZoavycyDsT9ptBvFEydOZOrUqUyfPp1zzjmnxm2++c1v8tBDD33Vtr9mzRqKioqAqmac2267jQsvvJCbbrqp3uN169aN/Px8KisrefbZZw94PS0tjZycHJ566imgqggtXry4ob9e61UwB2bdVfW9JvuuHi0EoUS2dj2Kh99byVn/eI+xf3yTO1/5jNSkeH47YRBzfjGOaVePZtJxfVUEpMWK6s1id68AhptZOvCsmQ0GbgHWA4nAZKoWs//t/vua2SRgEtBsl4CMZD2C8ePHk5+fz+jRo4GqG76PPvooL7/8MvHx8VxwwQVUVFRwzDHH8MYbb9CnT58D7hFcfvnlXHfdddxxxx2cfvrpZGZmMnjw4K+KS3WPPfYYP/jBD7j99tspKytj4sSJDBs2LDp/gJYowmafPec/y+dzX+KZzTn86+HtVPp2+ndvx02n9Ofbw3qQ0aFNMPlFosCqWnCa4EBmvwZ2u/ufqj13PHCju59e1765ubk+b968rz2Xn5/PgAEDohE1ZsXE33TWXfDG76qafSwEJ/6i6iqR8BVU4XamzlnNc4vXsqe0gqyObZgwvCdnDOtJv27tAg4vcnDMbL6759a3XTR7DXUBytx9m5mlAOOAP5hZD3dfZ1UN3mcCS6OVQeQANXQa2L63jBmL1vD47NV8un4nKQkhvj2sB98dlcXIrHT175dWL5pNQz2AKWYWoupexDR3f97M3ggXCQMWAVdHMUOTWbJkCRdffPHXnktKSmL27Nm17CFRUV+3z3CnAV85i0+Th3H/B4m8uOQ1issqGdwrjdvPHMyE4T3V1VNiSjR7DX0EjKjh+RMb8RjN5tPakCFDojoaONqaqokwqiJo/99bWsHTa7oxZV4uy4t20TZpA98ZmcH5o7IYktE+oOAiwWqxI4uTk5PZvHkznTp1ajbFoKVydzZv3kxycgufzbKObp8bd5bwyAereOTDL9m6p4whvdrzx7OH8q2hPUhNarH/DEQaRYv9F5CRkUFhYSEbN24MOkqrkJycTEZGRtAxDk0N7f/LN+zkwXdX8szCNZSWVzJuQFeuHNuHvJyO+gAhEtZiC0FCQgI5OTlBx5CmVtc9gGrt/0sShvCX1+DNz94hKT6Oc47M4Ipjc+irvv4iB2ixhUBiUATz+7+yPYt7FuXx8doddErdzo/GHc5FR2fRqW3TjVsUaWlUCKTlqOUegLvzxqdF/PnVZXy8dgc5nVO54ztDOHNELy3QLhIBFQJpOfa7B+C9j+WdZRv586vLWFywjayObfjTucM4c3hPTfEschBUCKTlqDZx4EfxQ/jtCxXM+3IOvdJTuOM7Qzj7yAwSVABEDpoKgTQfEawgN7fiMO7Kr+DDL7bQLW0Pt00YxHmjMkmKVxOQSEOpEEjzUM+N4IIte7j9hU945eMNdG6bxK9OH8gFR2XpHoBII1AhkOahlhvBu0vK+cdbn3P/rJWEzPjJyYfz/bF9SElUARBpLCoE0jzUcCP43wsLueOlT9mwo4Qzh/fk5lMH0L19Cx/9LNIMqRBI81DtRvDylOHc9Hw5C1YvZmhGe/5x4UiO7N0x6IQirZYKgTQbRelD+eP6JKbPL6RzW+eP5wzlnJEZxGndX5GoUiGQplFHj6DKSufxOau546VPKSmv4Kpv9OGHJxymqaBFmogKgURfHT2CVm/ew01Pf8QHX2xmzGGduP3MIeR0Tg04sEhsUSGQ6KuhR1Blr1E8/P4q7nzlM+LjjDu+M4TvjsrUjKAiAYjmUpXJwDtAUvg4093912aWA0wFOgILgIvdvTRaOaQZ2K9HUGH6kdzwvx8w78utHH9EF35/1hB6pqcEnVIkZkXziqAEONHdd5lZAvCumb0E/Bj4i7tPNbP7gCuAe6OYQ4IW7hFU8cU7/Gd7H256ci/JCaX8+bxhnDWil64CRAIWzaUqHdgVfpgQ/nLgROCC8PNTgFtRIWj1liUO4KdLSllcuJ3xA7tw+5mD6ZqmMQEizUFU7xGEF66fDxwG/B1YAWxz9/LwJoVAr2hmkGC5O//64Et+90I+bZPjuef8EZw+tIeuAkSakagWAnevAIabWTrwLDCgps1q2tfMJgGTALKysqKWURpBLV1Dt+8t4+anP+Klpes5sX9X7jxnqBaIEWmGmqTXkLtvM7O3gKOBdDOLD18VZABra9lnMjAZIDc3t8ZiIc1ALV1DFxds44dPLGDdtmJ+flp/vn9sHw0ME2mmojZ5u5l1CV8JYGYpwDggH3gTOCe82aXAjGhlkCawX9dQXzmLh95dyTn3vU9lJTx51WgmHddXRUCkGYvmFUEPYEr4PkEcMM3dnzezT4CpZnY7sBB4MIoZJNqqdQ31UAJ/+LQz933xCeMGdONP5w4lvU1i0AlFpB7R7DX0ETCihue/AGpedURannDX0LWLZnLrko68saoLv/xWf644Nkc3hEVaCI0slkPi7jy4qjN3fDCSbmnJPHX1CEZkdQg6logcBBUCabDS8kpueWYJTy8oZPzAbtx5zjDat9FEcSItjQqB1K2WrqFbd5dy1aPzmbNyCzeM68f1J/VTU5BIC6VCILWrpWvoio27uOLhuazdXszdE4czYbjGBIq0ZCoEUrsaZg19v7QPVz8yn4RQHE9ceZRWDhNpBVQIpHb7zRr62t5+XP3gHHI6p/LQZaPI7Ngm6IQi0ghUCKR24a6hlStn8ei6TH71Roix/Trx9wtHkqbVw0RaDRUCqdOebiO54c04Zn6ygYuOzuLWbw8iPhS1AekiEgAVAqnV5l0lXPbPuXy8dju/On0g3xuTrZ5BIq2QCoHUqGhHMRc+MJvVW/Zw/yW5nDSgW9CRRCRKVAjkAGu37eXCB2azYUcx//zeKI7p2znoSCISRSoEsayGwWKrN+/hggc+ZPueMh65Ik/dQ0VigApBrKphsNiK5IFceP9sissrePzKoxmS0T7olCLSBFQIYtV+g8U2Lnmd7y7YgbvzxJVHM6BHWtAJRaSJqBDEqmqDxSrjErhxbltCCfDY90dzWNe2QacTkSakQhCrwoPF1iyayU3z27MyeRDTrjyK3p1Sg04mIk0smktVZprZm2aWb2Yfm9n14edvNbM1ZrYo/HVatDJI3WaX9WX83FwKUgcz7erRKgIiMSqaVwTlwE/cfYGZtQPmm9mr4df+4u5/iuKxpR4LVm/lsn/OpVeHFB77/lF0S0sOOpKIBKTeQmBmXYExQE9gL7AUmOfulXXt5+7rgHXhn3eaWT6g+YqbgeUbdnL5w3PpmpbE41ceRdd2KgIisazWpiEzO8HMXgFeAE6lajH6gcAvgSVm9hszi6hriZllU7V+8ezwUz80s4/M7CEz07qGTWjNtr1c8tAcEkJxPHK5ioCI1H1FcBpwpbuv3v8FM4sHTgdOBp6u6wBm1ja8zQ3uvsPM7gVuAzz8/S7g8hr2mwRMAsjKyorol5FqahgstmV3KZc8OJtdJeVMu2o0WZ00jbSIgLl79N7cLAF4HnjF3f9cw+vZwPPuPriu98nNzfV58+ZFJWOrVMNgsd1dR3LBA7P5dN0OHrniKPJyNGJYpLUzs/nunlvfdvX2GjKz680szao8aGYLzGx8BPsZ8CCQX70ImFmPapudRdU9B2lM+w0WK//iHa5+dD5L12znfy4YqSIgIl8TSffRy919BzAe6AJ8D7gjgv3GABcDJ+7XVfSPZrbEzD4CTgB+1MDsUpt9g8UshIcSufvzbsxavon//s4QTh6oWURF5Osi6T66bwL604B/uvtii2BSend/t9q+1b14EPmkIcKDxXzlLB4o7MU9H7Xj5lP7c15uZtDJRKQZiqQQzDezmUAOcEt4TECdXUelGcjM455lHfjzR8uYdFwfrv5G36ATiUgzVWshMLN4dy8HrgCGA1+4+x4z60RV85A0Y49++CV/fnUZZ4/M4OZT+gcdR0SasbquCD40s0LgZeBld98G4O6bgc1NEU4a5t3lm/jVjKWc2L8rd5w9hLg4LS8pIrWrtRC4e66Z9aZqMNlfzawX8C7wEvC2u5c0UUY5CKs37+GHTyygX9d23HP+CBK00LyI1KPOs4S7f+nu97n7mcAxwH+AccAsM3uhKQJK5HaXlDPpkXlUVjqTLzmS1CRNLisi9Yv4TOHuZcAb4S/CVwjSTLg7P52+mGUbdvLw9/I0k6iIRCySAWWnm9lCM9tqZjvMbKeZ7XD3NU0RUCLzj7dW8OKS9dx8an+OO7xL0HFEpAWJ5Irgr8B3gCUezfkoJHL7zSP05qdF/GnmZ5wxrCdXju0TdDoRaWEiKQQFwFIVgWZiv3mE1kyYynVPlzCgexp/OHsoEYz1ExH5mkgKwc+AF83sbeCrnkI1TSInTaDaPEJeUcrLz08nITSByZccSUpiKOh0ItICRVIIfgfsApKBxOjGkXqF5xHyilLKiOelXYfx98tHktFBU0qLSMNEUgg6unu9s41KEwnPI/Tea8/y52VdOeP0CYzu2ynoVCLSgkUy2ui1SKadlqbzyo4sLvrsWPqMPJFLj8kOOo6ItHCRFIJrgZfNbG/17qPRDiY1K9iyhxunLWZYRntuP3Owbg6LyCGrt2nI3ds1RRCpX3lFJT96chEA/3PBSJITdHNYRA5dXYvXZ9e1Y3jFsozGDiS1u/etFcz7ciu3nTmYzI66OSwijaOupqE7zexpM7vEzAaZWVczyzKzE83sNuA9YEBtO5tZppm9aWb5ZvaxmV0ffr6jmb1qZsvD3zs08u/UKi0q2MZfX1/OhOE9OXOEZvcQkcZTayFw93OB/wccAfwdmAXMAL4PfAac6O6v1vHe5cBP3H0AcDRwrZkNBG4GXnf3fsDr4cdSh90l5Vw/dSHd05L57YTBQccRkVamznsE7v4J8IuGvLG7rwPWhX/eaWb5QC9gAnB8eLMpwFvATQ05Rqz4zX8+pmDLHqZOGk37lISg44hIK9Mkk9WH7zeMAGYD3cJFYl+x6NoUGVqql5asY9q8Qn5wfF/ycjoGHUdEWqGoT1hvZm2Bp4Eb3H1HpN0dzWwSMAkgKysregGbm2oTyq1PG8rNzyxhaEZ7bhh3eNDJRKSVimohMLMEqorAY+7+TPjpDWbWw93XmVkPoKimfd19MjAZIDc3NzYmvKs2oZyHErm3w+8pLe/J3RO10piIRE9di9ePrGtHd19Q1+tW9dH/QSB/vwnqngMuBe4If58RcdrWrvqEcuWlpK79kF9P+CU5nbXIjIhET11XBHeFvycDucBiwIChVLX1H1vPe48BLgaWmNmi8HM/p6oATDOzK4DVwLkNi94KVZtQrqQyhGeP4bujMoNOJSKtXF2L158AYGZTgUnuviT8eDBwY31v7O7vUlU4anLSwUeNAZl5lFzwLI8++Rjvlw/gTxecrykkRCTqIrlH0H9fEQBw96VmNjyKmWLa75e0Y8r2U3nkijw6pGrWbxGJvkgKQb6ZPQA8CjhwEZAf1VQxav6XW5nywZdcdkw2Y/tp3WERaRqRFILvAT8Arg8/fge4N2qJYlRZRSU/f2YJPdsn89NvHhF0HBGJIZHMPlpsZvcBL7r7Z02QKSY9MGsln23YyQOX5JKaFPXhHSIiX6m3c7qZnQEsAl4OPx5uZs9FO1gsWb15D3e/voxTBnVn3MBuQccRkRgTySilXwN5wDYAd18EZEcxU0xxd345YynxcXHcesagoOOISAyKpBCUu/v2qCeJUc9/tI53lm3kxvGH0719ctBxRCQGRdIYvdTMLgBCZtYPuA54P7qxYsP2vWX85j+fMDSjPRePzg46jojEqEiuCP4LGASUAI8D24EbohkqVvzx5U/ZsruE3581hFCcBo6JSDAi6TW0B/iFmf3e3Xc3QabWq9rMovMr+/HY7NV8/9gcBvdqH3QyEYlh9RYCMzsGeABoC2SZ2TDgKne/JtrhWpX9ZhZ9JOm39Gzflx+drOmlRSRYkTQN/QX4JrAZwN0XA8dFM1SrtN/Moj23zeO3EwZrzICIBC6iSe7dvWC/pyqikKV12zezqIUo8RCefazGDIhIsxDJx9GCcPOQm1kiVb2GNNfQwcrMwy+ZwVNPT2XGtj7cdd55QScSEQEiuyK4GriWqoXn1wDDw4/lID2/NZOfbRjHyeNP15gBEWk2Iuk1tAm4sAmytGp7Syu4/YVPGNJLYwZEpHmJZK6hPmb2HzPbaGZFZjbDzPpEsN9D4e2XVnvuVjNbY2aLwl+nHeov0FI8MOsLNuwo4dffHqgxAyLSrETSNPQ4MA3oAfQEngKeiGC/h4FTanj+L+4+PPz1YqRBW7KNO0u47+0VfHNQN3KzOwYdR0TkayIpBObuj7h7efhr3wI1dXL3d4Ath5ywFbj79WUUl1fys1P6Bx1FROQAkRSCN83sZjPLNrPeZvYz4AUz62hmDfl4+0Mz+yjcdNShAfu3KJ8X7eKJOQVckJdF3y5tg44jInKASLqPfjf8/ar9nr+cqiuDeu8XVHMvcFt4v9uAu8LvcwAzmwRMAsjKyjqIQzQvf3j5U1ISQlw/rl/QUUREahRJr6GcxjqYu2/Y97OZ3Q88X8e2k4HJALm5ufU2RTVHc1Zu4dVPNnDj+MPp3DYp6DgiIjWKpNfQuWbWLvzzL83sGTMb0ZCDmVmPag/PApbWtm1L5+78/sV8uqclc8WxB3PRJCLStCK5R/D/3H2nmR1L1ZxDU4D76tvJzJ4APgCOMLNCM7sC+KOZLTGzj4ATgB8dQvZm7YUl61hUsI0fjz+clMRQ0HFERGoVyT2CffMKfQu4191nmNmt9e3k7ufX8PSDB5GtxSopr+CPL39G/+7tOHtkRtBxRETqFMkVwRoz+1/gPOBFM0uKcL+Y9eiHq1m9ZQ83n9pfg8em9m6nAAANAUlEQVREpNmL5IrgPKoGhv3J3beF2/l/Gt1YLVB40ZldPUZzzxu7OPawznzj8C5BpxIRqVekK5Q9U+3xOmBdNEO1ONUWnUmyePoU38Itp12Oma4GRKT506oojaHaojNW6VyRsYZBPbX8pIi0DGrrbwzhRWcqiKOMePKOnxB0IhGRiKkQNIbMPFac+jh3lZ3L9MH/oMvAsUEnEhGJmJqGGsmti1JZmnQOb3/7hKCjiIgcFF0RNIKFq7cya/kmrvpGX9KSE4KOIyJyUFQIGsHf3/yc9DYJXHR076CjiIgcNBWCQ/Tx2u28ll/E5WNyaJukljYRaXlUCA7R39/8nHZJ8Vx6THbQUUREGkSF4BAs37CTl5au55JjetM+RfcGRKRlUiE4BP94awXJ8SEuH9NoSzaIiDQ5FYIG+nLzbmYsWsOFR2XRSYvOiEgLpkLQQPe+tYL4UByTjtOiMyLSsqkQNMCabXt5ekEhE0dl0jUtOeg4IiKHJGqFwMweMrMiM1ta7bmOZvaqmS0Pf+8QreNH0+S3V+AOV32jb9BRREQOWTSvCB6mah2D6m4GXnf3fsDr4cctStHOYp6YW8DZIzPolZ4SdBwRkUMWtULg7u8AW/Z7egJVax4T/n5mtI4fLQ/MWkl5RSU/OF5XAyLSOjT1PYJu4YVt9i1w07WJj39Ituwu5dEPv+SMYT3J7pwadBwRkUbRbOdEMLNJwCSArKysYMOEl6F8uSiLPaVtuPaEw4LNIyLSiJq6EGwwsx7uvi689nFRbRu6+2RgMkBubq43VcADhJeh9IpSzqoMUdj3Lvp1axdYHBGRxtbUTUPPAZeGf74UmNHExz944WUozStIoJxLehYEnUhEpFFFs/voE8AHwBFmVmhmVwB3ACeb2XLg5PDj5i17LB5KpJw4KuIS6D705KATiYg0qqg1Dbn7+bW8dFK0jhkVmXk8N+xePvvwJSacdR5HZOYFnUhEpFFpZHE9yisq+cPSNBb2vpwjcscFHUdEpNGpENTj1U82sHZ7Md8bkx10FBGRqFAhqMfD768io0MKJw3oFnQUEZGoUCGoQ/66HcxeuYVLRvcmFGdBxxERiQoVgjpMeX8VyQlxnJebGXQUEZGoUSGoxdbdpfx70RrOGtGL9DaJQccREYkaFYJaPDmvgOKySi1KLyKtngpBDSoqnUc++JKj+3Skf/e0oOOIiESVCkENXsvfwJpte7lMVwMiEgNUCGow5f1V9EpPYZy6jIpIDFAh2M9n63fy/orNXHR0b+JD+vOISOunM91+pnywiqT4OCaOUpdREYkNKgTVbN9TxrML1jBheE86pKrLqIjEBhWCaqbNK2BvWYW6jIpITFEhCKuodP714SrysjsyqGf7oOOIiDQZFYKwNz8tomDLXl0NiEjMCWTxejNbBewEKoByd88NIkd1Uz5YRY/2yYwfpC6jIhJbAikEYSe4+6YAj/+Vz4t2Mmv5Jn76zSNIUJdREYkxQRaC5qFgDitfnE5efC8mjtIKZCISe4IqBA7MNDMH/tfdJweSomAOPuXbnFBWyjcSEkjcOhraak1iEYktQbWDjHH3kcCpwLVmdtz+G5jZJDObZ2bzNm7cGJ0Uq2bh5aXEWyUJlMOqWdE5johIMxZIIXD3teHvRcCzwAEfw919srvnuntuly5dopOj97GUEk85cVgoEbLHRuU4IiLNWZMXAjNLNbN2+34GxgNLmzoHwJK4I7ig5OcsPfy/4NLnIFPNQiISe4K4R9ANeNbM9h3/cXd/OYAcTJ9fyMeh/vT5zjhITggigohI4Jq8ELj7F8Cwpj7u/orLKpixaC2nDO5OmoqAiMSwmO00/3p+Edv3lnHOkRlBRxERCVTMFoLp8wvo0T6ZY/p2DjqKiEigYrIQbNhRzNvLNnL2yAxCcRZ0HBGRQMVkIXh24RoqHc5Ws5CISOwVAndn+vxCcnt3IKdzatBxREQCF3OFYHHhdj4v2qWbxCIiYTFXCJ6aV0ByQhzfGtoj6CgiIs1CTBWC4rIKnlu8llMH96Cdxg6IiAAxVghe/WQDO4vL1SwkIlJNTBWCp+YX0is9hdF9OgUdRUSk2YiZQrB+ezHvLt/I2SN7EaexAyIiX4mZQvDMwkKNHRARqUFMFIJ9YwfysjvSu5PGDoiIVBcThWDB6m18sXE35+TqakBEZH8xUQimzy8kJSHEaUM0dkBEZH+tvhAUl1Xw/OK1nDqkO22TgliHR0SkeQukEJjZKWb2mZl9bmY3R/NYr3y8np0lGjsgIlKbINYsDgF/B04FBgLnm9nAaB1v+vxCMjqkcHSOxg6IiNQkiCuCPOBzd//C3UuBqcCEaBxo4yezGLLyQX542BaNHRARqUUQjea9gIJqjwuBoxr9KAVzSJ9+Nj8OlRGXPwMKsiAzr9EPIyLS0gVxRVDTR3M/YCOzSWY2z8zmbdy48eCPsmoWocoy4q2SuIoyWDWrAVFFRFq/IApBIZBZ7XEGsHb/jdx9srvnuntuly5dDv4o2WOJi08CC0EoEbLHNjiwiEhrFkTT0Fygn5nlAGuAicAFjX6UzDy49LmqK4HssWoWEhGpRZMXAncvN7MfAq8AIeAhd/84KgfLzFMBEBGpRyAjrNz9ReDFII4tIiJf1+pHFouISN1UCEREYpwKgYhIjFMhEBGJcSoEIiIxztwPGNTb7JjZRuDLBu7eGdjUiHGagjJHX0vLC8rcVFpa5rry9nb3ekfktohCcCjMbJ675wad42Aoc/S1tLygzE2lpWVujLxqGhIRiXEqBCIiMS4WCsHkoAM0gDJHX0vLC8rcVFpa5kPO2+rvEYiISN1i4YpARETq0KoLgZmdYmafmdnnZnZz0HnqY2YPmVmRmS0NOkskzCzTzN40s3wz+9jMrg86U33MLNnM5pjZ4nDm3wSdKVJmFjKzhWb2fNBZImFmq8xsiZktMrN5Qeepj5mlm9l0M/s0/P/06KAz1cXMjgj/bfd97TCzGxr0Xq21acjMQsAy4GSqFsOZC5zv7p8EGqwOZnYcsAv4l7sPDjpPfcysB9DD3ReYWTtgPnBmM/8bG5Dq7rvMLAF4F7je3T8MOFq9zOzHQC6Q5u6nB52nPma2Csh19xbRJ9/MpgCz3P0BM0sE2rj7tqBzRSJ8vlsDHOXuBz3mqjVfEeQBn7v7F+5eCkwFJgScqU7u/g6wJegckXL3de6+IPzzTiCfqjWpmy2vsiv8MCH81ew/DZlZBvAt4IGgs7RGZpYGHAc8CODupS2lCISdBKxoSBGA1l0IegEF1R4X0sxPUi2ZmWUDI4DZwSapX7iJZRFQBLzq7s0+M/BX4GdAZdBBDoIDM81svplNCjpMPfoAG4F/hpvfHjCz1KBDHYSJwBMN3bk1FwKr4blm/8mvJTKztsDTwA3uviPoPPVx9wp3H07Vetl5Ztasm+HM7HSgyN3nB53lII1x95HAqcC14abP5ioeGAnc6+4jgN1As7+vCBBuxjoDeKqh79GaC0EhkFntcQawNqAsrVa4nf1p4DF3fyboPAcjfOn/FnBKwFHqMwY4I9zmPhU40cweDTZS/dx9bfh7EfAsVc21zVUhUFjt6nA6VYWhJTgVWODuGxr6Bq25EMwF+plZTrhiTgSeCzhTqxK+8fogkO/ufw46TyTMrIuZpYd/TgHGAZ8Gm6pu7n6Lu2e4ezZV/x+/4e4XBRyrTmaWGu5AQLiJZTzQbHvDuft6oMDMjgg/dRLQbDs97Od8DqFZCAJas7gpuHu5mf0QeAUIAQ+5+8cBx6qTmT0BHA90NrNC4Nfu/mCwqeo0BrgYWBJucwf4eXhN6uaqBzAl3MsiDpjm7i2iO2YL0w14tuqzAvHA4+7+crCR6vVfwGPhD45fAN8LOE+9zKwNVT0jrzqk92mt3UdFRCQyrblpSEREIqBCICIS41QIRERinAqBiEiMUyEQEYlxKgQiIjFOhUCkBuEpia+p4/UUM3s7PB6htm1eM7MO0Uko0nhUCERqlg7UWgiAy4Fn3L2ijm0eqec9RJoFFQKRmt0B9A0v+HFnDa9fCMyAqnUZzOyd8LZLzWxseJvnqBr+L9KsaWSxSA3C02o/X9MCQeEpCFa7e/fw458Aye7+u3BTUZvw+gyY2XLgaHff3GThRQ5Sq51rSCSKOgPVFy2ZCzwUnon13+6+qNprRUBPQIVAmi01DYkcvL1A8r4H4ZXljqNqqcBHzOySatsmh7cXabZUCERqthNoV9ML7r4VCJlZMoCZ9aZq4Zj7qZqWe2T4eQO6A6uaIrBIQ6kQiNQg3Kb/Xvjmb003i2cCx4Z/Ph5YZGYLgbOBu8PPHwl86O7l0c4rcih0s1ikAcxsBPBjd7+4jm3uBp5z99ebLpnIwdMVgUgDuPtC4M26BpQBS1UEpCXQFYGISIzTFYGISIxTIRARiXEqBCIiMU6FQEQkxqkQiIjEuP8PZdb76ik9wp0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Instantiate the symbolic equations to the example.\n",
    "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",
    "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",
    "# Define simulation step size, time range, and total number of steps to be taken.\n",
    "h=0.2\n",
    "t_range = np.arange(0,7,h)\n",
    "num_steps = len(t_range)\n",
    "\n",
    "# Compute analytical trajectory\n",
    "xs_c = [v_t_example.subs(t, x) for x in t_range]\n",
    "\n",
    "# Compute approximate trajectories\n",
    "xs_ee = Explicit_Euler().simulate(lambda v_val: dv_dt_example.subs(v,v_val), \n",
    "                                      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.xlabel('t (s)')\n",
    "plt.ylabel('speed (m/s)')\n",
    "\n",
    "plt.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise: Implement the Midpoint Method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The midpoint method approximates the solution of a given IVP in the form of\n",
    "$$ \\dot{x} = f(x) $$\n",
    "by the following:\n",
    "$$ \\mathit{aux} = x(t) + f(x(t))\\frac{h}{2}$$\n",
    "$$ x(t+h) = x(t) + f(aux)h$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Midpoint:\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",
    "            # TODO: implement the following \n",
    "            x = ...\n",
    "            xs.append(x)\n",
    "        return xs\n",
    "\n",
    "# Run solver and plot\n",
    "xs_mid = Midpoint().simulate(lambda v_val: dv_dt_example.subs(v,v_val), \n",
    "                                v0_num, 0.0, h, num_steps)\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, xs_mid, 'o', markersize=3, label='v_mid')\n",
    "plt.xlabel('t (s)')\n",
    "plt.ylabel('speed (m/s)')\n",
    "\n",
    "plt.legend()\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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}