LKM.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The lumped kinetic model (LKM)\n",
    "\n",
    "Marek Szymanski"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.integrate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sparsegrad.base import *\n",
    "from sparsegrad.forward import nvalue\n",
    "\n",
    "N=1000\n",
    "\n",
    "def f_langmur(c, params):\n",
    "    return params.a*c/(1.+params.b*c)\n",
    "\n",
    "def internal_ode(c, q, c_dz, c_ddz, params, f):\n",
    "    qstar = f(c, params)\n",
    "    d=qstar-q\n",
    "    c_dt = -params.u*c_dz + params.D * c_ddz - params.k/params.epsilon*d\n",
    "    q_dt = params.k/(1-params.epsilon) * d\n",
    "    return c_dt, q_dt\n",
    "\n",
    "def cinj_prescribed(t, params):\n",
    "    return where(t<=params.tinj, params.cinj, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ode(t, x, params, f):\n",
    "    L=params.L\n",
    "    h=L/(N-1)\n",
    "    c_int = x[0:N-2]\n",
    "    q_int = x[N-2:2*N-4]\n",
    "    # BC0: Solve -params.D/params.u * c_dz[0] + c[0] - cinj_prescribed(t, params) for c[0]\n",
    "    v = -params.D/params.u\n",
    "    bc0 = (cinj_prescribed(t, params) - v*c_int[0]/h)/(1-v/h)\n",
    "    # BC1: \n",
    "    bc1 = c_int[-1]\n",
    "    c = stack(bc0, c_int, bc1)\n",
    "    h = L/(N-1)\n",
    "    c_dz = (c[1:]-c[:-1]) / h\n",
    "    c_dz_sym = 0.5*(c_dz[:-1]+c_dz[1:])\n",
    "    c_ddz = (c_dz[1:]-c_dz[:-1]) / h\n",
    "    assert np.allclose(nvalue(-params.D/params.u * c_dz[0] + c[0]), cinj_prescribed(t, params))\n",
    "    assert np.allclose(nvalue(c_dz[-1]), 0)\n",
    "    return stack(*internal_ode(c_int, q_int, c_dz_sym, c_ddz, params, f))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import namedtuple\n",
    "params_type = namedtuple('params', ['L','D', 'u', 'a', 'b', 'epsilon', 'cinj', 'tinj','k'])    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "params_table = params_type(L=1e-2, \n",
    "                           epsilon=0.4, \n",
    "                           u=1e-2/60, \n",
    "                           D=2e-8/60,\n",
    "                           cinj=1e-3/1e-3, \n",
    "                           a=1, \n",
    "                           #b=0.1e-3/1e-3,\n",
    "                           b=0,\n",
    "                           tinj=3*60, \n",
    "                           k=100/60)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sparsegrad.forward import seed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "cinit = 0 \n",
    "qinit = 0\n",
    "x0 = np.zeros(2*N-4)\n",
    "x0[:N-2] = cinit\n",
    "x0[N-2:2*N-4] = qinit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pylab as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = x0\n",
    "def jac(t, x):\n",
    "    return ode(t, seed(x), params_table, f_langmur).dvalue.tocsr()\n",
    "def func(t, x):\n",
    "    return ode(t,x,params_table,f_langmur)\n",
    "result = scipy.integrate.solve_ivp(func, [0, 8*60], x0, method='Radau', jac=jac, rtol=1e-6, atol=1e-8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fa170b14278>]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAHfVJREFUeJzt3WtwXPd53/Hvs4vFZXEHAd5JkbIoybRsyzYiKXGsUSMpIzm2lGkSjzhx2jSu2YzjxLdpIrep2rrtC8eZNM2M7Clju77Ulixf0rANY7q15VpOLZmQLNsiZdoQLRGgeAGIBYj74vL0xe5CIIjLAljw7P7P7zPDmb0cYJ8hsD/89znnPMfcHRERCUsi6gJERKT0FO4iIgFSuIuIBEjhLiISIIW7iEiAFO4iIgFSuIuIBEjhLiISIIW7iEiAqqJ64fb2dt+zZ09ULy8iUpGefvrpfnfvWGm7yMJ9z549dHV1RfXyIiIVycxeKmY7tWVERAKkcBcRCZDCXUQkQAp3EZEAKdxFRAK0Yrib2afN7IKZPbfE82Zmf2Vm3Wb2IzN7Y+nLFBGR1Shm5f4Z4J5lnr8X2Jf/dxD4xPrLEhGR9VjxOHd3/46Z7Vlmk/uBz3nuen1PmlmLmW1z97MlqlECMjI5Tc/AGGcy45y7NMHk9CxTM7NMTc8yPbv4JR+XvRBkMZeJNHvl5pUPYflHC4/Ne+qVx+Z/wcLn5n3Fst9jwetc/r1sma+7fLuEweamWq7b3MC17fWL1iZSipOYdgA98+735h+7ItzN7CC51T27d+8uwUtLpfjm8+f59D/8nCdPDTCzRIiv1XLZFvolgq/ZlOY9d7yKd3TuUsjLZa7qGarufgg4BNDZ2Rn4204AxrMzfOBLz/L14+fY3Zbm4O3X8prtTexoqWNbcx111UmqkwmqkkZVwq5KQM2/KHzhpi/23IJtco/5FY9d+b2u3Kao77/C6xRuTs/Mcu7SBD/qHeKrz/TyJ1/9MU+dGuDPf+v1JBIKeMkpRbifAXbNu78z/5jEXHZ6loOf7+K73f388T038O63XEsqGf0BWvP/gFz5t6QywnFzUy2v29nCb9+6m7/8Pz/jv3zzZ+xqS/OBu6+PujQpE6V4px0G/kn+qJnbgCH12wXg0Hde4Imf9fPRf/w63nPHdWUR7KExM95/1z5+/ebtfPzb3XRfGIm6JCkTxRwK+QjwPeAGM+s1s3eZ2e+b2e/nNzkCnAK6gb8G3rNh1UrFeLF/lL/6Vje/9tptvOMXdq38BbJmZsafvm0/NVVJHn68O+pypEwUc7TMgRWed+APSlaRBOGT3z0FwENv3x9xJfHQ3lDDb75pJ1946iU+fO+NbG6qjbokiZg+J0vJXZqY4mvPnOG+129ni0LmqvmdX7yGqRnnf/5IXVFRuMsG+NrTvYxlZ/jdX9oTdSmx8qqOBm7c2sjR585FXYqUAYW7lNzXj5/jhi2N3LSjOepSYueem7Zy7KUB+kcmoy5FIqZwl5IaGp+i68UMd756c9SlxNLt13fgDt//+UDUpUjEFO5SUk/8rI/pWVe4R+S1O5qpSyUV7qJwl9L61k8u0JpOcfOu1qhLiaVUMsGbrmnlyVMXoy5FIqZwl5J6+qUMt+7dRFKnwUemc08rJ88PM5adjroUiZDCXUpmYDTLSxfHuHl3S9SlxNprtjfjDs+fHY66FImQwl1K5oe9gwDcvEvhHqXXbG8C4MTLQxFXIlFSuEvJPHt6kITldupJdLY119KaTnH85UtRlyIRUrhLyTzbM8j1Wxqpr7mqk6RlATNj//Ymnj+rcI8zhbuUzPGXh7RqLxPXdTTwQt/oZfPjJV4U7lISmdEs/SNZrt/SGHUpAlzb0cDI5DR9wzpTNa4U7lIS3X25OeLXbW6IuBKB3JwZeOXnIvGjcJeSKFwkQuFeHq7tqAfghb7RiCuRqCjcpSS6L4xQm0qwo6Uu6lIE2NpUS10qySmt3GNL4S4l0X1hhFd1NOgCzWUikTB2t6XpGRiLuhSJiMJdSqL7wohaMmVmV1sdPQPjUZchEVG4y7qNZ2c4Mzg+txNPysPO1jS9mTEdDhlTCndZt55M7qP/NZvSEVci8+1srWM0O0NmbCrqUiQCCndZt0Jfd1ebwr2cFH4e6rvHk8Jd1m0u3FsV7uWk8PPozajvHkcKd1m3nsw4dakk7Q3VUZci8+xozR2WemZQK/c4UrjLuvUMjLGztQ4zHQZZTppqq6hLJTl/SSMI4kjhLuvWkxlXv70MmRlbm2s5d2ki6lIkAgp3WbfezBi7WnVmajna0lTD+SGFexwp3GVdhiemGJ6YZrvGDpSlrU1auceVwl3W5Xw+OLY210ZciSxmS1MtFy5N6kSmGFK4y7qczX/k39aslXs52tJUS3ZmVicyxZDCXdbllXDXyr0cFT5RnVPfPXaKCnczu8fMTppZt5k9uMjzu83scTP7gZn9yMzeWvpSpRwVQmNzU03ElchitjTlwv28+u6xs2K4m1kSeBi4F9gPHDCz/Qs2+1PgMXd/A/AA8PFSFyrl6ezQBO0N1dRUJaMuRRYxt3JXuMdOMSv3W4Budz/l7lngUeD+Bds40JS/3Qy8XLoSpZydGxqfWx1K+dncmPtEpZV7/FQVsc0OoGfe/V7g1gXb/DvgG2b2h0A9cFdJqpOyd3Zogp06xr1spZIJ2huqFe4xVKodqgeAz7j7TuCtwOfN7IrvbWYHzazLzLr6+vpK9NISpQvDk2zWyr2sbW6spW9YIwjipphwPwPsmnd/Z/6x+d4FPAbg7t8DaoH2hd/I3Q+5e6e7d3Z0dKytYikb0zOzZMaydDRoZ2o529RQTf9INuoy5CorJtyPAfvMbK+ZVZPbYXp4wTangTsBzOzV5MJdS/PADYxmcYf2RoV7OWtvqOHiqFbucbNiuLv7NPBe4CjwPLmjYo6b2UfM7L78Zh8C3m1mPwQeAX7XdUpc8PpGcoHRXq9Rv+VsU301/cNaucdNMTtUcfcjwJEFjz007/YJ4M2lLU3K3cX8R32t3Mtbe2MN41MzjGWnSVcX9ZaXAOgMVVmz/sLKXT33srYp/8nqovrusaJwlzV7JdzVlilnhT++hZ+XxIPCXdasfyRLTVWChhp91C9nm/J/fHXETLwo3GXN+kcmaW+o0eX1ylxh5X5RK/dYUbjLmvWPZNWSqQBthZ77qFbucaJwlzXrH57UztQKUJtK0lhTpZ57zCjcZc0KbRkpf+2NNeq5x4zCXdZkdtYZGM3S3qi2TCXYVF+tnnvMKNxlTYbGp5iedTbVa+VeCTY1VOs495hRuMuazB3jrrNTK0JbvebLxI3CXdakTycwVZS2+hSDY1No5FN8KNxlTQof8TXutzK0pquZnnWGJ6ejLkWuEoW7rInmylSWlnTuE9bg6FTElcjVonCXNekfmSSZMJrrUlGXIkVoTed+TgNj2qkaFwp3WZP+4Syb6qtJJDR6oBK05s9SzSjcY0PhLmtycVQnMFWS1kJbRuEeGwp3WZOB0ezczBIpf4W2TEY999hQuMuaDI5N0ZJWv71SNNWmSJjaMnGicJc1GRxXuFeSRMJoSVcr3GNE4S6rNjvrDI5l5/q4Uhla0ikyY2rLxIXCXVZteHKaWUeHQVaY1nQ1Gc10jw2Fu6zaUH7116KVe0Vp1co9VhTusmqD47nVX6t67hWlNV2tQyFjROEuq5aZW7kr3CtJa712qMaJwl1WrbD6a65TW6aStKRTTEzNMp6diboUuQoU7rJqQ+O5lbvaMpWlcHSTVu/xoHCXVSuc5aijZSqLwj1eFO6yaoPjWRprqqhK6tenkmgEQbzo3SmrNjQ2RbNaMhVHkyHjReEuqzY4PqWzUytQ4egmHQ4ZDwp3WbXMWFaHQVagwh/kAbVlYqGocDeze8zspJl1m9mDS2zzDjM7YWbHzeyLpS1TysnQ2JR2plagVDJBY03V3EloEraqlTYwsyTwMHA30AscM7PD7n5i3jb7gA8Db3b3jJlt3qiCJXpqy1Su5nRqbnyEhK2YlfstQLe7n3L3LPAocP+Cbd4NPOzuGQB3v1DaMqVcFCZCqi1TmXKTIbVyj4Niwn0H0DPvfm/+sfmuB643s38wsyfN7J7FvpGZHTSzLjPr6uvrW1vFEilNhKxsrelqBse1co+DUu1QrQL2AXcAB4C/NrOWhRu5+yF373T3zo6OjhK9tFxNhY/0astUpua6FINqy8RCMeF+Btg17/7O/GPz9QKH3X3K3X8O/JRc2EtgCh/p1ZapTJoMGR/FhPsxYJ+Z7TWzauAB4PCCbf4HuVU7ZtZOrk1zqoR1SpkofKRXuFemlnSKofEpZmc96lJkg60Y7u4+DbwXOAo8Dzzm7sfN7CNmdl9+s6PARTM7ATwO/Et3v7hRRUt0BudW7mrLVKLmuhSzDsMT01GXIhtsxUMhAdz9CHBkwWMPzbvtwAfz/yRghX5ti3aoVqTCvpLB8axGSAROZ6jKqhTCXUfLVKZXRhBop2roFO6yKpoIWdkK4a5j3cOnd6isyuDYFC31WrVXqsK+kiEd6x48hbusyuBYlhZdXq9iFfaVZEa1cg+dwl1WZXB8SodBVrDCvhKdpRo+hbusyuDYlA6DrGBVyQSNtVXaoRoDCndZlVxbRiv3StaSTuks1RhQuEvRZmedIbVlKp6Gh8WDwl2KVpgIqbZMZWuuS5FRWyZ4Cncp2tzoAbVlKlpLupohtWWCp3CXos2NHlBbpqK1plNqy8SAwl2KpomQYWipy02GnNFkyKAp3KVomggZhpZ0Ne4wPKHVe8gU7lI0TYQMg4aHxYPCXYqmiZBh0PCweFC4S9EyY1kaazURstK1zM1018o9ZHqXStF0AlMYCm01naUaNoW7FE0TIcMwdzUm9dyDpnCXomXGtHIPQVOddqjGgcJdipZry2jlXumSCaOptkptmcAp3KVomggZjtZ6DQ8LncJdilKYCNmqtkwQWjQ8LHgKdynK8ERuImSz2jJBaNbwsOAp3KUog+OaCBkSDQ8Ln8JdipLRRMigtNSldJHswCncpSgaGhaWlnQ1lyamNRkyYAp3KcqQxv0GpfBzvKTWTLAU7lIUTYQMi4aHhU/hLkUphIAmQoZBw8PCp3CXogyOTWkiZEA0PCx8Rb1TzeweMztpZt1m9uAy2/2GmbmZdZauRCkHmggZlhYNDwveiuFuZkngYeBeYD9wwMz2L7JdI/A+4KlSFynRy4xl56YJSuVr1dWYglfMyv0WoNvdT7l7FngUuH+R7f4D8FFgooT1SZkYHJtSvz0gjbUpzNSWCVkx4b4D6Jl3vzf/2BwzeyOwy93/roS1SRnRRMiwJBNGc53OUg3ZuveOmVkC+AvgQ0Vse9DMusysq6+vb70vLVdRRhMhg6PhYWErJtzPALvm3d+Zf6ygEbgJ+LaZvQjcBhxebKequx9y90537+zo6Fh71XJVaSJkmJrT1WrLBKyYcD8G7DOzvWZWDTwAHC486e5D7t7u7nvcfQ/wJHCfu3dtSMVy1Q1PTOOaCBmc1nRq7sxjCc+K4e7u08B7gaPA88Bj7n7czD5iZvdtdIESvcIJTGrLhCXXltHKPVRVxWzk7keAIwsee2iJbe9Yf1lSTgo73VrrFe4haUlX61DIgOl0Q1nR4NzoAbVlQtKSTjE8Mc30zGzUpcgGULjLigY1yz1IhTab+u5hUrjLigbyF3Vo0w7VoLTWa3hYyBTusqLBsSxm0KQdqkFp1vCwoCncZUUD+ROYkgmLuhQpIQ0PC5vCXVaUGZua+wgv4dDwsLAp3GVFmVFNhAxRS/7oJx3rHiaFu6woMzalcA9QY20VCdPRMqFSuMuKcit37UwNTaIwGVJtmSAp3GVZ7k5mLEubeu5BaklXqy0TKIW7LGt8aobJ6VnNcg9Ui4aHBUvhLssqzPtu01yZIGl4WLgU7rKsTP7sVK3cw6ThYeFSuMuyCqs69dzD1JJOMaRwD5LCXZZVmCujo2XC1FJXzfDkNFOaDBkchbssKzMX7lq5h6iwL0V99/Ao3GVZhR2qzRoaFqT2hhoALo4o3EOjcJdlZcayNNelqErqVyVEm/Lh3j8yGXElUmp6x8qycqMHtGoP1aaGXLtNK/fwKNxlWZnRrCZCBqy9Xiv3UCncZVmZMU2EDFlTXRWppNGvlXtwFO6yLI37DZuZsam+hotauQdH4S5Lcnf6R7K0NyrcQ7apoZqLo1q5h0bhLku6NDFNdmZ2ri8rYWpv0Mo9RAp3WVLhDa+Ve9g2NVSr5x4ghbssqfCGL5zoImFqb6ihf2QSd4+6FCkhhbssqbBy36S2TNDaG6qZnJ5lNDsTdSlSQgp3WVK/2jKxUPjjrb57WBTusqT+kSxm0KZDIYNWOEtVJzKFReEuS+ofmaQ1Xa25MoEr7FPpG1a4h0TvWllS/8gk7Q1atYdua3MtAOcvKdxDUlS4m9k9ZnbSzLrN7MFFnv+gmZ0wsx+Z2TfN7JrSlypX28WRrHamxkBbuppU0jh3aSLqUqSEVgx3M0sCDwP3AvuBA2a2f8FmPwA63f11wFeAPyt1oXL19Y9M0t6ocA9dImFsbqzlvMI9KMWs3G8But39lLtngUeB++dv4O6Pu/tY/u6TwM7SlilRyK3c1ZaJg81NNQr3wBQT7juAnnn3e/OPLeVdwN8v9oSZHTSzLjPr6uvrK75KueompmYYnpymQyv3WNjaVMu5IYV7SEq6Q9XM3gl0Ah9b7Hl3P+Tune7e2dHRUcqXlhLrnzuBSSv3ONjSVMsF7VANSjHhfgbYNe/+zvxjlzGzu4B/Ddzn7votqXCFIye25I+kkLBtaapleHKa0cnpqEuREikm3I8B+8xsr5lVAw8Ah+dvYGZvAP4ruWC/UPoy5WorfETfpnCPhS1Nufab+u7hWDHc3X0aeC9wFHgeeMzdj5vZR8zsvvxmHwMagC+b2bNmdniJbycV4uzQOJDrxUr4Cj9nHQ4ZjqpiNnL3I8CRBY89NO/2XSWuSyJ2/tIEtakEzXW6OHYcbG4qnMikcA+FzlCVRZ0dmmBbcx1mFnUpchXoLNXwKNxlUeeGJub6sBK+hpoqGmureHlwPOpSpEQU7rKoc5dyK3eJj52taXoGxlbeUCqCwl2uMDvrnL80wRbtTI2V3W119GS0cg+Fwl2uMDCWZWrGdRhkzOxqTdObGdPl9gKhcJcrFI5x36pwj5VdbWkmpmbp00U7gqBwlyucLYS72jKxsqstt4+lZ0CtmRAo3OUKhROY1JaJl12taQB6M9qpGgKFu1zh9MUxalMJTYSMmZ35cNcRM2FQuMsVXhoYY3dbWicwxUxddZL2hhq1ZQKhcJcrnL6YC3eJn11tdbw0MBp1GVICCne5jLtzemCM3W31UZciEXhVRwPdFxTuIVC4y2X6RiYZn5rhmk1aucfR9Vsa6B+ZJDOajboUWSeFu1zm9MXczjS1ZeJp35ZGAH56fjjiSmS9FO5ymZcK4a6VeyxdXwj3CyMRVyLrpXCXy7zQN0JVwuaOeZZ42d5cS0NNFd1auVc8hbtc5uS5Ya7tqKe6Sr8acWRmXLe5gZ+e18q90ukdLJc5eX6YG7Y2RV2GROj6LQ3quQdA4S5zRian6c2Mc+PWxqhLkQjt39bExdHs3BgKqUwKd5lz8lxutVbYqSbxdPPuVgB+cHow4kpkPRTuMqcQ7lq5x9v+bU1UVyV4tkfhXskU7jLn2Z4MLekUO1p0eb04q65KcNP2Jp55KRN1KbIOCneZ0/Vihs5r2kgkNDAs7m7Zu4lnewYZmZyOuhRZI4W7ANA/Msmp/lF+YU9r1KVIGXjLvnamZ52nTl2MuhRZI4W7AND14gAAnXvaIq5EysGbrmmlNpXg//60L+pSZI0U7gLAEz/rpy6V5KYdOsZdoDaV5PZ9HXz9uXPMzOqC2ZVI4S7MzDpHj5/jV27cTE1VMupypEy8/fXbuTA8yVM/V2umEinchWMvDtA/kuWtr90WdSlSRu569RYaa6r44lOnoy5F1kDhLjzW1UNdKsk/urEj6lKkjNRVJzlw627+/rlzuq5qBVK4x9zpi2P87bMvc+CW3aSrq6IuR8rM7715L9XJBP/x705EXYqsUlHhbmb3mNlJM+s2swcXeb7GzL6Uf/4pM9tT6kKl9Nyd/3TkBEkzDt5+bdTlSBna2lzLH925j6PHz/Plrp6oy5FVWDHczSwJPAzcC+wHDpjZ/gWbvQvIuPt1wH8GPlrqQqW03J2Pf/sFjh4/z4d+9Xq2NtdGXZKUqX/+lr388nXtfPhrP+ZLx07jrqNnKkExK/dbgG53P+XuWeBR4P4F29wPfDZ/+yvAnWam0xzL0NDYFI//5ALv/lwXHzt6kre/fjvvfotW7bK0VDLBx9/5Rm67dhN/8tUf8xuf+H98/nsvcvLcMBNTM1GXJ0sopsm6A5j/eawXuHWpbdx92syGgE1AfymKnO+xYz0ceuIU+deae/yKtYQvenPJr/HLtvdFH1/s/lq/75XPLfUai9eysIxiXn96ZpbRbO7N2FhbxYP33sjBt1yrcQOyoqbaFJ/9vVt45Pun+eQTp/g3f3t87rm6VJL6miS1qSQJM8zAgET+hpG7CIh+y17xR3fu4+2v376hr3FV96CZ2UHgIMDu3bvX9D1a66u5Yf5IWlv0ZuH1Fn3OlviapbZf+I3n/5ou/b1W3v7K+2v/vst/Te5ewoztLbXs29LIbde26Zh2WZVkwnjnbdfw27fu5oW+UU6cvcTpi6NcmphmdHKa8ewMDsy6487cbfzKhU3cNdelNvw1ign3M8Cuefd35h9bbJteM6sCmoErznxw90PAIYDOzs41/bTv3r+Fu/dvWcuXikgJFC7Fd93mhqhLkWUU03M/Buwzs71mVg08ABxesM1h4J/mb/8m8C3XXhcRkcisuHLP99DfCxwFksCn3f24mX0E6HL3w8CngM+bWTcwQO4PgIiIRKSonru7HwGOLHjsoXm3J4DfKm1pIiKyVjpDVUQkQAp3EZEAKdxFRAKkcBcRCZDCXUQkQBbV4ehm1ge8tMYvb2cDRhuUgOpaHdW1euVam+panfXUdY27r3jxhcjCfT3MrMvdO6OuYyHVtTqqa/XKtTbVtTpXoy61ZUREAqRwFxEJUKWG+6GoC1iC6lod1bV65Vqb6lqdDa+rInvuIiKyvEpduYuIyDIqLtxXulh3FMzs02Z2wcyei7qW+cxsl5k9bmYnzOy4mb0v6poAzKzWzL5vZj/M1/Xvo65pPjNLmtkPzOx/RV1LgZm9aGY/NrNnzawr6noKzKzFzL5iZj8xs+fN7BfLoKYb8v9PhX+XzOz9UdcFYGYfyP/OP2dmj5jZhl28uKLaMvmLdf8UuJvc5f6OAQfc/UTEdd0OjACfc/eboqxlPjPbBmxz92fMrBF4Gvj1Mvj/MqDe3UfMLAV8F3ifuz8ZZV0FZvZBoBNocve3RV0P5MId6HT3sjpm28w+Czzh7p/MX+8h7e6DUddVkM+MM8Ct7r7W82pKVcsOcr/r+9193MweA464+2c24vUqbeVezMW6rzp3/w65OfZlxd3Puvsz+dvDwPPkrncbKc8Zyd9N5f+VxSrDzHYCvwZ8Mupayp2ZNQO3k7ueA+6eLadgz7sTeCHqYJ+nCqjLX7EuDby8US9UaeG+2MW6Iw+rSmBme4A3AE9FW0lOvvXxLHAB+N/uXhZ1AX8J/DEwG3UhCzjwDTN7On8t4nKwF+gD/lu+jfVJM6uPuqgFHgAeiboIAHc/A/w5cBo4Cwy5+zc26vUqLdxlDcysAfgq8H53vxR1PQDuPuPuN5O7Ju8tZhZ5O8vM3gZccPeno65lEb/s7m8E7gX+IN8KjFoV8EbgE+7+BmAUKIv9YAD5NtF9wJejrgXAzFrJdRr2AtuBejN750a9XqWFezEX65Z58j3trwJfcPevRV3PQvmP8Y8D90RdC/Bm4L58f/tR4FfM7L9HW1JOftWHu18A/oZcizJqvUDvvE9dXyEX9uXiXuAZdz8fdSF5dwE/d/c+d58Cvgb80ka9WKWFezEX65a8/I7LTwHPu/tfRF1PgZl1mFlL/nYduR3kP4m2KnD3D7v7TnffQ+5361vuvmErq2KZWX1+hzj5tsevApEfmeXu54AeM7sh/9CdQKQ76xc4QJm0ZPJOA7eZWTr/3ryT3H6wDVHUNVTLxVIX6464LMzsEeAOoN3MeoF/6+6firYqILcS/R3gx/n+NsC/yl8TN0rbgM/mj2RIAI+5e9kcdliGtgB/k8sDqoAvuvvXoy1pzh8CX8gvtk4B/yzieoC5P4J3A/8i6loK3P0pM/sK8AwwDfyADTxTtaIOhRQRkeJUWltGRESKoHAXEQmQwl1EJEAKdxGRACncRUQCpHAXEQmQwl1EJEAKdxGRAP1/PtruTL6ddVsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(result.t/60, result.y[N-3,:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "88769499279848a784ecada0e63acc24",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=0, description='i', max=917), Output()), _dom_classes=('widget-interact'…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.plot(i=IntSlider(value=0, description='i', max=917))>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import ipywidgets\n",
    "\n",
    "def plot(i=ipywidgets.IntSlider(min=0,max=len(result.t)-1)):\n",
    "    print('time = ', result.t[i])\n",
    "    plt.plot(result.y[:N-2,i])\n",
    "    plt.show()\n",
    "\n",
    "ipywidgets.interact(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.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}