Grad worksheet for MP1

worksheets/Worksheets_Lab3/Worksheet_G_Lab3_2_2024.ipynb

@@ -0,0 +1,208 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "1533bd8e",
+   "metadata": {},
+   "source": [
+    "Names: "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8be05185",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# your import statements\n",
+    "import numpy as np"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "c17b4ce4",
+   "metadata": {},
+   "source": [
+    "Let's go back to the long hallway with the smoker.   \n",
+    "\n",
+    "Consider a long hallway (length $L = 20$m) in an office building. If we assume that any\n",
+    "cigarette smoke mixes across the width of the hallway and vertically\n",
+    "through the depth of the hallway much faster than it mixes along the\n",
+    "hallway, we can write the diffusion of cigarette smoke as an equation\n",
+    "$$\\frac {\\partial S} {\\partial t}\n",
+    "= \\frac {\\partial \\kappa \\partial S}{\\partial x^2} - \\gamma S + \\alpha(x)$$\n",
+    "where $S$ is the concentration of smoke, $\\kappa = 0.05 $m$^2 $s$^{-1}$ is the rate of diffusion of\n",
+    "smoke, $\\gamma = 1/(3600 {\\rm s})^{-1}$ is the rate that the smoke sticks to the ceiling, $\\alpha(x)$ is the sources of smoke, t is the time and x is\n",
+    "distance along the hallway.  No smoke leaves through either end of the hallway.  We will assume the smoker is at $x=L/3$, and puts out 0.005 su s$^{-1}$ (su = smoke units)."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "9eb3d006",
+   "metadata": {},
+   "source": [
+    "**Question 1**\n",
+    "\n",
+    "Under what conditions can you move $\\kappa$ outside the derivative?  Make that assumption here.  What is your new differential equation?"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5f110569",
+   "metadata": {},
+   "source": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "ace733e0",
+   "metadata": {},
+   "source": [
+    "**Question 2**\n",
+    "\n",
+    "Using a centre-difference scheme, separating your hallway into $N=15$ divisions (so $N+1=16$ grid points) and noting that the 0th and $N=16$th grid points are boundary points we can use the following function to create a matrix that solves your equations in Question 1, assuming steady state.  That is writing your systems as:\n",
+    "$$ {\\rm matrix} \\cdot \\vec S = \\vec \\alpha $$\n",
+    "where $\\vec S$ is your vector of $S$ values at the grid points and $\\vec \\alpha$ is your vector of sources, all zero except at index 5."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2a4c62bd-7372-4dea-b8f3-08a32073b465",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def create_matrix_alpha(L, kappa, gamma, alpha_o):\n",
+    "    N = 15\n",
+    "    dx = L/N\n",
+    "    matrix = np.zeros((N+1, N+1))\n",
+    "    alpha = np.zeros(N+1)\n",
+    "    matrix[0, 0] = 1\n",
+    "    matrix[0, 1] = -1\n",
+    "    for i in range(1, N):\n",
+    "        matrix[i, i-1] = kappa/dx**2\n",
+    "        matrix[i, i] = -(2*kappa/dx**2 + gamma)\n",
+    "        matrix[i, i+1] = kappa/dx**2\n",
+    "    matrix[N, N] = 1\n",
+    "    matrix[N, N-1] = -1\n",
+    "    alpha[int(N/3)] = -alpha_o\n",
+    "\n",
+    "    return matrix, alpha"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6bfe0890-3961-4989-9cf9-d799739e9a2f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "matrix, alpha = create_matrix_alpha(L=20, kappa=0.05, gamma=1/3600, alpha_o=0.005)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "79181ef3-4a40-4a24-a90c-10a8d5151b1e",
+   "metadata": {},
+   "source": [
+    "Find the steady state solution."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c6c0e067-c35e-4df2-935a-7c40eebe2657",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "6bcb3435",
+   "metadata": {},
+   "source": [
+    "**Question 3**\n",
+    "\n",
+    "Now consider the case where the smoker stops.  Starting from the above steady state, we want to consider how the system changes in time. So we have:\n",
+    "$$ \\frac {\\partial \\vec S}{\\partial t} = matrix \\cdot \\vec S$$\n",
+    "\n",
+    "Discretize this equation using a forward Euler method for the time stepping.  (I suggest you continue to use the vector notation but add subscripts to mark the time at which variables are evaluated, e.g. $\\vec S_i$ would be the smoke at time step $i$.)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "1abcc6cf",
+   "metadata": {},
+   "source": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "305884ee",
+   "metadata": {},
+   "source": [
+    "**Question 4**\n",
+    "\n",
+    "Solve the system and plot the solution at various times.  \n",
+    "Hints:\n",
+    "   * Enforce the boundary conditions ($S[0] = S[1]$ and $S[N] = S[N-1]$ at each time step.\n",
+    "   * Try a time step of 5 s and run for, say 15 time steps\n",
+    "   * To multiple a matrix by a vector, use np.matmul(matrix, vector)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "263f8bac-3041-4f39-8361-b2d9411bf57e",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3ffd3078",
+   "metadata": {},
+   "source": [
+    "**Question 5 (optional)**\n",
+    "\n",
+    "If you have time, try solving the system using the midpoint method."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "039dbe49-25a3-4fb0-bbad-16dfbf84a9f3",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "577152a7",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "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.10.13"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}