{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Ordinary Differential Equations\n", "\n", "Differential equations are systems whose variables change with time. Mathematically, differentials are functions that contain derivatives of itself. When these systems are integrated, they provide analytical functions that are dependent with time. Integration of differentials can be performed by calculus, or numerically. \n", "\n", "INSERT image here)\n", "\n", "Ordinary differential equations (ODES) are equations that have a single dependent variables. Partial differential equations are equations that are dependent on 2 or more variables. Furthermore, the order of a differential equation is characterized by the highest derivative of an independent variable. \n", "\n", "Biological cell growth is a first order differential that describes the exponential increase of biomass concentration (X) accumulation over time (assuming $ \\mu _g $ is constant ). \n", "\n", "\n", "$$ \\frac{dX}{dt} = \\mu_g X $$\n", "\n", "The Navier-Stokes equations in contrast, are 2nd order partial differentials describing particle position of a fluid as a function of the derivatives of speed, sheer, time and, pressure. \n", "\n", "$$ \\rho \\big( \\frac{du}{dt} +u \\frac{du}{dx} + v \\frac{dv}{dy} + w \\frac{dw}{dz} \\big) = - \\frac{dP}{dx} + \\mu \\big( \\frac{d^2u}{dx^2} + \\frac{d^2v}{dy^2} + \\frac{d^2w}{dz^2} \\big) $$ \n", "\n", "The scope of this tutorial focuses on first order, ordinary differential equations. \n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "When ODE's are solved, some have an analytical solution. For example, solving the Biomass growth equation gives this analytical solution.\n", "\n", "$$ X = X_0 e^{\\mu _g (t-t_0)} $$\n", "\n", "Where $ X_0 $ and $ t_0 $ are initial conditions of the cell concentration and lag time. \n", "\n", "If you have worked with numerical methods in Excel, you may have come across solving differentials by using this form:\n", "\n", "\n", "$$ \\frac{dX}{dt} = \\mu_g X$$\n", "\n", "$$ \\frac{X_{i+1}-X_i}{t_{i+1} - t{i}} = \\mu_g X_i $$\n", "\n", "$$ X_{i+1} = X_{i}+ \\mu _g X_i \\Delta t $$ \n", "\n", "\n", "This form of discretizing is the basis of numerical methods for ODE's and are called *single step* or *Runge Kutta* methods.\n", "\n", "$$ Value_{new} = Value_{old} + slope*stepsize $$ \n", "\n", "$$ y_{i+1} = y_i + \\phi h$$\n", "\n", "Well studied methods have minimized the error in the $ \\phi$ term to provide highly accurate ODE integrators. The following sections show you the difference between an analytical solution and a basic ODE solver. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Comparing Analytical solutions with Numerical Solver\n", "\n", "Using the Cellular biomass equation as a test, we observe the efficacy of the simplest Runge kutta method.\n", "$$ X_{i+1} = X_{i}+ \\mu _g X_i \\Delta t $$ \n", "and compare it with the known analytical solution. \n", "$$ X = X_0 e^{\\mu _g (t-t_0)} $$\n", "\n", "As we know that the smaller the stepsize, the more accurate the solver, we observe the difference between a solution presented with a stepsize of 10 and 100. The following code illustrates this difference and is very useful for simple implementations such as use in Excel." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd4VFX6wPHvm0ISWui9BAELkIIG6SgiggVQVCyIYCG6\niooFRXGVVdlFZRWR/a0LqytqhCgqoIg0QRQEBEQQEBEltCT0QICQdn5/3Js4SWYmkzKZSfJ+nifP\nzC1z77mT5Lz3nirGGJRSSlVdAb5OgFJKKd/SQKCUUlWcBgKllKriNBAopVQVp4FAKaWqOA0ESilV\nxWkgUGVKRIaLyBJfp6O4RGSiiHxQws/2FpGdZZCGPSJyZWmPU+CYRkTalfCzZXJdyv9pIPBTInK7\niGwQkTQRSRKRRSLSy9fpKooxJt4Yc5U3zyEiK0XkuIiEePM8bs6fL3M1xnxrjLnAF2kpS+V1XSIS\nYZ8rqISfv1xEcuz/jVMislNE7iqwT75rEZEn7P+jjqVNf2WkgcAPichjwFTg70BjoBXwL2CwL9NV\nlJL+YxfzHBFAb8Dg59+H8qqDxpiaQG3gUWCmiDgNWiLyLDAWuMwYs60c01hhaCDwMyISDrwAPGiM\n+dQYc9oYk2mM+cIY86S9T4iITBWRg/bP1Ny7Y/tuab+IPCkih+y7oOtF5BoR+VVEjonIMw7nmygi\nc0Ukwb672iQi0Q7bx4vIbnvbdhG5wWHbKBFZLSKvi8hRYKK97juHfYyI3C8iu0TkhIj8S0TE3hYo\nIv8UkSMi8oeIjPHgTvFOYC3wLjCywHf3rn38hXZ614lIW4ftb4jIPhE5KSIbRaS3i9/BQhF5qMC6\nLSJyg4isslf9ZN+R3pL7nTvs21JEPhWRwyJyVESm2+vbisjX9rojIhIvInXcXKvj+a+xv/9TInJA\nRJ5w2DZaRH6zf7cLRKSZi2OsFJF7HZbzflceXtdF9jFOiMg2ERnssM3td19A7rlO2OfqLiIBIvKs\niCTaf7fv2f8LbhnLl8AxIMrJNb8E3Av0Mcb8WtTxqixjjP740Q8wEMgCgtzs8wJWZtgIaAisAV60\nt11uf/45IBgYDRwGPgRqAR2Bs0Abe/+JQCZwk73/E8AfQLC9/WagGdZNwy3AaaCpvW2Ufa6HgCAg\nzF73nUNaDfAFUAfryeYwMNDedj+wHWgB1AWW2fu7u/bfgAeAS+x0N3bY9i5wFLjUTk88MMdh+x1A\nfXvb40AyEOrwPXxgvx8GrHP4XLR93GoO19TOYfvlwH77fSDwE/A6UAMIBXrZ29oB/YEQ+/e2Cpjq\ncJw9wJUurjsJ6G2/rwtcbL+/AjgCXGwf901gVYHvv539fiVwr8M2Z78rV9cVbH/3zwDV7POeAi7w\n5LsvcC0RBX/PwN328c8DagKfAu+7+LxjugKwngxzgM4FrmUusAto5ev/a3//8XkC9KfALwSGA8lF\n7LMbuMZheQCwx35/OVZGH2gv17L/Kbo67L8RuN5+PxFY67AtwDHTcXLuzcAQ+/0oYG+B7c4yl14O\nyx8B4+33XwP3OWy7smAGUeDYvbAy/wb28i/Aow7b3wX+67B8DfCLm+/xOBDt8D3kBoJQe1t7e3kK\n8H8FrslVhtkdK9i5DGYOn7se+NFheQ+uA8Fe4D6gdoH1bwOvOCzXtL+jiIJppXSBoDdW4Axw2D4b\nmFjc7x7ngWA58IDD8gX2dRT6Hu105QAngHNANjC2wD4GOAm86a3/1cr0o0VD/uco0KCI4pFmQKLD\ncqK9Lu8Yxphs+/1Z+zXFYftZrAwj177cN8aYHGB/7vFE5E4R2WwXB5wAOgENnH3WjWSH92cczt2s\nwOeLOtZIYIkx5oi9/CEFiofcnCu3wnCHiKTa1xJO/msBwBiTDiQAd4hIAHAb8H4RacvVEkg0xmQV\n3CAijUVkjl20cxL4wNn5XbgRK3NNFJFvRKS7vT7f34IxJg3rb6i5h8f1VDNgn/33kSuxwHlcfvce\nHr/g33QQVh2ZMweNMXWw6gimYT2hFHQrcJOI/K0Y6aiSNBD4n++x7nKud7PPQaC1w3Ire11Jtcx9\nY2d8LYCDItIamAmMAerb/3g/A+Lw2dIMX5tkn6tQOgoSkTCsIpvLRCRZRJKxKgmjHes03Hy+N/Ck\nfYy69rWkkv9aHM3CejrrB5wxxnzvwfWAFcxauQjkf8f6viKNMbWxiqpcnT8fY8wPxpghWMWB87Ce\nrKDA34KI1MAq/jrg5DCngeoOy008ObfDeVrafx+5Wrk4T1Gc/c04+5vOIv8NTOEDGXMOeAqIFJGC\n/zO/Yj1lPiAi40uQzipDA4GfMcakYpXv/0usSt7qIhIsIleLyCv2brOBZ0WkoYg0sPcvURt42yUi\nMtTOvMZiBaK1WGXcBquoA7Ga6HUqxXkK+gh4RESa25WmT7nZ93qsIoAOQIz9cxHwLVYFclFqYWUs\nh4EgEXkO627SKTvjzwH+SeGngRSssmxn1mMFuMkiUkNEQkWkp0Ma0oBUEWkOjPMg3YhINbH6Z4Qb\nYzKxijxy78xnA3eJSIxYDQb+jlW/scfJoTYDQ+2/qXbAPcW4rnVYd/lP2n+PlwODgDmeXEMBh+30\nO55rNvCoiLQRkZr2dSQ4e7IqyBiTgfV7es7Jtm1YwWCciIwtQVqrBA0EfsgY80/gMeBZrH+afVh3\n5fPsXV4CNgBbgK3AJntdSc3Hqgg+DowAhhqrpdJ2rH+w77EyiUhgdSnOU9BMYAnWdfwIfImVWWc7\n2Xck8D9jzF5jTHLuDzAdGF5EURrAYuArrLvERCCdooui3sO65oJBdiIwyy4uG+a4wS6SG4RVMbwX\nq5jtFnvz37AqdVOBhVgVop4aAeyxi5Tux3pawRizDPgr8AlWAGqLVSTizOtABtbvchZWha6n15Vh\nX9fVWJXT/wfcaYz5pRjXkHusM8AkYLV9rm7AO1gBdxVWY4V0rEYInnoH60lskJPz/YRVj/a8iNxf\n3PRWBWJXrKgqSkQmYlUQ3uEHabkaeMsY07rIncuBiNwJxBlj/L4jn1KloU8EymdEJMxuHx9kF5U8\nD3zm63QBiEh1rGaqM3ydFqW8TQOB8iXBKi45jlU0tAMn5bzlTUQGYBXJpWC1TFKqUtOiIaWUquL0\niUAppao4rw8SVhYaNGhgIiIifJ0MpZSqUDZu3HjEGNOwqP0qRCCIiIhgw4YNvk6GUkpVKCKSWPRe\nWjSklFJVngYCpZSq4jQQKKVUFVch6gicyczMZP/+/aSnp/s6KUp5RWhoKC1atCA4ONjXSVGVXIUN\nBPv376dWrVpEREQg4tEAjkpVGMYYjh49yv79+2nTpo2vk6MquQpbNJSenk79+vU1CKhKSUSoX7++\nPvFWZfHxEBEBAQHWa3zBMQLLToV9IgA0CKhKTf++q7D4eIiLgzNnrOXERGsZYPjwMj9dhX0iUEqp\nSmvChD+DQK4zZ6z1XqCBoBREhMcffzxvecqUKUycOLFc07BhwwYefvjhEn328ssvL1FHveeee45l\ny5YBMHXqVM44/MHWrFn07IS//PIL3bt3JyQkhClTpuTb9tVXX3HBBRfQrl07Jk+e7PTz8+bNY/v2\n7U7T4w0Fr1Epr9u7t3jrS6lKBIK3vtnNmt1H8q1bs/sIb32zu1THDQkJ4dNPP+XIkSNF7+wFWVlZ\nxMbGMm3atHI97wsvvMCVV14JlCyTrFevHtOmTeOJJ57Itz47O5sHH3yQRYsWsX37dmbPnp0vw89V\nMBA4pscbNBCocteqVfHWl1KVCARRLcIZ8+GPecFgze4jjPnwR6JahJfquEFBQcTFxfH6668X2jZq\n1Cjmzp2bt5x7p7xy5Uouu+wyhgwZwnnnncf48eOJj4/n0ksvJTIykt27reB0+PBhbrzxRrp06UKX\nLl1YvdqaGGzixImMGDGCnj17MmLECFauXMl1110HQFpaGnfddReRkZFERUXxySefAPCXv/yF2NhY\nOnbsyPPPP+/2mn744QeGDh0KwPz58wkLCyMjI4P09HTOO++8fNc2bdo0Dh48SN++fenbt2/eMSZM\nmEB0dDTdunUjJaXwlLONGjWiS5cuhZpFrl+/nnbt2nHeeedRrVo1br31VubPn59vnzVr1rBgwQLG\njRtHTEwMu3fvzvddR0RE8PTTTxMTE0NsbCybNm1iwIABtG3blrfeeivvOK+++ipdunQhKioq7zs5\nffo01157LdHR0XTq1ImEhASn17hkyRK6d+/OxRdfzM0330xaWlreuZ988kkiIyO59NJL+e233wD4\n+OOP6dSpE9HR0fTp08ft968U4LwIqHp1mDTJK6er0JXFuf72+Ta2Hzzpdp9GtUK48+31NK4dQsrJ\nc7RrVJM3lu3ijWW7nO7foVltnh/UschzP/jgg0RFRfHkk096nN6ffvqJHTt2UK9ePc477zzuvfde\n1q9fzxtvvMGbb77J1KlTeeSRR3j00Ufp1asXe/fuZcCAAezYsQOA7du389133xEWFsbKlSvzjvvi\niy8SHh7O1q1bATh+/DgAkyZNol69emRnZ9OvXz+2bNlCVFSU07R17tyZzZs3A/Dtt9/SqVMnfvjh\nB7KysujatWu+fR9++GFee+01VqxYQYMGDQArM+3WrRuTJk3iySefZObMmTz77LMefS8HDhygZcs/\n569v0aIF69aty7dPjx49GDx4MNdddx033XST0+O0atWKzZs38+ijjzJq1ChWr15Neno6nTp14v77\n72fJkiXs2rWL9evXY4xh8ODBrFq1isOHD9OsWTMWLlwIQGpqKuHh4fmu8ciRI7z00kssW7aMGjVq\n8PLLL/Paa6/x3HPWNAq53/97773H2LFj+eKLL3jhhRdYvHgxzZs358SJEx59F6qKy71JatIEUlKs\nJ4FJk7xSUQxeDgQisgc4hTUHbZYxJlZE6gEJQASwBxhmjDnuzXQAhIcF07h2CAdOpNO8TijhYWXT\nSad27drceeedTJs2jbCwMI8+06VLF5o2bQpA27ZtueqqqwCIjIxkxYoVACxbtixf8cfJkyfz7jwH\nDx7s9FzLli1jzpw/5xKvW7cuAB999BEzZswgKyuLpKQktm/f7jIQBAUF0bZtW3bs2MH69et57LHH\nWLVqFdnZ2fTu3bvIa6tWrVreE8oll1zC0qVLi/xMWRs8eDBgfZ9paWnUqlWLWrVqERISwokTJ1iy\nZAlLliyhc+fOgPUktWvXLnr37s3jjz/OU089xXXXXef0eteuXcv27dvp2dOajz4jI4Pu3bvnbb/t\nttvyXh999FEAevbsyahRoxg2bFje05ZSbiUkWE1Gf/8dyqH1WHk8EfQ1xjgWoo8HlhtjJovIeHv5\nqdKcwJM799zioIevaMcH6/byyJXt6dG2QWlOm2fs2LFcfPHF3HXXXXnrgoKCyMnJASAnJ4eMjIy8\nbSEhIXnvAwIC8pYDAgLIysrK+8zatWsJDQ0tdL4aNWp4nLY//viDKVOm8MMPP1C3bl1GjRpVZNv0\nPn36sGjRIoKDg7nyyisZNWoU2dnZvPrqq0WeLzg4OK/ZY2BgYN71eKJ58+bs2/fnfPL79++nefPm\nHn8+l+P3WfC7zsrKwhjD008/zX333Vfos5s2beLLL7/k2WefpV+/fnl3+rmMMfTv35/Zs2c7Pbdj\nk8/c92+99Rbr1q1j4cKFXHLJJWzcuJH69esX+7pUFXHkCCxdCk88US5BAHxTRzAEmGW/nwVc7+0T\n5gaB6bd35rGrLmD67Z3z1RmUVr169Rg2bBhvv/123rqIiAg2btwIwIIFC8jMzCzWMa+66irefPPN\nvOXc4hp3+vfvz7/+9a+85ePHj3Py5Elq1KhBeHg4KSkpLFq0qMjj9O7dm6lTp9K9e3caNmzI0aNH\n2blzJ506dSq0b61atTh16pSHV+Vely5d2LVrF3/88QcZGRnMmTMn7+6+LM85YMAA3nnnnbwnrAMH\nDnDo0CEOHjxI9erVueOOOxg3bhybNm0qdL5u3bqxevXqvPL/06dP8+uvv+YdOyEhIe8190lh9+7d\ndO3alRdeeIGGDRvmC3ZKFfLpp5CdDbfcUm6n9PYTgQGWiUg28B9jzAygsTEmyd6eDDR29kERiQPi\nwCrzLY0t+1OZfnvnvCeAHm0bMP32zmzZn1pmTwWPP/4406dPz1sePXo0Q4YMITo6moEDBxbrLh5g\n2rRpefUPWVlZ9OnTJ19lpzPPPvssDz74IJ06dSIwMJDnn3+eoUOH0rlzZy688EJatmyZV6ThTteu\nXUlJScmr2IyKiiI5OdlpB6e4uDgGDhxIs2bN8oq1ipKcnExsbCwnT54kICCAqVOnsn37dmrXrs30\n6dMZMGAA2dnZ3H333XTsWPhp79Zbb2X06NFMmzYtX4W8p6666ip27NiRl1HXrFmTDz74gN9++41x\n48YREBBAcHAw//73v51e47vvvsttt93GuXPnAHjppZc4//zzASv4RkVFERISkvfUMG7cOHbt2oUx\nhn79+hEdHV3sNKsqJCEB2reHmJhyO6VX5ywWkebGmAMi0ghYCjwELDDG1HHY57gxpq6748TGxpqC\n7d137NjBRRdd5I1kK1UiuRMo5VaclwX9O69ikpOheXOr1dALL5T6cCKy0RgTW9R+Xi0aMsYcsF8P\nAZ8BlwIpItLUTmRT4JA306CUUhXG3LmQk1OuxULgxUAgIjVEpFbue+Aq4GdgATDS3m0kMN/5EZSq\nWPbs2VOmTwOqCkpIgI4drZ9y5M06gsbAZ3a5chDwoTHmKxH5AfhIRO4BEoFhXkyDUkpVDPv2wXff\nwYsvlvupvRYIjDG/A4VqxYwxR4F+3jqvUkpVSB9/bL2Wc7EQVJEhJpRSyu8lJEDnzlaLoXKmgUAp\npXzt999h/Xq49VafnF4DQSnoMNTlNwz1sWPH6N+/P+3bt6d///554yg52rNnDx9++GHecmm+G0+s\nXLmSNWvWeO34qgr56CPrdZhvqkyrTCCI3xpPxNQIAv4WQMTUCOK3ln7aNx2GuvyGoZ48eTL9+vVj\n165d9OvXz+lcBQUDgbe/Gw0EqswkJEDXrtb4Qj5QJQJB/NZ44j6PIzE1EYMhMTWRuM/jSh0MdBjq\n8huGev78+YwcabU6HjlyJPPmzSt03PHjx/Ptt98SExPD66+/nu+7mThxIiNHjqR37960bt2aTz/9\nNG/I6IEDB+YNAbJx40Yuu+wyLrnkEgYMGEBSktUJftq0aXTo0IGoqChuvfVW9uzZw1tvvcXrr79O\nTEwM3377bZG/s+7du9O+fXtmzpwJQFJSEn369CEmJoZOnTrx7bffuv3dqEpq507YvNlnxUJQSYah\nHvvVWDYnux6LZ+3+tZzLPpdv3ZnMM9wz/x5mbpzp9DMxTWKYOnBqkefWYajLZxjqlJSUvBFbmzRp\n4jTATJ48mSlTpvDFF18A5PtuwBrzZ8WKFWzfvp3u3bvzySef8Morr3DDDTewcOFCrr32Wh566CHm\nz59Pw4YNSUhIYMKECbzzzjtMnjyZP/74I28E0zp16nD//fdTs2bNvCeb22+/3eXvbMuWLaxdu5bT\np0/TuXNnrr32WmbPns2AAQOYMGEC2dnZOvlNVZWQYA0ud/PNPktCpQgERSkYBIpaXxw6DHV+5TEM\ntYiUaGL3q6++muDgYCIjI8nOzmbgwIGA9b3v2bOHnTt38vPPP9O/f3/AKqrK/T1FRUUxfPhwrr/+\neq6/3vk4ie5+Z0OGDCEsLIywsDD69u3L+vXr6dKlC3fffTeZmZlcf/31xJTj2DLKjyQkQK9e1tAS\nPlIpAkFRd+4RUyNITE0stL51eGtWjlpZ6vPrMNR/8tYw1I0bNyYpKYmmTZuSlJREo0aNPD5uLsfv\n2TGdjsNTd+zYke+//77QZxcuXMiqVav4/PPPmTRpUt5TlyN3v7OCgUtE6NOnD6tWrWLhwoWMGjWK\nxx57jDvvvLPY16UqsJ9/hu3bwWHUYF+oEnUEk/pNonpw9XzrqgdXZ1K/spn2TYeh9v4w1IMHD2bW\nLGv08lmzZjFkyJAyT8sFF1zA4cOH8wJBZmYm27ZtIycnh3379tG3b19efvllUlNT8ya8cTyfu9/Z\n/PnzSU9P5+jRo6xcuZIuXbqQmJhI48aNGT16NPfee2/esNeqCpkzBwIC4MYbfZqMKhEIhkcOZ8ag\nGbQOb40gtA5vzYxBMxgeWXbTvj3++OP5Wg+NHj2ab775hujoaL7//vsSDUO9YcMGoqKi6NChQ5FD\nUIM1DPXx48fz5sddsWIF0dHRecNQ33777SUehjoyMtLtMNSOlcVFSU5OpkWLFrz22mu89NJLtGjR\ngpMnTxIUFJQ3DPVFF13EsGHD8oahHj9+PEuXLqV9+/YsW7aM8ePHFzpuVFQUgYGBREdHO63AL0q1\natWYO3cuTz31FNHR0cTExLBmzRqys7O54447iIyMpHPnzjz88MPUqVOHQYMG8dlnn+VVFrv7nUVF\nRdG3b1+6devGX//6V5o1a8bKlSvzfj8JCQk88sgjxU6zqsCMsYqF+vaFxk5H4y83Xh2GuqzoMNSq\nIps4cWK+SuXi0L/zSmzTJrjkEpg5E+691yun8IthqJVSSrkwZw4EBYEfzGNdKSqLlfJn5d3bXFUA\nxli9ifv3h3r1fJ0afSJQSqlyt24dJCb6tBOZIw0ESilV3ubMgWrVwEnrN1/QQKCUUuUpJ8eae+Dq\nqyE83NepATQQKKVU+fruOzh40G+KhUADQanNmzcPEeGXX34p8TEKDlDnzN///vd8yz169CjRuSZO\nnFho6GelVDmaMwfCwsAeisUfVJ1AEB9vDfEaEGC9xpd+GGqA2bNn06tXL2bPnl0mx3OlYCDQ4Y+V\nqoCysmDuXCsIFDF3hzeGznelagSC+HiIi7Nq6Y2xXuPiSh0M0tLS+O6773j77bfzBntbuXIll19+\nOTfddBMXXnghw4cPJ7fT3gsvvECXLl3o1KkTcXFxFOzM9/XXX+cb0Gzp0qXccMMNjB8/nrNnzxIT\nE8Pw4VZvaMcJYF5++WUiIyOJjo7O63E7c+ZMunTpQnR0NDfeeKOObKmUP1i5Eg4fLnJeYm8Nne9K\n5ehHMHasNZ63K2vXwrkCI42eOQP33GP16nMmJgamuh/Mbv78+QwcOJDzzz+f+vXr540t9OOPP7Jt\n2zaaNWtGz549Wb16Nb169WLMmDE899xzAIwYMYIvvviCQYMG5R2vb9++PPDAAxw+fJiGDRvyv//9\nj7vvvptBgwYxffp0p+MNLVq0iPnz57Nu3TqqV6/OsWPHABg6dCijR48GrKEn3n77bR566CG316OU\n8rI5c6wngWuucbvbhOUTOJOZ/+btTOYZJiyfUKZD4+SqGk8EBYNAUes9NHv2bG61K3xuvfXWvOKh\nSy+9lBYtWhAQEEBMTAx79uwBYMWKFXTt2pXIyEi+/vprtm3blu94IsKIESP44IMPOHHiBN9//z1X\nX3212zQsW7aMu+66i+rVrUH16tmdU37++Wd69+5NZGQk8fHxhc6llCpnGRnw6adWk9Eihqzfm7q3\nWOtLq3I8ERRx505EhFUcVFDr1tajWgkcO3aMr7/+mq1btyIiZGdnIyJce+21+YaZzh2KOT09nQce\neIANGzbQsmVLJk6c6HQ46LvuuotBgwYRGhrKzTffTFBQyX5Fo0aNYt68eURHR/Puu+8WmqRFKVXO\nli2D48eLLBYCaBXeyunQ+a3CW3kjZVXkiWDSJKiefxhqqle31pfQ3LlzGTFiBImJiezZs4d9+/bR\npk0bl9MN5mb6DRo0IC0tzWUroWbNmtGsWTNeeumlfPMbBAcHOx3Kun///vzvf//LqwPILRo6deoU\nTZs2JTMzk/gyqhhXSpVCQgLUqQP2RFTuvNj3xULrynLo/IKqRiAYPhxmzLCeAESs1xkzrPUlNHv2\nbG644YZ862688UaXrYfq1KnD6NGj6dSpEwMGDKBLly5ukjucli1b5ht1Mi4uLm+WLEcDBw5k8ODB\nxMbGEhMTk9c09MUXX6Rr16707NmTCy+8sKSXqZQqC+np8NlncMMN4FBi4Eqbum0AaBDWwGtD5zvS\nYaj90JgxY+jcuTP33HOPr5OifKwy/51XKfPmWUHgq69gwIAid3988eNM/2E6h8cdpnZI7RKf1tNh\nqCtHHUElcskll1CjRg3++c9/+jopSqmykpAADRrAFVcUuasxhnk759GvTb9SBYHi0EDgZ3KboCql\nKonTp2HBAhgxAoKDi9x966Gt/H78d8b3LDwLn7dU6DqCilCspVRJ6d93JbFwodVvyYNOZBFTI4h+\nKxqAzJzizXNeGl4PBCISKCI/isgX9nI9EVkqIrvs17olOW5oaChHjx7VfxZVKRljOHr0KKGhob5O\niiqthARo0gTsOcCdcexJnGvc0nFeHVbCkdcri0XkMSAWqG2MuU5EXgGOGWMmi8h4oK4x5il3x3BW\nWZyZmcn+/fudtsVXqjIIDQ2lRYsWBHtQnKD81MmT0KiRNaTNtGkud4uYGuG030Dr8NbsGbunxKf3\ni8piEWkBXAtMAh6zVw8BLrffzwJWAm4DgTPBwcG0adOm9IlUSilvWbDAGsGgiGKh8u5JXJC3i4am\nAk8COQ7rGhtjkuz3yUBjZx8UkTgR2SAiGw4fPuzlZCqllBckJEDLltC9u9vdXPUY9lZP4oK8FghE\n5DrgkDHGZTMYY5VLOS2bMsbMMMbEGmNiGzZs6K1kKqWUdxw/DosXw7Bh1vD3bkzqN4mwoPzjD3mz\nJ3FB3nwi6AkMFpE9wBzgChH5AEgRkaYA9ushL6ZBKaV847PPIDPTo7GFhkcOZ/AFgwHKpSdxQV6r\nIzDGPA08DSAilwNPGGPuEJFXgZHAZPt1vrfSoJRSPpOQAOedB7FF1tUC8Nux37i0+aWsu3edlxNW\nmC/6EUwG+ovILuBKe1kppSqPw4dh+XLraUCkyN13H9vNxqSN3NKx6KcHbyiXnsXGmJVYrYMwxhwF\n+pXHeZVSyic++QSysz3qRDZh+YS8pqNBAb4Z7EGHmFBKqbKWkAAXXghRUS53ye1E5jgT2dPLn6Z+\n9frlVjeiC0MHAAAgAElEQVSQq0IPMaGUUn4nKQm++abIYiF301GWNw0ESilVlj7+GIzx+05kjjQQ\nKKVUWUpIsIqEiphHwtedyBxpIFBKqbKydy+sWeNR34FJ/SYRGpR/UMHy7ETmSAOBUkqVlY8+sl49\n7ER2eevLAd90InOkrYaUUqqsJCTAJZdA27ZF7pqVk8WPyT8y9KKhfDLsk3JInGv6RKCUUmVh927Y\nsAFuvdWj3b/+42tSTqf45AmgILdPBCLSHbgD6A00Bc4CPwMLgQ+MMaleT6FSSlUECQnW67Bhbndz\n7EQmCCfPnSyHxLnnMhCIyCLgINZYQJOwBocLBc4H+gLzReQ1Y8yC8kioUkr5tYQEa7jpVq5b/RTs\nRGYwPPjlgwQHBvv0ycDlDGUi0sAYc8Tthz3Ypyw4m6FMKaX8xo4d0KEDvPEGPPywy928NROZK57O\nUOayjsBdBi8iq4vaRymlqoyEBKsX8U03ud3NnzqROSppZXH593hQSil/ZIwVCPr0gWbN3O7qT53I\nHJU0EHh3xnullKootm6FX37xqLXQpH6TCA4IzrfOV53IHLmrLB7qahMQ5mKbUkpVLQkJEBgIN95Y\n5K63d7qdxxY/Rmp6KhnZGbQKb8WkfpN83oTUXfPRQW62fVHWCVFKqQrHGJgzB664AjyYW/27vd9x\n6PQhZl0/izuj7yyHBHrGXSCYCXxvXDUrUkqpqm7jRvj9d5jgfujogn0HsnKyyimBnnFXRzAC2Cgi\nc0RklIg0Ka9EKaVUhZCQAMHBcMMNLnfJ7TuQ22zUYHho0UPEb40vr1QWyWU/grwdRC4ErgYGAOHA\nCuArYLUxJtvrKUT7ESil/FBODkREWENOf+G6tLy8+w44KnU/glzGmF+MMa8bYwYCVwDfATcD60qf\nTKWUqqDWroV9+4psLeSvfQccFTn6qIjUK7BqLbDUGJPpnSQppVQFkJAAISEweLDb3VqFt3L6RODr\nvgOOPOlHsAk4DPwK7LLf7xGRTSJyiTcTp5RSfik725p74JproHZtt7tO6jeJoID899z+0HfAkSeB\nYClwjTGmgTGmPlZ9wULgAeD/vJk4pZTyS99+C8nJHnUiG3LBEAIlkBrBNXw+AY0rnkxM080YMzp3\nwRizRESmGGPiRCTEi2lTSin/lJAA1avDtdcWueucn+dwLvscy+9cTs9WPcshccXnyRNBkog8JSKt\n7Z8ngRQRCQRyvJw+pZTyL1lZMHcuDBoENWq43C1+azwRUyMY/floggOC2XNiT/mlsZg8eSK4HXge\nmIc1xtBqe10g4H4GBqWUqmy+/hqOHHFbLFRw3oHMnEzivogDwa+KhHIV2Y/AH2g/AqWU37jnHuuJ\nICUFQkOd7uLLvgOOSt2PQERmikiki201RORuEfG/0KaUUt6SkQGffgpDhrgMAlAx+g44clc09C/g\nr3Yw+Bmr2Wgo0B6oDbwD+E8faaWU8rYlS+DECbjlFre7VYS+A45cBgJjzGZgmIjUBGL5c/L6HcaY\nnUUdWERCgVVAiH2eucaY5+0OaglABLAHGGaMOV7K61BKKe9LSIC6daF/f7e7vdj3RUbOG4lxmLrF\n3/oOOPJkiIk0Y8xKY8xsY8w8T4KA7RxwhTEmGogBBopIN2A8sNwY0x5Ybi8rpZR/O3sW5s+HoUOh\nWjW3u9aoVgODoWH1hn7bd8CRJ62GSsQevjrNXgy2fwwwBLjcXj8LWAk85a10KKVUmVi0CE6dclss\n5DjcdKAE8s+r/smI6BHlmMiSKelUlR4RkUAR2QwcwhqfaB3Q2BiTZO+SDDT2ZhqUUqpMJCRYk8/0\n7et0c8HhprNNNvcvvN+vhpt2pchA4KrlkCeMMdnGmBigBXCpiHQqsN3gYv5jEYkTkQ0isuHw4cMl\nTYJSSpXe6dPWUNM33QRBzgtSJiyfkNdvINeZzDNMWO5+0hp/4MkTwf+JyHoReUBEwktyEmPMCax5\nDAZi9UpuCmC/HnLxmRnGmFhjTGxDD6aAU0opr/n8czhzxm2xUEVrMurIk8ri3sBwoCXWjGUfioj7\nKnNARBqKSB37fRjQH/gFWACMtHcbCcwvYdqVUqp8JCRAs2bQq5fLXVw1DfXXJqOOPKojMMbsAp7F\nqtS9DJgmIr+IyFA3H2sKrBCRLcAPWHUEXwCTgf4isgu40l5WSin/lJpqVRTffDMEBrrc7a99/lpo\nnT83GXXkycQ0UcBdwLVYQ1IPMsZsEpFmwPfAp84+Z4zZAnR2sv4o0K80iVZKqXIzfz6cO1dkJ7K0\nDKuRZJMaTUg5nUKr8FZM6jfJb5uMOvKk+eibwH+BZ4wxZ3NXGmMOisizXkuZUkr5g4QEaN0aunVz\nujl+azzPLH+Gval7CQkMYcqAKRUi83fkSdHQZ8aY9x2DgIg8AmCMed9rKVNKKV87dswaVmLYMBAp\ntDm3yWhuhfC57HPEfR5XIZqMOvIkENzpZN2oMk6HUkr5n08/teYfcFEsVJGbjDpyWTQkIrdhzTvQ\nRkQWOGyqBRzzdsKUUsrnEhKgXTu4+GKnmytyk1FH7uoI1gBJQAPgnw7rTwFbvJkopZTyuUOHrElo\nnn7aabEQQMvwlk4z/YrQZNSRy6IhY0yiPdhcd2PMNw4/m4wxWeWZSKWUKneffAI5OW5bC93e6fZC\n6ypKk1FH7oqGvjPG9BKRU+QfBkKwRoeo7fXUKaWUr8yZAx06QKdOhTY5Di4nCPXC6nHs7LEK1WTU\nkbv5CHrZr7XKLzlKKeVj8fHw1FNw4ACEh8OHH8LwPzP2gvMRGwxns87y/tD3K1wAyOXJoHNtRSTE\nfn+5iDycO3SEUkpVKvHxEBdnBQGwehXHxVnrbZWlpZAjT5qPfgJki0g7YAbWmEMfejVVSinlCxMm\nWIPLOTpzxlpvqywthRx5Eghy7MrhG4A3jTHjsMYRUkqpymWvi8zcYX3L8JZOd6loLYUceRIIMu0+\nBSOBL+x1wd5LklJK+Ujz5s7Xt/ozk68sLYUceRII7gK6A5OMMX+ISBtAh5ZQSlU+rVsXXle9Okya\nRPzWeFpPbc3k1ZMRhPph9SvEfMSeKHLQOWPMduBhh+U/gJe9mSillCp3n30Gq1fDDTfApk1WcVCr\nVlYQiKLStRRyJNZskW52EOkJTARaYwWO3H4E53k9dbbY2FizYcOG8jqdUqqqOXTI6i/QsiWsXQvB\n+Uu/I6ZG5M1F7Kh1eGv2jN1TToksPhHZaIyJLWo/T4ahfht4FNgIZJc2YUop5VeMgfvus5qKrlhR\nKAhA5Wwp5MiTQJBqjFnk9ZQopZQvvP8+zJsHr7wCHTs63aVF7RbsO7mv0PqK3FLIkSeVxStE5FUR\n6S4iF+f+eD1lSinlbfv2wcMPW3MRP/ZYoc3xW+OJmBrhNAhU9JZCjjx5IuhqvzqWMxngirJPjlJK\nlRNj4J57rPkG3n230HzEBYeSABAEg6F1eOsKOaaQK560GupbHglRSqly9e9/w9Kl1mvbtoU2OxtK\nIjcI+HMFcUl4MtZQYxF5W0QW2csdROQe7ydNKaW85LffYNw4GDDAqih2orJXEDvypI7gXWAx0Mxe\n/hUY660EKaWUV2Vnw8iRUK0avP22y0lnXFUEV5YKYkeeBIIGxpiPgBwAe9whbUaqlKqYpkyBNWtg\n+nTXQ0oAwzoOK7SuMlUQO/Kksvi0iNTHnpxGRLoBqV5NlVJKecPWrfDcczB0KNxeeMwgsCqJn1n+\nDHtT91aKSWc84UkgeAxYALQVkdVAQ+Amr6ZKKaXKWkYG3Hkn1KkDb73ltEioMk4644kii4aMMZuA\ny4AewH1AR2OMTl6vlKpYXnwRNm+GGTOgYUOnu1TGSWc84W7O4qEuNp0vIhhjPvVSmpRSqmytXw//\n+IdVSTxkiMvdqlJLIUfuioYG2a+NsJ4GvraX+wJrAA0ESin/d/asVSTUrBm88YbbXRvWaMih04cK\nra+MLYUcuZu8/i4AEVkCdDDGJNnLTbGalCqllP97+mnYudPqPBYe7nSX+K3xPL3saadBoLK2FHLk\nSWVxy9wgYEsBKnd4VEpVDitWWE8BDz4IV17pdJeqNJSEK54EguUishiYbS/fAiwr6kMi0hJ4D2iM\n1fR0hjHmDRGpByQAEcAeYJgx5njxk66UUm6cPAl33QXt2sHLrufSqkpDSbjiSauhMcBbQLT9M8MY\n85AHx84CHjfGdAC6AQ+KSAdgPLDcGNMeWG4vK6VU2XrsMWt00ffegxo1XO5WVSuIHXnyRIAx5jPg\ns+Ic2C5OSrLfnxKRHUBzYAhwub3bLGAl8FRxjq2UUm4tXGgNHzF+PHTv7nSX+K3xTFg+AYPzWRor\newWxI48CQWmJSATQGVgHNHaoc0jGKjpy9pk4IA6gVauq8wtRSpXS0aNw770QFQUTJzrdxVm9gKOq\nUEHsyJOxhkpFRGoCnwBjjTEnHbcZa8Jkp+HYGDPDGBNrjIlt6KLzh1JKFfLAA1YweO89CAlxuouz\neoFcrcNbM2PQjEpfQezIq08EIhKMFQTiHTqgpYhIU2NMkt0UtXB7LaWUKok5c+Cjj2DSJIiOdrmb\nq/J/QapMBbEjdz2Lt+L8bl2wbuaj3B1YRARr4vsdxpjXHDYtAEYCk+3X+cVNtFJKFXLwoPU00LUr\nPPmk212b1GxCUlpSofVVqV7AkbsngutKeeyewAhgq4hsttc9gxUAPrInt0kECo/1qpRSxWGMVS+Q\nng6zZkFQ4awtt3I4d1TRgqpavYAjdz2LE0tzYGPMd+Dk27b0K82xlVIqn7ffhkWLrM5jF1xQaLOz\nUUUDJZA6oXUq/RDTnnBXNHSKP4uGcjN0w59FQ7W9nDallCraH3/Ao49C374wZozTXZxVDmebbGpW\nq8mRJ4+URyr9mrsnglrlmRCllCq2nByr97AI/O9/EOC8IaR2GnPPo+ajItJLRHIHoWsgIm28myyl\nlPLAG2/AN99Yr61bF9ocvzWeiKkR2mmsCEU2HxWR54FY4ALgf0A14AOsymCllPKNHTuskUUHDYJR\nowpt1k5jnvPkieAGYDBwGsAYcxDQYiOllO9kZlpzDNSsac045mTaSe005jlPOpRlGGOMiOROXu96\n9CallCoP//gHbNhgdR5r0sTpLtppzHOePBF8JCL/AeqIyGisIahnejdZSinlwqZN1vzDt90GN99c\naLPWCxRfkU8ExpgpItIfOIlVT/CcMWap11OmlFIFpadbRUING8L06YU2a71AybjrR9AOa6TQ1XbG\nv9Re30tE2hpjdpdXIpVSCoDnnoNt2+DLL6FevUKbi6oXqMqdxtxxVzQ0FespoKBUe5tSSpWf776D\nKVMgLg6uvjrfptzioMRU5wMi5NYLaBBwzl3RUGNjzNaCK40xW+35BZRSqnykpcHIkRARYQUDB0UV\nB4HWCxTFXSCo42ZbWFknRCmlXBo3zhpKYuVKqJW/9bq74iDQegFPuCsa2mC3EspHRO4FNnovSUop\n5WDxYnjrLWs8oT59Cm12N0yE9hfwjLsngrHAZyIynD8z/lisnsU3eDthSinF8eNwzz1w0UXWZDMO\nippzuHV4a+0v4CF3g86lAD1EpC/QyV690BjzdbmkTCmlHn4YkpNh3jwIDc1brc1Ey5Yn/QhWACvK\nIS1KKfWnTz+FDz6A55+H2Nh8m7SZaNny6pzFSilVIikpcN99cPHFMGFC3urc4qCimomq4tFAoJTy\nL8ZYQeDUKXjvPQgOBrSZqDd5NB+BUkp5XXy81U8gIADmz4ehQ6Fjx7zN2kzUezQQKKV8Lz7e6jGc\n6FDkM38+xMcX2WsYtJloaWnRkFLK9yZMgDMF7vbPnCFt3CPEjTnr9klAm4mWngYCpZTv7XXeKax6\n0lHOZLr+mBYHlQ0tGlJK+U5GBjz5pFVB7MTecNcf1eKgsqNPBEop3/j9d7j1VvjhB7jiCvj+ezh7\nNm/z6WB4pp/zj1b24qC3vtlNVItwerRtkLduze4jbNmfyv2XtS3z8+kTgVKq/M2eDTEx8Ouv8PHH\nsHw5zJwJrVuTA+wJh9GDYHZU4Y9WheKgqBbhjPnwR9bsPgJYQWDMhz8S1cLNI1Ip6BOBUqr8nD5t\nDRvxzjvQowd8+CG0bg1AfBRMGAuJqa4/XlV6Dfdo24C/39CJu9/9gZHdI/h4436m39453xNCWdJA\noJQqH5s3W0VBv/5qtRKaOBGCrCzIk85ilb04KFfauSxmrvqd/377O+mZOfxn1e88fEU7rwUB0ECg\nlPI2Y+Bf/4LHH4f69WHZMqtOgKKHjMhVFYqDzmVl8+G6vUz/+jeOns7g0oi67ExJY2T31nywbi/d\n2tbXJwKlVAV09Kg1jPT8+XDNNfDuu9bE83j2FACVvzgoO8cwf/MBXlv6K/uPn6X7efW5JrIJry/b\nxb/vuJgebRvQrW19xnz4o9eKh7wWCETkHeA64JAxppO9rh6QAEQAe4Bhxpjj3kqDUsqHVq2C4cOt\nAeReew0eecQaPsJW1JARULmLg4wxfP3LIV5dvJNfkk/RsVlt/n5DJL3bN+A/q37Pl+n3aNuA6bd3\nZsv+VK8EAm+2GnoXGFhg3XhguTGmPbDcXlZKVSbZ2fC3v0HfvtYcAt9/b80uZgcBT4aMgMpdHLRh\nzzGG/ed77pm1gfTMbN68rTOfj+lFn/MbIiLcf1nbQhl+j7YNvNJ0FLz4RGCMWeVkkvshwOX2+1nA\nSuApb6VBKVXO9u+3ngJWrYIRI6y6AYc5hqt6cdAvySeZsngny3YcomGtEF66vhO3dGlJcKBvW/KX\ndx1BY2NMkv0+GWjsakcRiQPiAFq10qFllfJ78+fD3XfDuXMwaxbceWfepuJUClfG3sL7j5/htaW/\n8tmPB6gZEsS4ARdwV88Iqlfzj2pan6XCGGNExHm/cmv7DGAGQGxsrMv9lFI+lp4O48bB9OnQuTPM\nmQPnn5+3uSo/BRxNO8f0Fb8Rv3YvCIzufR5/uawtdWtU83XS8invQJAiIk2NMUki0hQ4VM7nV0qV\npV9+sfoG/PQTjB0LkydDSAjg+VMAVL5K4bRzWfz329+Zuep3zmZmc/MlLXnkyvY0qxPm66Q5Vd6B\nYAEwEphsv84v5/MrpcqCMVZT0DFjoHp1+PxzuO66vM2ePgVA5aoULtgXYGDHJjwx4HzaNapV9Id9\nyJvNR2djVQw3EJH9wPNYAeAjEbkHSASGeev8SikvOXkS7r/fGi+ob19rgvlmzYDiPQVA5SkOys4x\nLPjpAP9cYvUF6HZePf478EI6t6rr66R5xJuthm5zscnFeIJKKb/3ww9WUVBiIrz0EowfD4GBQPGf\nAipDpbAxhhU7D/HKV1ZfgA5NazPr7kj6tG+AiPg6eR7zjyprpZR/y8mxOoU9/bR19//NN9CzZ94T\nwN7UvQRIANkmu8hDVZangI2Jx3h50U7W7zlG6/rVmXZbZ66LbEpAQMUJALk0ECil3EtJgZEjYfFi\na0L5//4X6tYt9ARQVBCoLE8BO5NP8erinSzbkUKDmiG8eH0nboltSbWgijuqvwYCpZRrS5daHcNS\nU+Hf/4b77iP+5w+ZMMvzegCoHE8B+4+f4fWlu/j0x/3UrOZ/fQFKo+JfgVKq7GVmwl//Ci+/DB06\n8MX/jWXM3skkvvAXBMHgWdeeivYU4GxmsK9+TuK/3/7Blv2pIHBvrzY8cHk7v+sLUBoaCJRS+f3x\nB9x2G6xbB3FxzLmnK/cseyivCKioIBAogeSYHFqFt6pwTwG5M4NNv70z0S3q8Pz8bczdtB8Bbo5t\nwdgrz/fbvgCloYFAqaouPt6aKGbvXmu+gLQ0CAnh2ykPMSJwAYmLZnh8qIr2BFBQj7YNeHFIR+6d\ntQFjDGczc4iNqMs/boikfWP/7gtQGhoIlKrK4uMhLg7O2E0+jxwhW+DRy9KZnjbd4yIgqNj1AL8f\nTmPxthSWbE/mx70n8tbfHNuCV2+K9mHKyocGAqWqKmOsWcPO5G/3H2jgsbXwZvfKWQ8AVvv/rQdS\nWbwtmSXbUth1KA2AyObh3HxJC5ZsT+HO7q2JX7eXNbuPeHWaSH+ggUCpqub33+H9962flBSnu7Ry\nM4E8kFdhXJGeAjKzc/jhj2NW5r89haTUdAIDhEsj6jG8ayv6d2xC4tHTjPnwx7yZwbp7eWYwf6GB\nQKmq4Phx+PhjeO89WL0agDXtQjg/DBqcLbz73nDXh6pImf/ZjGxW7TrM4m3JLN9xiNSzmYQEBdDn\n/IY8ftUF9LuwUb7WP5//dLBcZwbzFxoIlKqsMjLgq6+sO/8FCyAjg12Ng3mnH3wYCXvrnOO2LTDz\nc6iR+efHTgfDM04GgqkoRUAnzmSwbMchlmxLZtWuw6Rn5hAeFky/CxtxVccm9Dm/gcu2/85mAOvR\ntkGlDgKggUCpysUY2LDBuvOfMweOHOFojQA+6JzDB1GwoVkmOIyAMDvKev37cqs4aG+4FQRy11eU\nIqCDJ86yxC7yWffHMbJzDE3DQ7kltiVXdWzCpW3q+XwWMH+mgUCpyiAx0WoB9N57sHMn54Jg/vnw\nfn/4ql0OWYGuPzo76s+M35E/Z/7GGH47lMbibcks3pbC1gNWpUa7RjW5/7LzuKpDE6JahFeogd98\nSQOBUhXVyZMwd65V9LNyJQCrWsH7g+DjDpBawn5P/loElJNj+HHfCZZst1r6/HHkNAAxLevw1MAL\nuapjY9o2rOnjVFZMGgiUqkiysqzxf957j6zPPiHoXCa/1oMP+sL7UbCnhMPf+7oIyNnQDmt2H+HH\nvSfo1DycJduSWbo9hUOnzhEUIHRvW5+7e7Xhqg6NaVw7tFzTWhlpIFDK3xkDmzfD++9z9r13CDua\nytEwSIiC96JgXQvylft7yteZv6OCQzvMWPU7/165m8AAOJuZQ/VqgVx+QUOu6tCEvhc2Ijws2Gdp\nrYw0ECjlrw4c4McpT1B9zidckJxJRiB8dYEw60r4sj1kluC/158yf7Da9v9++DSHTp6jR9v63Pn2\nenKMIcdArdAgru7UhKs6NKFX+waEBrup6FClooFAKV9xHOOnVSuYNAmGDGHNG0+Q/d4sev6aTmdg\nTQu4/1r4qCMcr+75kA+5/CXzP5meyY6DJ9mRdJLt9s+vKWlkZOUAUC0wgLo1gjl8KoOhFzfnlRuj\nCNKWPuVCA4FSvhAfT9a9dxOUnmEtJyaSNeIOsgR65MDvdeDFy+CDKPitfvEP78vM3xjD/uNn2Z5k\nZ/oHrUx///E/e67Vq1GNDk1rM6pHBBc1rUWHpuEcOpXOI3M28/AV7fhg3V7W7zlW6dvv+wsNBEp5\nmeN0jvVC69L4lOHrqcdpnJ5/vyADZ4Oh3x2wpiXFLvf3ReafnpnNb4fS8jL73Mz/VHqWlSaBNg1q\nEN2yDrdd2ooOzWrToWltGtUKyde0c83uIzwyZ3Ner95uVWRoB3+hgUApL4jfGs/Exc9QffdeYpLh\noRSISoHo5GM0cjO3e40MWNPK8/OUNvN31Vpny/7UQr1sj6adY0fSKbYnpbL94El2JJ3it8NpZOdY\nxVXVqwVyYZNaDIlpxkVNrQz/gia1PJrBa8v+1Co5tIO/EGOKX+ZY3mJjY82GDRt8nQylCsm9208/\nkEiv47XolJRN2/1niE6Gi45AsFX8TXog/NwIfmoCPzWGZ76FJqcLH29POLR51PX5ggOCqR1Sm2Nn\nj5XJxC9rdh/Jd+e9ZvcRxsT/yDPXXkhocGDenf6OpJOknDyX97mm4aF5mf1FTWvToVltWterXiEn\nbq/MRGSjMSa2qP30iUApD+Rm+ElHE+maFk7HpGzO25dGdDKsS4HGpwFOAbC/lpXhf3E+bGlsvd9V\nD7IdGr0cqe75GD/eKPLJys7h0KlzhAQFMqJba+6dtYHzGtZgR9IpAkV44uMtAAQFCO0a1aRn2wZ0\naGZl+hc1rU29SjRNo9JAoKqiAq11vrv/Gu4I/dIqww+rB8Cxs8c4P7sOHZKyaLP3FNHJMD8FLjoM\n1XKs4QzSA2FbI1jY/s8Mf0tjOFa96CR4c4yfjKwcUk6mk5SaTlLqWZJTrffJqekknUwnOfUsh0+d\nI6dAYcDPB07SvE4oV3VsQgf7Lr9do5qEBGmzzcpOi4ZUpZWvktbO4AesO8rMz6F6gTvxv/WBA+EQ\nnWyX5adA07Q/9zlQyyrSyc3wf2oMv9bPf5dfErkZfv2w+mTlGFLPHae1myKf9MzsPzP2k2f/zOAd\nXo+knSv0uZohQTQND6VJeKj9Gpa3fOhkOpMX/cId3ayJWLSCtvLwtGhIA4GqMJxl7MfOHnP6/ujZ\no1Ymawx10qFxmlV8M/cjaOimsvZcIGxrmD/D39IYjtYoXdodM/zctAbRkGd6/o2JV97Pmt1HeDB+\nExOuvYim4WHW3fyJs/YdfG5Gf5bjZzILHTs8LDh/Jl87rECmH0qtUOc9cZ3WEWhrnUpDA4HyS99N\nfoCIV2bQ7Hg2B+oE8OKAMP570RnPM3YMGAhPhyZ25t44zfX7xqchJLvodOUAUX+BnQ1wO1Jn8Qhg\naFS9Bbdf+CQx9QeRejaT1LOZnDybyR9HTrN53wnCw4I5djrD6ezA9WpUo0nt0EJ3883s5SbhoR61\nynGlOK2GVMWjgUCVSnHuvt29bxXeimvaX8OXu76kx7eJTitIRw+C2ZFW5t74tJ2RpxV+75jJO8vc\nswQO1YDkmpBSE1KcvH//U2ieVvizRbXW8UQgtTBGyOEUgaYBdbLupGZ230L7Va8WSHhYMLVDgzmV\nnsnB1HSimodzdWTTfBl+49qhOqyCKhUNBH6uNBmtY+bqWTHJMe7eFsbEpek0P5HDgToB/O2qMN7u\ncIb6Rd19eygoG2pmWO3ga2ZYmX3NjPzrpiyBeumFP5slkBUAoW4y95Sadqbu5v2xMDBFjEjgakau\n0YOcj8mfj/11CLUIECGbU4QFNKJz7b8Q02AQtcOCqR0WlJfJh4cFUzvMfg211tcKDaZakJXI3GKY\nO8T/RSwAAAiLSURBVLq24gMtm1de4NeBQEQGAm8AgcB/jTGT3e1f3EDw1je7afjVJLr+5728Igh3\nGV9pM9riZ8Cni53RloazzO9MEPy1L3wT4TzTdrbOcX3BdZ4Uv7higFd6Os/kj3qQuRfXbVvyt9aZ\n0E/4JKY253JOEhZQGwkQzmSlUjM4nAARTmWcoEnNFjzT82+M6jyCGtUCSz3hiZbNq/Lgt4FARAKB\nX4H+wH7gB+A2Y8x2V58pbiD45Om7GfjP/5Xsri+XgQADgTkQ6OTV3bbc12t3wl+/hbCsPw+bHgjT\nL4XvW0K1bAjOtl4df4JznKwrxX71z1gRt7iyBNKqwelq1mtaNet7zHtfxPqC6757B1qeLHyesiiW\ncc0qp68bUo+AAOHY2WPUqlaH0xnZZHOKYIcK2/KkZfOqPPhzIOgOTDTGDLCXnwYwxvzD1WeKGwj2\n1wuixfHCt6hZAkm1PMzU/aDE7FwgZAZARmD+n8xAJ+vc7Hf/BufD1uQAQ25znblnBFKice5dKVWx\njM1Z6xt3T3YFm2HqnbiqSvy5Z3FzYJ/D8n6ga8GdRCQOiANo1aoYg68AzZwEAbAy96XnQXYAZMuf\nrzlSeF1Rr558Zu5H4KxUIweIud99Bp8VQJllwlfvgojUwuv3hsMXF5TNOTxRsBPV/joBPHNFNeZE\nnfO4yK60PWt1TBulCvPbnsXGmBnADLCeCIrz2YN1A50+ESSGwz3Xl036PLE33HUGvLVJ+aXjmX6e\nD2fgqDh3357UtySm7uXr7s1ZPXEyEZHDaQU8sPsIvcqxOMTZeXq0baBBQFVpvggEB4CWDsst7HVl\nZt19d1LXSR1BURlfWSsqAy5pRlvcyuzZUccICyxcaT3bTeV5Wdx9e0IzYaV8zxeB4AegvYi0wQoA\ntwK3l+UJDg+cwFdQqNWQu4zPG62G3GXA7oYR8LaWwH/tH6WU8lXz0WuAqViNWd4xxkxyt39l7Eeg\nlFLe5s+VxRhjvgS+9MW5lVJK5aczQyulVBWngUAppao4DQRKKVXFaSBQSqkqrkKMPioih4HEEn68\nAXCkDJPjC3oNvlfR0w96Df6gvNPf2hjTsKidKkQgKA0R2eBJ8yl/ptfgexU9/aDX4A/8Nf1aNKSU\nUlWcBgKllKriqkIgmOHrBJQBvQbfq+jpB70Gf+CX6a/0dQRKKaXcqwpPBEoppdzQQKCUUlVcpQ4E\nIjJQRHaKyG8iMt7X6SkOEWkpIitEZLuIbBORR3ydppISkUAR+VFEvvB1WkpCROqIyFwR+UVEdtjT\nrVYYIvKo/Tf0s4jMFpFQX6epKCLyjogcEpGfHdbVE5GlIrLLfq3ryzQWxcU1vGr/HW0Rkc9EpI4v\n05ir0gYCEQkE/gVcDXQAbhORDr5NVbFkAY8bYzoA3YAHK1j6HT0C7PB1IkrhDeArY8yFQDQV6FpE\npDnwMBBrjOmENfT7rb5NlUfeBQYWWDceWG6MaQ8st5f/v737C62yjuM4/v7EKtKi6EbSXSgjlBDT\nogiFLrIoShxFXVkpeScaFRZIF3VTCUV0ERVRlNAockkZUSlUFKERSiYW9EdFp7ONQhsJs9qni99v\nNZZzO8dtv53zfF8w9pxnh+d8nnGe8z3P7znn+5vKXuf/+7AdmG97AfADsGGyQ51O0xYC4FrgJ9v7\nbZ8C3gLaC2caM9vdtnfn5T7Si8+ssqlqJ6kVuI0GnQdH0sXA9cCrALZP2T5eNlXNWoALJLUA04Cj\nhfOMyvbnwG/DVrcDm/LyJmASJ56t3en2wfY223/lmztJMzQW18yFYBZweMjtLhrwhRRA0mxgEfBV\n2SR1eQ54BBgoHaROc4Be4LU8vPWKpOmlQ42V7SPAM8AhoBs4YXtb2VR1m2G7Oy8fA2aUDDMO7gM+\nLB0CmrsQNAVJFwLvAA/Y/r10nlpIWgb02N5VOstZaAGuAl60vQj4g6k/JPGvPI7eTipoM4Hpku4u\nm+rsOX3uvWE/+y7pUdLwb0fpLNDcheAIaXreQa15XcOQdC6pCHTY3lI6Tx2WAMslHSQNzd0g6Y2y\nkWrWBXTZHjwb6yQVhkZxI3DAdq/tP4EtwOLCmer1i6TLAPLvnsJ56iJpFbAMWOEp8kWuZi4EXwOX\nS5oj6TzSBbKthTONmSSRxqW/t/1s6Tz1sL3Bdqvt2aT//ye2G+rdqO1jwGFJc/OqpcB3BSPV6hBw\nnaRp+Tm1lAa62D3MVmBlXl4JvFcwS10k3UIaKl1u+2TpPIOathDkCzJrgY9JT/y3be8rm6omS4B7\nSO+iv8k/t5YOVVHrgA5J3wILgScL5xmzfCbTCewG9pKO+SnZ5mAoSW8CO4C5krokrQY2AjdJ+pF0\nprOxZMbRjLAPzwMXAdvzMf1S0ZBZtJgIIYSKa9ozghBCCGMThSCEECouCkEIIVRcFIIQQqi4KAQh\nhFBxUQhC5eXuomvy8kxJnRPwGI9LWj/e2w1hPEQhCAEuAdYA2D5q+87CeUKYVFEIQkhfTGrLX/DZ\nPNg/XtIqSe/m3vcHJa2V9FBuPrdT0qX5fm2SPpK0S9IXkuaN8DhXSPpM0n5J90/WzoUwmigEIaQm\ncj/bXgg8POxv84E7gGuAJ4CTufncDuDefJ+XgXW2rwbWAy+M8DjzgJtJLdIfy72kQiiupXSAEKa4\nT/N8EH2STgDv5/V7gQW5O+xiYHNq5QPA+SNs6wPb/UC/pB5SG+WuiYsewthEIQjhzPqHLA8MuT1A\nOn7OAY7ns4latvU3cfyFKSKGhkKAPlIjsJrlOSIOSLoLUtdYSVfm5dslPTV+MUOYGFEIQuXZ/hX4\nMl8kfrqOTawAVkvaA+zjvylR24CGmkwoVFN0Hw1hguRJeB603Vs6SwhnEoUghBAqLoaGQgih4qIQ\nhBBCxUUhCCGEiotCEEIIFReFIIQQKi4KQQghVNw/GQNvO2ksQ0MAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt \n", "import numpy as np\n", "from scipy import integrate\n", "%matplotlib inline\n", "\n", "y0 = 0.1 \n", "ug = 0.5 \n", "tot = 12.5 \n", "stepsize = 10\n", "\n", "\n", "## using simple RK method\n", "t = np.linspace(0,tot,stepsize) # creating a time vector\n", "t2 = np.linspace(0,tot,stepsize*10) # creating a time vector with many more timesteps\n", "dt = t[1]-t[0] # establishing dt\n", "dt2 = t2[1]-t2[0] # establishing dt2\n", "y = np.zeros(len(t)) # creating a Cell density vector \n", "y1 = np.zeros(len(t2)) # creating a Cell density vector \n", "\n", "y[0] = y0\n", "y1[0] = y0\n", "##### Creating the differential integration with 10 stepsize and the analytical solution\n", "for i in range (1, len(t)):\n", " y[i] = y[i-1] + ug*y[i-1]*dt\n", " X = np.exp(ug*(t))*y0 #analytical solution\n", "##### Creating the differential integration with 100 stepsize\n", "for i in range (1, len(t2)):\n", " y1[i] = y1[i-1] + ug*y1[i-1]*dt2\n", " \n", "# plotting both in comparison\n", "\n", "plt.xlabel('time,h')\n", "plt.ylabel('Cell density (g/L)')\n", "plt.title ('Comparing Analytical solution to RK')\n", "plt.plot(t,y,'-x',label='Numerical with 10 timesteps')\n", "plt.plot(t2,y1,'-og',label='Numerical with 100 timesteps')\n", "plt.plot(t,X,'-or',label = 'Analytical')\n", "plt.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see from the graph above, that earlier values both numerical solutions are very close to the actual solution. As the stepsizes increase, the value approach that of the analytical solution though there is still an offset. In a more complicated system, this offset may be much larger than that presented here. \n", "\n", "More advanced ODE solvers are used to eliminate this offset (e.g ode45). We will not learn to code these more advanced sovlers and rather, we will learn how to implement them into the equations we want to analyze." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example of a single ODE: \n", "\n", "Suppose Ace Chemical Engineering (ACE) has a 100 L tank that has 50 g/L of NaCl. An inlet stream flows at 15 L/min and has a concentration of 2.5 g/L. ACE wants to empty the tank at 15L/min while keeping the inlet flow rate constant. How many minutes will it take for the tank to reach 30 g/L? \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The concentration balance with no reaction is written out as:\n", "\n", "$$ V\\frac{dC}{dt} = Q_{in}C_{in} - Q_{out}C_{out} $$ Where $ C_{out} = C $ \n", "\n", "The known values are substituted in and the formula is rearranged to isolate $\\frac{dC}{dt} $\n", "\n", "$$ \\frac{dC}{dt} = \\frac{15*2.5}{100} -\\frac{15*C}{100}$$\n", "\n", "This can now be put into OdeInt to obtain a solution.\n", "\n", "OdeInt has the following syntax: \n", "\n", "> integrate.odeint(function,inital_condition,time_vector)\n", "\n", "Where the \"function\" arguement is an actual python function returning us the change in the value desired. In this case the function returns us the change in concentration at timestep $t_i$. odeint returns an array of values that is equal in size to the time vector." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEKCAYAAAAIO8L1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmYlXX9//Hn+8w+wzrDgOyLosi+jIogZRJGVmjm1tcF\n/VZ8y9LM0sj6ZnblL0rLvHL7khuaZe6YLaaAu5FsguyggCDLsO+znffvj3MPDnhm5jAz59xnhtfj\nus51r+e+39xcM6+5t8/H3B0REZEjRcIuQERE0pMCQkRE4lJAiIhIXAoIERGJSwEhIiJxKSBERCSu\npAWEmT1oZlvM7L0a8wrN7CUzWxkM29dY9iMzW2Vmy83sc8mqS0REEpPMM4iHgfFHzJsMzHD3vsCM\nYBoz6w9cAgwIvnOPmWUksTYREalH0gLC3V8Dth8x+1xgWjA+DTivxvzH3b3M3T8AVgGnJqs2ERGp\nX2aK99fJ3TcG45uATsF4V+DfNdZbH8z7BDObBEwCKCgoGNGvX78klSoi0jLNnTt3q7sX17deqgPi\nEHd3Mzvqdj7cfSowFaCkpMTnzJnT5LWJiLRkZrY2kfVS/RTTZjPrDBAMtwTzNwDda6zXLZgnIiIh\nSXVAPA9MDMYnAtNrzL/EzHLMrDfQF/hPimsTEZEaknaJycz+DJwJdDCz9cDNwBTgCTP7GrAWuAjA\n3Reb2RPAEqAS+La7VyWrNhERqV/SAsLdv1rLorG1rH8rcGuy6hGRpldRUcH69es5ePBg2KVIHLm5\nuXTr1o2srKwGfT+0m9Qi0vytX7+e1q1b06tXL8ws7HKkBndn27ZtrF+/nt69ezdoG2pqQ0Qa7ODB\ngxQVFSkc0pCZUVRU1KizOwWEiDSKwiF9Nfb/RgEhIiJxKSBERJrAmjVr+NOf/nTU39u5cyf33HPP\noemPPvqICy64oClLazAFhIhIE6grICorK2v93pEB0aVLF5566qkmr68hFBAi0qw98sgjDB48mCFD\nhnD55ZezZs0azjrrLAYPHszYsWNZt24dAFdeeSXXXnsto0aNok+fPof9Ev7Vr37FoEGDGDJkCJMn\nTwZg9erVjB8/nhEjRjBmzBiWLVtW53YmT57M66+/ztChQ7njjjt4+OGHmTBhAmeddRZjx45l7969\njB07luHDhzNo0CCmT59+6HurV69m6NCh3HDDDaxZs4aBAwcCsYcArrrqKgYNGsSwYcOYNWsWAA8/\n/DDnn38+48ePp2/fvtx4441JObZ6zFVEmsQtf13Mko92N+k2+3dpw81fGlDr8sWLF/OLX/yCt956\niw4dOrB9+3YmTpx46PPggw9y7bXX8txzzwGwceNG3njjDZYtW8aECRO44IIL+Mc//sH06dOZPXs2\n+fn5bN8ea4R60qRJ3HffffTt25fZs2dz9dVXM3PmzFq3M2XKFG6//XZeeOEFIPZLfN68eSxcuJDC\nwkIqKyt59tlnadOmDVu3bmXkyJFMmDCBKVOm8N5777FgwQIgdiZS7e6778bMWLRoEcuWLePss89m\nxYoVACxYsID58+eTk5PDSSedxDXXXEP37jVbLGo8BYSINFszZ87kwgsvpEOHDgAUFhby9ttv88wz\nzwBw+eWXH/bX9XnnnUckEqF///5s3rwZgJdffpmrrrqK/Pz8Q9vYu3cvb731FhdeeOGh75aVldW5\nnXjGjRtHYWEhEHsv4aabbuK1114jEomwYcOGOr8L8MYbb3DNNdcA0K9fP3r27HkoIMaOHUvbtm0B\n6N+/P2vXrlVAiEh6qusv/XSRk5NzaNy99sako9Eo7dq1O/RXfUO3U1BQcGj8scceo7S0lLlz55KV\nlUWvXr0a9Y5CzRoyMjLqvM/RULoHISLN1llnncWTTz7Jtm3bANi+fTujRo3i8ccfB2K/lMeMGVPn\nNsaNG8dDDz3E/v37D22jTZs29O7dmyeffBKIhcC7775b53Zat27Nnj17al2+a9cuOnbsSFZWFrNm\nzWLt2rX1fm/MmDE89thjAKxYsYJ169Zx0kkn1VlHU1JAiEizNWDAAH784x/z6U9/miFDhnD99dfz\n+9//noceeojBgwfz6KOPcuedd9a5jfHjxzNhwgRKSkoYOnQot99+OxALlwceeIAhQ4YwYMCAQzeV\nazN48GAyMjIYMmQId9xxxyeWX3rppcyZM4dBgwbxyCOPUN3ZWVFREaNHj2bgwIHccMMNh33n6quv\nJhqNMmjQIC6++GIefvjhw84cks3qOj1Kd+owSCRcS5cu5eSTTw67DKlDvP8jM5vr7iX1fVdnECIi\nEpcCQkRE4lJAiEijNOfL1C1dY/9vFBAi0mC5ubls27ZNIZGGqvuDyM3NbfA29B6EiDRYt27dWL9+\nPaWlpWGXInFU9yjXUAoIEWmwrKysBvdWJulPl5hERCQuBYSIiMSlgBARkbgUECIiEpcCQkRE4lJA\niIhIXAoIERGJSwEhIiJxKSBERCQuBYSIiMSlgBARkbgUECIiEpcCQkRE4lJAiIhIXKEEhJl9z8wW\nm9l7ZvZnM8s1s0Ize8nMVgbD9mHUJiIiMSkPCDPrClwLlLj7QCADuASYDMxw977AjGBaRERCEtYl\npkwgz8wygXzgI+BcYFqwfBpwXki1iYgIIQSEu28AbgfWARuBXe7+L6CTu28MVtsEdIr3fTObZGZz\nzGyOujkUEUmeMC4xtSd2ttAb6AIUmNllNdfxWA/ocXtBd/ep7l7i7iXFxcVJr1dE5FgVxiWmzwIf\nuHupu1cAzwCjgM1m1hkgGG4JoTYREQmEERDrgJFmlm9mBowFlgLPAxODdSYC00OoTUREApmp3qG7\nzzazp4B5QCUwH5gKtAKeMLOvAWuBi1Jdm4iIfCzlAQHg7jcDNx8xu4zY2YSIiKQBvUktIiJxKSBE\nRCQuBYSIiMSlgBARkbgUECIiEpcCQkRE4lJAiIhIXAoIERGJSwEhIiJxKSBERCQuBYSIiMSlgBAR\nkbgUECIiEpcCQkRE4lJAiIhIXAoIERGJSwEhIiJxNSggzGxdUxciIiLppaFnENakVYiISNppaEB4\nk1YhIiJpJ7O2BWZ2fW2LgFbJKUdERNJFrQEBtK5j2Z1NXYiIiKSXugJiJfCiu29LVTEiIpI+6gqI\n7sCTZpYFzAD+AfzH3XX/QUTkGFDrTWp3/5W7nwWcA7wL/Dcwz8z+ZGZXmFmnVBUpIiKpV9cZBADu\nvgd4NvhgZv2BzwOPAJ9LanUiIhKaegPCzIbHmf0culEtItKi1RsQwD3AcGAhsUdcBwKLgbZm9i13\n/1cS6xMRkZAk8qLcR8Awdy9x9xHAMOB9YBzw62QWJyIi4UkkIE5098XVE+6+BOjn7u8nrywREQlb\nIpeYFpvZvcDjwfTFwBIzywEqklaZiIiEKpEziCuBVcB1wef9YF4F8JlkFSYiIuFK5DHXA8Bvgs+R\n9jZkp2bWDrif2A1vJ/aOxXLgL0AvYA1wkbvvaMj2RUSk8cLqMOhO4J/u3g8YAiwFJgMz3L0vsTe3\nJ4dUm4iIEEJAmFlb4FPAAwDuXu7uO4FzgWnBatOA81Jdm4iIfCyMM4jeQCnwkJnNN7P7zawA6OTu\nG4N1NgFxm/Iws0lmNsfM5pSWlqaoZBGRY89RB4SZ/T8z+6GZFTVwn5nEXry7192HAfs44nJS0CBg\n3EYB3X1q8E5GSXFxcQNLEBGR+jTkDOI/QCVwRwP3uR5Y7+6zg+mniAXGZjPrDBAMtzRw+yIi0gQS\neQ/iMO7+XGN26O6bzOxDMzvJ3ZcDY4ElwWciMCUYTm/MfkREpHESaayvGPgGscdPD63v7v/diP1e\nAzxmZtnE3qu4itjZzBNm9jVgLXBRI7YvIiKNlMgZxHTgdeBloKopduruC4CSOIvGNsX2RUSk8RIJ\niHx3/2HSKxERkbSSyE3qF8zsnKRXIiIiaSWRgPgusZA4aGZ7gs/uZBcmIiLhSqQtptapKERERNJL\nQo+5mtkEYs1jALzi7i8kryQREUkH9V5iMrMpxC4zVb+r8F0z+2WyCxMRkXAlcgZxDjDU3aMAZjYN\nmA/8KJmFiYhIuBJtaqNdjfG2yShERETSSyJnEL8E5pvZLMCI3YtQXw0iIi1cIk8x/dnMXgFOCWb9\n0N03JbUqEREJXa2XmMysXzAcDnQmaIUV6BLMExGRFqyuM4jrgUnE74vagbOSUpGIiKSFWgPC3ScF\no59394M1l5lZblKrEhGR0CXyFNNbCc5LubhdzomISJOo9QzCzI4DugJ5ZjaM2BNMAG2A/BTUVq8V\nm/awt6ySVjlH3e+RiIjUo67frJ8DrgS6Ab+tMX8PcFMSa0pYeVWUe19ZxQ2f6xd2KSIiLU5d9yCm\nAdPM7Cvu/nQKa0pYu7ws7n/9Ay49rSdd2uWFXY6ISIuSyHsQT5vZF4ABQG6N+T9PZmGJ6NQ2l4PA\nbS8u546Lh4ZdjohIi5JIY333ARcT60fagAuBnkmuKyHZGRG+fkZvnp2/gXc/3Bl2OSIiLUoiTzGN\ncvcrgB3ufgtwOnBicstK3LfOPJ4OrbK59W9LcddzTSIiTSWRgKh+B2K/mXUBKoi9WZ0WWudm8b1x\nJ/KfNdt5cfHmsMsREWkxEgmIv5pZO+A2YB6wBvhTMos6WheXdKdvx1b88h9LOVhRFXY5IiItQp0B\nYWYRYIa77wyeZOoJ9HP3n6akugRlZkT46Zf6s3bbfv7w2vthlyMi0iLUGRBBJ0F315guc/ddSa+q\nAcb0LeYLgzpz16xVfLh9f9jliIg0e4lcYpphZl8xM6t/1XD95IsnkxExfv7CkrBLERFp9hIJiP8B\nngTKzGy3me0xs91JrqtBOrfN49qxfXlpyWZmLtMNaxGRxqg3INy9tbtH3D3b3dsE021SUVxD/Pfo\n3hxfXMDPnl+iG9YiIo2QyItyMxKZly6yMyP8/NyBrNu+n/teXR12OSIizVZdPcrlmlkh0MHM2ptZ\nYfDpRayV17Q1+oQOfGlIF+6ZtZpVW/aEXY6ISLNU1xnE/wBzgX7BsPozHbgr+aU1zs1f6k9+TgY3\nPrWQqqjesBYROVq1BoS73+nuvYEfuHsfd+8dfIa4e9oHRIdWOfz0i/2Zt24nj769JuxyRESanURa\nc/29mY0CetVc390fSWJdTeLLw7oyfcFH/PrF5Xy2fye6tU+Lfo5ERJqFRG5SPwrcDpwBnBJ8Shq7\nYzPLMLP5ZvZCMF1oZi+Z2cpg2L4J9sGtXx4IwE3PvqfG/EREjkIi70GUAKPd/Wp3vyb4XNsE+/4u\nsLTG9GRizXr0BWYE043WrX0+Pxzfj9dWlPL0vA1NsUkRkWNCIgHxHnBcU+7UzLoBXwDurzH7XGBa\nMD4NOK+p9nf5yJ6c0qs9tzy/mA07DzTVZkVEWrREAqIDsMTMXjSz56s/jdzv74AbgWiNeZ3cfWMw\nvgno1Mh9HBKJGL+5cChRd77/xAKieqpJRKRe9d6kBn7WlDs0sy8CW9x9rpmdGW8dd3czi/tb3Mwm\nAZMAevTokfB+exTlc/OXBnDj0wt58M0P+PqYPkdfvIjIMSSRpjZeJdYHRFYw/g6xfiEaajQwwczW\nAI8DZ5nZH4HNZtYZIBhuqaWeqe5e4u4lxcXFR7XjC0u6Ma5/J379z+Us36QX6ERE6pLIU0zfAJ4C\n/i+Y1RV4rqE7dPcfuXs3d+8FXALMdPfLgOeBicFqE4m9kNekzIxfnj+INnmZXPeXBZRVqq0mEZHa\nJHIP4tvE/urfDeDuK4GOSahlCjDOzFYCnw2mm1yHVjlMOX8wSzfu5rZ/Lk/GLkREWoRE7kGUuXt5\ndXcQZpYJNMldXnd/BXglGN8GjG2K7dbns/07ccXpPbn/jQ84rU8R4/o32f1wEZEWI5EziFfN7CYg\nz8zGEesb4q/JLSv5bjrnZAZ2bcMPnnyX9TvUA52IyJESCYjJQCmwiFgDfn8HfpLMolIhNyuDu746\nnKqoc82f51NRFa3/SyIix5BEAiIPeNDdL3T3C4AHg3nNXq8OBUz5yiDmr9vJbS/qfoSISE0J9UnN\n4YGQB7ycnHJS74uDu3DZyB5Mfe19Xly8KexyRETSRiIBkevue6sngvEW1SzqT77QnyHd2nL9Xxaw\ncrPejxARgcQCYp+ZDa+eMLMRQItq0Cg3K4P7Lh9BXnYGkx6dy64DFWGXJCISukQC4jrgSTN73cze\nAP4CfCe5ZaVe57Z53HPpCD7cvp/rHp+vXuhE5JiXSFMb7xDrdvRbwDeBk919brILC8OpvQu5ecIA\nZi0v5Y6XVoRdjohIqBJ5UQ5inQT1CtYfbmbNoke5hrjstB4s3rCLu2at4oSOrThvWNewSxIRCUW9\nARH0KHc8sACobrzIgRYZEGbGz88dyJpt+7jxqYV0bpvLaX2Kwi5LRCTlrL5uOM1sKdDf07C/zpKS\nEp8zZ05Str1rfwXn3/smW/eW88zVozi+uFVS9iMikmpmNtfd6+06OpQe5ZqDtvlZPHTlqWRGjKse\neodte8vCLklEJKXC6lGuWehRlM/9E0vYvPsgX39kDvvLK8MuSUQkZVLeo1xzM6xHe+68ZBhXPzaX\nb/5xHvdfUUJ2ZiK5KiLSvCXao9wyoHXwWRrMO2aMH3gcU84fzGsrSvneXxboHQkROSYk0qPcRcB/\ngAuBi4DZZnZBsgtLNxed0p0fn3Myf1u0kZ88t4g0vGcvItKkErnE9GPgFHffAmBmxcQa63sqmYWl\no298qg+7DlRw16xVtMnLYvL4flR3pCQi0tIkEhCR6nAIbCOxm9st0vfPPpFdByr4v1ffJysS4ftn\nn6iQEJEWKZGA+KeZvQj8OZi+GPhH8kpKb2bGLRMGUBmNctesVTjOD84+SSEhIi1OvQHh7jeY2fnA\nGcGsqe7+bHLLSm+RiHHreYMA4+5Zq4k63Pg5hYSItCy1BoSZnQB0cvc33f0Z4Jlg/hlmdry7r05V\nkekoFhIDMYN7X1mNO/xwvEJCRFqOuu4l/A7YHWf+rmDZMS8SMX5x7kAuPa0H9726mpufX0xUj8CK\nSAtR1yWmTu6+6MiZ7r7IzHolraJmJhIxfnHeQApyMpn62vvs2F/Bby4copfpRKTZqysg2tWxLK+O\nZcccM+Omc06msCCbKf9Yxq4DFdx32XDysxNtTV1EJP3U9WfuHDP7xpEzzezrQIvsMKixvvnp4/nV\nVwbxxspS/usPs9mxrzzskkREGqyuP3GvA541s0v5OBBKgGzgy8kurLm6+JQetM3L5trH5/Ple97k\ngStPUVPhItIs1XoG4e6b3X0UcAuwJvjc4u6nu/um1JTXPI0feBx//sZp7DlYyfn3vMXbq7eFXZKI\nyFFLpLG+We7+++AzMxVFtQQjehby3LdHU9w6h8sfmM0Tcz4MuyQRkaOiR22SqHthPk9/axSnH1/E\njU8t5Na/LaGyKhp2WSIiCVFAJFnbvCwevPIUrji9J394/QMue2A2pXvUO52IpD8FRApkZUT4+bkD\n+e1FQ5i/bidf+v0bzFu3I+yyRETqpIBIofOHd+OZq0eRlWlc/H9v8+jba9SvhIikLQVEig3o0pYX\nvjOGM07owP9OX8y3/jiPnfv1voSIpJ+UB4SZdTezWWa2xMwWm9l3g/mFZvaSma0Mhu1TXVuqtM3P\n4oGJp3DTOf2YsWwz43/3uh6FFZG0E8YZRCXwfXfvD4wEvm1m/YHJwAx37wvMCKZbrEjEmPSp43nm\nW6PJy87gv+7/N7e/uJwKPeUkImki5QHh7hvdfV4wvgdYCnQFzgWmBatNA85LdW1hGNStLS9ccwYX\nDO/GXbNW8eV73mTpxniN6IqIpFao9yCCVmGHAbOJtR67MVi0CehUy3cmmdkcM5tTWlqakjqTrSAn\nk9suHMK9lw5n066DTLjrDX738grKK3U2ISLhCS0gzKwV8DRwnbsf9iezxx7tift4j7tPdfcSdy8p\nLi5OQaWp8/lBnfnX9z7NOYM687uXVzLhrjd4b8OusMsSkWNUKAFhZlnEwuGxoLc6gM1m1jlY3hnY\nEkZtYSssyObOS4bxhytK2L6vnHPvfpNfvLCEvWWVYZcmIseYMJ5iMuABYKm7/7bGoueBicH4RGB6\nqmtLJ+P6d+Kl732ai0q68cCbHzD2N6/w13c/0nsTIpIylupfOGZ2BvA6sAiovsh+E7H7EE8APYC1\nwEXuvr2ubZWUlPicOXOSWG16mL9uB/87/T3e27Cb0ScUccuEgZzQUU2Ii0jDmNlcdy+pd73m/Bfp\nsRIQAFVR50+z13Lbi8vZX17Fpaf14NqxfSlqlRN2aSLSzCQaEHqTupnIiBiXn96LmT84k4tP6c4f\nZ6/jzNte4e5ZqzhYURV2eSLSAikgmpkOrXK49cuDePG6MZzWp5DbXlzOWbe/wpNzPlRT4iLSpBQQ\nzdQJHVtz/8RTeHzSSDq0zuGGpxby2d++ytNz1ysoRKRJKCCauZF9ipj+7dFMvXwE+dmZfP/Jdxl3\nx2s8M09BISKNo5vULUg06vxryWZ+9/IKlm3aQ6+ifL4+pg8XjOhGblZG2OWJSJrQU0zHsGjUeXHx\nJu59dTUL1++iqCCbK07vxeWn96SwIDvs8kQkZAoIwd2Z/cF2pr72PjOXbSE3K8KFI7ozcVQvvUch\ncgxLNCAyU1GMhMPMGNmniJF9ili5eQ9/eP19/vLOhzz677Wc3qeIy0b25OwBncjK0K0oEfkknUEc\nY7buLeOJOR/yp9nrWL/jAMWtc/jqKd25+NQedG2XF3Z5IpICusQkdaqKOq+u2MIf/72OWctj7SKe\n3qeIrwzvxviBx1GQo5NLkZZKASEJ+3D7fp6Zt4Gn561n3fb95Gdn8PmBnfnK8K6M7FNEJGJhlygi\nTUgBIUfN3ZmzdgfPzFvPC+9uZE9ZJR1b5/D5gcdxzqDOlPQqJENhIdLsKSCkUQ5WVPHSks38beFG\nZi3fQllllOLWOYwfEAuLU3srLESaKwWENJl9ZZXMXLaFvy+KhcXBiijt8rM488RiPtOvI58+sZh2\n+Xq/QqS5UEBIUuwvr2TWslJmLNvMK8tL2b6vnIyIMaJHez7TryOf6VfMSZ1aE+sXSkTSkQJCkq4q\n6ry7ficzl25h5rItLNkY61q8Q6tsTj++A6OPL2L0CR3oXpgfcqUiUpMCQlJu464DvL5yK2+t2sqb\nq7dRuqcMgO6FeYzq04GRxxcyokch3QvzdIYhEiIFhITK3Vm1ZS9vrd7Gm6u28u/3t7H7YCUAxa1z\nGNGjPSW92jO8Z3sGdGlDTqYaExRJFTW1IaEyM/p2ak3fTq2ZOKoXVVFnxeY9zFm7g3lrdzBn7Xb+\nuXgTANmZEQZ3bcvArm0ZFAyPLy4gU02AiIRKZxASmi27DzJ37Q7mrt3B/A93suSj3RwIuk/NzYrQ\nv3ObQ4Fxcuc2nNCxlZotF2kCusQkzU5V1Hm/dC+LNuxi0YZdvLdhF4s/2s3+8lhoRAx6FRXQt1Mr\nTurUmhOPa81JnVrTq0OBGhwUOQq6xCTNTkbk48tS5w/vBsRC44Ot+1i+aQ/LN+9hxaY9rNiyh5eW\nbCYa/G2TlWH0Kiqgd4fYp1eHAnoVFdCnuICOrXN0Q1ykgRQQktYyIsYJHVtxQsdWfIHOh+YfrKhi\ndeleVm7ey/LNe1i1ZS8fbN3HKytKKa/8uKvV/OwMehYV0LtDPj2LCujaLo+u7fPoFgzzs/UjIFIb\n/XRIs5SblcGALm0Z0KXtYfOros5HOw+wZts+1mzdxwdb97Nm2z6WbYyddVRUHX5JtX1+Fl3b58WC\no11+MJ5LcetcOrXJobh1jp6wkmOWAkJalIyI0b0wn+6F+YzpW3zYsqqoU7qnjA0797N+xwE27DwQ\nG+44wOrSfby2Yuuhm+Q1tc/PolObXIpb59CpTSw4OgYB0qFVDu0LsikqyKZNbpZavpUWRQEhx4yM\niHFc21yOa5vLiJ6fXO7u7NhfwcZdB9iyu4wtew6yueZw90FWbdnLlj1lVEU/+XBHRsRon59NYUEW\nhQXZFBXk0L4gi8KCHIoKsmlfkE27vCza5GXRNi+LNrmZtMnL0g12SVsKCJGAmVFYkE1hQTYDutS+\nXjTqbN9fzubdB9m2t5zt+8rZtq+cHcFw+74ytu8rZ9mm3WzfV87OAxXU9bBgfnYGbXKzaJOXGQTH\n4SHSOjeL/JwMWuVkkp+dSUFOBgXZmRTkBOM5mRRkZ6p1XWlyCgiRoxSJGB1axS4vJaKyKsrOAxVs\n31fO7gMV7D5Ywa4DFew+UBkMD5+3cddBlm/ew+4DFewpq6wzXGrKzYocCo787CBQcjLJzYyQm5VB\nblZsmJeVQU71dGbGoWV5WbHxnGC93MwM8rJjy3IyM8jOjJCVYWRnRPRk2DFCASGSZJkZkaMKlJqq\nos6+8kr2l1Wxr7ySfWWV7Curig3LY+P7yyvZW1bJ/vJgflkl+4LxXQcq2FJRxcGKKg5UVHGwIsrB\niirKajzp1RBZGUZWRiQIjQjZh8bt0LysjAg5h8aN7MyMQwFTvU5mxMiIWDCMkJlhh8/7xDpGZiRy\naDq2fjCdUft6GRHDDCJ2+HjEjAwzLMLH4zXWixjHdBgqIETSWEbEYpeccrOadLvRqFNeFf1EcBys\nHq+s4mB5VWwYLKuoilJR5ZRXRimvilJRGaWiKjZeXumx8cPmRdlXVhmsG1teFiyvXrfKnaqof+Lp\nsnRTMywyLDYeCcIktsxiyyIfB0+k5niN9SNB4FRPV88zYts3AztsefWyj+cVZGcy9Yp633NrNAWE\nyDEoEjFyI7FLSu3CLobYAwJRh8polKqoUxl1qqqCYdQPnx+NhU3N6cqq+OtVRp3KqihRh6g70ah/\nPB5vOt56Ncar3HGPBeyh8eC7VdHqf8cR4zW24w5ObHvusXW8xr+/evzQetEa60ehiti/JWKffNou\nGRQQIhK62F/mkBHROyfpJO2erzOz8Wa23MxWmdnksOsRETlWpVVAmFkGcDfweaA/8FUz6x9uVSIi\nx6a0CgjgVGCVu7/v7uXA48C5IdckInJMSreA6Ap8WGN6fTDvEDObZGZzzGxOaWlpSosTETmWpFtA\n1Mvdp7orj3E5AAAGB0lEQVR7ibuXFBcX1/8FERFpkHQLiA1A9xrT3YJ5IiKSYukWEO8Afc2st5ll\nA5cAz4dck4jIMSmt3oNw90oz+w7wIpABPOjui0MuS0TkmJRWAQHg7n8H/h52HSIix7p0u8QkIiJp\nQgEhIiJxKSBERCQuBYSIiMSlgBARkbjME+3PMA2Z2R5gedh11KEDsDXsIuqg+hpH9TVcOtcGLb++\nnu5eb1MUafeY61Fa7u7J71apgcxsjuprONXXOOlcXzrXBqqvmi4xiYhIXAoIERGJq7kHxNSwC6iH\n6msc1dc46VxfOtcGqg9o5jepRUQkeZr7GYSIiCSJAkJEROJqtgFhZuPNbLmZrTKzyWHXcyQzW2Nm\ni8xsgZnNSYN6HjSzLWb2Xo15hWb2kpmtDIbt06y+n5nZhuAYLjCzc0KqrbuZzTKzJWa22My+G8xP\ni+NXR33pcvxyzew/ZvZuUN8twfx0OX611ZcWxy+oJcPM5pvZC8F0So5ds7wHYWYZwApgHLF+q98B\nvuruS0ItrAYzWwOUuHtavGxjZp8C9gKPuPvAYN6vge3uPiUI2fbu/sM0qu9nwF53vz2MmmrU1hno\n7O7zzKw1MBc4D7iSNDh+ddR3Eelx/AwocPe9ZpYFvAF8Fzif9Dh+tdU3njQ4fgBmdj1QArRx9y+m\n6me3uZ5BnAqscvf33b0ceBw4N+Sa0pq7vwZsP2L2ucC0YHwasV8qoailvrTg7hvdfV4wvgdYCnQl\nTY5fHfWlBY/ZG0xmBR8nfY5fbfWlBTPrBnwBuL/G7JQcu+YaEF2BD2tMryeNfiACDrxsZnPNbFLY\nxdSik7tvDMY3AZ3CLKYW15jZwuASVGiXwKqZWS9gGDCbNDx+R9QHaXL8gkskC4AtwEvunlbHr5b6\nID2O3++AG4FojXkpOXbNNSCagzPcfSjweeDbwSWUtOWxa41p81dT4F6gDzAU2Aj8JsxizKwV8DRw\nnbvvrrksHY5fnPrS5vi5e1Xw89ANONXMBh6xPNTjV0t9oR8/M/sisMXd59a2TjKPXXMNiA1A9xrT\n3YJ5acPdNwTDLcCzxC6LpZvNwfXr6uvYW0Ku5zDuvjn4wY0CfyDEYxhcm34aeMzdnwlmp83xi1df\nOh2/au6+E5hF7Pp+2hy/ajXrS5PjNxqYENzTfBw4y8z+SIqOXXMNiHeAvmbW28yygUuA50Ou6RAz\nKwhuFmJmBcDZwHt1fysUzwMTg/GJwPQQa/mE6h+AwJcJ6RgGNzEfAJa6+29rLEqL41dbfWl0/IrN\nrF0wnkfs4ZJlpM/xi1tfOhw/d/+Ru3dz917Efs/NdPfLSNWxc/dm+QHOIfYk02rgx2HXc0RtfYB3\ng8/idKgP+DOx0+QKYvdsvgYUATOAlcDLQGGa1fcosAhYGPxAdA6ptjOIncIvBBYEn3PS5fjVUV+6\nHL/BwPygjveAnwbz0+X41VZfWhy/GnWeCbyQymPXLB9zFRGR5Guul5hERCTJFBAiIhKXAkJEROJS\nQIiISFwKCBERiSsz7AJEwmJm1Y8KAhwHVAGlwfR+dx8VUl0/B15z95fD2L9INT3mKkL6tBwrkk50\niUkkDjPbGwzPNLNXzWy6mb1vZlPM7NKg/4BFZnZ8sF6xmT1tZu8En9FxtnmlmT0XtN+/xsy+Y2bX\nB+38/9vMCoP1HjazC4LxNWZ2i5nNC/bXL5XHQY5tCgiR+g0BvgmcDFwOnOjupxJrfvmaYJ07gTvc\n/RTgKxzeNHNNA4n1g3AKcCuxS1nDgLeBK2r5zlZ3H06s8bgfNP6fI5IY3YMQqd87HjStbGargX8F\n8xcBnwnGPwv0jzWLBEAbM2vlH/czUG2Wx/ps2GNmu4C/1tjW4Fr2X9044Fxi4SKSEgoIkfqV1RiP\n1piO8vHPUAQY6e4Hm2BbtX2nqo51RJqcLjGJNI1/8fHlJsxsaDA81cweCa0qkUZQQIg0jWuBkqD3\nsSXE7lkA9AAOhFeWSMPpMVeRJDKz24BH3X1h2LWIHC0FhIiIxKVLTCIiEpcCQkRE4lJAiIhIXAoI\nERGJSwEhIiJxKSBERCSu/w9Z/xReKR/dQQAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "############Constants for the Differential equation\n", "V = 100 # L \n", "Q = 15 #flow in and out\n", "\n", "########### Intial conditions and vectors for time\n", "time_scale =40\n", "t = np.linspace(0,time_scale,1000) # Time vector is implemented\n", "C0 = 50 # inital condition of concentration in g/L\n", "\n", "###########function for the differential equation\n", "def dcdt(C,t):\n", " return (Q/V*2.5-Q/V*C)\n", "\n", "########### Implementing the ODE solver and storing the solution as a vector\n", "sol = integrate.odeint(dcdt,C0,t)\n", "\n", "\n", "# print(t)\n", "plt.plot(t,sol,label=('concentration'))\n", "plt.axis([0 ,42,-5,100])\n", "plt.legend()\n", "plt.xlabel('Time,min')\n", "plt.ylabel('Concentration, g/L')\n" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(array([40, 41, 42, 43]), array([0, 0, 0, 0]))\n", "[ 39.85584786 39.63216078 39.40981316 39.18879699]\n", "It will take approximately 1.6 minutes to reach 40 g/L\n" ] } ], "source": [ "#rudimentarily, the time that it will take to reach 40 g/L is presented here.\n", "h = np.where((sol>=39) & (sol<=40))\n", "print (h)\n", "print(sol[h])\n", "print(\"It will take approximately {} minutes to reach 40 g/L\".format(h[0][0]*time_scale/len(t)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "From this example, we can see that the concentration drops following an exponential curve. This can be easily shown through an analytical solution in the form of $$ C = C_0 *\\exp{(t-t_0)} $$ However, this is a simplistic case that ignores the change of volume of the tank. This is largely in part due to the constant volume operation that ACE is employing. Note: the axis of the graph above is set to reflect the answer in the next part" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example of a system of ODEs:\n", "\n", "\n", "Now suppose ACE wants to begin draining this tank while further diluting the system. To accomplish this, they introduce another inlet stream of pure water at __10 L/min__ (total of __25 L/min__ in) and open the outlet stream such that it flows at __30L/min__. Now, how long will it take for the tank to drain? What will the concentration profile look like? \n", "\n", "Now this is an interesting case as the volume is no longer constant and must be taken into account. Analytically, this can be done through linear algebra and eigenvectors/eigenvalues. \n", "\n", "The volume portion is rather easy as\n", "$$ \\frac{dV}{dt} = \\sum{Q_{in}}-\\sum{Q_{out}} $$ \n", "From the given information, the change in volume per minute is simply $$ 25 - 30 = -5\\frac{L}{min}$$\n", "Thus in 20 minutes, the tank will be empty. \n", "\n", "\n", "But now, the concentration equation is more difficult. Starting from the concentration balance again, ommiting the stream of pure water (Cin = 0 ) and rearranging\n", "$$ \\frac{dC}{dt} = \\frac{15*2.5}{V} -\\frac{15*C}{V}$$\n", "\n", "To answer this, we will take advantage of ```odeint's``` ability to take in arrays of ODEs to make a system of ODE's\n", "\n", "\n", "From before, \n", "> integrate.odeint(function,inital_condition,time_vector)\n", "\n", "This time, our function arguements will take both Volume and concentration and return both differentials back. Keep in mind that the function will need these values inside array objects and so instead of inputting variables, array positions are used. \n", "\n", "In the following code for the differentials, the input arguement y is a vector containing the __concentration__ and __volume__. In essence, ```y[0]``` is the concentration, ```y[1]``` is the Volume.\n", "\n", "\n", "now, the differentials returned are in the order of $\\frac{dC}{dt}$ and $\\frac{dV}{dt}$ i.e. __dydt[0] = change in concentration__ and __dydt[1] = change in volume__ " ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def tank_Concentration(y,t):\n", " Q_in1 = 15\n", " Q_in2 = 10\n", " Q_out = 30\n", " Cin_1 = 2.5\n", " Cin_2 = 0\n", " \n", " \"\"\"This function outputs the following values in this order:\n", " 1. the change in concentration C per unit time t\n", " 2. the change in volume V per unit time t\"\"\"\n", " dydt = np.zeros_like(y)\n", " dydt[0] = (Q_in1 * Cin_1)/y[1] + (Q_in2 * Cin_2)/y[1] - (Q_out*y[0])/y[1]\n", " dydt[1] = Q_in1 + Q_in2 - Q_out\n", " return dydt\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The condition are then implemented into the same form as the first example using odeint and a solution vector is obtained then plotted.\n", "\n", "Note that the solution vector (sol) is in the order of the differential function order (concentration and then volume). \n" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEKCAYAAAAIO8L1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4VHX2x/H3SSMECBAIIRQNIEoLNaCCiFIUFRGQsrpL\ns6Crov7YRVG3uLu6i4p9bahIWXeliWCnKigIhCZSBOm9905yfn/cG4yYhEkykzszOa/H+8zcO+V+\ncn3CyW3nK6qKMcYYc74IrwMYY4wJTlYgjDHG5MgKhDHGmBxZgTDGGJMjKxDGGGNyZAXCGGNMjgJW\nIERkhIjsFpEfsi1LEJFpIrLWfSyf7bXHROQnEflRRK4PVC5jjDG+CeQexEig43nLhgAzVLU2MMOd\nR0TqAb8B6rufeV1EIgOYzRhjzAUErECo6mxg/3mLbwFGuc9HAV2yLf9AVU+p6gbgJ6BFoLIZY4y5\nsKgiXl+Squ5wn+8EktznVYHvsr1vq7vsV0RkADAAoEJJmqWkpEBchcCkNcaYMLRo0aK9qpp4ofcV\ndYE4R1VVRPLd50NVhwPDAdJqlNP024/DnZMhuaHfMxpjTDgSkU2+vK+or2LaJSLJAO7jbnf5NqB6\ntvdVc5flrXwKlCwP43rDiYN+jmqMMcVbUReIKUBf93lfYHK25b8RkRIiUgOoDSy44LdFREGPUXBo\nK3x0H1jjQWOM8ZtAXub6P2AecJmIbBWRO4GhQAcRWQu0d+dR1RXAOGAl8AVwv6pm+LSiiy6HDv+A\nHz+Fb18OwE9ijDHFk4Ryu++0tDRNT0939hzG94NVU6Dvx5ByldfRjCkWzpw5w9atWzl58qTXUUwO\nYmNjqVatGtHR0b9YLiKLVDXtQp/37CS1X4lA51dh1w8wvj/cOwfKVPY6lTFhb+vWrZQpU4aUlBRE\nxOs4JhtVZd++fWzdupUaNWoU6DvCp9VGbDz0HAOnj8KEOyDjrNeJjAl7J0+epEKFClYcgpCIUKFC\nhULt3YVPgQBIqgedXoJN38KMv3mdxphiwYpD8Crs/5vwKhAAjXpB2h0w9xVY9YnXaYwxJmSFX4EA\n6DgUqjSBj34P+9Z5ncYYUwxs3LiR//73v/n+3MGDB3n99dfPzW/fvp3u3bv7M1qBhWeBiCrh3B8h\nETCuL5w54XUiY0yYy6tAnD2b+znR8wtElSpVmDBhgt/zFUR4FgiA8hdDt7dh13L49A92E50xYWr0\n6NE0bNiQRo0a0bt3bzZu3Ejbtm1p2LAh7dq1Y/PmzQD069ePBx98kJYtW1KzZs1f/CP8zDPPkJqa\nSqNGjRgyZAgA69ato2PHjjRr1ozWrVuzevXqPL9nyJAhzJkzh8aNG/Piiy8ycuRIOnfuTNu2bWnX\nrh1Hjx6lXbt2NG3alNTUVCZPnnzuc+vWraNx48YMHjyYjRs30qBBA8C5CKB///6kpqbSpEkTZs2a\nBcDIkSPp1q0bHTt2pHbt2jzyyCMB2bbhcZlrbi69Dq5+BGY/C9Uvh2Z9L/wZY0yB/O3jFazcftiv\n31mvSjx/vbl+rq+vWLGCp556irlz51KxYkX2799P3759z00jRozgwQcf5KOPPgJgx44dfPPNN6xe\nvZrOnTvTvXt3Pv/8cyZPnsz8+fOJi4tj/36nCfWAAQN48803qV27NvPnz+e+++5j5syZuX7P0KFD\nGTZsGJ984pz7HDlyJIsXL+b7778nISGBs2fPMmnSJOLj49m7dy9XXHEFnTt3ZujQofzwww8sXboU\ncPZEsrz22muICMuXL2f16tVcd911rFmzBoClS5eyZMkSSpQowWWXXcbAgQOpXj17x6LCC+8CAXDN\nENi6ED4bDMmNoEpjrxMZY/xk5syZ9OjRg4oVKwKQkJDAvHnz+PDDDwHo3bv3L/667tKlCxEREdSr\nV49du3YBMH36dPr3709cXNy57zh69Chz586lR48e5z576tSpPL8nJx06dCAhIQFw7kt4/PHHmT17\nNhEREWzbti3PzwJ88803DBw4EIA6depw8cUXnysQ7dq1o2zZsgDUq1ePTZs2WYHIt4hIuPUdeOtq\nGNcH7vnaafBnjPGrvP7SDxYlSpQ49zyvLhKZmZmUK1fu3F/1Bf2eUqVKnXv+/vvvs2fPHhYtWkR0\ndDQpKSmFukche4bIyMg8z3MUVPieg8iuVEXoMRIOb4NJv4fMTK8TGWP8oG3btowfP559+/YBsH//\nflq2bMkHH3wAOP8ot27dOs/v6NChA++99x7Hjx8/9x3x8fHUqFGD8ePHA04RWLZsWZ7fU6ZMGY4c\nOZLr64cOHaJSpUpER0cza9YsNm3adMHPtW7dmvfffx+ANWvWsHnzZi677LI8c/hT8SgQANVbwHVP\nw5rP4duXvE5jjPGD+vXr88QTT9CmTRsaNWrEoEGDePXVV3nvvfdo2LAhY8aM4eWX827i2bFjRzp3\n7kxaWhqNGzdm2LBhgFNc3n33XRo1akT9+vXPnVTOTcOGDYmMjKRRo0a8+OKLv3r9t7/9Lenp6aSm\npjJ69Gjq1KkDQIUKFWjVqhUNGjRg8ODBv/jMfffdR2ZmJqmpqfTq1YuRI0f+Ys8h0MKjWZ+vVGFC\nf1g5GfpMhhpXBy6cMcXAqlWrqFu3rtcxTB5y+n/ka7O+4rMHAT839atwidOv6fCOC3/GGGOKqeJV\nIABKlIGeo+H0MWdvIuOM14mMMSYoFb8CAVCpLtz8CmyeB9Of9DqNMcYEpeJZIAAa9oDmd8G8f8PK\nKV6nMcaYoFN8CwTA9f+Eqs1g8v3W1M8YY85TvAtEVlO/iCgY2xtOH/c6kTHGBI3iXSAAylWHW9+G\n3Svh00HW1M+YMJW9CZ7xjRUIgEvaQ5tHYdn/YNFIr9MYY0xQsAKRpc0jUKstfP4IbF/idRpjjA+G\nDBnCa6+9dm7+ySef5LnnnmPw4ME0aNCA1NRUxo4d+6vPjRw5kgceeODcfKdOnfjqq68AKF26NIMH\nD6Z+/fq0b9+eBQsWcM0111CzZk2mTHEuaMnIyGDw4ME0b96chg0b8tZbbwX2B/VI+Dfr81VEJHTL\n1tRvwNcQl+B1KmNCx+dDYOdy/35n5VS4YWiuL/fq1YuHH36Y+++/H4Bx48bx6KOPMnXqVJYtW8be\nvXtp3rw5V1/te9eEY8eO0bZtW5577jm6du3Kn/70J6ZNm8bKlSvp27cvnTt35t1336Vs2bIsXLiQ\nU6dO0apVK6677jpq1KhR6B85mBRoD0JEhvk7SFAoVQF6jnLusJ50rzX1MybINWnShN27d7N9+3aW\nLVtG+fLlWbp0KbfddhuRkZEkJSXRpk0bFi5c6PN3xsTE0LFjRwBSU1Np06YN0dHRpKamnhurYerU\nqYwePZrGjRtz+eWXs2/fPtauXRuIH9FTBd2D6An80Z9Bgka1NOfy188HwzcvwNXh+WMa43d5/KUf\nSD169GDChAns3LmTXr16sWHDhgt+JioqisxsfwBmb7sdHR2NiAAQERFxrjleRETEuZbaqsqrr77K\n9ddf788fJegU9ByE+DVFsGlxNzS4FWY9Deu/8jqNMSYPvXr14oMPPmDChAn06NGD1q1bM3bsWDIy\nMtizZw+zZ8+mRYsWv/hMSkoKS5cuJTMzky1btrBgwYJ8rfP666/njTfe4MwZp1XPmjVrOHbsmN9+\npmCR6x6EiOR2AF4I9wIh4rTi2PkDTLgT7p0D8VW8TmWMyUH9+vU5cuQIVatWJTk5ma5duzJv3jwa\nNWqEiPDss89SuXLlXwzl2apVK2rUqEG9evWoW7cuTZs2zdc677rrLjZu3EjTpk1RVRITE88NaxpO\ncm33LSIbACXnYqCqWjOQwXyR73bf+bXnRxh+rXOirN8nEBkduHUZE4Ks3XfwC0i7b1Wtoao13cfz\nJ8+LQ5FIvAw6vwJbvoNpf/U6jTHGFKl8n4MQkWQRKbohjbyW2h1aDIDvXoMV4bcLaYwxuSnISeox\nwOqwvdQ1J9c9DVXTYPIDsDf8LmUzpjBCeVTKcFfY/zf5LhCq2h6oCbxXqDWHkqgY5/6IqBi3qV/4\nXa1gTEHExsayb98+KxJBSFXZt28fsbGxBf4On+6DEJGrgNqq+p6IVATKqOqKAq81FJWtBre+A2O6\nwSf/B13fcq52MqYYq1atGlu3bmXPnj1eRzE5iI2NpVq1agX+/AULhIj8FUgDLsPZa4gB/gO0KuhK\nReT/gLtwrpJaDvQH4oCxQAqwEeipqgcKuo6AqNUWrnkMvvonVL8cmt/pdSJjPBUdHR127SXMz3w5\nxNQV6AwcA1DV7UCZgq5QRKoCDwJpqtoAiAR+AwwBZqhqbWCGOx98rh7sdH/9YghsW+x1GmOMCRhf\nCsRpdQ4wKoCIlPLDeqOAkiIShbPnsB24BRjlvj4K6OKH9fhfRAR0extKJ8G4vnB8v9eJjDEmIHwp\nEONE5C2gnIjcDUwH3i7oClV1GzAM2AzsAA6p6lQgSVV3uG/bCSTl9HkRGSAi6SKS7tlxz7gEZyS6\nIzvgwwHW1M8YE5YuWCBUdRgwAZiIcx7iL6r6akFXKCLlcfYWagBVgFIi8rvz1nlujyWHPMNVNU1V\n0xITEwsao/CqNYOO/4KfpsGc573LYYwxAeLTVUyqOk1E5me9X0QSVLWgx1baAxtUdY/7XR8CLYFd\nIpKsqjtEJBnYXcDvLzrN74It852mftXSoNa1Xicyxhi/ueAehIjcIyI7ge+BdGCR+1hQm4ErRCRO\nnJ667YBVwBSgr/uevsDkQqyjaIjAzS87LTkm3gmHtnmdyBhj/MaXcxB/BBqoakq23kwF7sWkqvNx\nDlktxrnENQIYDgwFOojIWpy9DG+ay+dXTCnoOQbOnoLx/eDsaa8TGWOMX/hSINYBx/25UlX9q6rW\nUdUGqtpbVU+p6j5VbaeqtVW1fSEOYRW9xEuh86uwdQFM+4vXaYwxxi98OQfxGDDXPQdxKmuhqj4Y\nsFShqEE32LIA5r8B1Vs488YYE8J8KRBvATNxDgfZ9Zx56fB32LYIpgyEpAbOnoUxxoQoXwpEtKoO\nCniScBAVAz1GwlutYVxvuGsGlCjtdSpjjCkQX85BfO7enJYsIglZU8CThaqyVeHWd53R6D55GKzL\npTEmRPmyB3Gb+/hYtmWK0/Lb5KTWtXDtEzDrKaepX4u7vU5kjDH5dsECoarWqrEgWv/Buarpi8eg\nSlPnzmtjjAkhvrT77pPTclUd7f84YSQiwhkz4q02ML4v3DPb6eFkjDEhwpdzEM2zTa2BJ3Haf5sL\niUtwRqI7ugs+vNua+hljQoovh5gGZp8XkXLABwFLFG6qNoWOQ+HTQTD7ObjmUa8TGWOMT/I9JjXO\nwEF2XiI/0u6Ahr3gq3/BTzO8TmOMMT7x5RzEx/zcejsCqAeMC2SosCMCnV6Encth4l1w7xxnjGtj\njAlivlzmOizb87PAJlXdGqA84Surqd/wa5yR6Pp/7txYZ4wxQcqXAYO+zjZ9a8WhECpeArf8G7al\nw9Q/eZ3GGGPylOsehIgcIedR3QRn0Lf4gKUKZ/W7wJb74bvXnKZ+qd29TmSMMTnKtUCoapmiDFKs\ndPib29TvQaepX6U6Xicyxphf8ekqJhFpJCIPuFPDQIcKe5HR0OM9iImDcX3g1FGvExljzK/4MuTo\nQ8D7QCV3el9EBub9KXNB8VWcpn771sLHD1pTP2NM0PFlD+JO4HJV/Yuq/gW4ArDuc/5Qs43T1O+H\nibDgba/TGGPML/hSIATIyDaf4S4z/nDVILi0I3z5OGxZ6HUaY4w5x5cC8R4wX0SeFJEnge+AdwOa\nqjiJiICub0J8MozvB8f2eZ3IGGMA3+6DeAG4A9jvTv1V9aVABytWSpaHnqPh2G748C7IzLjwZ4wx\nJsByLRAi8pmI/E5ESqvqIlV9xZ2WFGXAYqNKE7jhWVg3E75+1us0xhiT5x7EW8BNwAYRGSciXUXE\nekMEUrN+0Oh2+PoZWDvd6zTGmGIu1wKhqpNV9TbgYmAi0AfYLCLviUiHogpYrIjATc9DUn3nUNPB\nLV4nMsYUY76cgziuqmNVtStwHdAY+CLgyYqrmDjnfERmhjMS3dlTXicyxhRTvtwolyQiA0XkW+Aj\n4EugacCTFWcVasEtrzntOL583Os0xphiKq+T1HeLyExgMVAbGKyqNVV1iKouK7KExVW9znDlA7Dw\nHfh+vNdpjDHFUF7jQVwJ/AuYoao2mLIX2j/p7EV8/CBUbgCV6nqdyBhTjOR1kvoOVZ1mxcFDkdHQ\n/T2IKQ1je8OpI14nMsYUIwUZk9oUpfhk6D4C9q+DKQOtqZ8xpshYgQgFNVpD2z/Dikkw/y2v0xhj\niom8RpRLyOuDqrrf/3FMrlo9DFsWwNQnoGpTZzQ6Y4wJoLz2IBYB6e7jHmANsNZ9vqgwKxWRciIy\nQURWi8gqEblSRBJEZJqIrHUfyxdmHWEnIgK6vgHxVd2mfnu9TmSMCXN5naSuoao1genAzapaUVUr\nAJ2AqYVc78vAF6paB2gErAKG4FwxVRuY4c6b7EqWh15jnOIw8U5r6meMCShfzkFcoaqfZc2o6udA\ny4KuUETKAlfjtgxX1dOqehC4BRjlvm0U0KWg6whryY3gpmGw/iv4aqjXaYwxYcyXArFdRP4kIinu\n9ASwvRDrrIFzmOo9EVkiIu+ISCkgSVV3uO/ZCSTl9GERGSAi6SKSvmfPnkLECGFN+0Dj38HsZ2Ht\nNK/TGGPClC8F4jYgEZjkTpXcZQUVhdOq4w1VbQIc47zDSaqqQI7Xc6rqcFVNU9W0xMTEQsQIcTcN\ng6RU+PBuOLDJ6zTGmDDkS7O+/ar6kKo2caeHCnkF01Zgq6rOd+cn4BSMXSKSDOA+7i7EOsJfdEno\nOcqa+hljAsaXZn2XishwEZkqIjOzpoKuUFV3AltE5DJ3UTtgJTAF6Osu6wtMLug6io0KtaDLG7B9\nCXxh5/SNMf6VVy+mLOOBN4F3AH9dNjMQeN8dgGg90B+nWI0TkTuBTUBPP60rvNXtBC0fhLmvQPUr\noFEvrxMZY8KELwXirKq+4c+VqupSIC2Hl9r5cz3FRru/uk39HoLKqZBUz+tExpgw4MtJ6o9F5D4R\nSXZvZku40F3WpohFRjn9mkqUgXG94eRhrxMZY8KALwWiLzAYmItzB3XWHdYmmJSpDD3eg/0bYMoD\n1tTPGFNovlzFVCOHqWZRhDP5lHIVtPsLrJwM3/n1qKAxphjy5RwEItIAqAfEZi1T1dGBCmUKodVD\nTlO/aX92mvpddIXXiYwxIcqXy1z/CrzqTtcCzwKdA5zLFJQIdHkdylZ3mvodLaZ3mxtjCs2XcxDd\nca4u2qmq/XGa65UNaCpTOCXLOU39Thywpn7GmALzpUCccIcdPSsi8Th3OFcPbCxTaJVT4abnYcPX\nMOufXqcxxoQgX85BpItIOeBtnCuYjgLzAprK+EeT38Hm72DOMGeAoUuv9zqRMSaEiObjckgRSQHi\nVfX7QAXKj7S0NE1Ptytu83TmBLzbAQ5uhntmQ/kUrxMZYzwmIotUNaeblX8hX2NSq+rGYCkOxkfR\nJaHnGKc37rg+cOak14mMMSEiXwXChKiEGtD1TdixDL541Os0xpgQYQWiuKhzI7R6GBaNhKX/8zqN\nMSYEFKhAiMhmfwcxRaDtnyGlNXzyf7BrhddpjDFBrqB7EOLXFKZoREbBre9CbFkY2xtOHvI6kTEm\niBW0QARFJ7jVO49w8ozdBJYvZZKcpn4HNsLk+62pnzEmV7neByEig3J7CSgdmDj5cyYjk/SNB7iq\ndkWvo4SWi1tCh7/B1D/BvNeg5QNeJzLGBKG89iDK5DKVBl4OfLQLE2DOWus1VCBXPgB1b4Zpf4FN\ndt+jMebX8rqTei3wparuK6ow+RUXE8XXa/bw2I11vY4SekTgltdg1zVOU79750DpSl6nMsYEkbz2\nIKoD40Vkjog8KSKXi0hQnZwuExvF6p1H2H3Ybv4qkNiyzk10Jw/BhDsg46zXiYwxQSTXAqGqz6hq\nW+BGYBlwB7BYRP4rIn1EJKmoQuamdKyzAzRn7V6Pk4Swyg2g0wuwcQ7MetrrNMaYIOLLiHJHVHWS\nqt6jqk2Ap4BEwPMBg0pGR1KxdIydhyisxrdD077wzQuw+jOv0xhjgoQvAwY1zT7hjCr3EXBTwNP5\n4KpLKjJn7V4yM+1yzUK54VlIbgST7nXGtTbGFHu+3AfxOvAdMByn5fc8YDzwo4hcF8BsPrn60kT2\nHTvNyh2HvY4S2qJjoedo59Iwa+pnjMG3ArEdaKKqaaraDGgCrAc64Aw/6qmseyBm22GmwiufAl2H\nw87v4fPBXqcxxnjMlwJxqaqea9yjqiuBOqq6PnCxfFepTCx1k+OZs8ZOVPvFZR3hqkGweDQsed/r\nNMYYD/lSIFaIyBsi0sadXgdWikgJ4EyA8/nk6ksrkr5pP0dOBkWc0HftE05Tv08Hwc7lXqcxxnjE\nlwLRD/gJeNid1rvLzgDXBipYfrSrk8SZDGW27UX4R2QUdB8BJcs75yOsqZ8xxZIvl7meUNXnVbWr\nOw1T1eOqmqmqR4si5IU0vagc5eKimbFql9dRwkfpStBjpDNU6Uf3WVM/Y4qhsBgwKCoygraXVWLW\nj7s5m5HpdZzwcdEV0OHvsPoTmPuq12mMMUUsLAoEQLu6SRw4fobFmw96HSW8XHEf1LsFpj8JG7/1\nOo0xpgiFTYG4+tKKREeKHWbyNxHo/G9nXOsJ/eGIbV9jiot8FwgR+aeIPCoiFQIRqKDKxEZzeY0K\nTLcC4X+x8c5NdCcPW1M/Y4qRguxBLADOAi8WZsUiEikiS0TkE3c+QUSmicha97F8fr+zfd1KrNtz\njA17jxUmmslJUn24+SXY9A3M/IfXaYwxRSDfBUJVP3KvaupTyHU/BKzKNj8EmKGqtYEZ7ny+tKvr\nNJi1w0wB0ug30Kw/fPsSrP7U6zTGmADzpVlfoog8LiLDRWRE1lSYlYpINZxmf+9kW3wLMMp9Pgro\nkt/vrZ4Qx2VJZZi20gpEwHQcCsmNYdLvYX9Q3ExvjAkQX/YgJgNlgenAp9mmwngJeATIfk1qkqru\ncJ/vBAo03sT19ZNYuHE/e4+eKmREk6NzTf0ExvaBMye8TmSMCRBfCkScqj6qquNUdWLWVNAVikgn\nYLeqLsrtPaqqQI53ZonIABFJF5H0PXt+3aDvhtRkMhWmrrC9iIApfzF0Gw67lsNnf/Q6jTEmQHwp\nEJ+IyI1+XGcroLOIbAQ+ANqKyH+AXSKSDOA+7s7pw6o63O0sm5aYmPir1+tULkONiqX4bPmOHD5t\n/ObS66H1H2HJf2DxGK/TGGMCwJcC8RBOkTgpIkfcqcCDL6jqY6paTVVTgN8AM1X1d8AUoK/7tr44\nh7byTUS4MbUy89bvY/+x0wWNaXxx7eNQo42zF7Hje6/TGGP8zJdeTGVUNUJVY93nZVQ1PgBZhgId\nRGQt0N6dL5AbGiSTkalMXbHTb+FMDiIi3aZ+CTCuN5ywu9iNCSc+XeYqIp1FZJg7dfLXylX1K1Xt\n5D7fp6rtVLW2qrZX1f0F/d76VeK5KCGOz36wAhFwpSpCz1FwaKs19TMmzPhymetQnMNMK93pIRH5\nV6CDFYZzmCmZuT/t5eBxO8wUcNVbwHVPwY+fwrcve53GGOMnvuxB3Ah0UNURqjoC6IhzD0NQuzG1\nMmczlS/tMFPRuPxeqN8VZvwNNn7jdRpjjB/4eid1uWzPywYiiL+lVi1LSoU4Ji/d7nWU4kEEOr8K\nCbVgfH84YoXZmFDnS4H4F7BEREaKyChgEfB0YGMVnojQpUlV5q3fx45DdjNXkShRBnqNgdNHnSKR\nYUPAGhPKfLmK6X/AFcCHwETgSlUdG+hg/tClcVVUYYrtRRSdSnXh5pdh81zncJMxJmTlWiBEpI77\n2BRIBra6UxV3WdBLqViKJheVY9KSbV5HKV4a9oS0O51R6FZ97HUaY0wBReXx2iBgAPB8Dq8p0DYg\nifysS+Oq/HXKClbvPEydyoG4fcPkqOO/YPsS59LXSvWgQi2vExlj8inXPQhVHeA+vUFVr80+4VzZ\nFBI6NUwmMkL4aIkdZipSUSWc+yMiImFcHzh93OtExph88uUk9VwflwWlCqVL0ObSRCYv3UZGpt3E\nVaTKXQTd3oFdK5x2HHYTnTEhJa9zEJVFpBlQUkSaiEhTd7oGiCuyhH7QvVk1dhw6yey1v+7+agKs\ndnto8wgsfR8Wj/Y6jTEmH/I6B3E90A+oBryQbfkR4PEAZvK79nWTSCgVw9gFW7j2skpexyl+2jwK\nWxfCZ4MhuRFUaex1ImOMD/I6BzHKPd/Q77xzEJ1V9cMizFhoMVER3Nq0KtNX7WLPERtIqMhFRDqH\nmkpVdM5HnDjgdSJjjA98uQ9ioojcJCKPiMhfsqaiCOdPvZpX52ym8uHirV5HKZ5KVYAeo+Dwdme4\n0szMC3/GGOMpX5r1vQn0AgYCAvQALg5wLr+7pFIZ0i4uz9iFW1A7WeqN6s3h+qdhzefw7UtepzHG\nXIAvVzG1VNU+wAFV/RtwJXBpYGMFRq/m1Vm/9xgLNhS4k7gprBYDoMGtMPMfsGG212mMMXnwpUCc\ndB+Pi0gV4AzOndUh56aGyZSJjeL9+Zu9jlJ8icDNr0CFS2DCHXDYhoY1Jlj5UiA+FpFywHPAYmAj\n8N9AhgqUuJgoejSrzmfLd7Dr8MkLf8AERonS0HOMc/Pc+H7W1M+YIJVngRCRCGCGqh5U1Yk45x7q\nqGrInaTO0ufKi8lQtb0Ir1WqA51fgS3fwfQnvU5jjMlBngVCVTOB17LNn1LVQwFPFUApFUtxzaWJ\n/Hf+Zk6ftStpPJXaHZrfDfP+DSsne53GGHMeXw4xzRCRW0VEAp6miPRrVYO9R0/x2XI7/u2565+G\nqmnw0f2w9yev0xhjsvGlQNwDjAdOichhETkiIocDnCugWl9SkZoVSzFy7kavo5ioEtBjJERGW1M/\nY4KMLzfKlVHVCFWNUdV4dz6k+2ZHRAh9W6awdMtB0jfaJa+eK1cdbn0Hdq+ETwdZUz9jgoQvN8rN\n8GVZqOmZ5PzBAAAT6UlEQVSRVo3ycdG8+fU6r6MYgEvawTVDYNn/YNFIr9MYY8i7m2usiCQAFUWk\nvIgkuFMKULWoAgZKXEwU/VrWYPqq3fy484jXcQzA1Y9ArXbw+SPOYEPGGE/ltQdxD7AIqOM+Zk2T\ngX8HPlrg9bnyYkpGR/LWbNuLCAoREdDtbShVyTkfcdwO/xnjpby6ub6sqjWAP6pqTVWt4U6NVDUs\nCkT5UjH8pkV1pizdzraDJ7yOY8Bp6tdztHOH9aR7ramfMR7y5ST1qyLSUkRuF5E+WVNRhCsKd7Wu\nCcDbs9d7nMScU62ZM6b12i/hmxcu/H5jTED4cpJ6DDAMuApo7k5pAc5VZKqWK8mtTavx3/mb2XHI\n9iKCRvO7oEF3mPU0rP/K6zTGFEu+3AeRBrRS1ftUdaA7PRjoYEVpYLtLUJTXZtmNWkFDBG5+GSrU\nhgl3OuNIGGOKlC8F4gegcqCDeKla+Th6Na/O2IVb2HrAbtQKGiVKQ68xcOaENfUzxgO+FIiKwEoR\n+VJEpmRNgQ5W1O6/9hJEhFdn2F5EUEm8DG55FbbMh2kh2yPSmJAU5cN7ngx0iGCQXLYkt7e4iDHf\nbeKeNjWpmVja60gmS4NbYcsC+O51qN4C6nf1OpExxYIvVzF9jTMGRLT7fCHOuBBh5/5rL6FkdCRD\nP1/tdRRzvg7/gGotYPIDsHet12mMKRZ8uYrpbmAC8Ja7qCrwUUFXKCLVRWSWiKwUkRUi8pC7PEFE\nponIWvexfEHXUVCJZUrw+2tqMXXlLr5bv6+oV2/yEhXjNPWLKgFje8PpY14nMibs+XIO4n6gFXAY\nQFXXApUKsc6zwB9UtR5wBXC/iNQDhuAMTlQbmOHOF7k7r6pBlbKxPPXpSjIzrWlcUClb1Wnqt2c1\nfPJ/1tTPmADzpUCcUtXTWTMiEgUU+DdTVXeo6mL3+RFgFc5eyS3AKPdto4AuBV1HYcRGR/JIxzr8\nsO0wk5Zs8yKCyUuttnDt4/D9WEgf4XUaY8KaLwXiaxF5HCgpIh1wxob42B8rdxv/NQHmA0mqmjWC\nz04gKZfPDBCRdBFJ37Nnjz9i/ErnRlVoVL0cQ79YzeGTdmll0Gn9R7ikA3wxBLaF5ekwY4KCLwVi\nCLAHWI7TwO8z4E+FXbGIlAYmAg+r6i8GIFJVJZe9FFUdrqppqpqWmJhY2Bg5iogQ/nFLffYePcUL\nU9cEZB2mECIioNtwKJ0E4/paUz9jAsSXAlESGKGqPVS1OzDCXVZgIhKNUxzeV9UP3cW7RCTZfT0Z\n2F2YdRRWw2rl6HPFxYyet5Hvtx70MorJSVwC9BwFR3fChwOsqZ8xAeDTmNT8siCUBKYXdIXu2Nbv\nAqtUNXsntilAX/d5X5y24p76w/WXUaF0CZ6Y9AMZdsI6+FR1m/r9NA3mPO91GmPCji8FIlZVj2bN\nuM/jCrHOVkBvoK2ILHWnG4GhQAcRWQu0d+c9FR8bzZ871WP5tkOM+GaD13FMTtLuhNSeTlO/dTO9\nTmNMWPHlTupjItI068ojEWkGFLjtqap+A0guL7cr6PcGys0Nk5mydDvPTf2Ra+skckmlMl5HMtmJ\nwM0vwc7lMPEuuGc2lK3mdSpjwoIvexAPA+NFZI6IfAOMBR4IbKzgISL8s1sD4mIi+cO4ZZzNsGPd\nQSemlNPU7+wpp6nf2dMX/Igx5sJ8abWxEGfY0d8D9wJ1VXVRoIMFk0plYnmqSwOWbT3EG1/Z8KRB\nqWJtuOU12LoQpv3Z6zTGhAVf9iDAGSSoIdAUuC2cRpTzVaeGVejUMJmXZqxl0Sa7rDIo1e8CV9wH\n89+EHyZ6ncaYkFfsR5TLj392S6VKuVgG/ncJB47ZYYyg1OHvUP1ymPIg7LF7WIwpDBtRLh/iY6N5\n7fam7Dl6isETlqHWCyj4REa7Tf1iYVxvOHX0gh8xxuTMRpTLp4bVyvH4jXWZvmo3r9v5iOAUXwW6\nvwt718AnD1tTP2MKyJfLXLNGlFsAnMpaqKqdA5YqyPVrmcKSzQd57ssfqV2pNNfVt/oZdGpe4zT1\nm/mUc8ipxd1eJzIm5NiIcgUgIjzbvSGb9h3j4bFLmfj7ltRNjvc6ljnfVX+ALQvhi8egSlOo1szr\nRMaEFF9HlFsNlHGnVe6yYi02OpK3+6QRHxvNXaPS2X3kpNeRzPkiIqDrmxCfDOOtqZ8x+eXLVUw9\ngQVAD6AnMF9Eugc6WCioFB/L233SOHD8NH3eXcCh49YaPOjEJUCPUXB0l3OndWaG14mMCRm+nKR+\nAmiuqn1VtQ/QArA7kVyp1coyvHca6/cc445RCzl++qzXkcz5qjaFG56BdTNg9nNepzEmZPhSICJU\nNXvr7X0+fq7YuKp2RV65rTFLNh/g3v8s5uQZ+ys16DTrD41ug6+Gwk8FbkZsTLHiyz/0X4jIlyLS\nT0T6AZ8Cnwc2Vujp2CCZod0aMnvNHu4alW57EsFGBG56ASrVg4l3w8EtXicyJuj5cpJ6MPAWTquN\nhsBwVX0k0MFCUc/m1RnWoxFz1+2l74gFHLHhSoNLTBz0HA0ZZ6ypnzE+yLVAiMglItIKQFU/VNVB\nqjoI2CMitYosYYjp3qwar9zWhCWbD3L72/Pt6qZgU/ES6PI6bEuHqU94ncaYoJbXHsRLwOEclh9y\nXzO56NSwCsP7NOOn3Ufp8u9vWbUjp81oPFOvM1z5ACwYDssneJ3GmKCVV4FIUtXl5y90l6UELFGY\naFsnifH3XkmGKt3fmMuMVbu8jmSya/8kXHSl09Rv92qv0xgTlPIqEOXyeK1kHq8ZV4OqZZl8/1XU\nSCzFnaPSeeaL1TbgULCIjIbu7znnJcb1saZ+xuQgrwKRLiK/amAjIncBxWrAoMKoXDaWCfe25LYW\n1Xnjq3Xc9vZ37DhU4BFbjT/FJ0P3EbBvLXz8oDX1M+Y8klvLahFJAiYBp/m5IKQBMUBXVd1ZJAnz\nkJaWpunp6V7H8NlHS7bx+KTlREUIf+5Uj+7NqiGS2/DcpsjMeR5m/B1ueA4uH+B1GmMCTkQWqeoF\nx/XJtUBk+6JrgQbu7ApVnemHfH4RagUCYMPeYzwyYRkLNx7gmssSebprKlXL2RE7T2Vmwge3wU8z\noP/nUL2514mMCSi/FYhgFooFAiAzUxk9byPPfPEjinLP1bW4t00tSsZEeh2t+DpxAN5qA5ln4Z7Z\nUKqi14mMCRhfC4S1zPBARITQr1UNpg26mvZ1k3h5xlraPf8VHy7eSkZm6BbskFayvHMT3bG91tTP\nGJcVCA9VKx/Hv29vyrh7rqR8qRgGjVtGhxe+ZtKSrXa1kxeqNIYbn4P1s+DrZ7xOY4znrEAEgRY1\nEvj4gat447dNiYmK4P/GLqPdC1/z7jcbOHTC2nUUqaZ9oPFv4etnYa019TPFm52DCDKZmcrUlbt4\ne856Fm06QFxMJF2bVKVX8+qkVi1rVz0VhdPH4d0OcHibcz6i3EVeJzLGr+wkdRj4YdshRs3dyORl\n2zl9NpOaiaW4pVFVOjeuQo2KpbyOF972rYPh10CFS+COLyCqhNeJjPEbKxBh5NDxM3z2ww4+WrKN\n+RucYTNrJpaiXZ1KXFunEs1TEoiOtKOFfrfqYxj7O2h+F9z0vNdpjPEbKxBhavvBE3y5YiczV+9m\n/vr9nM7IJC4mkiYXlaN5SgItUhJoclF5u2TWX6b+Cea+Ct3egYY9vE5jjF9YgSgGjp06yzc/7WXu\nT3tZsPEAq3ceRhUiBGomlqZecjz1qsRTLzmeWpVKkxwfS0SEncPIl4yzMOpm2LEU7p4Jlep6nciY\nQrMCUQwdOnGGxZsOsGTzAVbuOMKqHYfZdvDnvk8xURFcnBDHxRVKkVIhjirlSpIUH0ul+BIklXEe\nY6Ntz+NXjuyEN1tDbFkYMAtKlPE6kTGFYgXCAHDw+GlW7jjMxr3H2bTvGBv2HmPTvuNs3HeMU2d/\nfa9FfGwU5eJiKFsy+ucpznksExtFyehIZ4qJJDbambKWxUZHEBUZQVSEEBUpREYI0RERREYKURE/\nz4fkXsyGOTC6M9S7xekCa1eTmRDma4GIKoow+SEiHYGXgUjgHVUd6nGkkFYuLoaWtSrS8rwxAFWV\ng8fPsOvISXYdPsXuwyfZfcR5PHjiDIfcafuhExx2n5/J8M8fEyIQFSGICOLOC+I+cm452efdf4+z\nvy7um37+nF/i5ap3xO08sOI/PL+6POMibwrsyozJQ5nYaKYPahPw9QRVgRCRSOA1oAOwFVgoIlNU\ndaW3ycKPiFC+VAzlS8VQp/KF36+qnDqbyckzGZw4k8GJ087jyTMZnDj98/KzmZmczVAyMpUzmUpG\nRiZnM5Wzmc4y57VMzmQqmarg/Ieqoueeg6Lnum9n7eWe/1rWPNneG0hbdQArN23k4SOjKHlxczaX\nanDhDxkTAEV1KDioCgTQAvhJVdcDiMgHwC2AFQiPici5Q0p5jSQV9k68D8PbcN+ef0C32VA60etE\nxgRMsF08XxXYkm1+q7vsHBEZICLpIpK+Z8+eIg1nDCXLQc8xcGI/TLzTmvqZsBZsBeKCVHW4qqap\nalpiov31ZjyQ3BBuHAYbvoav/uV1GmMCJtgKxDagerb5au4yY4JL097Q5Hcw+zlYM9XrNMYERLAV\niIVAbRGpISIxwG+AKR5nMiZnNw6Dyqnw4d1wYJPXaYzxu6AqEKp6FngA+BJYBYxT1RXepjImF9El\nnUGGVGF8Xzh7yutExvhVUBUIAFX9TFUvVdVaqvq013mMyVNCTej6BmxfAl8M8TqNMX4VdAXCmJBT\n5yZo9RCkj4BlY71OY4zfWIEwxh/a/gUuvgo+fgh22W07JjxYgTDGHyKjoPsIiI2Hcb3h5GGvExlT\naFYgjPGXMklOI7/9G2DKAxRJ/w9jAsgKhDH+lNIK2j8JKyfDd294ncaYQgnpdt8icgT40esceagI\n7PU6RB4sX+FYvoIL5mwQ/vkuVtULtqIItmZ9+fWjLz3NvSIi6Zav4Cxf4QRzvmDOBpYvix1iMsYY\nkyMrEMYYY3IU6gViuNcBLsDyFY7lK5xgzhfM2cDyASF+ktoYY0zghPoehDHGmACxAmGMMSZHIVsg\nRKSjiPwoIj+JSNC10RSRjSKyXESWikh6EOQZISK7ReSHbMsSRGSaiKx1H8sHWb4nRWSbuw2XisiN\nHmWrLiKzRGSliKwQkYfc5UGx/fLIFyzbL1ZEFojIMjff39zlwbL9cssXFNvPzRIpIktE5BN3vki2\nXUiegxCRSGAN0AFn3OqFwG2qGjRd0kRkI5CmqkFxs42IXA0cBUaragN32bPAflUd6hbZ8qr6aBDl\nexI4qqrDvMiULVsykKyqi0WkDLAI6AL0Iwi2Xx75ehIc20+AUqp6VESigW+Ah4BuBMf2yy1fR4Jg\n+wGIyCAgDYhX1U5F9bsbqnsQLYCfVHW9qp4GPgBu8ThTUFPV2cD+8xbfAoxyn4/C+UfFE7nkCwqq\nukNVF7vPj+AMZlWVINl+eeQLCuo46s5Gu5MSPNsvt3xBQUSqATcB72RbXCTbLlQLRFVgS7b5rQTR\nL4RLgekiskhEBngdJhdJqrrDfb4TSPIyTC4Gisj37iEozw6BZRGRFKAJMJ8g3H7n5YMg2X7uIZKl\nwG5gmqoG1fbLJR8Ex/Z7CXgEyMy2rEi2XagWiFBwlao2Bm4A7ncPoQQtdY41Bs1fTa43gJpAY2AH\n8LyXYUSkNDAReFhVf9HPOxi2Xw75gmb7qWqG+/tQDWghIg3Oe93T7ZdLPs+3n4h0Anar6qLc3hPI\nbReqBWIbUD3bfDV3WdBQ1W3u425gEs5hsWCzyz1+nXUce7fHeX5BVXe5v7iZwNt4uA3dY9MTgfdV\n9UN3cdBsv5zyBdP2y6KqB4FZOMf3g2b7ZcmeL0i2Xyugs3tO8wOgrYj8hyLadqFaIBYCtUWkhojE\nAL8Bpnic6RwRKeWeLERESgHXAT/k/SlPTAH6us/7ApM9zPIrWb8Arq54tA3dk5jvAqtU9YVsLwXF\n9sstXxBtv0QRKec+L4lzcclqgmf75ZgvGLafqj6mqtVUNQXn37mZqvo7imrbqWpITsCNOFcyrQOe\n8DrPedlqAsvcaUUw5AP+h7ObfAbnnM2dQAVgBrAWmA4kBFm+McBy4Hv3FyLZo2xX4ezCfw8sdacb\ng2X75ZEvWLZfQ2CJm+MH4C/u8mDZfrnlC4rtly3nNcAnRbntQvIyV2OMMYEXqoeYjDHGBJgVCGOM\nMTmyAmGMMSZHViCMMcbkyAqEMcaYHEV5HcAYr4hI1qWCAJWBDGCPO39cVVt6lOvvwGxVne7F+o3J\nYpe5GkPwdI41JpjYISZjciAiR93Ha0TkaxGZLCLrRWSoiPzWHT9guYjUct+XKCITRWShO7XK4Tv7\nichHbv/+jSLygIgMcvv8fyciCe77RopId/f5RhH5m4gsdtdXpyi3gynerEAYc2GNgHuBukBv4FJV\nbYHTfnmg+56XgRdVtTlwK79szZxdA5xxEJoDT+McymoCzAP65PKZvaraFKd53B8L/+MY4xs7B2HM\nhS1Ut7WyiKwDprrLlwPXus/bA/WctkgAxItIaf15nIEss9QZs+GIiBwCPs72XQ1zWX9Wc8BFOMXF\nmCJhBcKYCzuV7XlmtvlMfv4digCuUNWTfviu3D6Tkcd7jPE7O8RkjH9M5efDTYhIY/exhYiM9iyV\nMYVgBcIY/3gQSHNHH1uJc84C4CLghHexjCk4u8zVmAASkeeAMar6vddZjMkvKxDGGGNyZIeYjDHG\n5MgKhDHGmBxZgTDGGJMjKxDGGGNyZAXCGGNMjqxAGGOMydH/A5/yGyzMjb6IAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "y0 = [50,100] #initial conditions array\n", "time_scale =40 \n", "t = np.linspace(0,time_scale,1000) # minutes \n", "sol = integrate.odeint(tank_Concentration,y0,t)\n", "lineObjects = plt.plot(t,sol)\n", "\n", "\n", "plt.axis([0 ,42,-5,100])\n", "plt.legend(iter(lineObjects), ('concentration', 'volume'))\n", "plt.xlabel('Time,min')\n", "plt.ylabel('Concentration, g/L and Volume,L')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see that the concentration of the tank still follows a similar curve and that the volume of the tank steadily decreases. Even with the axis in both axis the same, we can check to see if there is a difference in thechange in concentration. " ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(array([ 19, 20, 300, 301, 302, 303, 304]), array([0, 0, 1, 1, 1, 1, 1]))\n", "[ 39.87975966 39.39989287 39.93993994 39.73973974 39.53953954\n", " 39.33933934 39.13913914]\n", "It will take approximately 0.76 minutes to reach 40 g/L\n" ] } ], "source": [ "\n", "#rudimentarily, the time that it will take to reach 40 g/L is presented here.\n", "h = np.where((sol>=39) & (sol<=40))\n", "print (h)\n", "print(sol[h])\n", "print(\"It will take approximately {} minutes to reach 40 g/L\".format(h[0][0]*time_scale/len(t)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is sensible as there is more water diluting the feed stream." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Going Further :\n", "\n", "- change both the C_in concentrations to 60, what happens to the concentration?\n", "- change both the C_in concentrations to 60, with the net flow rate equal to 0; what happens to the concentration?\n", "- change both the C_in concentrations to 60, with the net flow rate equal to -10; what happens to the concentration?\n", "- replicate the results in example 1 by varying the parameters of the flow conditions." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "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.6.1" } }, "nbformat": 4, "nbformat_minor": 2 }