{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# 06 Differentiation – Student Activity\n",
    "See *Computational Physics* (Landau, Páez, Bordeianu), Chapter 5.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "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{ep} 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": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def D_fd(y, t, h):\n",
    "    \"\"\"Forward difference\"\"\"\n",
    "    return (y(t + h) - y(t))/h\n",
    "\n",
    "def D_cd(y, t, h):\n",
    "    \"\"\"Central difference\"\"\"\n",
    "    # implement\n",
    "\n",
    "def D_ep(y, t, h):\n",
    "    \"\"\"Extended difference\"\"\"\n",
    "    # implement"
   ]
  },
  {
   "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)}}$ (compared to the analystical value – 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": [
    "#### Function definitions "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "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",
    "def y2(t):\n",
    "    # IMPLEMENT ME (remove `pass`)\n",
    "    pass\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": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Evaluate the finite difference derivatives\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": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "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": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "y2(t_values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Calculation of the **absolute error**: subtract the two arrays that you got previously:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# IMPLEMENT ME"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Calculate the relative error $E$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "def error(Dxx, y, y2, t, h):\n",
    "    \"\"\"Relative error.\"\"\"\n",
    "    # INCOMPLETE ... replace None with appropriate terms\n",
    "    return (Dxx(y, t, h) - None)/None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "error(D_fd, y, y2, 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": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "h_values = 10**(np.arange(-15, -1, 0.1))\n",
    "abs_errors = np.abs(error(D_fd, y, y2, 0.1, h_values))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# INCOMPLETE, replace None with appropriate terms\n",
    "plt.loglog(None, None, label=\"t=0.1\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the three different $t$ values together in one plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# INCOMPLETE: replace None values by appropriate terms\n",
    "for t in t_values:\n",
    "    abs_errors = None\n",
    "    plt.loglog(None, None, label=r\"$t={}$\".format(t))\n",
    "ax = plt.gca()\n",
    "ax.legend(loc=\"best\")\n",
    "ax.set_xlabel(r\"$h$\")\n",
    "ax.set_ylabel(r\"$|E|$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* error behavior depends on $t$ and on cancellation of errors (e.g. for $t=1$\n",
    "* algorithmic error decreases for decreasing $h$ until the round of error starts dominating"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}