differentiation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Activity 11: Differentiation: Finite Difference Operators\n",
    "See *Computational Physics* (Landau, Páez, Bordeianu), Chapter 5.5\n",
    "\n",
    "The background is stated in notebook [09-differentiation-theory.ipynb](https://nbviewer.jupyter.org/github/ASU-CompMethodsPhysics-PHY494/PHY494-resources/blob/master/09_differentiation/09-differentiation-theory.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Implementation of Finite Difference Algorithms in Python\n",
    "\n",
    "*Forward difference* and *central difference* (discussed in class):\n",
    "\n",
    "\\begin{align}\n",
    "D_\\text{fd} y(t) &\\equiv \\frac{y(t+h) - y(t)}{h} \\\\\n",
    "D_\\text{cd} y(t) &\\equiv \\frac{y\\Big(t + \\frac{h}{2}\\Big) - y\\Big(t - \\frac{h}{2}\\Big)}{h}\n",
    "\\end{align}\n",
    "\n",
    "and also the *extended difference algorithm*\n",
    "\n",
    "\\begin{align}\n",
    "D_\\text{ed} y(t) &\\equiv \\frac{4 D_\\text{cd}y(t, h/2) - D_\\text{cd}y(t, h)}{3} \\\\\n",
    "  &= \\frac{8\\big(y(t+h/4) - y(t-h/4)\\big) - \\big(y(t+h/2) - y(t-h/2)\\big)}{3h}\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Implement the three finite difference operators `D_fd`, `D_cd`, and `D_ed` in the module **`diffops.py`**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from importlib import reload"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import diffops\n",
    "reload(diffops)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Test your implementations\n",
    "Test function: $y(t) = \\cos t$\n",
    "1. What is the analytical derivative $\\frac{d\\cos(t)}{dt}$?\n",
    "1. Calculate the derivative of $y(t) = \\cos t$ at $t=0.1, 1, 100$.\n",
    "1. Print derivative and relative error $E = \\frac{D y(t) - y^{(1)}(t)}{y^{(1)}(t)}$ (finite difference value $D y(t)$ compared to the analystical value $y^{(1)}(t)$– use numpy functions for \"exact\" values) as function of $h$.\n",
    "1. Reduce $h$ until you reach machine precision, $h \\approx \\epsilon_m$\n",
    "1. Plot $\\log_{10} |E(h)|$ against $\\log_{10} h$.\n",
    "\n",
    "Try to do the above for all three algorithms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Test function definitions \n",
    "\n",
    "For analyzing your finite difference operators, define the test function $y(t) = \\cos t$ and the analytical derivative function $y^{(1)}(t) = \\frac{dy}{dt}$. We also set up the $t$ values that we want to analyze:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "# test function\n",
    "y = np.cos\n",
    "\n",
    "# Analytical derivative (use np.cos, np.sin, np.exp or whatever else you need)\n",
    "# replace NotImplemented\n",
    "def y1(t):\n",
    "    return NotImplemented\n",
    "\n",
    "t_values = np.array([0.1, 1, 100], dtype=np.float64)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use numpy functions for everything because then you can operate on all `t_values` at once."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Forward difference "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Evaluate the finite difference derivatives\n",
    "\n",
    "Calculate $D_\\text{fd}\\cos(t)$ for $h=0.1$ for our `t_values`.\n",
    "\n",
    "Note that we pass *a function* `y` to the forward difference function `D_fd` and we can also pass a whole array of `t_values`!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "diffops.D_fd(y, t_values, 0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Evaluate the exact derivatives\n",
    "Compute the exact derivatives (again, operate on all $t$ together... start thinking in numpy arrays!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "y1(t_values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Calculation of the **absolute error**: subtract the two arrays that you got previously:\n",
    "\n",
    "$$\n",
    "\\Delta = D y(t) - y^{(1)}(t)\n",
    "$$\n",
    "\n",
    "(finite difference value $D y(t)$ compared to the analystical value $y^{(1)}(t)$– use numpy functions for \"exact\" values)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Calculate the relative error $E$\n",
    "\n",
    "$$\n",
    "E = \\frac{\\Delta}{y^{(1)}(t)} = \\frac{D y(t) - y^{(1)}(t)}{y^{(1)}(t)}\n",
    "$$\n",
    "\n",
    "Implement a function `diffops.error(Dxx, y, y1, t, h)` that calculates the error $E$ for function $y(t)$ with finite difference operator $D$ and analytical derivative $y^{(1)}$ for all values $t$ at step size $h$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Open `diffops.py` in your editor and complete the `error()` function in the file.\n",
    "\n",
    "Look at the template code for `error()`: Note that we pass again a general function for the difference operator so that we can use `error()` with `D_fd()`, `D_cd()` and `D_ep()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import diffops\n",
    "reload(diffops)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "diffops.error(diffops.D_fd, y, y1, t_values, 0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Plot $|E|$\n",
    "Plot $\\log_{10} |E(h)|$ against $\\log_{10} h$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "plt.matplotlib.style.use(\"ggplot\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are looking for the error over a range of $h$ values, $10^{-15} \\le h \\le 10^{-1}$ and at $t=0.1$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h_values = 10**(np.linspace(-15, -1, 141))\n",
    "\n",
    "# compute your relative errors for all h_values and t=0.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot in a log-log plot $|E(h)|$ (for $t=0.1$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the three different $t$ values together in one plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# INCOMPLETE: replace NotImplemented values by appropriate terms\n",
    "for t in NotImplemented:\n",
    "    rel_errors = NotImplemented\n",
    "    plt.loglog(NotImplemented, NotImplemented, label=r\"$t={}$\".format(t))\n",
    "    \n",
    "ax = plt.gca()\n",
    "ax.legend(loc=\"best\")\n",
    "ax.set_xlabel(r\"$h$\")\n",
    "ax.set_ylabel(r\"$|E|$\")\n",
    "ax.set_title(r\"Forward difference $D_\\mathrm{fd}$\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Conclusions?\n",
    "\n",
    "1. ...\n",
    "2. ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Central difference \n",
    "\n",
    "For the central difference algorithm, plot in a log-log plot $|E(h)|$ for all $t$ values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Extended difference \n",
    "\n",
    "For the extended difference algorithm, plot in a log-log plot $|E(h)|$ for all $t$ values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "celltoolbar": "Slideshow",
  "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}