{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2f0879a1-70d3-4f20-8e6c-9909d6dcee58",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.integrate import odeint\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pylab as pl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "afc846e9-2846-498f-9c8b-24d1b3923f1e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Figure size 1200x1200 with 0 Axes>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 1200x1200 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize = (12,12))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "77e10b6e-7708-4e3a-9e4e-0fceceed0ca1",
   "metadata": {},
   "outputs": [],
   "source": [
    "omega = 2\n",
    "omegasq = omega**2\n",
    "a, b, c, d = 0, 1, omegasq, -0.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "149f2096-f875-4d54-97ed-2f248de87585",
   "metadata": {},
   "outputs": [],
   "source": [
    "def dx_dt(xx,t):\n",
    "    x = xx[0]\n",
    "    y = xx[1]\n",
    "    term1 = a*x + b*y\n",
    "    term2 = -omegasq*np.sin(x) + d*y\n",
    "    return [term1, term2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "930bbc07-44eb-439c-8ec6-4ab4a5e96458",
   "metadata": {},
   "outputs": [],
   "source": [
    "maxtime = 20\n",
    "xmin = -8\n",
    "xmax = 8\n",
    "ts = np.linspace(0, maxtime, 1000)\n",
    "ic = np.linspace(xmin, xmax, 8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "621cb3f0-f1c7-4a78-8cc3-2e10da1fe553",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-8.0, 8.0)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 132340 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for r in ic:\n",
    "    for s in ic:\n",
    "        x0 = [r,s]\n",
    "        xs = odeint(dx_dt,x0,ts)\n",
    "        plt.plot(xs[:,0],xs[:,1],\"b\")\n",
    "\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.xlim(xmin,xmax)\n",
    "plt.ylim(xmin,xmax)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb6d7420-6149-4d90-8fe6-83b354ec8519",
   "metadata": {},
   "source": [
    "---\n",
    "# Competing Species Example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "12418800-37ab-422b-b2c1-f76a67e4b700",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.style.use('ggplot')\n",
    "\n",
    "import sympy as sm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "4854aa7b-832f-4866-b849-f9b869a3790e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Paramters for \"mutual exclusion\n",
    "beta = 1\n",
    "delta = 2\n",
    "gamma = 3\n",
    "b = 1\n",
    "d = 2\n",
    "c = 3\n",
    "\n",
    "def dz_dt(z,t):\n",
    "   return [z[0]*(beta - delta*z[0] - gamma*z[1]), z[1]*(b - d*z[1] - c*z[0])]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "75328322-db3f-4573-a9b1-dc30d4975c42",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find critical points\n",
    "\n",
    "#define symbolic exprsssion\n",
    "x, y = sm.symbols('x,y')\n",
    "P = x*(beta - delta*x - gamma*y)\n",
    "Q = y*(b - d*y - c*x)\n",
    "\n",
    "#Set P(x,y) = 0 and Q(x,y) = 0\n",
    "Peqn = sm.Eq(P,0)\n",
    "Qeqn = sm.Eq(Q,0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "ed7258d8-1908-43bd-89e4-0048e9ff5a1f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(0, 0), (0, 1/2), (1/5, 1/5), (1/2, 0)]\n"
     ]
    }
   ],
   "source": [
    "#solve our equations for populations of species of x and y\n",
    "criticalpoints = sm.solve((Peqn,Qeqn),x,y)\n",
    "print(criticalpoints)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "9ec4806e-3199-49f0-9217-e7f20f903266",
   "metadata": {},
   "outputs": [],
   "source": [
    "#calculate a vectorfield\n",
    "X,Y = np.mgrid[0:0.5:20j,0:0.5:20j]\n",
    "u = X*(beta - delta*X - gamma*Y)\n",
    "v = Y*(b-d*Y-c*X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "092f82d4-5aae-4056-b9ee-2453b651c988",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 200468 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for i in range(4):\n",
    "    plt.plot(criticalpoints[i][0],criticalpoints[i][1],\"yo\")\n",
    "    \n",
    "#we can color our critical points algorithmically\n",
    "plt.quiver(X,Y,u,v,color = \"red\")\n",
    "\n",
    "#plot trajectories\n",
    "ts = np.linspace(0,10.500)\n",
    "ic = np.linspace(0,0.5,10)\n",
    "for r in ic:\n",
    "    for s in ic:\n",
    "        z0 = [r,s]\n",
    "        zs = odeint(dz_dt,z0,ts)\n",
    "        plt.plot(zs[:,0],zs[:,1], \"royalblue\",linewidth = 0.75)\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "024812c6-e894-4ede-8594-99da95b5b92c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1. 1.]\n",
      "[-1.  -0.5]\n"
     ]
    },
    {
     "ename": "AttributeError",
     "evalue": "'EigResult' object has no attribute 'eignvalues'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[41], line 9\u001b[0m\n\u001b[1;32m      6\u001b[0m \u001b[38;5;28mprint\u001b[39m(np\u001b[38;5;241m.\u001b[39mlinalg\u001b[38;5;241m.\u001b[39meig(J)\u001b[38;5;241m.\u001b[39meigenvalues)\n\u001b[1;32m      8\u001b[0m pointcolor \u001b[38;5;241m=\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mro\u001b[39m\u001b[38;5;124m'\u001b[39m\n\u001b[0;32m----> 9\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m(\u001b[38;5;28meval\u001b[39m\u001b[38;5;241m.\u001b[39meigenvalues[\u001b[38;5;241m0\u001b[39m]\u001b[38;5;241m<\u001b[39m\u001b[38;5;241m0\u001b[39m \u001b[38;5;129;01mand\u001b[39;00m \u001b[38;5;28;43meval\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43meignvalues\u001b[49m[\u001b[38;5;241m1\u001b[39m]\u001b[38;5;241m<\u001b[39m\u001b[38;5;241m0\u001b[39m):\n\u001b[1;32m     10\u001b[0m     pointcolor \u001b[38;5;241m=\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124myo\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m     12\u001b[0m plt\u001b[38;5;241m.\u001b[39mplot(criticalpoints[k][\u001b[38;5;241m0\u001b[39m],criticalpoints[k][\u001b[38;5;241m1\u001b[39m],\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124myo\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n",
      "\u001b[0;31mAttributeError\u001b[0m: 'EigResult' object has no attribute 'eignvalues'"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 12372 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for k in range(4):\n",
    "    xx = float(criticalpoints[k][0])\n",
    "    yy = float(criticalpoints[k][1])\n",
    "    J = np.array([[beta-2*delta*xx - gamma*yy, -gamma*xx], [-c*yy, b -2*d*yy-c*xx]])\n",
    "    eval = np.linalg.eig(J)\n",
    "    print(np.linalg.eig(J).eigenvalues)\n",
    "\n",
    "    pointcolor = 'ro'\n",
    "    if(eval.eigenvalues[0]<0 and eval.eignvalues[1]<0):\n",
    "        pointcolor = \"yo\"\n",
    "    \n",
    "    plt.plot(criticalpoints[k][0],criticalpoints[k][1],\"yo\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "841551b2-00d3-41ed-a2d7-3a93424f44bc",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95808e94-ad18-48c4-a566-d4dd6e89336a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "221bd086-1ca7-40ae-b656-ce645c86f409",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "b9baea92-ed70-4183-9b65-28c3166449f4",
   "metadata": {},
   "source": [
    "---\n",
    "# Animation fun time!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "f2daadfc-ce2b-4a33-b882-22c7a14c51c5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Exothermic Reactor with Tc = 305.0 K')"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 35264 characters]",
      "text/plain": [
       "<Figure size 800x800 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "#import seaborn as sns\n",
    "from scipy.integrate import odeint\n",
    "\n",
    "from matplotlib import animation\n",
    "from IPython.display import HTML\n",
    "\n",
    "# https://github.com/jckantor/CBE30338\n",
    "# THIS FILE IS ORIGINALLY BY Jeff Kantor at University of Notre Dame\n",
    "\n",
    "# mode parameteres\n",
    "Ea = 72750 # activation energy J/gmol\n",
    "R = 8.314 # gas constant J/gmol/K\n",
    "k0 = 7.2e10 # Arrhenius rate constant 1/min\n",
    "V = 100.0 # Volume [L]\n",
    "rho = 1000.0 # Density [g/L]\n",
    "Cp = 0.239 # Heat capacity [J/g/K]\n",
    "dHr = -5.0e4 # Enthalpy of reaction [J/mol]\n",
    "UA = 5.0e4 # Heat transfer [J/min/K]\n",
    "q = 100.0 # Flowrate [L/min]\n",
    "cAi = 1.0 # Inlet feed concentration [mol/L]\n",
    "Ti = 350.0 # Inlet feed temperature [K]\n",
    "cA0 = 0.5; # Initial concentration [mol/L]\n",
    "T0 = 350.0; # Initial temperature [K]\n",
    "Tc = 305.0 # Coolant temperature [K]\n",
    "\n",
    "# Arrhenius rate expression\n",
    "def k(T):\n",
    "    return k0*np.exp(-Ea/R/T)\n",
    "\n",
    "def deriv(y,t):\n",
    "    cA,T = y\n",
    "    dcA = (q/V)*(cAi - cA) - k(T)*cA\n",
    "    dT = (q/V)*(Ti - T) + (-dHr/rho/Cp)*k(T)*cA + (UA/V/rho/Cp)*(Tc-T)\n",
    "    return [dcA,dT]\n",
    "\n",
    "# create a set of initial conditions\n",
    "ICs = [[cA0,T0] for cA0 in [0] for T0 in np.linspace(295,480,19)]\n",
    "ICs += [[cA0,T0] for cA0 in np.linspace(0,1,21) for T0 in [290]]\n",
    "ICs += [[cA0,T0] for cA0 in [1] for T0 in np.linspace(295,475,18)]\n",
    "\n",
    "# create background figure and axes\n",
    "#sns.set(font_scale=1.5)\n",
    "fig, ax = plt.subplots(figsize=(8,8))\n",
    "ax.set_xlim((0,1))\n",
    "ax.set_ylim((290,450))\n",
    "ax.set_xlabel('Concentration [gmol/liter]')\n",
    "ax.set_ylabel('Temperature [K]')\n",
    "ax.set_title('Exothermic Reactor with Tc = {0:.1f} K'.format(Tc))\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
