exercises.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introductory Numerical Methods for Simulating Batch Reactors\n",
    "This notebook comprises the tutorials and exercises we will use for this lecture."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preliminary Code\n",
    "The cell below defines helper functions to assist with the exercises. Please run this cell be before proceeding below to the exercises."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "def plot_solution(ts, *ys, labels = (), marker_size = 3):\n",
    "    if len(labels) == 0:\n",
    "        labels = (\"y\" + str(i) + \"(t)\" for i in range(1, len(ys)+1))\n",
    "    f = plt.figure()\n",
    "    plt.plot(ts, ys[0], marker = 'o', markersize = marker_size, label = labels[0])\n",
    "    plt.xlabel(\"time (t)\")\n",
    "    for i in range(1, len(ys)):\n",
    "        plt.plot(ts, ys[i], marker = 'o', markersize = marker_size, label = labels[i])\n",
    "    plt.legend(loc='best')\n",
    "    # return f"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tutorial 1: Simple ODE Simulation via Explicit Euler\n",
    "Solve the following ODE:\n",
    "$$\\begin{gathered}\n",
    "\\frac{d y(t)}{dt} = e^{-t} \\\\\n",
    "y(0) = -1\n",
    "\\end{gathered}$$\n",
    "for $t \\in [0, 1]$ using $\\Delta t=0.1$. Plot your result against the analytical solution (you also need to determine this).\n",
    "\n",
    "Recall that the update rule for explicit Euler is:\n",
    "$$\\tilde{y}(t_{k+1}) = \\tilde{y}(t_k) + f(y(t_k), t_k) \\Delta t$$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that let's compute the approximate solution using Euler's method. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define parameters\n",
    "y0 = -1\n",
    "dt = 0.1 # TODO try changing this\n",
    "ts = np.arange(0, 1.0001, dt)\n",
    "num_steps = len(ts) - 1\n",
    "\n",
    "# Evaluate the analytical function for each t in ts\n",
    "y = -np.exp(-ts)\n",
    "\n",
    "# Define storage vector for y(t)\n",
    "y_euler = np.zeros(len(ts))\n",
    "\n",
    "# Fill in y_euler using explicit Euler\n",
    "y_euler[0] = y0\n",
    "for k in range(num_steps):\n",
    "    y_euler[k + 1] = y_euler[k] + np.exp(-ts[k]) * dt\n",
    "\n",
    "# Plot the results\n",
    "plot_solution(ts, y, y_euler, labels = (\"Exact\", \"Explicit Euler\"))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tutorial 2: Simulate Coupled ODEs\n",
    "Simulate the following system of ODEs:\n",
    "$$\\begin{gathered}\n",
    "\\frac{d y_1(t)}{dt} = -5y_1(t) + 5y_2(t) \\\\\n",
    "\\frac{d y_2(t)}{dt} = 14y_1(t) - 2y_2(t) \\\\\n",
    "y_1(0)=y_2(0)=1 \n",
    "\\end{gathered}$$\n",
    "for $t \\in [0, 1]$ with $\\Delta t = 0.1$. Recall that the vectorized Euler update rule is as follows:\n",
    "$$\\begin{gathered}\n",
    "\\frac{d \\textbf{y}(t)}{dt} = \\textbf{f}(\\textbf{y}(t), t)\\\\\n",
    "\\textbf{y}(0) = \\boldsymbol{y_0}\n",
    "\\end{gathered}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define parameters\n",
    "y0 = np.ones(2)\n",
    "dt = 0.1\n",
    "ts = np.arange(0, 1.0001, dt)\n",
    "num_steps = len(ts) - 1\n",
    "\n",
    "# Define storage vector for y(t) --> index as y[k, i]\n",
    "y_euler = np.zeros((len(ts), 2))\n",
    "\n",
    "# Approximate y(t) using explicit euler and fill in the storage vector\n",
    "y_euler[0, :] = y0\n",
    "for k in range(num_steps):\n",
    "    fk = np.array([-5 * y_euler[k, 0] + 5 * y_euler[k, 1], 14 * y_euler[k, 0] - 2 * y_euler[k, 1]])\n",
    "    y_euler[k+1, :] = y_euler[k, :] + fk * dt\n",
    "\n",
    "# Plot the results\n",
    "plot_solution(ts, y_euler[:, 0],  y_euler[:, 1], labels = (\"y1(t)\", \"y2(t)\"))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tutorial 3: Stability of Explicit Euler\n",
    "Simulate the following:\n",
    "$$\\begin{gathered}\n",
    "\\frac{d y(t)}{dt} = -20y(t) \\\\\n",
    "y(0) = 1\n",
    "\\end{gathered}$$\n",
    "for $t \\in [0, 1]$. Let's compare with the analytical solution and see how the stability is affected by our choice of $\\Delta t$ (it should be stable only for $\\Delta t < 0.1$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define parameters\n",
    "y0 = 1\n",
    "dt = 0.11 # TODO adjust the stepsize to see what happens\n",
    "ts = np.arange(0, 1.0001, dt)\n",
    "num_steps = len(ts) - 1\n",
    "\n",
    "# Define storage vectors for y(t)\n",
    "y_euler = np.zeros(len(ts))\n",
    "\n",
    "# Define the analytical solution of y(t) for each t in ts\n",
    "y = np.exp(-20 * ts)\n",
    "\n",
    "# Do explicit euler and fill y_euler\n",
    "y_euler[0] = y0\n",
    "for k in range(num_steps):\n",
    "    y_euler[k + 1] = y_euler[k] + -20 * y_euler[k] * dt\n",
    "\n",
    "# Plot the solution \n",
    "plot_solution(ts, y, y_euler, labels = (\"Exact\", \"Explicit Euler\"))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 1: Simulate Reaction Network in Batch Reactor\n",
    "Simulate the response of the reaction system:\n",
    "$$\\begin{gathered}\n",
    "\\text{A} \\rightleftharpoons \\text{B} \\\\\n",
    "\\text{A} \\rightleftharpoons \\text{C}\n",
    "\\end{gathered}$$\n",
    "for $t \\in [0, 1]$ using implicit Euler. The ODEs are given by:\n",
    "$$\\begin{aligned}\n",
    "&\\frac{d \\textbf{c}(t)}{dt} = \\gamma \\textbf{r}(\\textbf{c}, t) \\\\\n",
    "&\\textbf{c}(0) = c_{i, 0}\n",
    "\\end{aligned}$$\n",
    "where we have:\n",
    "$$\\textbf{r}(\\textbf{c}, t) = \\textbf{k}^T(t)\\textbf{c}(t)$$\n",
    "omitting $\\beta$ since $\\beta=1$ and we have:\n",
    "$$k_{ij}(t) = A_{ij} \\exp\\Bigg(\\frac{-E_{a,ij}}{RT(t)}\\Bigg)$$\n",
    "The needed parameters are given below.\n",
    "\n",
    "**Hint:** Use `np.dot` for dot products and matrix multiplication. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the parameters\n",
    "R = 1.987\n",
    "A = np.array([[3.6362e6, 190.6879], [-2.5212e16, 0], [0, -8.7409e24]])\n",
    "Ea = np.array([[10000., 5000.], [25000., 0], [0., 40000.]])\n",
    "gamma = np.array([[-1, -1], [1, 0], [0, 1]])\n",
    "c0 = np.array([1., 0., 0.])\n",
    "dt = 0.1\n",
    "ts = np.arange(0, 1.0001, dt)\n",
    "num_steps = len(ts) - 1\n",
    "Ts = np.ones(len(ts)) * 325.\n",
    "Ts[Ts < 0.5] = 333\n",
    "\n",
    "# Define collection array c_euler\n",
    "c_euler = np.zeros((len(ts), 3)) # time by cs\n",
    "\n",
    "# Implement explicit Euler and fill c\n",
    "# TODO implement the explicit Euler simulation here --> filling in c_euler\n",
    "\n",
    "#  Plot the solution\n",
    "plot_solution(ts, c_euler[:, 0], c_euler[:, 1], c_euler[:, 2], labels = (\"cA(t)\", \"cB(t)\", \"cC(t)\"))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use a shared notebook [here](https://colab.research.google.com/drive/10qk9W7t1_t-3VG_Az_f553GrvBw1X3xF?usp=sharing) to share code with the class. Please DO NOT open until instructed to do so. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "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.9.7 (default, Sep 16 2021, 16:59:28) [MSC v.1916 64 bit (AMD64)]"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "8d689450f2d38e0e8f630e8a92d7c9c51deb0f64dc8e967a759d71fbd476c3a7"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}