11. Differential_Equations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 11. Differential Equations\n",
    "\n",
    "[Data Science Playlist on YouTube](https://www.youtube.com/watch?v=HReAo38LoM4&list=PLLBUgWXdTBDg1Qgmwt4jKtVn9BWh5-zgy)\n",
    "[![Python Data Science](https://apmonitor.com/che263/uploads/Begin_Python/DataScience11.png)](https://www.youtube.com/watch?v=HReAo38LoM4&list=PLLBUgWXdTBDg1Qgmwt4jKtVn9BWh5-zgy \"Python Data Science\")\n",
    "\n",
    "Specific types of equations with differential terms arise from fundamental relationships such as conservation of mass, energy, and momentum. For example, the accumulation of mass $\\frac{dm}{dt}$ in a control volume is equal to the mass inlet $\\dot m_{in}$ minus mass outlet $\\dot m_{out}$ from that volume.\n",
    "\n",
    "$\\frac{dm}{dt} = \\dot m_{in} - \\dot m_{out}$\n",
    "\n",
    "Dynamic models can either be regressed (identified) from data or developed without data with fundamental relationships. Even fundamental relationships may have unknown or uncertain parameters. A combined approach for dynamic system modeling is to use fundamental physics-based relationships with data science. This approach uses the best features of both methods. It creates a model that aligns to measured values and extrapolates to regions where there is no or little data.\n",
    "\n",
    "![exercise](https://apmonitor.com/che263/uploads/Begin_Python/exercise.png)\n",
    "\n",
    "The first exercise is to [solve differential equations](https://apmonitor.com/pdc/index.php/Main/SolveDifferentialEquations) with several examples using `odeint`. The same examples are also [solved with Gekko](https://apmonitor.com/pdc/index.php/Main/PythonDifferentialEquations). For simulation, the two give equivalent results but Gekko is built to use differential equations in optimization or combined with machine learning. The `odeint` function has a specific purpose to solve ODEs. The function `odeint` solves the differential equation requires 3 inputs.\n",
    "\n",
    "    y = odeint(model, y0, t)\n",
    "\n",
    "1. `model` Function name that returns derivative values at requested `y` and `t` values as `dydt = model(y,t)`\n",
    "2. `y0` Initial conditions of the differential states\n",
    "3. `t` Time points where the solution is reported \n",
    "\n",
    "![idea](https://apmonitor.com/che263/uploads/Begin_Python/idea.png)\n",
    "\n",
    "### Solve Differential Equation\n",
    "\n",
    "We solve the differential equation with initial condition $y(0) = 5$:\n",
    "\n",
    "$ k \\, \\frac{dy}{dt} = -y$\n",
    "\n",
    "where $k=10$. The solution of `y` is reported from an initial time `0` to final time `20` and with a plot of the result for $y(t)$ versus $t$. Notice how the equation is rearranged to return just the derivative value as `dydt = -(1.0/k) * y` from the function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEGCAYAAABvtY4XAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3xUdb7/8dcnhYSQACEJoRNCU1ooASlSLKtgXxU76LqIfXV1r9e7d39Xd+/u3dW1rLjqguBd29qwXLsgoILSQgfpNdSEEkoggSTf3x8zeHPZJATMmZOZvJ+PxzwymXOO349nzrw5+Z7vfI855xARkcgT5XcBIiLiDQW8iEiEUsCLiEQoBbyISIRSwIuIRKgYvwsoLzU11WVkZPhdhohI2FiwYMFu51xaRctqVcBnZGSQk5PjdxkiImHDzDZXtkxdNCIiEUoBLyISoRTwIiIRSgEvIhKhFPAiIhHK01E0ZrYJOAiUAiXOuWwv2xMRkf8VimGS5zjndoegHRERKSfsu2hKSst4/qt1LM4t8LsUEZFaxeuAd8AUM1tgZmMrWsHMxppZjpnl5Ofnn3IDh4+V8trszTz49mKKjpX+2HpFRCKG1wE/yDnXGxgB3G1mQ05cwTk3wTmX7ZzLTkur8Nu2VWoYH8vjV2exPr+QJ75YXQMli4hEBk8D3jm3PfgzD3gf6OdFO2d3TGVU/7ZM+nYjczfs8aIJEZGw41nAm1kDM0s6/hy4AFjuVXsPjziDNk0S+NXkJRQWl3jVjIhI2PDyDD4dmGVmS4B5wCfOuc+9aqxBXAxPjMxi674j/NenK71qRkQkbHg2TNI5twHI8uq/X5G+GU0Yc3Y7Xpy5kQu6NmNop1Pv0xcRiRRhP0zyRA9e0JkOTRP518lL2X/kmN/liIj4JuICPj42midHZpF/qJjfffS93+WIiPgm4gIeIKt1Y+4e1p53F25lyoqdfpcjIuKLiAx4gHvO7UiX5g35t/eWkX+w2O9yRERCLmIDvl5MFH+5rieHikt4aPISnHN+lyQiElIRG/AAndKT+PVFZzJjdT6vzqn0toUiIhEpogMeYPSAtgzrnMYfPlnJ2l0H/S5HRCRkIj7gzYzHr+5Bg7gY7ntzMcUlmpBMROqGiA94gKZJ8Tx+VQ++33GAp6as8bscEZGQqBMBD3B+l3RuOKsNE2Zu4Lt1uv+IiES+OhPwAL+5+EzapTbggbeXUHD4qN/liIh4qk4FfEK9GJ65the7DxXz7+8v19BJEYlodSrgAbq3asSDF3Tmk2U7eHN+rt/liIh4ps4FPMDtQzIZ3DGVRz9cweqdGjopIpGpTgZ8VJTx1DU9SYqP5Z5/LOTIUQ2dFJHIUycDHiAtKY6nr81iXf4hfvvRCr/LERGpcXU24AEGd0zjzqHteXN+Lh8u2e53OSIiNapOBzzAL3/SiT5tk/n1e8vYvKfQ73JERGpMnQ/42OgonrmuJ1EG976xiKMlZX6XJCJSI+p8wAO0Sk7gzyOzWLp1P499vsrvckREaoQCPujCrs24eUBbJs3aqLtAiUhEUMCX8+uLz6R7y0Y8+M4S9ceLSNhTwJcTFxPN8zf2JsqMO19bSNExjY8XkfClgD9B6yYJPH1tFt/vOMCjH2p8vIiELwV8Bc49I527zwmMj38nR/PViEh4UsBX4oGfdGZg+xR+88Fyvt9+wO9yREROmQK+EtFRxjPX9aJR/Vjuen0BB4qO+V2SiMgpUcBXIS0pjudu7E3uviM89M5SzR8vImFFAX8SfTOa8PDwM/h8xU4mfLPB73JERKpNAV8NYwa34+LuzXns81XMXJvvdzkiItWigK8GM+Pxq3vQsWkS976xiNy9h/0uSUTkpBTw1dQgLoYJo/tQVua47ZUcDh8t8bskEZEqKeBPQduUBoy7vherdx3kocm66CoitZsC/hQN69yUf7mwMx8v3aGLriJSqyngT8OdQ9tzUfdmuugqIrWaAv40mBl/vjrrh4uuW/booquI1D4K+NN04kXXQ8W66CoitYvnAW9m0Wa2yMw+9rqtUGub0oDnbuzNuvxD3PfGIkrLdNFVRGqPUJzB3wesDEE7vhjcMY1HLu3CtFV5PK7b/YlILeJpwJtZK+BiYKKX7fht9IAMburfhvHfbND0wiJSa3h9Bv8X4CGgrLIVzGysmeWYWU5+fviOSHnk0q4M6pDCr99fxvxNe/0uR0TEu4A3s0uAPOfcgqrWc85NcM5lO+ey09LSvCrHc7HRUTx/Qx9aJSdw+6sLNJ2BiPjOyzP4QcBlZrYJeBM418xe87A93zVKiGXizdmUlJYx5mWNrBERf3kW8M65f3POtXLOZQDXAdOdczd51V5t0T4tkedv7KORNSLiO42D98DZHVN59LKuTFuVx28/WqE5a0TEFzGhaMQ59xXwVSjaqi1G9W9L7t7DTPhmA22aJDBmcKbfJYlIHROSgK+rHh5+Brl7D/OHT1fSKrk+w7s197skEalD1EXjoago4+lre9KzdWPue3Mxi7bs87skEalDFPAei4+NZuLobNIbxjPm5RxNTCYiIaOAD4GUxDj+/rO+lDrHLX+fR8Hho36XJCJ1gAI+RDLTEpkwKpute48w9pUFFB0r9bskEYlwCvgQ6teuCU9ck8W8TXu5/83FGiMvIp5SwIfYZVkt+H+XdOHzFTt55MPlGiMvIp7RMEkf/PzsduQdLGL81xtomhTPL87r6HdJIhKBFPA+eXj4Gew+eJSnpq4hNTGOG85q43dJIhJhFPA+MTP+dFV39hYW85sPlpGSWI8LuzbzuywRiSDqg/dRbHQUz93Ymx6tGnPvG4uYt1HzyItIzVHA+yyhXgwv3dKXVsn1GfPyfFbuOOB3SSISIRTwtUCTBvV45dZ+NIiLYdSkeWzcXeh3SSISARTwtUSr5ARe/flZOOe4aeJcthUc8bskEQlzCvhapEPTRF6+tR8Hio4xauJc8g8W+12SiIQxBXwt061lI/77lr7s2F/E6Jfmsf/wMb9LEpEwpYCvhbIzmjB+VB/W5x3iZ3+fR6Hu7Soip0EBX0sN6ZTGuOt7sji3gLGv5mhyMhE5ZQr4Wmx4t+Y8fnUW367bw12vL+RoSZnfJYlIGFHA13JX92nF76/oxvRVedzzj4UcK1XIi0j1KODDwE392/Lby7oy5ftd3PfmIkoU8iJSDZqLJkzcPDCDY6Vl/P6TlURHLeHpa7KIida/zyJSOQV8GBkzOJOSMsefPltFTJTxxMgsoqPM77JEpJZSwIeZO4a2p6S0jCemrCE6ynj8qh5EKeRFpAIK+DB0z7kdOVbqeGbaWqLN+OOV3RXyIvJPFPBh6v7zO+KcY9z0dZSUOR6/uoe6a0Tk/1DAhykz44ELOhMdFcXTX66hpKyMJ0fqwquI/C8FfJi77/yOxEQbf/5iNSWljr9c15NYhbyIoICPCHef04F60VH84dOVlJSV8ez1vakXo5AXqeuUAhHitiGZPHJpF75YsYs7X1tAcYnmrhGp6xTwEeRng9rxn1d0Y9qqPG57ZQFHjirkReoyBXyEGdW/LY9d1Z2Za/MZ/dJcDhRpPnmRukoBH4Gu7duGcdf1YtGWAm54cQ57DunOUCJ1kQI+Ql2a1YIXR2ezdtchrhk/m537i/wuSURCTAEfwc45oymv3NqPXQeKufpv37Fpd6HfJYlICCngI9xZmSm8cVt/CotLGDl+Nqt3HvS7JBEJEQV8HdC9VSPevn0AUQbXjJ/Ngs37/C5JRELAs4A3s3gzm2dmS8xshZn91qu25OQ6picx+Y6BJCfEcuPEOUxbucvvkkTEY16ewRcD5zrnsoCewHAz6+9he3ISrZskMPnOgXRKT2Lsqwt4e36u3yWJiIc8C3gXcCj4a2zw4bxqT6onNTGON27rz6AOqTz07lL+On0tzultEYlEnvbBm1m0mS0G8oCpzrm5Fawz1sxyzCwnPz/fy3IkqEFcDBNHZ/PTXi15YsoaHvlwBaVlCnmRSONpwDvnSp1zPYFWQD8z61bBOhOcc9nOuey0tDQvy5Fy6sVE8eTILMYOyeSV2Zu5942FFB3T1AYikaTK2STNLB64BBgMtACOAMuBT5xzK6rbiHOuwMy+AoYHt5daICrK+PVFZ5KWGMcfPl1J/sG5TBiVTXKDen6XJiI1oNIzeDN7FPgWGADMBcYDbwMlwJ/MbKqZ9ahi+zQzaxx8Xh84H1hVc6VLTbltSCbPXt+LJVv3c+UL+kKUSKSo6gx+vnPu0UqWPWVmTYE2VWzfHHjZzKIJ/EPytnPu49MrU7x2aVYLmjeK57ZXcrjyhe94cXQ2fdom+12WiPwIlZ7BO+c+ATCzkScuM7ORzrk851xOFdsvdc71cs71cM51c879rmZKFq9kZzThvbsG0TA+hutfnMMnS3f4XZKI/AjVucj6b9V8TSJAu9QGvHfXIHq0bMTd/1jI+K/XaxilSJiqtIvGzEYAFwEtzWxcuUUNCfTDS4Rq0qAer405i1+9s4Q/fraKTXsK+e1l3XQbQJEwU1Uf/HZgAXBZ8OdxB4FfelmU+C8+Nppx1/WibUoCz81Yz4b8Qv52Ux+NsBEJI3ayP7/NLNY5F5LbAmVnZ7ucnEq79cUnHyzaxkPvLqVZw3gm3ZxNx/Qkv0sSkSAzW+Ccy65oWVXDJD8ys0srWZZpZr8zs1trqkipva7o1ZI3x/bn8NFSfvr8d8xYled3SSJSDVV1qt5G4AtOK81svpl9amYzzGwjgTHxC5xzL4WkSvFd7zbJfHjPINo0SeDnL89n4swNuvgqUstV2gfvnNsJPGRmucAsIJ7AN1nXOOcOh6g+qUVaNK7P5DsH8MBbS/j9JytZtfMgv7+iG/Gx0X6XJiIVqM6wiHTgHQIXVpsRCHmpoxLqxfD8jb35xXkdmbxgK9eOn82O/TokRGqjkwa8c+43QEdgEnALsNbM/svM2ntcm9RSUVHGAz/pxPhRfVifX8ilz85i7oY9fpclIieo1sBmF+hs3Rl8lADJwGQze9zD2qSWu7BrMz64eyAN68dy48S5vPzdJvXLi9QiJw14M/uFmS0AHicw+Vh359ydQB/gKo/rk1quQ9MkPrh7EMM6p/HIhyv41TtLNe2wSC1R5XTBQanAlc65zeVfdM6Vmdkl3pQl4aRhfCwTRmUzbvpa/vLlWlbvOsALN/ahdZMEv0sTqdOq0wf/HyeGe7llK2u+JAlHUVHG/ed3YuLobDbvOcwlz87Sjb1FfKbJRaRGnd8lnU/uHUyr5Pr8/OUcHv98FSWlZX6XJVInKeClxrVJSeDdOwdyfb/WPP/Vem6aNJe8g0V+lyVS5yjgxRPxsdH88coePDEyi8W5BVw8TkMpRUJNAS+eurpPKz64exCJcTHcMHEuf52+ltIyDaUUCQUFvHjujGYN+fCeQVzcvTlPTFnDqElzyTugLhsRryngJSSS4mN55rqePH5VDxZu2ceIZ2by1WrNSiniJQW8hIyZcU3f1nx0z9mkJcVxy3/P54+fruRoiUbZiHhBAS8h1zE98O3Xm/q3Yfw3Gxg5fjab9xT6XZZIxFHAiy/iY6P5/RXdeeHG3mzIP8RFz8zk7ZxczWUjUoMU8OKrEd2b8/n9Q+jeqhEPTV7KXa8vZF/hUb/LEokICnjxXcvG9Xl9TH8eHnEGX67cxfBnvmHW2t1+lyUS9hTwUitERxl3DG3P+3cFxszfNGku//nx95qZUuRHUMBLrdKtZSM+vncwowe0ZdKsjVz211ks27rf77JEwpICXmqd+vWi+d3l3fj7z/qy/8gxfvr8tzw9dQ3HNGmZyClRwEutNaxzU6bcP5RLs1rwzLS1/PT5b1m986DfZYmEDQW81GqNEmJ5+tqe/O2mPuwoKOLSZ2fxwlfrNZ+NSDUo4CUsDO/WjCm/HMJ5Zzblsc9XcdUL37Fml87mRaqigJewkZIYx/M39uaZ63qyeU8hl4ybxbhpa9U3L1IJBbyEFTPj8p4tmfrAUC7s1oynpq7h0mc10kakIgp4CUupiXE8e30vXhydzd7Co1z+3Cz++NlKjZsXKUcBL2HtJ13SmfrAUK7Jbs34rzcw4pmZfLdO34IVAQW8RIBG9WP501U9eH3MWZQ5xw0T5/LAW4vZc6jY79JEfKWAl4gxqEMqX9w/hHvO6cBHS7dz3lNf8/Z8zVApdZdnAW9mrc1shpmtNLMVZnafV22JHBcfG82vLuzMJ78YTMemiTz07lKunTCHdXkaUil1j5dn8CXAg865M4H+wN1m1sXD9kR+0Ck9ibfGDuCxq7qzeudBRjwzk8c+X8XhoyV+lyYSMp4FvHNuh3NuYfD5QWAl0NKr9kROFBVlXNu3DdMeDEx38MJX6zn/ya/5bNkOddtInRCSPngzywB6AXMrWDbWzHLMLCc/Pz8U5Ugdk5oYx1PX9GTyHQNolFCPO19fyOiX5rE+/5DfpYl4yrw+kzGzROBr4A/OufeqWjc7O9vl5OR4Wo/UbSWlZbw2ZzNPTl1D0bFSxgzO5J5zOtAgLsbv0kROi5ktcM5lV7TM0zN4M4sF3gVeP1m4i4RCTHQUtwxqx/QHh3FZVkte+Go95z75Fe8t3EqZJjCTCOPlKBoDJgErnXNPedWOyOlIS4rjyWuyePfOgTRrGM8Dby/hyhe+Y+GWfX6XJlJjvDyDHwSMAs41s8XBx0Uetidyyvq0Teb9uwbxxMgsthcc4crnv+OXby1m5/4iv0sT+dE874M/FeqDFz8VFpfw/FfreHHmRqItcI/Y24a0I6Ge+uel9vKtD14knDSIi+FfLjyDaQ8MZVjnNJ7+cg3D/vwVb83fohuMSFhSwIucoHWTBF64qQ+T7xhAy+T6/Ou7y7h43Ey+XqNhvBJeFPAilcjOaMJ7dw7kuRt6c/hoKTe/NI9Rk+by/fYDfpcmUi0KeJEqmBkX92jO1AeG8P8u6cLSrfu5+NmZ3PfmIjbvKfS7PJEq6SKryCnYf/gY479Zz0vfbqSk1HF9vzbce24HmjaM97s0qaOqusiqgBc5DXkHihg3fS1vzsslNjqKnw3K4Pah7WlUP9bv0qSOUcCLeGTT7kKemrqGD5dsp2F8DGOHZHLLoHYkauoDCREFvIjHVmzfz9NT1/DlyjySE2K5fWh7Rg9oqzH04jkFvEiILM4t4Ompa/h6TT6pifW4Y2h7burflvjYaL9LkwilgBcJsQWb9/L01LXMWrebpklx3D60PTf0a0P9egp6qVkKeBGfzNmwh798uYY5G/aSmliPMYMzGdW/raYnlhqjgBfx2fxNexk3bS0z1+6mcUIsY85ux+iBGTSM16gb+XEU8CK1xKIt+/jr9HVMW5VHw/gYbh6YwS0DM0hJjPO7NAlTCniRWmb5tv38dfo6vvh+J3ExUVyb3ZrbhmTSKjnB79IkzCjgRWqpdXmHmPDNet5ftI0yB5dnteD2oe3p3CzJ79IkTCjgRWq57QVHmDRrI2/M28Lho6Wcd0ZTxgzOpH9mEwI3RxOpmAJeJEzsKzzKK7M388rsTewpPEq3lg25bXAmF3VvTmy05gaUf6aAFwkzRcdKeW/hNibO2sCG/EJaNIrn1rPbcW3f1iRp5I2Uo4AXCVNlZY4Zq/OY8M0G5m7cS2JcDCOzW3HzgAwyUhv4XZ7UAgp4kQiwdGsBL83ayCfLdlBS5ji3c1NuGZTB2R1S1U9fhyngRSJI3oEiXpu7hX/M3czuQ0fp2DSRWwZlcEXPlvqGbB2kgBeJQMUlpXy0ZAf//e1GVmw/QFJcDFf1acVN/dvSoWmi3+VJiCjgRSKYc46FW/bx6uzNfLpsJ0dLyxjYPoVR/dvyky7pxGj0TURTwIvUEbsPFfPW/Fz+MXcL2wqO0KxhPNf0bc21fVvTsnF9v8sTDyjgReqY0jLHjFV5vDpnM9+szceAYZ2bcn2/NpzTOU1n9RFEAS9Sh+XuPczbObm8nZPLrgPFpDeM49rs1ozMbk3rJpr7Jtwp4EWEktIyZqzO5x9zN/PVmnycg0EdUrgmuzUXdm2mu06FKQW8iPwf2wqO8O6CrbyzIJfcvUdIio/h8p4tGNmnNT1aNdK4+jCigBeRCpWVOeZs3MM7OVv5dNkOikvK6JSeyJW9W/HTXi1Jbxjvd4lyEgp4ETmpA0XH+GjJdt5dsJWFWwqIMhjUIZWrerfigq7pJNTTl6hqIwW8iJySjbsLeX/hVt5btI2t+47QoF40w7s154peLRjYPpXoKHXh1BYKeBE5LWVljvmb9vLuwq18tmwnB4tLSE2M49Ks5lzesyVZ6q/3nQJeRH60omOlzFiVx/8s3s70VXkcLS0jIyWBy7JacGlWCzqm6y5UflDAi0iN2n/kGF+s2Mn/LN7Gd+v34Bx0Tk/i4h7NuaRHczLTNBdOqCjgRcQzeQeL+GzZTj5eup35m/YB0KV5Qy7u0ZyLuzfXvPUeU8CLSEjs3F/EJ8t28PHS7SzaUgDAGc2SGNGtORd1b6ZuHA8o4EUk5LYVHOHz5Tv5fPkOcjbvwzlon9aAEd2aM7xbM7q2aKgLtDXAl4A3s5eAS4A851y36myjgBeJTHkHivhixU4+W76TORv2UOagRaN4LujajAu6pNO3XRPdVPw0+RXwQ4BDwCsKeBE5bm/hUaat3MWU73fxzZp8ikvKaFQ/lvPOaMr5XdIZ3DFVNxY/Bb510ZhZBvCxAl5EKnL4aAkz1+5myopdTFu1i4LDx4iNNvpnpnDuGU05/8x0zXh5ErU64M1sLDAWoE2bNn02b97sWT0iUnuVlJaxYPM+pq/K48uVu1ifXwhAx6aJnHdmOud0TqN322R15ZygVgd8eTqDF5HjNu0u5MuVu5i+Ko95G/dSUuZIiothcKdUhnVuyrBOaTTVZGhVBrxmDxKRWikjtQFjBmcyZnAmB4uO8e263cxYlc+M1Xl8umwnAF1bNGRIpzSGdEyjT9tk6sXo7L48ncGLSFhxzrFyx0FmrM7j69X5LNyyj5IyR0K9aAZkpgQCv1MaGSkJdWIYpl+jaN4AhgGpwC7gEefcpKq2UcCLyKk6WHSM2ev38M3afL5Zs5stew8D0LJxfc7ukMqgjqkMbJ9CamKcz5V6Q190EpE6Y9PuQmauzWfWut18t34PB4tKgMA3as/ukMqgDqn0bdeExLjI6KFWwItInVRa5li2bT/frtvNt+t2k7NpH0dLy4iOMrq3bMSA9ikMyEwhOyM5bG9oooAXEQGOHC1l4ZZ9zF6/h9kb9rAkt4CSMkdstJHVqjH92jXhrMwU+rRNDpszfAW8iEgFCotLmL9pL7M37GHuhr0s27af0jJHdJTRrUVDzspMoV9GE7IzkmmcUM/vciukgBcRqYbC4hIWbtnH3A17mbdxL4tzCzhaWgZAp/REsjOa0Dcjmey2TWiVXL9WjNJRwIuInIaiY6Uszi1gweZ9zN+0lwWb9nGwOHDRtlnDePq0TaZXm8b0aZtM1xaNfBmHry86iYichvjYaPpnptA/MwUIXLRds+sgOZv2Mm/TPhZu3scny3YAUC8mih4tG/0Q+j1bJ9Oskb/ftNUZvIjIj7DrQBELN+9jweZ9LNyyj+XbDvzQrdO8UTw9WzemZ+vG9GqTTPeWjahfL7pG29cZvIiIR9IbxjOie3NGdG8OQHFJKd9vP8CiLQUszg08PlsemFohyqBTehJZrRrTo3Ujslo1pnOzJM8mUFPAi4jUoLiYaHq1SaZXm+QfXttzqJjFuQUsyS1gydb9TPl+J2/l5AKBrp2erRrz5tj+REXV7EVbBbyIiMdSEuM478x0zjszHQjMp5O79whLthawdGsBB4tKajzcQQEvIhJyZkablATapCRwaVYLz9rR3JoiIhFKAS8iEqEU8CIiEUoBLyISoRTwIiIRSgEvIhKhFPAiIhFKAS8iEqFq1WRjZpYPbD7NzVOB3TVYTk1RXadGdZ0a1XVqIrGuts65tIoW1KqA/zHMLKeyGdX8pLpOjeo6Narr1NS1utRFIyISoRTwIiIRKpICfoLfBVRCdZ0a1XVqVNepqVN1RUwfvIiI/F+RdAYvIiLlKOBFRCJU2AW8mQ03s9Vmts7MHq5geZyZvRVcPtfMMkJQU2szm2FmK81shZndV8E6w8xsv5ktDj7+w+u6gu1uMrNlwTb/6Y7mFjAuuL+WmlnvENTUudx+WGxmB8zs/hPWCcn+MrOXzCzPzJaXe62JmU01s7XBn8mVbHtzcJ21ZnZzCOr6s5mtCr5P75tZ40q2rfI996CuR81sW7n36qJKtq3ys+tBXW+Vq2mTmS2uZFsv91eF2RCyY8w5FzYPIBpYD2QC9YAlQJcT1rkL+Fvw+XXAWyGoqznQO/g8CVhTQV3DgI992GebgNQqll8EfAYY0B+Y68N7upPAlzVCvr+AIUBvYHm51x4HHg4+fxh4rILtmgAbgj+Tg8+TPa7rAiAm+PyxiuqqznvuQV2PAr+qxvtc5We3pus6YfmTwH/4sL8qzIZQHWPhdgbfD1jnnNvgnDsKvAlcfsI6lwMvB59PBs4zs5q/2WE5zrkdzrmFwecHgZVASy/brEGXA6+4gDlAYzNrHsL2zwPWO+dO9xvMP4pz7htg7wkvlz+GXgauqGDTC4Gpzrm9zrl9wFRguJd1OeemOOdKgr/OAVrVVHs/pq5qqs5n15O6gp//a4A3aqq96qoiG0JyjIVbwLcEcsv9vpV/DtIf1gl+GPYDKSGpDgh2CfUC5laweICZLTGzz8ysa4hKcsAUM1tgZmMrWF6dfeql66j8g+fH/gJId87tgMAHFGhawTp+77dbCfzlVZGTvedeuCfYdfRSJd0Nfu6vwcAu59zaSpaHZH+dkA0hOcbCLeArOhM/cZxnddbxhJklAu8C9zvnDpyweCGBbogs4Fngg1DUBAxyzvUGRgB3m9mQE5b7ub/qAZcB71Sw2K/9VV1+7rd/B0qA1ytZ5WTveU17AWgP9AR2EOgOOZFv+wu4nqrP3j3fXyfJhko3q+C1U9pn4RbwW4HW5X5vBWyvbB0ziwEacXp/Up4SM4sl8Aa+7px778TlzrkDzrlDweefArFmlup1Xc657cGfecD7BP5ULq86+9QrI4CFzrldJy7wa38F7VxeKVIAAAL4SURBVDreTRX8mVfBOr7st+CFtkuAG12wo/ZE1XjPa5RzbpdzrtQ5Vwa8WEl7fu2vGOBK4K3K1vF6f1WSDSE5xsIt4OcDHc2sXfDs7zrgwxPW+RA4frX5amB6ZR+EmhLs45sErHTOPVXJOs2OXwsws34E9v0ej+tqYGZJx58TuEi3/ITVPgRGW0B/YP/xPx1DoNIzKz/2Vznlj6Gbgf+pYJ0vgAvMLDnYJXFB8DXPmNlw4F+By5xzhytZpzrveU3XVf6azU8raa86n10vnA+scs5trWih1/urimwIzTHmxZVjLx8ERn2sIXBF/t+Dr/2OwEEPEE/gT/51wDwgMwQ1nU3gT6elwOLg4yLgDuCO4Dr3ACsIjB6YAwwMQV2ZwfaWBNs+vr/K12XAc8H9uQzIDtH7mEAgsBuVey3k+4vAPzA7gGMEzph+TuCazTRgbfBnk+C62cDEctveGjzO1gE/C0Fd6wj0yR4/xo6PFmsBfFrVe+5xXa8Gj52lBIKr+Yl1BX//p8+ul3UFX//78WOq3Lqh3F+VZUNIjjFNVSAiEqHCrYtGRESqSQEvIhKhFPAiIhFKAS8iEqEU8CIiEUoBL3WWmTU2s7uCz1uY2WS/axKpSRomKXVWcG6Qj51z3XwuRcQTMX4XIOKjPwHtg/OErwXOdM51M7NbCMzuFw10IzC3Sj1gFFAMXOSc22tm7Ql8SSwNOAzc5pxbFfr/DZGKqYtG6rKHCUxV3BP4lxOWdQNuIDAvyR+Aw865XsBsYHRwnQnAvc65PsCvgOdDUrVINekMXqRiM1xg/u6DZrYf+Cj4+jKgR3B2wIHAO+VuNxAX+jJFKqeAF6lYcbnnZeV+LyPwuYkCCoJn/yK1krpopC47SOA2aqfMBeb03mhmI+GHe9tm1WRxIj+WAl7qLOfcHuDb4I2a/3wa/4kbgZ+b2fGZCGvsFnQiNUHDJEVEIpTO4EVEIpQCXkQkQingRUQilAJeRCRCKeBFRCKUAl5EJEIp4EVEItT/B1TSECwHw+CzAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.integrate import odeint\n",
    "\n",
    "# function that returns dy/dt\n",
    "def model(y,t):\n",
    "    k = 10.0\n",
    "    dydt = -(1.0/k) * y\n",
    "    return dydt\n",
    "\n",
    "y0 = 5                 # initial condition\n",
    "t = np.linspace(0,20)  # time points\n",
    "y = odeint(model,y0,t) # solve ODE\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.plot(t,y)\n",
    "plt.xlabel('time'); plt.ylabel('y(t)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![gekko](https://apmonitor.com/che263/uploads/Begin_Python/gekko.png)\n",
    "\n",
    "### Solve Differential Equation with Gekko\n",
    "\n",
    "[Python Gekko](https://gekko.readthedocs.io/en/latest/) solves the same differential equation problem. Gekko is built for large-scale problems. [Additional tutorials on Gekko](https://apmonitor.com/wiki/index.php/Main/GekkoPythonOptimization) show how to solve other types of equations and optimization problems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEGCAYAAABvtY4XAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3gU173/8fd3VRBNCJAQIIpML6LLmGrcwOAYN4K7cbvB3U5ix3FycxMn9+aXOI57i7FxcNwNLnHBBVyopkj0ZkTvkmhCFPXz+0OLry6WZGG0O7urz+t59tFqZ8bz9ezuh6MzZ86Ycw4REYk8Pq8LEBGRwFDAi4hEKAW8iEiEUsCLiEQoBbyISISK9rqAihITE11qaqrXZYiIhI3MzMy9zrmkypaFVMCnpqaSkZHhdRkiImHDzLZWtUxdNCIiEUoBLyISoRTwIiIRSgEvIhKhFPAiIhEqoKNozGwLkA+UAiXOufRA7k9ERP5XMIZJnu2c2xuE/YiISAVh30VTWuZ49usNLNt+0OtSRERCSqAD3gGfm1mmmU2sbAUzm2hmGWaWkZube9I7OFJUwqvfbOWXby/jWFHpqdYrIhIxAh3wQ51z/YExwB1mduaJKzjnJjnn0p1z6UlJlV5tW634uBgeHt+HTblHeOjTdbVQsohIZAhowDvndvl/5gDvAQMDsZ+hnRK5YUgqU+ZvYd4GdfeLiEAAA97MGppZ4+PPgVHAqkDt79eju9EhqSH3TV1O3rHiQO1GRCRsBLIFnwzMNbPlwCLgY+fcp4HaWf3YKB67vC85+YX88YPVgdqNiEjYCNgwSefcJqBPoP77lenTNoE7z+7EE19kMbJHMmN6tQrm7kVEQkrYD5M80Z3ndKJXShN++95KcvILvC5HRMQzERfwMVE+HruiD0eLSvntuytxznldkoiIJyIu4AE6tWjMr0d3Y+baHN7O2O51OSIinojIgAe4YUgqQzo2548frmHz3iNelyMiEnQRG/A+n/HI5X2IifJxz5tLKSop87okEZGgitiAB2jVpD4PjevFih15PDZzvdfliIgEVUQHPMDotFZceXpb/jFrI/M36ipXEak7Ij7gAX4/tgenNW/IL99azsGjRV6XIyISFHUi4BvERvPkVf3Yd6SQB97R0EkRqRvqRMADpKU04b5RXfl09R7eWqyhkyIS+epMwAP8bHgHhnYqHzq5Mfew1+WIiARUnQp4n8949PK+xMX4uPuNpRSW6AYhIhK56lTAAyTHx/G3n/Zh9a5D/GW6bhAiIpGrzgU8wMgeydw09DSmzN/Cp6t2e12OiEhA1MmAB3hgTDf6tGnCr6atYPv+o16XIyJS6+pswMdG+3j66v4A3Pn6Ek1lICIRp84GPEDbZg14+Ke9Wb4jTzfsFpGIU6cDHsqnMrh+cHsmz93MjDXZXpcjIlJr6nzAA/z2J91JS4nnvqnL2XFA/fEiEhkU8EC96Cievqo/pWWOu95YSnGp+uNFJPwp4P1SExvy13G9WLrtoMbHi0hEUMBXcGHv1twwJJWX5m3moxW7vC5HROSUKOBP8NsLujOgfVPun7aCrOx8r8sREfnRFPAniI328czV/WkQG8Wtr2ZyuLDE65JERH4UBXwlWjaJ46mr+rN57xF+PW2F5o8XkbCkgK/C4I7NuX90Nz5euZvJczd7XY6IyElTwFfjljM7cH7PZP7yyToWbd7vdTkiIidFAV8NM+Ph8X1o16wBd7y+hJxDBV6XJCJSYwr4HxAfF8Nz1/bncEEJt722RDcJEZGwoYCvgW4t43l4fG8ytx7gwQ/WeF2OiEiNRHtdQLi4sHdr1uw6xLNfb6Rn63iuHdTe65JERKqlFvxJuHdUV87umsSDH6zWSVcRCXkK+JMQ5TOeuKof7Zo14PbXMtl18JjXJYmIVEkBf5Li42KYNCGdwuIyJr6SQUGxTrqKSGhSwP8InVo04vEr+7J61yEeeEdXuopIaFLA/0jndk/m3pFdeH/ZLl6Ys8nrckREvkcBfwruOLsTP+nVir98so6Zut2fiISYgAe8mUWZ2VIz+yjQ+wo2M+Pv4/vQK6UJd7+5lDW7DnldkojId4LRgr8HWBuE/XiifmwUL0xIJz4uhv94eTE5+ZrOQERCQ0AD3szaAD8BXgzkfryWHB/Hi9enc+BoMRP/lamRNSISEgLdgn8cuB+o8i7WZjbRzDLMLCM3NzfA5QROWkoTHruiL8u2H+RXmkNeREJAwALezC4EcpxzmdWt55yb5JxLd86lJyUlBaqcoBid1pL7R3flw+W7eOKLLK/LEZE6LpBz0QwFLjKzC4A4IN7MXnXOXRvAfXruthEd2ZBzmMdnZtEhqREX9WntdUkiUkcFrAXvnPuNc66Ncy4VuBL4MtLDHcpH1vzlsl6cntqU+6YuZ/EWzVkjIt7QOPgAqBcdxaTr0klJqM/P/pXBxtzDXpckInVQUALeOfe1c+7CYOwrVDRtGMuUG08nyowb/rmI3PxCr0sSkTpGLfgAat+8IS9en05ufiH/8fJijhaVeF2SiNQhCvgA69euKU9e2Y8VO/O4+41llJZp+KSIBIcCPghG9WzJg2N7MnNtNn/8cLXGyItIUOiWfUFy/ZBUdhw4ygtzNtO2aQN+dmYHr0sSkQingA+i34zpzq6DBfx5+lpaxNfj4r4pXpckIhFMAR9EPp/xyOV92HekkHvfXk6T+jGc1bWF12WJSIRSH3yQxcVEMWlCOl2SG3Pbq0tYuu2A1yWJSIRSwHsgPi6GKTedTlLjetw0ZTEbcnQhlIjUPgW8R1o0juOVmwcS5fMxYfJCducd87okEYkwCngPtW/ekCk3nk5+QQkTJi/i4NEir0sSkQiigPdYWkoTJk1IZ+v+o9w0RVe7ikjtUcCHgMEdm/PkleU3C7nllUwKS3RHKBE5dQr4EDE6rRUPjevNnKy93PX6UopLq7wJlohIjSjgQ8j49Lb88aKefL4mm19NXU6Z5q0RkVOgC51CzPVDUjlcWMLDn31Lw3rR/M8laZiZ12WJSBhSwIegO87uxOHCEp77eiMN60XzmzHdFPIictIU8CHq/vO7cqSwhEmzN9GoXjR3n9vZ65JEJMwo4EOUmfHg2J4cKSzl0RnrqR8TpRkoReSkKOBDmM9nPDSuFwXFpfx5+lp8PuPmYad5XZaIhAkFfIiLjvLx+JV9KXOO//5oDdE+4/ohqV6XJSJhQMMkw0BMlI8nr+rHqB7J/OGD1byyYKvXJYlIGFDAh4mYKB9PX92f87q34L/eX8XrC7d5XZKIhDgFfBiJjfbxzDX9OadbC3773kreWqyQF5GqKeDDTL3oKJ69pj8juiTxwLsreTtju9cliUiIUsCHobiYKJ6/bgDDOiVy/7QV6q4RkUop4MNUXEwUL0xI/6675uX5W7wuSURCjAI+jMXFRPGPawd8N7rmhdmbvC5JREKIAj7MHT/x+pPerfjz9LU889UGr0sSkRChC50iQEyUjyeu6EtslI+HP/uWopIyfn5eZ01QJlLHKeAjRHSUj7+P70O0z3jiiywKS8r49eiuCnmROkwBH0GifMZD43oTG+3jH7M2kl9QzH9fnIbPp5AXqYsU8BHG5zP+55I0GsfF8I9ZGzlcWMLfx/chJkqnW0TqGgV8BDIzHhjTjfj60fzt0285XFDCM9f0Jy4myuvSRCSI1KyLYLef1Yn/viSNL7/N4fqXFpFfUOx1SSISRAr4CHfdoPY8fkVfMrce4OoXFrL/SJHXJYlIkCjg64CL+6bw/HUDWJ+dz+XPf8Oug8e8LklEgkABX0ec2z2Zl28aSHZeAZc9O5/12flelyQiARawgDezODNbZGbLzWy1mf0xUPuSmhnUoTlv3TKYMuf46XPzWbR5v9cliUgABbIFXwic45zrA/QFRpvZoADuT2qgR+t43rltCImN63Ht5IV8umqP1yWJSIAELOBducP+X2P8Dxeo/UnNtW3WgGm3DqFHq3hufy2TV3ULQJGIFNA+eDOLMrNlQA4wwzm3sJJ1JppZhpll5ObmBrIcqaBZw1he/9kZnNW1Bb97fxWPzliPc/r3VySSBDTgnXOlzrm+QBtgoJmlVbLOJOdcunMuPSkpKZDlyAkaxEYz6boBXJ7ehie/yOJX01ZQVFLmdVkiUkuqvZLVzOKAC4HhQGvgGLAK+Ng5t7qmO3HOHTSzr4HR/u0lRERH+XhoXG9aJ9Tn8ZlZ7Dp4jOeuHUCT+jFelyYip6jKFryZPQjMAwYDC4HngbeBEuCvZjbDzHpXs32SmSX4n9cHzgPW1V7pUlvMjJ+f14VHxvdh8Zb9jHtuPtv3H/W6LBE5RdW14Bc75x6sYtmjZtYCaFfN9q2Al80sivJ/SN52zn3048qUYBg3oA2tE+pzyysZXPrsPF68/nT6tk3wuiwR+ZGqbME75z4GMLPxJy4zs/HOuRznXEY1269wzvVzzvV2zqU55/5UOyVLIA3u2Jx3bx9C/dgorpz0jYZRioSxmpxk/U0NX5MI0alFY967fSjdWsZz22uZTJq9USNsRMJQlV00ZjYGuABIMbMnKyyKp7wfXiJYYqN6vPGzQdw7dRn/b/o61mcf5s+XplEvWlMOi4SL6lrwu4BMoMD/8/jjA+D8wJcmXqsfG8XTV/XnnnM7My1zB1e/sJC9hwu9LktEash+6E9vM4txzgVlIvH09HSXkVFlt7546KMVu7j37eUkNqrHi9en071VvNcliQhgZpnOufTKllU3TPJDMxtbxbIOZvYnM7uptoqU0HZh79ZMvXUwJWVljHtuPp+v1slXkVBXXRfNzyi/wGmtmS02s+lm9pWZbaZ8THymc+6loFQpIaF3mwQ+uHMYnVs04pZXM3n6yyydfBUJYVWeZHXO7QHuN7PtwFwgjvIrWdc753QVTB2VHB/HW7cM5tfvrODvn69n1c5D/P3yPjSqp9v7ioSamgyTTAamAr8AWlIe8lKHxcVE8fgVffndT7rz+Zo9XPrMPDbvPeJ1WSJygh8MeOfc74DOwGTgBiDLzP6fmXUMcG0SwsyM/xjegVdvPoO9hwu56Om5fLE22+uyRKSCGs0m6co7Wvf4HyVAU2Camf0tgLVJGBjSKZEP7xpGu2YNuPnlDJ6YmUVZmfrlRULBDwa8md1tZpnA3yiffKyXc+42YAAwLsD1SRho07QB79w2hEv7pfDYzPVMfCWTvGNBGVkrItWoSQs+EbjMOXe+c27q8THxzrkyyqcSFiEuJopHL+/DH8b24Otvcxj71FxW7czzuiyROq0mffC/d85Vek8359za2i9JwpWZcePQ03jrlkEUlZRx2XPzeXPRNg2lFPFIQO/oJHXTgPbN+PjuYZxxWjMeeHcl901dwbGiUq/LEqlzFPASEM0b1WPKjQO559zOvLt0B5c+O49NuYd/eEMRqTUKeAmYKJ/xi5FdmHLjQLIPFXDR0/P497KdXpclUmco4CXgRnRJ4uO7h9OtZWPueXMZ909bztEizTgtEmgKeAmK1gn1eXPiIO48uxNTM3dw0dPzWLfnkNdliUQ0BbwETXSUj/vO78qrN59B3rFiLn56Hq8t3KpRNiIBooCXoBvaKZHpdw/njA7N+c/3VnHH60t0YZRIACjgxRNJjesx5YbT+c2Ybny+Opsxj89mwaZ9XpclElEU8OIZn8+4ZURH3rltCPViorjqhQX89ZN1FJWUeV2aSERQwIvn+rRN4KO7hnHl6W35x6yNXPbcPDbkaMy8yKlSwEtIaFgvmr9c1pvnrxvAzgPHuPCpObyyQCdgRU6FAl5Cyvk9W/LZz89k4GnN+a/3V3HjlMVkHyrwuiyRsKSAl5DTIj6OKTeczoNje7Bg0z5GPTabfy/bqda8yElSwEtI8vmMG4aexvS7h9MhqSH3vLmM219bwr7DhV6XJhI2FPAS0jokNWLarUN4YEw3vlibw6jHZvPZ6j1elyUSFhTwEvKifMatIzry4V3DaNkkjlteyeTnby7lwJEir0sTCWkKeAkbXVs25v07hnLPuZ35aMVuRj42i49X7FbfvEgVFPASVmKifPxiZBc+vGsYrZrU547Xl3Drq5nkaKSNyPco4CUsdW8Vz3u3l/fNf/1tLuc9OoupGdvVmhepQAEvYSs6ysetIzryyT3D6dqyMb+atoIJLy1i276jXpcmEhIU8BL2OiQ14q2Jg/nTxT1ZsvUAox6fxXNfb6S4VHPaSN2mgJeI4PMZEwanMvPeEYzoksRDn65j7FNzWbLtgNeliXhGAS8RpVWT+jx/XTqTrhtA3rFixj03n9+9v5JDBZpvXuoeBbxEpFE9WzLjlyO4YUgqry/cxnmPzNJ0B1LnBCzgzaytmX1lZmvNbLWZ3ROofYlUplG9aP4wtifv3zGU5Pg47nlzGVe/sJCs7HyvSxMJikC24EuAe51z3YFBwB1m1iOA+xOpVO82Cbx/x1D+55I01uw+xJgn5vCX6Ws5UljidWkiARWwgHfO7XbOLfE/zwfWAimB2p9IdaJ8xrWD2vPlvSMY178Nz8/exLmPzOKjFbvUbSMRKyh98GaWCvQDFlaybKKZZZhZRm5ubjDKkTqseaN6PPTT3rxz2xCaNYzlzteXcvULC1m7+5DXpYnUOgt068XMGgGzgD87596tbt309HSXkZER0HpEjistc7y+cCuPzFjPoWPFXHNGe345sgtNG8Z6XZpIjZlZpnMuvbJlAW3Bm1kM8A7w2g+Fu0iwRfmM6wan8vV9ZzFhcCqvL9rGWX//minzNusiKYkIgRxFY8BkYK1z7tFA7UfkVCU0iOXBi3ryyT3D6ZXShAc/XMMFT8xh9np1GUp4C2QLfihwHXCOmS3zPy4I4P5ETkmX5Ma8cvNAJl03gMKSMia8tIgJLy3i2z0aVinhKeB98CdDffASKgpLSnnlm6089eUG8guKuTy9Lb8c2YUW8XFelybyf1TXB6+AF6nGwaNFPPXlBv71zRZionxMPLMDE8/sQIPYaK9LEwE8PMkqEu4SGsTyXxf2YMYvRnBW1yQen5nFiIe/5tUFW3UiVkKeAl6kBlITG/LsNQN457bBtG/WgN+9v4pRj83moxW7KCsLnb+CRSpSwIuchAHtmzH11sFMvj6d2Cgfd76+lIufmcecLI24kdCjgBc5SWbGud2TmX7PcB4Z34f9R4q4bvIirnlxgeafl5Cik6wip6iwpJTXFmzjma82sO9IEed0a8EvR3YhLaWJ16VJHaBRNCJBcKSwhJe/2cLzszaRd6yY83sm84uRXejWMt7r0iSCKeBFguhQQTEvzd3M5DmbOVxUwoW9W3P3OZ3onNzY69IkAingRTxw8GgRk2ZvYsr8LRwrLuWCXq2465xOatFLrVLAi3ho3+FCJs/dzMvzt3CkqJQxaS2565zO9GitoJdTp4AXCQEHjxbx0tzN/HPeFvILSxjZI5m7zulE7zYJXpcmYUwBLxJC8o4VM2XeFibP3cShghKGdUrk9rM7MrhDc8onYRWpOQW8SAjKLyjmtYXbeHHOZvYeLqRv2wTuOLsT53Zrgc+noJeaUcCLhLCC4lKmZu7g+Vkb2XHgGF2SG3HriI6M7dOamChdiyjVU8CLhIGS0jI+WrGbZ7/ewPrsw7RqEseNQ1O5amA7GsfFeF2ehCgFvEgYKStzzFqfy6TZm/hm0z4a14vmqjPacePQVFo1qe91eRJiFPAiYWrFjoO8MGczH6/Yhc+MsX1ac/Ow0zQNgnxHAS8S5rbvP8rkuZuZmrGdI0WlDExtxk3DUhnZoyVROiFbpyngRSJE3rFipmZs55/ztrDz4DHaNK3PDUNSufz0tsSrn75OUsCLRJiS0jJmrs3mpblbWLRlPw1jo7i0fwoTBqfSRXPe1CkKeJEItmpnHv+ct4UPV+yiqKSMwR2aM2Fwe0b2SCZawywjngJepA7Yf6SItxZv59UFW9l58Bgt4+O45ox2XDGwLS0ax3ldngSIAl6kDiktc3y5Lod/fbOFOVl7ifYZI3skc80Z7RnSsbmuko0w1QV8dLCLEZHAivIH+sgeyWzee4Q3Fm1jasZ2Plm1h/bNG3DVwHaMH9CG5o3qeV2qBJha8CJ1QEFxKZ+t3sNrC7axaMt+YqN8jOyZzBXpbRnWKVGt+jCmLhoR+U5Wdj6vL9rGe0t3cvBoMSkJ9Rmf3obx6W1JSdCVsuFGAS8i31NQXMqMNdm8nbGdOVl7MYPhnZO4PL0N53VPJi4myusSpQYU8CJSre37jzI1cwfTMrazK6+A+LhoxvZpzbgBbejXNkHz1IcwBbyI1EhpmeObjfuYlrmdT1fvoaC4jA5JDRnXvw2X9U/RZGchSAEvIictv6CY6St3My1zB4u3HMAMBndoziV9Uxjdq6WmRggRCngROSVb9h7h/WU7eX/pTrbsO0pstI+R3ZO5pF8KI7okERutK2a9ooAXkVrhnGP5jjzeX7qTD5fvYt+RIhIaxDAmrSVje7fmjA7NNbtlkCngRaTWFZeWMXfDXv69dCcz1mRzpKiUxEb1uLB3K8b2aUW/tk01vj4IFPAiElDHikr56tscPly+iy/W5VBUUkZKQn3GpLXkgt6tNBIngBTwIhI0+QXFzFybzYfLdzMnK5fiUkfrJnGM6dWKC3q1VMu+lingRcQTeceK+WJtNtNX7mb2+r0UlZbRMj6O0WktGdUzmYGpzTSl8SlSwIuI5/ILivlibQ4fr9zN7PW5FJaU0bRBDOd1T+b8ni0Z1jlRV8/+CJ4EvJm9BFwI5Djn0mqyjQJepG44WlTCrG9z+Wz1Hr5Yl0N+QQkNYqM4q2sS53VP5uyuLWjaMNbrMsOCVwF/JnAY+JcCXkSqUlRSxjeb9vHZ6j3MXJNNTn4hUT4jvX1TRvZI5rzuyaQmNvS6zJDlWReNmaUCHyngRaQmysocK3fmMXNtNjPWZLNuTz4AHZMacq6/ZZ+e2pQY9dt/J6QD3swmAhMB2rVrN2Dr1q0Bq0dEwsv2/Uf5Ym02M9fmsHDzPopLHY3jojmzcxJnd2vBWV2TSKzjNy4J6YCvSC14EanK4cIS5mbt5at1OXz1bQ45+YWYQa+UJozoksSILkn0bZtQ50blKOBFJKKUlTnW7D7El+tymL0+lyXbDlDmID4umuGdkzizSyLDOyfRug7cwET3ZBWRiOLzGWkpTUhLacLd53Ym72gx8zbuZda3ucxan8vHK3cD5X33wzsnMaxTIoM6NqdRvboVeYEcRfMGcBaQCGQDf3DOTa5uG7XgReRUOef4NjufuVl7mZO1l4Wb91FQXEa0z+jfrilDOyUytFNz+rRNiIiTtbrQSUTqrILiUpZsPcCcDXuZk5XL6l2HcA4axEYx8LRmDOnYnCEdE+nRKj4sp1BQwIuI+B08WsSCTfuYt2Ef8zfuZWPuEQASGsQwMLUZgzo0Z1CH5nRr2TgsAl998CIifgkNYhmd1orRaa0A2JNXwPyNe1mwaR8LNu3n8zXZ/vX+N/AHntaM7q3iw26ue7XgRUQq2HnwGAs37fsu8LftPwpA43rR9G/flIGnNWPgac3o3aYJ9aK9nztHXTQiIj/SroPHWLxlP4s2lz+ycg4DEBvto0+bJgxo34z09k0Z0L6pJ/PnKOBFRGrJ/iNFLN6yn8Wb95Ox9QCrd+VRXFqeox2TGpLevhn92yfQv11TOiY1Cng/vgJeRCRACopLWb79IBlbD5Dpf+QdKwbKL7zq264p/dom0L99U/q2TaBJ/Zha3b9OsoqIBEhcTBRndGjOGR2aA+VX2W7ae4Ql2w6wdNtBlm47wJNfZnG8Ld0hqSF92yTQt10Cfdok0L1VPLHRgRmPrxa8iEiA5RcUs3x7Hsu2H2DZ9jyWbT/I3sOFAMRG+ejbNoE3Jw76Ud05asGLiHiocVwMwzonMqxzIlB+te2uvAKWbz/Isu0HyS8oDkhfvQJeRCTIzIyUhPqkJNTngl6tAraf8J+IQUREKqWAFxGJUAp4EZEIpYAXEYlQCngRkQilgBcRiVAKeBGRCKWAFxGJUCE1VYGZ5QJbf+TmicDeWiyntqiuk6O6To7qOjmRWFd751xSZQtCKuBPhZllVDUfg5dU18lRXSdHdZ2culaXumhERCKUAl5EJEJFUsBP8rqAKqiuk6O6To7qOjl1qq6I6YMXEZH/K5Ja8CIiUoECXkQkQoVdwJvZaDP71sw2mNkDlSyvZ2Zv+ZcvNLPUINTU1sy+MrO1ZrbazO6pZJ2zzCzPzJb5H78PdF3+/W4xs5X+fX7vfohW7kn/8VphZv2DUFPXCsdhmZkdMrOfn7BOUI6Xmb1kZjlmtqrCa83MbIaZZfl/Nq1i2+v962SZ2fVBqOthM1vnf5/eM7OEKrat9j0PQF0PmtnOCu/VBVVsW+13NwB1vVWhpi1mtqyKbQN5vCrNhqB9xpxzYfMAooCNQAcgFlgO9DhhnduBf/ifXwm8FYS6WgH9/c8bA+srqess4CMPjtkWILGa5RcAnwAGDAIWevCe7qH8Yo2gHy/gTKA/sKrCa38DHvA/fwB4qJLtmgGb/D+b+p83DXBdo4Bo//OHKqurJu95AOp6ELivBu9ztd/d2q7rhOWPAL/34HhVmg3B+oyFWwt+ILDBObfJOVcEvAlcfMI6FwMv+59PA841s9q/2WEFzrndzrkl/uf5wFogJZD7rEUXA/9y5RYACWYWuHuIfd+5wEbn3I+9gvmUOOdmA/tPeLniZ+hl4JJKNj0fmOGc2++cOwDMAEYHsi7n3OfOuRL/rwuANrW1v1Opq4Zq8t0NSF3+7//lwBu1tb+aqiYbgvIZC7eATwG2V/h9B98P0u/W8X8Z8oDmQakO8HcJ9QMWVrJ4sJktN7NPzKxnkEpywOdmlmlmEytZXpNjGkhXUvUXz4vjBZDsnNsN5V9QoEUl63h93G6i/C+vyvzQex4Id/q7jl6qorvBy+M1HMh2zmVVsTwox+uEbAjKZyzcAr6ylviJ4zxrsk5AmFkj4B3g5865QycsXkJ5N0Qf4Cng/WDUBAx1zvUHxgB3mNmZJyz38njFAhcBUytZ7NXxqikvj9t/AiXAa1Ws8kPveW17DugI9AV2U94dciLPjhdwFdW33gN+vH4gG06tu9QAAANVSURBVKrcrJLXTuqYhVvA7wDaVvi9DbCrqnXMLBpowo/7k/KkmFkM5W/ga865d09c7pw75Jw77H8+HYgxs8RA1+Wc2+X/mQO8R/mfyhXV5JgGyhhgiXMu+8QFXh0vv+zj3VT+nzmVrOPJcfOfaLsQuMb5O2pPVIP3vFY557Kdc6XOuTLghSr259XxigYuA96qap1AH68qsiEon7FwC/jFQGczO83f+rsS+OCEdT4Ajp9t/inwZVVfhNri7+ObDKx1zj1axTotj58LMLOBlB/7fQGuq6GZNT7+nPKTdKtOWO0DYIKVGwTkHf/TMQiqbFl5cbwqqPgZuh74dyXrfAaMMrOm/i6JUf7XAsbMRgO/Bi5yzh2tYp2avOe1XVfFczaXVrG/mnx3A+E8YJ1zbkdlCwN9vKrJhuB8xgJx5jiQD8pHfayn/Iz8f/pf+xPlH3qAOMr/5N8ALAI6BKGmYZT/6bQCWOZ/XADcCtzqX+dOYDXlowcWAEOCUFcH//6W+/d9/HhVrMuAZ/zHcyWQHqT3sQHlgd2kwmtBP16U/wOzGyimvMV0M+XnbL4Asvw/m/nXTQderLDtTf7P2QbgxiDUtYHyPtnjn7Hjo8VaA9Ore88DXNcr/s/OCsqDq9WJdfl//953N5B1+V+fcvwzVWHdYB6vqrIhKJ8xTVUgIhKhwq2LRkREakgBLyISoRTwIiIRSgEvIhKhFPAiIhFKAS91lpklmNnt/uetzWya1zWJ1CYNk5Q6yz83yEfOuTSPSxEJiGivCxDx0F+Bjv55wrOA7s65NDO7gfLZ/aKANMrnVokFrgMKgQucc/vNrCPlF4klAUeBnznn1gX/f0OkcuqikbrsAcqnKu4L/OqEZWnA1ZTPS/Jn4Khzrh/wDTDBv84k4C7n3ADgPuDZoFQtUkNqwYtU7itXPn93vpnlAR/6X18J9PbPDjgEmFrhdgP1gl+mSNUU8CKVK6zwvKzC72WUf298wEF/618kJKmLRuqyfMpvo3bSXPmc3pvNbDx8d2/bPrVZnMipUsBLneWc2wfM89+o+eEf8Z+4BrjZzI7PRFhrt6ATqQ0aJikiEqHUghcRiVAKeBGRCKWAFxGJUAp4EZEIpYAXEYlQCngRkQilgBcRiVD/H262Dre49HjqAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from gekko import GEKKO\n",
    "\n",
    "m = GEKKO(remote=False)    # GEKKO model\n",
    "m.time = np.linspace(0,20) # time points\n",
    "y = m.Var(5.0); k = 10.0   # GEKKO variable and constant\n",
    "m.Equation(k*y.dt()+y==0)  # GEKKO equation\n",
    "\n",
    "m.options.IMODE = 4        # Dynamic simulation\n",
    "m.solve(disp=False)        # Solve\n",
    "\n",
    "plt.plot(m.time,y)\n",
    "plt.xlabel('time'); plt.ylabel('y(t)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![expert](https://apmonitor.com/che263/uploads/Begin_Python/expert.png)\n",
    "\n",
    "### Differential Equation Activity\n",
    "\n",
    "Solve the differential equation with initial condition $y(0) = 10$:\n",
    "\n",
    "$ k \\, \\frac{dy}{dt} = -y$\n",
    "\n",
    "Compare the five solutions of `y` from time `0` to `20` with `k=[1,2,5,10,20]`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![idea](https://apmonitor.com/che263/uploads/Begin_Python/idea.png)\n",
    "\n",
    "### Symbolic Solution\n",
    "\n",
    "Compact differential equation problems may have an analytic solution that can be expressed symbolically. A symbolic math package in Python is `sympy`. Sympy determines the analytic solution as $y(x)=C_1 \\, \\exp{\\left(-\\frac{x}{k}\\right)}$. With the initial condition $y(0)=5$, the constant $C_1$ is equal to 5."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import display\n",
    "import sympy as sym\n",
    "from sympy.abc import x, k\n",
    "y = sym.Function('y')\n",
    "ans = sym.dsolve(sym.Derivative(y(x), x) + y(x)/k, y(x))\n",
    "display(ans)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![idea](https://apmonitor.com/che263/uploads/Begin_Python/idea.png)\n",
    "\n",
    "### Solve Differential Equation with Input `u`\n",
    "\n",
    "Differential equations can also have an input (feature) that changes from an external source (exogenous input) such as actively changed by a measurement sensor, a person (manually), or selected by a computer.\n",
    "\n",
    "![expert](https://apmonitor.com/che263/uploads/Begin_Python/expert.png)\n",
    "\n",
    "Calculate the response `y(t)` when the input `u` changes from `0` to `2` at `t=5`.\n",
    "\n",
    "$2 \\frac{dy(t)}{dt} + y(t) = u(t)$\n",
    "\n",
    "The initial condition is `y(0)=1` and the solution should be calculated until `t=15`. **Tip**: The expression `y(t)` does not mean `y` multiplied by `t`. It indicates that `y` changes with time and is written as a function of time. There are additional examples for [ODEINT](https://apmonitor.com/pdc/index.php/Main/SolveDifferentialEquations) and [GEKKO](https://apmonitor.com/pdc/index.php/Main/PythonDifferentialEquations) if you need help."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### TCLab Activity\n",
    "\n",
    "![expert](https://apmonitor.com/che263/uploads/Begin_Python/expert.png)\n",
    "\n",
    "### Data Collection\n",
    "\n",
    "![connections](https://apmonitor.com/che263/uploads/Begin_Python/connections.png)\n",
    "\n",
    "Turn on heater 1 to 100% and record $T_1$ every 5 seconds for 3 minutes. The data should include a total of 37 data points for each temperature sensor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import tclab\n",
    "import time\n",
    "# collect data for 3 minutes, every 5 sec\n",
    "n = 37\n",
    "tm = np.linspace(0,180,n)\n",
    "t1s = np.empty(n); t2s = np.empty(n)\n",
    "with tclab.TCLab() as lab:\n",
    "    lab.Q1(100); lab.Q2(0)\n",
    "    print('Time T1  T2')\n",
    "    for i in range(n):\n",
    "        t1s[i] = lab.T1; t2s[i] = lab.T2\n",
    "        print(tm[i],t1s[i],t2s[i])\n",
    "        time.sleep(5.0)\n",
    "# put into dataframe\n",
    "data = pd.DataFrame(np.column_stack((tm,t1s,t2s)),\\\n",
    "                    columns=['Time','T1','T2'])\n",
    "data.to_csv('11-data.csv',index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![expert](https://apmonitor.com/che263/uploads/Begin_Python/expert.png)\n",
    "\n",
    "### Solve Differential Equation\n",
    "\n",
    "Use the parameters `a`, `b`, and `c` from [10. Solve Equations](https://github.com/APMonitor/data_science/blob/master/10.%20Solve_Equations.ipynb) or use the following default values.\n",
    "\n",
    "| Parameter | Value |\n",
    "|------|------|\n",
    "| a    | 78.6 |\n",
    "| b    | -50.3 |\n",
    "| c    | -0.003677 |\n",
    "\n",
    "Solve the ordinary differential equation (ODE) with these values.\n",
    "\n",
    "$\\frac{dT_1}{dt} = c (T_1-a)$\n",
    "\n",
    "The initial condition for $T_1$ is $a + b$. Show the ODE solution at time points between the initial time `0` and the final time `180` sec. Plot the measured $T_1$ on the same plot as the ODE predicted temperature. Add appropriate labels to the plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}