Lecture-17-Explicit-Finite-Difference.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Lecture 17:  Numerical Solutions to the Diffusion Equation (Explicit Methods)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "## Reading and Reference\n",
    "\n",
    "* Numerical Recipes, W. Press, Cambridge University Press, 1986\n",
    "* Numerical Methods, R. Hornbeck, Quantum Publishers, 1975\n",
    "* B. Gustafsson, Fundamentals of Scientific Computing, Springer, 2011\n",
    "* S. Farlow, Partial Differential Equations for Scientists and Engineers, Dover, 1993 \n",
    "\n",
    "## What to Learn?\n",
    "\n",
    "* How to set up the finite difference equations\n",
    "* How to set up a data structure to organize the solution\n",
    "* How to apply the initial and boundary conditions to permit computation of the solution\n",
    "\n",
    "## What to do?\n",
    "\n",
    "*  Write the solver code for three example problems\n",
    "*  Visualize the results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "### Introduction\n",
    "----\n",
    "\n",
    "In each of the next three lectures a different numerical technique will be studied and implemented to solve the diffusion equation.  Each technique is built from mathematics that we've studied in previous lectures.  These or similar methods are used in numerical software and an in-depth understanding of the basic presented will give the student a foundation for more advanced methods that may be encountered in commercial codes.\n",
    "\n",
    "The first technique we will study is the Explicit Finite Difference method.  This is one of three common [finite difference methods](https://en.wikipedia.org/wiki/Finite_difference_method) that are easy to program and built from the definition of Taylor's polynomial."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Taylor Series and Derivatives\n",
    "----"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Taylor's approximation of a first and second derivatives are defined as central differences and are second order accurate:\n",
    "\n",
    "$$\n",
    "- 2 h \\left. \\frac{d}{d \\xi_{1}} f{\\left (\\xi_{1} \\right )} \\right|_{\\substack{ \\xi_{1}=c }} - f{\\left (c - h \\right )} + f{\\left (c + h \\right )} = 0\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "$$\n",
    "h^{2} \\left. \\frac{d^{2}}{d \\xi_{1}^{2}}  f{\\left (\\xi_{1} \\right )} \\right|_{\\substack{ \\xi_{1}=c }} + 2 f{\\left (c \\right )} - f{\\left (c - h \\right )} - f{\\left (c + h \\right )} = 0\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Dividing Space and Time\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The diffusion equation is a partial differential equation in two independent variables, space and time.  In order that we may use our Taylor's series approximations for the time and space derivatives we need to descretize the domain of the problem.  One easy way to visualize time and space descretization is to use a _grid_ construction.  The figure below shows one way that the time and space variables can be represented."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "![](./images/grid.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Each cell in the grid holds the value of the dependent parameter at a particular value of time and space.  We indicate this using the following symbols:\n",
    "\n",
    "$$\n",
    "u_{i,j}\n",
    "$$\n",
    "\n",
    "The $i$ and $j$ are the _indices_ of the grid and reference a particular location where the value of the dependent parameter is stored.  How much memory is required to store these values depends on the type of data and the size of the grid.  In addition, the _grid spacing_ must be specified for each of the independent variables, in this case we need both a $\\delta x$ and a $\\delta t$.  For example, a difference in time might be represented as:\n",
    "\n",
    "$$\n",
    "u(i,j) - u(i,j+1) = c(x,t) - c(x,t+\\delta t)\n",
    "$$\n",
    "\n",
    "Typically, the grid is a _uniform_ grid in each independent variable - meaning that the distance in the independent variable between grid points is the same in any one variable.  There are cases where non-uniform grids may be desirable.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Finite Difference Form for the Diffusion Equation\n",
    "----"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The finite difference representation of the diffusion equation with $u_{i,j}$ as the stored value that represents $c(x,t)$ is:\n",
    "\n",
    "$$\n",
    "\\frac{u_{i,\\, j+1} - u_{i,\\, j}}{\\Delta t} = D \\frac{u_{i - 1,\\, j} - 2 u_{i,\\, j} + u_{i + 1,\\, j}}{\\Delta x^2}\n",
    "$$\n",
    "\n",
    "Note:\n",
    "* A _forward_ difference is used on the LHS\n",
    "* A _central_ difference is used on the RHS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The general steps required to solve a finite difference equation are as follows:\n",
    "\n",
    "* Identify the geometry of the system\n",
    "* Define the initial conditions and boundary conditions\n",
    "* Write the difference equation in terms of the unknowns\n",
    "* Compute the unknowns subject to any stability requirments and boundary conditions\n",
    "* Store and update the results in the grid\n",
    "* Write the results to disk\n",
    "* Visualize the results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The following parameters are needed to solve the finite difference equation above:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "* `numberOfPoints` - the number of grid points within our computational domain\n",
    "* `numberOfIterations` - the number of timesteps within our computational domain\n",
    "* `lengthOfDomain` - the physical length of our computational domain\n",
    "* `dx` - the distance between successive grid points in our domain\n",
    "* `xPoints` - a linear space divided into numberOfPoints points\n",
    "* `initialCondition` - our starting distribution of solute (i.e. $c(x,0)$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "In addition to the simulation parameters, we start with an initial seed of concentration data. Unlike our other analytical strategies there are no coefficients to compute, no functions to fit.  The data could be generated from a function or from measurements.  Here we choose a `sin` function as our initial condition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "numberOfPoints = 100\n",
    "numberOfIterations = 1000\n",
    "lengthOfDomain = 1.0\n",
    "dx = lengthOfDomain/numberOfPoints\n",
    "xPoints = np.linspace(0.0, lengthOfDomain, numberOfPoints)\n",
    "initialCondition = np.sin(xPoints*np.pi/lengthOfDomain)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "def plotIC():\n",
    "    fig, ax = plt.subplots(figsize=(10,5))\n",
    "    ax.plot(xPoints, initialCondition, 'ro')\n",
    "    ax.set_xlabel(r'Distance $x$')\n",
    "    ax.set_ylabel(r'Concentration $c(x,t)$')\n",
    "    ax.set_title(r'Initial Conditions')\n",
    "    plt.show()\n",
    "    return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "plotIC()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Set up other simulation parameters\n",
    "----"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We now set:\n",
    "\n",
    "* `diffusionCoefficient` - the diffusion coefficient\n",
    "* `dt` - the discrete time step (the formulaic choice of `dt` is needed to satisfy stability constraints)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "diffusionCoefficient = 10.0\n",
    "dt = dx**2/(4*diffusionCoefficient)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Our choice of timestep is restricted:\n",
    "\n",
    "$$\n",
    "dt \\leq \\frac{\\Delta x^2}{2 \\, D}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Set up the data structure\n",
    "----"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "There are potentially three strategies for storing the results of numerical computations.  One is to store ALL the data, another is to store SOME of the data, and the last is to store NONE of the data except the very last computation.  Each strategy has advantages and disadvantages although all strategies may seem equally difficult to implement.  In this lecture we will design design a data structure that stores all the data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Let us create a `numpy` array that has one dimension equal to the number of points in our grid and another dimension that is equal to the number of iterations we wish to compute:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "arrayWithAllTheData = np.zeros((numberOfPoints,numberOfIterations), dtype='float32')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "You can think of the 2D array as having one space axis (the first index, we will use `i` for this one) and one time axis (the second index, we will use `j` for this one).\n",
    "\n",
    "We will set our initial conditions by assigning the `initialCondition` array to the first row of the `arrayWithAllTheData`.  Note the slicing in the first index so that we can copy the contents of the `initialCondition` into the whole first row with a single assignment statement."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "arrayWithAllTheData[:,0] = initialCondition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "With our initial conditions in place we need to develop the computational steps to advance our solution in time.  The PDE we are solving (with a constant diffusion coefficient) is:\n",
    "\n",
    "$$\n",
    "\\frac{\\partial c(x,t)}{\\partial t} = D \\frac{\\partial^2 c(x,t)}{\\partial x^2}\n",
    "$$\n",
    "\n",
    "we transform this into:\n",
    "\n",
    "$$\n",
    "\\frac{u_{i,\\, j+1} - u_{i,\\, j}}{\\Delta t} = D \\frac{u_{i - 1,\\, j} - 2 u_{i,\\, j} + u_{i + 1,\\, j}}{\\Delta x^2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "so that we can algebraically solve for $u_{i+1,\\, j}$:\n",
    "\n",
    "$$\n",
    "u_{i,\\, j+1} = \\frac{D \\Delta t}{\\Delta x^2} \\left( u_{i - 1,\\, j} - 2 u_{i,\\, j} + u_{i + 1,\\, j} \\right) + u_{i,\\, j}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "From the expression above you can see that all the terms on the RHS of the expression are at the index $j$ (the last iteration) and all the terms on the LHS are for the $j+1$ index (the next iteration).  This scheme defines a simple method (with a restrictive timestep) for integrating a PDE.  Re-examine the figure below in comparison to the finite difference scheme:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "![](./images/grid.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Write the Solver and Compute the Results\n",
    "----\n",
    "\n",
    "Add your code below.  As a reminder, the algebraic equation we want to solve is:\n",
    "\n",
    "$$\n",
    "u_{i,\\, j+1} = \\frac{D \\Delta t}{\\Delta x^2} \\left( u_{i - 1,\\, j} - 2 u_{i,\\, j} + u_{i + 1,\\, j} \\right) + u_{i,\\, j}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "To make all of this work we proceed as follows:\n",
    "\n",
    "* Compute the prefactor $D \\Delta t/ \\Delta x^2$.\n",
    "* Apply the boundary conditions in the $j$th row of the array.\n",
    "* Using the $j$ row of the array (plus the boundary conditions), fill in row $j+1$ of the array with values corresponding to the new time $t + \\Delta t$ according to the equation above.\n",
    "* Advance the index and repeat until all rows are filled.\n",
    "* Visualize the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Note the counting for j in this loop.  You may wish to print out\n",
    "# the values of i,j to help build this operation.\n",
    "for j in range(1,numberOfIterations):\n",
    "    for i in range(1,numberOfPoints-1):\n",
    "        arrayWithAllTheData[i,j] = 0 # What should you put here?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "### DIY:  Sketch the algorithm up to this point and for the cell below.\n",
    "----\n",
    "\n",
    "Doing this will help you visualize the operations and it will increase your ability to make modifications in the future and devise new more compact ways to integrate this PDE.\n",
    "\n",
    "If you've sketched the algorithm as advised above then you see that in our development of this solution we implicitly set the boundary conditions.  We initialize `arrayWithAllTheData` with `np.zeros` and then compute on all the interior rows/columns.  This creates a condition where all the boundary cells are set to zero and their values remain untouched throughout the computation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Plot the results\n",
    "----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "from ipywidgets import interact, fixed\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def plotArray(xPoints, dataArray, rowID=0):\n",
    "    \"\"\"\n",
    "    This function in conjunction with interact() permits\n",
    "    inspection of the contents of an array, row by row.  This\n",
    "    is useful for some small tasks such as examining the results\n",
    "    of a PDE solution.\n",
    "    \"\"\"\n",
    "    x = xPoints\n",
    "    y = dataArray[:,rowID]\n",
    "    fig = plt.figure(figsize=(7,4))\n",
    "    axes = fig.add_axes([0.1, 0.1, 0.8, 0.8])\n",
    "    axes.set_ylim(0,1)\n",
    "    axes.plot(x, y, 'ro', label=r\"$c(x,t)$\")\n",
    "    axes.legend()\n",
    "    axes.grid(False)\n",
    "    plt.show()\n",
    "    return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "interact(plotArray, \n",
    "         xPoints=fixed(xPoints), \n",
    "         dataArray=fixed(arrayWithAllTheData),\n",
    "         rowID=(0,numberOfIterations-1,1), );"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "### DIY:  Compute a solution where you change the boundary conditions on the LHS to be $c(x=L,t) = 1.0$.\n",
    "\n",
    "----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Your solver code goes here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "### DIY:  Vectorize the above solution method.\n",
    "----\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Your solver code goes here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "## Homework\n",
    "----\n",
    "\n",
    "Problem 1.  Solve the following problem:\n",
    "\n",
    "$$\n",
    "\\frac{\\partial c(x,t)}{\\partial t} = D \\frac{\\partial^2 c(x,t)}{\\partial x^2}\n",
    "$$\n",
    "\n",
    "with the initial condition \n",
    "\n",
    "$$c(x,t=0) = \\sin \\pi x,$$\n",
    "\n",
    "over the domain \n",
    "\n",
    "$$( x \\, \\vert \\, 0 \\le x \\le 1.0 ).$$\n",
    "\n",
    "with zero flux boundary conditions.\n",
    "\n",
    "----\n",
    "Problem 2.  Solve the following problem:\n",
    "\n",
    "$$\n",
    "\\frac{\\partial c(x,t)}{\\partial t} = \\frac{\\partial}{\\partial x} D(c) \\frac{\\partial c(x,t)}{\\partial x}\n",
    "$$\n",
    "\n",
    "with the initial condition:\n",
    "\n",
    "$$\n",
    "c(x \\leq 0.5 , t=0) = 1.0\\\\\n",
    "c(x \\gt 0.5 , t=0) = 0.0\n",
    "$$\n",
    "\n",
    "with D having the dependence:\n",
    "\n",
    "$$\n",
    "D(c) = D_0\\cdot c \\cdot(1-c)\n",
    "$$\n",
    "\n",
    "over the domain:\n",
    "\n",
    "$$\n",
    "( x \\, \\vert \\, 0 \\le x \\le 1.0 )\n",
    "$$\n",
    "\n",
    "with zero flux boundary conditions."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}