03_finite_differences.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "<table>\n",
    " <tr align=left><td><img align=left src=\"https://i.creativecommons.org/l/by/4.0/88x31.png\">\n",
    " <td>Text provided under a Creative Commons Attribution license, CC-BY. All code is made available under the FSF-approved MIT license. (c) Kyle T. Mandli</td>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Review:  Finite Differences"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "Finite differences are expressions that approximate derivatives of a function evaluated at a set of points, often called a *stencil*.  These expressions can come in many different flavors including types of stencils, order of accuracy, and order of derivatives.  In this lecture we will review the process of derivation, error analysis and application of finite differences."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Derivation of Finite Differences\n",
    "\n",
    "The general approach to deriving finite differences should be familiar for at least the first order differences.  Consider three different ways to define a derivative at a point $x_i$\n",
    "$$\n",
    "    u'(x_i) = \\lim_{\\Delta x \\rightarrow 0} \\left \\{ \\begin{aligned} \n",
    "        &\\frac{u(x_i + \\Delta x) - u(x_i)}{\\Delta x} & \\equiv D_+ u(x_i)\\\\\n",
    "        &\\frac{u(x_i + \\Delta x) - u(x_i - \\Delta_x)}{2 \\Delta x} & \\equiv D_0 u(x_i)\\\\\n",
    "        &\\frac{u(x_i) - u(x_i - \\Delta_x)}{\\Delta x} & \\equiv D_- u(x_i).\n",
    "    \\end{aligned} \\right .\n",
    "$$\n",
    "\n",
    "![Approximations to $u'(x)$](./images/fd_basic.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "If instead of allowing $\\Delta x \\rightarrow 0$ we come up with an approximation to the slope $u'(x_i)$ and hence our definitions of derivatives can directly be seen as approximations to derivatives when $\\Delta x$ is perhaps small but non-zero.\n",
    "\n",
    "For the rest of the review we will delve into a more systematic way to derive these approximations as well as find higher order accurate approximations, higher order derivative approximations, and understand the error associated with the approximations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Interpolating Polynomials\n",
    "\n",
    "One was to derive finite difference approximations is by finding an interpolating polynomial through the given stencil and differentiating that directly.  Given $N+1$ points $(x_0,u(x_0)), (x_1,u(x_1)), \\ldots, (x_{N},u(x_{N}))$ assuming the $x_i$ are all unique, the interpolating polynomial $P_N(x)$ can be written as\n",
    "\n",
    "$$\n",
    "    P_N(x) = \\sum^{N}_{i=0} u(x_i) \\ell_i(x)\n",
    "$$\n",
    "\n",
    "where \n",
    "\n",
    "$$\n",
    "    \\ell_i(x) = \\prod^{N}_{j=0, j \\neq i} \\frac{x - x_j}{x_i - x_j} = \\frac{x - x_0}{x_i - x_0} \\frac{x - x_1}{x_i - x_1} \\cdots \\frac{x - x_{i-1}}{x_i - x_{i-1}}\\frac{x - x_{i+1}}{x_i - x_{i+1}} \\cdots \\frac{x - x_{N}}{x_i - x_{N}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "Note that $\\ell_i(x_i) = 1$ and $\\forall j\\neq i, ~~ \\ell_i(x_j) = 0$.  \n",
    "\n",
    "Since we know how to differentiate a polynomial we should be able to then compute the given finite difference approximation given these data points."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Example: 2-Point Stencil\n",
    "\n",
    "Say we have two points to form the approximation to the derivative with.  The interpolating polynomial through two points is a linear function with the form\n",
    "$$\n",
    "    P_1(x) = u(x_0) \\frac{x - x_1}{x_0 - x_1} + u(x_1) \\frac{x - x_0}{x_1 - x_0}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Differentiating $P_1(x)$ leads to\n",
    "$$\n",
    "    P'_1(x) = u(x_0) \\frac{1}{x_0 - x_1} + u(x_1) \\frac{1}{x_1 - x_0}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "If we allow the spacing between $x_0$ and $x_1$ to be $\\Delta x = x_1 - x_0$ we can then write this as\n",
    "$$\n",
    "    P'_1(x) = \\frac{u(x_1) - u(x_0)}{\\Delta x}\n",
    "$$\n",
    "which is the general form of $D_-u(x)$ and $D_+u(x)$ above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "If we extend this to have three points we have the interpolating polynomial\n",
    "$$\n",
    "    P_2(x) = u(x_0) \\frac{x - x_1}{x_0 - x_1} \\frac{x - x_2}{x_0 - x_2} + u(x_1) \\frac{x - x_0}{x_1 - x_0} \\frac{x - x_2}{x_1 - x_2} + u(x_2) \\frac{x - x_0}{x_2 - x_0} \\frac{x - x_1}{x_2 - x_1}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Differentiating this leads to\n",
    "$$\\begin{aligned}\n",
    "    P'_2(x) &= u(x_0) \\left( \\frac{1}{x_0 - x_1} \\frac{x - x_2}{x_0 - x_2} + \\frac{x - x_1}{x_0 - x_1} \\frac{1}{x_0 - x_2}\\right )+ u(x_1) \\left ( \\frac{1}{x_1 - x_0} \\frac{x - x_2}{x_1 - x_2} + \\frac{x - x_0}{x_1 - x_0} \\frac{1}{x_1 - x_2} \\right )+ u(x_2)\\left ( \\frac{1}{x_2 - x_0} \\frac{x - x_1}{x_2 - x_1} + \\frac{x - x_0}{x_2 - x_0} \\frac{1}{x_2 - x_1} \\right ) \\\\\n",
    "    &= u(x_0) \\left(\\frac{x - x_2}{2 \\Delta x^2} + \\frac{x - x_1}{2 \\Delta x^2} \\right )+ u(x_1) \\left ( \\frac{x - x_2}{-\\Delta x^2} + \\frac{x - x_0}{-\\Delta x^2} \\right )+ u(x_2)\\left ( \\frac{x - x_1}{2\\Delta x^2} + \\frac{x - x_0}{2 \\Delta x^2} \\right ) \\\\\n",
    "    &=\\frac{u(x_0)}{2\\Delta x^2} (2x - x_2 - x_1)+ \\frac{u(x_1)}{-\\Delta x^2} ( 2x - x_2 - x_0)+ \\frac{u(x_2)}{2\\Delta x^2}( 2x - x_1 - x_0) \\\\\n",
    "    &=\\frac{u(x_0)}{2\\Delta x^2} (2x - x_2 - x_1)+ \\frac{u(x_1)}{-\\Delta x^2} ( 2x - x_2 - x_0)+ \\frac{u(x_2)}{2\\Delta x^2}( 2x - x_1 - x_0).\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "If we now evaluate the derivative at $x_1$, assuming this is the central point, we have\n",
    "$$\\begin{aligned}\n",
    "    P'_2(x_1) &= \\frac{u(x_0)}{2\\Delta x^2} (x_1 - x_2)+ \\frac{u(x_1)}{-\\Delta x^2} ( x_1 - x_2 + x_1 - x_0)+ \\frac{u(x_2)}{\\Delta x^2}( x_1 - x_0) \\\\\n",
    "    &= \\frac{u(x_0)}{2\\Delta x^2} (-\\Delta x)+ \\frac{u(x_1)}{-\\Delta x^2} ( -\\Delta x + \\Delta x)+ \\frac{u(x_2)}{\\Delta x^2}( 2\\Delta x) \\\\\n",
    "    &= \\frac{u(x_2) - u(x_0)}{2 \\Delta x}\n",
    "\\end{aligned}$$\n",
    "giving us the third approximation from above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Taylor-Series Methods\n",
    "\n",
    "Another way to derive finite difference approximations can be computed by using the Taylor series and the method of undetermined coefficients.\n",
    "\n",
    "$$u(x) = u(x_n) + (x - x_n) u'(x_n) + \\frac{(x - x_n)^2}{2!} u''(x_n) + \\frac{(x - x_n)^3}{3!} u'''(x_n) + \\mathcal{O}((x - x_n)^4)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Say we want to derive the second order accurate, first derivative approximation that just did, this requires the values $(x_{n+1}, u(x_{n+1})$ and $(x_{n-1}, u(x_{n-1})$.  We can express these values via our Taylor series approximation above as\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    u(x_{n+1}) &= u(x_n) + (x_{n+1} - x_n) u'(x_n) + \\frac{(x_{n+1} - x_n)^2}{2!} u''(x_n) + \\frac{(x_{n+1} - x_n)^3}{3!} u'''(x_n) + \\mathcal{O}((x_{n+1} - x_n)^4) \\\\\n",
    "    &= u(x_n) + \\Delta x u'(x_n) + \\frac{\\Delta x^2}{2!} u''(x_n) + \\frac{\\Delta x^3}{3!} u'''(x_n) + \\mathcal{O}(\\Delta x^4) \n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "and \n",
    "\n",
    "$$\\begin{aligned}\n",
    "    u(x_{n-1}) &= u(x_n) + (x_{n-1} - x_n) u'(x_n) + \\frac{(x_{n-1} - x_n)^2}{2!} u''(x_n) + \\frac{(x_{n-1} - x_n)^3}{3!} u'''(x_n) + \\mathcal{O}((x_{n-1} - x_n)^4) \\\\\n",
    "&= u(x_n) - \\Delta x u'(x_n) + \\frac{\\Delta x^2}{2!} u''(x_n) - \\frac{\\Delta x^3}{3!} f'''(x_n) + \\mathcal{O}(\\Delta x^4) \n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Now to find out how to combine these into an expression for the derivative we assume our approximation looks like\n",
    "\n",
    "$$u'(x_n) + R(x_n) = A u(x_{n+1}) + B u(x_n) + C u(x_{n-1})$$\n",
    "\n",
    "where $R(x_n)$ is our error.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Plugging in the Taylor series approximations we find\n",
    "\n",
    "$$u'(x_n) + R(x_n) = A \\left ( u(x_n) + \\Delta x u'(x_n) + \\frac{\\Delta x^2}{2!} u''(x_n) + \\frac{\\Delta x^3}{3!} u'''(x_n) + \\mathcal{O}(\\Delta x^4)\\right ) + B u(x_n) + C \\left ( u(x_n) - \\Delta x u'(x_n) + \\frac{\\Delta x^2}{2!} u''(x_n) - \\frac{\\Delta x^3}{3!} u'''(x_n) + \\mathcal{O}(\\Delta x^4) \\right )$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Since we want $R(x_n) = \\mathcal{O}(\\Delta x^2)$ we want all terms lower than this to disappear except for those multiplying $u'(x_n)$ as those should sum to 1 to give us our approximation.  Collecting the terms with common $\\Delta x^n$ we get a series of expressions for the coefficients $A$, $B$, and $C$ based on the fact we want an approximation to $u'(x_n)$.  The $n=0$ terms collected are $A + B + C$ and are set to 0 as we want the $u(x_n)$ term to disappear\n",
    "\n",
    "$$\\Delta x^0: ~~~~ A + B + C = 0$$\n",
    "\n",
    "$$\\Delta x^1: ~~~~ A \\Delta x - C \\Delta x = 1 $$\n",
    "\n",
    "$$\\Delta x^2: ~~~~ A \\frac{\\Delta x^2}{2} + C \\frac{\\Delta x^2}{2} = 0 $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "This last equation $\\Rightarrow A = -C$, using this in the second equation gives $A = \\frac{1}{2 \\Delta x}$ and $C = -\\frac{1}{2 \\Delta x}$.  The first equation then leads to $B = 0$.  Putting this altogether then gives us our previous expression including an estimate for the error:\n",
    "\n",
    "$$u'(x_n) + R(x_n) = \\frac{u(x_{n+1}) - u(x_{n-1})}{2 \\Delta x} + \\frac{1}{2 \\Delta x} \\frac{\\Delta x^3}{3!} u'''(x_n) + \\mathcal{O}(\\Delta x^4) + \\frac{1}{2 \\Delta x} \\frac{\\Delta x^3}{3!} u'''(x_n) + \\mathcal{O}(\\Delta x^4) $$\n",
    "\n",
    "$$R(x_n) = \\frac{\\Delta x^2}{3!} u'''(x_n) + \\mathcal{O}(\\Delta x^3) = \\mathcal{O}(\\Delta x^2)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example: First Order Derivatives"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "f = lambda x: numpy.sin(x)\n",
    "f_prime = lambda x: numpy.cos(x)\n",
    "\n",
    "# Use uniform discretization\n",
    "x = numpy.linspace(-2 * numpy.pi, 2 * numpy.pi, 1000)\n",
    "N = 20\n",
    "x_hat = numpy.linspace(-2 * numpy.pi, 2 * numpy.pi, N)\n",
    "delta_x = x_hat[1] - x_hat[0]\n",
    "\n",
    "# Compute forward difference using a loop\n",
    "f_prime_hat = numpy.empty(x_hat.shape)\n",
    "for i in xrange(N - 1):\n",
    "    f_prime_hat[i] = (f(x_hat[i+1]) - f(x_hat[i])) / delta_x\n",
    "f_prime_hat[-1] = (f(x_hat[i]) - f(x_hat[i-1])) / delta_x\n",
    "\n",
    "# Vector based calculation\n",
    "# f_prime_hat[:-1] = (f(x_hat[1:]) - f(x_hat[:-1])) / (delta_x)\n",
    "\n",
    "# Use first-order differences for points at edge of domain\n",
    "f_prime_hat[-1] = (f(x_hat[-1]) - f(x_hat[-2])) / delta_x  # Backward Difference at x_N\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "axes.plot(x, f_prime(x), 'k')\n",
    "axes.plot(x_hat + 0.5 * delta_x, f_prime_hat, 'ro')\n",
    "axes.set_xlim((x[0], x[-1]))\n",
    "axes.set_ylim((-1.1, 1.1))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Example: Second Order Derivative\n",
    "\n",
    "Using our Taylor series approach lets derive the second order accurate second derivative formula.  Again we will use the same points and the Taylor series centered at $x = x_n$ so we end up with the same expression as before:\n",
    "\n",
    "$$u''(x_n) + R(x_n) = A \\left ( u(x_n) + \\Delta x u'(x_n) + \\frac{\\Delta x^2}{2!} u''(x_n) + \\frac{\\Delta x^3}{3!} u'''(x_n) + \\frac{\\Delta x^4}{4!} u^{(4)}(x_n) + \\mathcal{O}(\\Delta x^5)\\right ) + B u(x_n) + C \\left ( u(x_n) - \\Delta x u'(x_n) + \\frac{\\Delta x^2}{2!} u''(x_n) - \\frac{\\Delta x^3}{3!} u'''(x_n) + \\frac{\\Delta x^4}{4!} u^{(4)}(x_n) + \\mathcal{O}(\\Delta x^5) \\right )$$\n",
    "\n",
    "except this time we want to leave $u''(x_n)$ on the right hand side.  Doing the same trick as before we have the following expressions:\n",
    "\n",
    "$$\\Delta x^0: ~~~~ A + B + C = 0$$\n",
    "\n",
    "$$\\Delta x^1: ~~~~ A \\Delta x - C \\Delta x = 0$$\n",
    "\n",
    "$$\\Delta x^2: ~~~~ A \\frac{\\Delta x^2}{2} + C \\frac{\\Delta x^2}{2} = 1$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The second equation implies $A = C$ which combined with the third implies\n",
    "\n",
    "$$A = C = \\frac{1}{\\Delta x^2}$$\n",
    "\n",
    "Finally the first equation gives\n",
    "\n",
    "$$B = -\\frac{2}{\\Delta x^2}$$\n",
    "\n",
    "leading to the final expression\n",
    "\n",
    "$$u''(x_n) + R(x_n) = \\frac{u(x_{n+1}) - 2 u(x_n) + u(x_{n-1})}{\\Delta x^2} + \\frac{1}{\\Delta x^2} \\left(\\frac{\\Delta x^3}{3!} u'''(x_n) + \\frac{\\Delta x^4}{4!} u^{(4)}(x_n)  - \\frac{\\Delta x^3}{3!} u'''(x_n) + \\frac{\\Delta x^4}{4!} u^{(4)}(x_n)  \\right) + \\mathcal{O}(\\Delta x^5)$$\n",
    "\n",
    "with\n",
    "\n",
    "$$R(x_n) = \\frac{\\Delta x^2}{12} u^{(4)}(x_n) + \\mathcal{O}(\\Delta x^3)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "f = lambda x: numpy.sin(x)\n",
    "f_dubl_prime = lambda x: -numpy.sin(x)\n",
    "\n",
    "# Use uniform discretization\n",
    "x = numpy.linspace(-2 * numpy.pi, 2 * numpy.pi, 1000)\n",
    "N = 10\n",
    "x_hat = numpy.linspace(-2 * numpy.pi, 2 * numpy.pi, N)\n",
    "delta_x = x_hat[1] - x_hat[0]\n",
    "\n",
    "# Compute derivative\n",
    "f_dubl_prime_hat = numpy.empty(x_hat.shape)\n",
    "f_dubl_prime_hat[1:-1] = (f(x_hat[2:]) -2.0 * f(x_hat[1:-1]) + f(x_hat[:-2])) / (delta_x**2)\n",
    "\n",
    "# Use first-order differences for points at edge of domain\n",
    "f_dubl_prime_hat[0] = (2.0 * f(x_hat[0]) - 5.0 * f(x_hat[1]) + 4.0 * f(x_hat[2]) - f(x_hat[3])) / delta_x**2\n",
    "f_dubl_prime_hat[-1] = (2.0 * f(x_hat[-1]) - 5.0 * f(x_hat[-2]) + 4.0 * f(x_hat[-3]) - f(x_hat[-4])) / delta_x**2\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "axes.plot(x, f_dubl_prime(x), 'k')\n",
    "axes.plot(x_hat, f_dubl_prime_hat, 'ro')\n",
    "axes.set_xlim((x[0], x[-1]))\n",
    "axes.set_ylim((-1.1, 1.1))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### General Derivation\n",
    "\n",
    "For a general finite difference approximation located at $\\bar{x}$ to the $k$th derivative with the arbitrary stencil $N \\geq k + 1$ points $x_1, \\ldots, x_N$ we can use some generalizations of the above method.  Note that although it is common that $\\bar{x}$ is one of the stencil points this is not necessary.  We also assume that $u(x)$ is sufficiently smooth so that our Taylor series are valid."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "At each stencil point we have the approximation\n",
    "$$\n",
    "    u(x_i) = u(\\bar{x}) + (x_i - \\bar{x})u'(\\bar{x}) + \\cdots + \\frac{1}{k!}(x_i - \\bar{x})^k u^{(k)}(\\bar{x}) + \\cdots.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Following our methodology above we want to find the linear combination of these Taylor series expansions such that\n",
    "$$\n",
    "    u^{(k)}(\\bar{x}) + \\mathcal{O}(\\Delta x^p) = a_1 u(x_1) + a_2 u(x_2) + a_3 u(x_3) + \\cdots + a_n u(x_n).\n",
    "$$\n",
    "Note that $\\Delta x$ can vary in general and the asymptotic behavior of the method will be characterized by some sort of average distance or sometimes the maximum distance between the stencil points."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Generalizing the approach above with the method of undetermined coefficients we want to eliminate the pieces of the above approximation that are in front of the derivatives less than order $k$.  The condition for this is\n",
    "$$\n",
    "    \\frac{1}{(i - 1)!} \\sum^N_{j=1} a_j (x_j - \\bar{x})^{(i-1)} = \\left \\{ \\begin{aligned}\n",
    "        1 & & \\text{if}~i - 1 = k, \\\\\n",
    "        0 & & \\text{otherwise}\n",
    "    \\end{aligned} \\right .\n",
    "$$\n",
    "for $i=1, \\ldots, N$.  Assuming the $x_j$ are distinct we can write the system of equations in a Vandermonde system which will have a unique solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "import scipy.special\n",
    "\n",
    "def finite_difference(k, x_bar, x):\n",
    "    \"\"\"Compute the finite difference stencil for the kth derivative\"\"\"\n",
    "    \n",
    "    N = x.shape[0]\n",
    "    A = numpy.ones((N, N))\n",
    "    x_row = x - x_bar\n",
    "    for i in range(1, N):\n",
    "        A[i, :] = x_row ** i / scipy.special.factorial(i)\n",
    "    b = numpy.zeros(N)\n",
    "    b[k] = 1.0\n",
    "    c = numpy.linalg.solve(A, b)\n",
    "    return c\n",
    "\n",
    "print finite_difference(2, 0.0, numpy.asarray([-1.0, 0.0, 1.0]))\n",
    "print finite_difference(1, 0.0, numpy.asarray([-1.0, 0.0, 1.0]))\n",
    "print finite_difference(1, -2.0, numpy.asarray([-2.0, -1.0, 0.0, 1.0, 2.0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Error Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Polynomial View\n",
    "\n",
    "Given $N + 1$ points we can form an interpolant $P_N(x)$ of degree $N$ where\n",
    "\n",
    "$$u(x) = P_N(x) + R_N(x)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We know from Lagrange's Theorem that the remainder term looks like\n",
    "\n",
    "$$R_N(x) = (x - x_0)(x - x_1)\\cdots (x - x_{N})(x - x_{N+1}) \\frac{u^{(N+1)}(c)}{(N+1)!}$$\n",
    "\n",
    "noting that we need to require that $u(x) \\in C^{N+1}$ on the interval of interest.  Taking the derivative of the interpolant $P_N(x)$ (in terms of Newton polynomials) then leads to \n",
    "\n",
    "$$P_N'(x) = [u(x_0), u(x_1)] + ((x - x_1) + (x - x_0)) [u(x_0), u(x_1), u(x_2)] + \\cdots + \\left(\\sum^{N-1}_{i=0}\\left( \\prod^{N-1}_{j=0,~j\\neq i} (x - x_j) \\right )\\right ) [u(x_0), u(x_1), \\ldots, u(x_N)]$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Similarly we can find the derivative of the remainder term $R_N(x)$ as\n",
    "\n",
    "$$R_N'(x) = \\left(\\sum^{N}_{i=0} \\left( \\prod^{N}_{j=0,~j\\neq i} (x - x_j) \\right )\\right ) \\frac{u^{(N+1)}(c)}{(N+1)!}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Now if we consider the approximation of the derivative evaluated at one of our data points $(x_k, y_k)$ these expressions simplify such that\n",
    "\n",
    "$$u'(x_k) = P_N'(x_k) + R_N'(x_k)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "If we let $\\Delta x = \\max_i |x_k - x_i|$ we then know that the remainder term will be $\\mathcal{O}(\\Delta x^N)$ as $\\Delta x \\rightarrow 0$ thus showing that this approach converges and we can find arbitrarily high order approximations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Truncation Error\n",
    "\n",
    "If we are using a Taylor series approach we can also look at the dominate term left over from in the Taylor series to find the *truncation error*.\n",
    "\n",
    "As an example lets again consider the first derivative approximations above, we need the Taylor expansions\n",
    "$$\n",
    "    u(\\bar{x} + \\Delta x) = u(\\bar{x}) + \\Delta x u'(\\bar{x}) + \\frac{1}{2} \\Delta x^2 u''(\\bar{x}) + \\frac{1}{3!} \\Delta x^3 u'''(\\bar{x}) + \\mathcal{O}(\\Delta x^4)\n",
    "$$\n",
    "and\n",
    "$$\n",
    "    u(\\bar{x} - \\Delta x) = u(\\bar{x}) - \\Delta x u'(\\bar{x}) + \\frac{1}{2} \\Delta x^2 u''(\\bar{x}) - \\frac{1}{3!} \\Delta x^3 u'''(\\bar{x}) + \\mathcal{O}(\\Delta x^4).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Plugging these into our expressions we have\n",
    "$$\\begin{aligned}\n",
    "    D_+ u(\\bar{x}) &= \\frac{u(\\bar{x} + \\Delta x) - u(\\bar{x})}{\\Delta x} \\\\\n",
    "    &= \\frac{\\Delta x u'(\\bar{x}) + \\frac{1}{2} \\Delta x^2 u''(\\bar{x}) + \\frac{1}{3!} \\Delta x^3 u'''(\\bar{x}) + \\mathcal{O}(\\Delta x^4)}{\\Delta x} \\\\\n",
    "    &= u'(\\bar{x}) + \\frac{1}{2} \\Delta x u''(\\bar{x}) + \\frac{1}{3!} \\Delta x^2 u'''(\\bar{x}) + \\mathcal{O}(\\Delta x^3).\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "If we now difference $D_+ u(\\bar{x}) - u'(\\bar{x})$ we get the truncation error\n",
    "$$\n",
    "    \\frac{1}{2} \\Delta x u''(\\bar{x}) + \\frac{1}{3!} \\Delta x^2 u'''(\\bar{x}) + \\mathcal{O}(\\Delta x^3)\n",
    "$$\n",
    "so the error for $D_+$ goes as $\\mathcal{O}(\\Delta x)$ and is controlled by $u''(\\bar{x})$.  Note that this approximation is dependent on $\\Delta x$ as the derivatives evaluated at $\\bar{x}$ are constants."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Similarly for the centered approximation we have\n",
    "$$\n",
    "    D_0 u(\\bar{x}) - u'(\\bar{x}) = \\frac{1}{6} \\Delta x^2 u'''(\\bar{x}) + \\mathcal{O}(\\Delta x^4).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Computing Order of Accuracy Graphically\n",
    "\n",
    "Model the error as\n",
    "$$\\begin{aligned}\n",
    "    e(\\Delta x) &= C \\Delta x^n \\\\\n",
    "    \\log e(\\Delta x) &= \\log C + n \\log \\Delta x\n",
    "\\end{aligned}$$\n",
    "\n",
    "Slope of line is $n$ when computing this!  We can also match the first point by solving for $C$:\n",
    "\n",
    "$$C = e^{\\log e(\\Delta x) - n \\log \\Delta x}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "f = lambda x: numpy.sin(x) + x**2 + 3.0 * x**3\n",
    "f_prime = lambda x: numpy.cos(x) + 2.0 * x + 9.0 * x**2\n",
    "\n",
    "# Compute the error as a function of delta_x\n",
    "delta_x = []\n",
    "error = []\n",
    "# for N in xrange(2, 101):\n",
    "for N in xrange(50, 1000, 50):\n",
    "    x_hat = numpy.linspace(-2 * numpy.pi, 2 * numpy.pi, N)\n",
    "    delta_x.append(x_hat[1] - x_hat[0])\n",
    "\n",
    "    # Compute forward difference\n",
    "    f_prime_hat = numpy.empty(x_hat.shape)\n",
    "    f_prime_hat[:-1] = (f(x_hat[1:]) - f(x_hat[:-1])) / (delta_x[-1])\n",
    "\n",
    "    # Use first-order differences for points at edge of domain\n",
    "    f_prime_hat[-1] = (f(x_hat[-1]) - f(x_hat[-2])) / delta_x[-1]  # Backward Difference at x_N\n",
    "    \n",
    "    error.append(numpy.linalg.norm(numpy.abs(f_prime(x_hat + delta_x[-1]) - f_prime_hat), ord=numpy.infty))\n",
    "    \n",
    "error = numpy.array(error)\n",
    "delta_x = numpy.array(delta_x)\n",
    "    \n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "axes.loglog(delta_x, error, 'ko', label=\"Approx. Derivative\")\n",
    "\n",
    "order_C = lambda delta_x, error, order: numpy.exp(numpy.log(error) - order * numpy.log(delta_x))\n",
    "axes.loglog(delta_x, order_C(delta_x[0], error[0], 1.0) * delta_x**1.0, 'r--', label=\"1st Order\")\n",
    "axes.loglog(delta_x, order_C(delta_x[0], error[0], 2.0) * delta_x**2.0, 'b--', label=\"2nd Order\")\n",
    "axes.legend(loc=4)\n",
    "axes.set_title(\"Convergence of 1st Order Differences\")\n",
    "axes.set_xlabel(\"$\\Delta x$\")\n",
    "axes.set_ylabel(\"$|f'(x) - \\hat{f}'(x)|$\")\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "marked": false
    },
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "f = lambda x: numpy.sin(x) + x**2 + 3.0 * x**3\n",
    "f_prime = lambda x: numpy.cos(x) + 2.0 * x + 9.0 * x**2\n",
    "\n",
    "# Compute the error as a function of delta_x\n",
    "delta_x = []\n",
    "error = []\n",
    "# for N in xrange(2, 101):\n",
    "for N in xrange(50, 1000, 50):\n",
    "    x_hat = numpy.linspace(-2 * numpy.pi, 2 * numpy.pi, N + 1)\n",
    "    delta_x.append(x_hat[1] - x_hat[0])\n",
    "\n",
    "    # Compute derivative\n",
    "    f_prime_hat = numpy.empty(x_hat.shape)\n",
    "    f_prime_hat[1:-1] = (f(x_hat[2:]) - f(x_hat[:-2])) / (2 * delta_x[-1])\n",
    "\n",
    "    # Use first-order differences for points at edge of domain\n",
    "#     f_prime_hat[0] = (f(x_hat[1]) - f(x_hat[0])) / delta_x[-1]  \n",
    "#     f_prime_hat[-1] = (f(x_hat[-1]) - f(x_hat[-2])) / delta_x[-1]\n",
    "    # Use second-order differences for points at edge of domain\n",
    "    f_prime_hat[0] =  (-3.0 * f(x_hat[0])  +  4.0 * f(x_hat[1])  + - f(x_hat[2]))  / (2.0 * delta_x[-1])\n",
    "    f_prime_hat[-1] = ( 3.0 * f(x_hat[-1]) + -4.0 * f(x_hat[-2]) +   f(x_hat[-3])) / (2.0 * delta_x[-1])\n",
    "    \n",
    "    error.append(numpy.linalg.norm(numpy.abs(f_prime(x_hat) - f_prime_hat), ord=numpy.infty))\n",
    "    \n",
    "error = numpy.array(error)\n",
    "delta_x = numpy.array(delta_x)\n",
    "    \n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "axes.loglog(delta_x, error, \"ro\", label=\"Approx. Derivative\")\n",
    "\n",
    "order_C = lambda delta_x, error, order: numpy.exp(numpy.log(error) - order * numpy.log(delta_x))\n",
    "axes.loglog(delta_x, order_C(delta_x[0], error[0], 1.0) * delta_x**1.0, 'b--', label=\"1st Order\")\n",
    "axes.loglog(delta_x, order_C(delta_x[0], error[0], 2.0) * delta_x**2.0, 'r--', label=\"2nd Order\")\n",
    "axes.legend(loc=4)\n",
    "axes.set_title(\"Convergence of 2nd Order Differences\")\n",
    "axes.set_xlabel(\"$\\Delta x$\")\n",
    "axes.set_ylabel(\"$|f'(x) - \\hat{f}'(x)|$\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "f = lambda x: numpy.sin(x) + x**2 + 3.0 * x**3\n",
    "f_dubl_prime = lambda x: -numpy.sin(x) + 2.0 + 18.0 * x\n",
    "\n",
    "# Compute the error as a function of delta_x\n",
    "delta_x = []\n",
    "error = []\n",
    "# for N in xrange(2, 101):\n",
    "for N in xrange(50, 1000, 50):\n",
    "    x_hat = numpy.linspace(-2 * numpy.pi, 2 * numpy.pi, N)\n",
    "    delta_x.append(x_hat[1] - x_hat[0])\n",
    "\n",
    "    # Compute derivative\n",
    "    f_dubl_prime_hat = numpy.empty(x_hat.shape)\n",
    "    f_dubl_prime_hat[1:-1] = (f(x_hat[2:]) -2.0 * f(x_hat[1:-1]) + f(x_hat[:-2])) / (delta_x[-1]**2)\n",
    "\n",
    "    # Use second-order differences for points at edge of domain\n",
    "    f_dubl_prime_hat[0] = (2.0 * f(x_hat[0]) - 5.0 * f(x_hat[1]) + 4.0 * f(x_hat[2]) - f(x_hat[3])) / delta_x[-1]**2\n",
    "    f_dubl_prime_hat[-1] = (2.0 * f(x_hat[-1]) - 5.0 * f(x_hat[-2]) + 4.0 * f(x_hat[-3]) - f(x_hat[-4])) / delta_x[-1]**2\n",
    "    \n",
    "    error.append(numpy.linalg.norm(numpy.abs(f_dubl_prime(x_hat) - f_dubl_prime_hat), ord=numpy.infty))\n",
    "    \n",
    "error = numpy.array(error)\n",
    "delta_x = numpy.array(delta_x)\n",
    "    \n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "# axes.plot(delta_x, error)\n",
    "axes.loglog(delta_x, error, \"ko\", label=\"Approx. Derivative\")\n",
    "order_C = lambda delta_x, error, order: numpy.exp(numpy.log(error) - order * numpy.log(delta_x))\n",
    "axes.loglog(delta_x, order_C(delta_x[2], error[2], 1.0) * delta_x**1.0, 'b--', label=\"1st Order\")\n",
    "axes.loglog(delta_x, order_C(delta_x[2], error[2], 2.0) * delta_x**2.0, 'r--', label=\"2nd Order\")\n",
    "axes.legend(loc=4)\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "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.11"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}