08_quadrature.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "<table>\n",
    " <tr align=left><td><img align=left src=\"./images/CC-BY.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",
    "from __future__ import print_function\n",
    "import numpy\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Numerical Quadrature\n",
    "\n",
    "**Goal:** Evaluate integrals\n",
    "\n",
    "$$ \\int^b_a f(x) dx$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Many integrals do not have closed form solutions\n",
    "$$ \n",
    "   \\int^b_a \\sqrt{1 + \\cos^2 x} dx\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Solution to ordinary differential equations\n",
    "   \n",
    "   $$\\frac{\\text{d}^2 u}{\\text{d}t^2} = f\\left(u, \\frac{\\text{d} u}{\\text{d}t}, t \\right)$$\n",
    "   \n",
    "   Defining $v = \\frac{\\text{d} u}{\\text{d}t}$ then leads to\n",
    "\n",
    "   $$\\begin{bmatrix}\n",
    "   \\frac{\\text{d} v}{\\text{d}t} \\\\ \\frac{\\text{d} u}{\\text{d}t} \\end{bmatrix} = \\begin{bmatrix} f(u, v, t) \\\\ v \\end{bmatrix}$$\n",
    "   \n",
    "   which can be solved by integration\n",
    "   \n",
    "   $$\\begin{bmatrix}\n",
    "   v \\\\ u \\end{bmatrix} = \\begin{bmatrix} v(t_0) + \\int^t_{t_0} f(u, v, \\hat{t}) d\\hat{t} \\\\ u(t_0) + \\int^t_{t_0} v d\\hat{t} \\end{bmatrix}$$  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Solving partial differential equations\n",
    "$$\n",
    "    u_t = \\nabla^2 u\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Basics of Quadrature\n",
    "\n",
    "We want to approximate an integral $I$ with some approximation $I_N$ such that\n",
    "$$\n",
    "    I = \\int^b_a f(x) dx \\approx I_N = \\sum^{N}_{i=1} w_i f(x_i)\n",
    "$$\n",
    "where the $x_i$ are the *quadrature points* or *nodes* and the $w_i$ are the *weights*.  Usually a particular quadrature rule specifies the points $x_i$ resulting in a particular set of weights $w_i$.\n",
    "\n",
    "Convergence requires that\n",
    "$$\n",
    "    \\lim_{N \\rightarrow \\infty} I_N = I.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Riemann Sums\n",
    "\n",
    "Given $f(x)$ and a partition of the interval $[a,b]$ with $\\{x_i\\}^N_{i=0}$ and $a = x_0 < x_1 < \\ldots < x_N = b$ and $x^*_i \\in [x_i, x_{i+1}]$ we define the Riemann integral as\n",
    "\n",
    "$$\\int^b_a f(x) dx = \\lim_{N\\rightarrow \\infty} \\sum^{N-1}_{i=0} f(x_i^*) (x_{i+1} - x_i)$$\n",
    "\n",
    "This is a general definition and leads to a number of quadrature approaches based on how we pick $x_i^* \\in [x_i, x_{i+1}]$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Midpoint Rule\n",
    "\n",
    "Choose $x_i^*$ such that\n",
    "\n",
    "$$x_i^* = \\frac{x_{i+1} + x_i}{2}$$\n",
    "\n",
    "so that\n",
    "\n",
    "$$I[f] = \\int^b_a f(x) dx \\approx \\sum^{N-1}_{i=0} f\\left(\\frac{x_{i+1} + x_i}{2} \\right ) (x_{i+1} - x_i) = I_N[f]$$\n",
    "\n",
    "over $\\Delta x_i = x_{i+1} - x_i$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example:  Integrate using midpoint rule\n",
    "\n",
    "Calculate and illustrate the midpoint rule.  Note that we are computing the cummulative integral here:\n",
    "\n",
    "$$\n",
    "    \\int^x_0 sin(\\hat{x}) d\\hat{x} = \\left . -\\cos \\hat{x} \\right|^x_0 = 1 - \\cos x\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Note that this calculates the cummulative integral from 0.0\n",
    "\n",
    "f = lambda x: numpy.sin(x)\n",
    "I = lambda x: 1.0 - numpy.cos(x)\n",
    "x = numpy.linspace(0.0, 2.0 * numpy.pi, 100)\n",
    "\n",
    "num_partitions = 10\n",
    "x_hat = numpy.linspace(0.0, 2.0 * numpy.pi, num_partitions + 1)\n",
    "x_star = 0.5 * (x_hat[1:] + x_hat[:-1])\n",
    "delta_x = x_hat[1] - x_hat[0]\n",
    "\n",
    "fig = plt.figure()\n",
    "fig.subplots_adjust(hspace=.5)\n",
    "axes = fig.add_subplot(2, 1, 1)\n",
    "\n",
    "axes.plot(x, numpy.zeros(x.shape), 'k--')\n",
    "axes.plot(x, f(x), 'b')\n",
    "\n",
    "for i in range(num_partitions):\n",
    "    axes.plot([x_hat[i], x_hat[i]], [0.0, f(x_star[i])], 'k--')\n",
    "    axes.plot([x_hat[i + 1], x_hat[i + 1]], [0.0, f(x_star[i])], 'k--')\n",
    "    axes.plot([x_hat[i], x_hat[i + 1]], [f(x_star[i]), f(x_star[i])], 'k--')\n",
    "    \n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"$f(x)$\")\n",
    "axes.set_title(\"Partition and $f(x)$\")\n",
    "axes.set_xlim((0.0, 2.0 * numpy.pi))\n",
    "axes.set_ylim((-1.1, 1.1))\n",
    "\n",
    "I_hat = numpy.zeros(x_star.shape)\n",
    "I_hat[0] = f(x_star[0]) * delta_x\n",
    "for i in range(1, num_partitions):\n",
    "    I_hat[i] = I_hat[i - 1] + f(x_star[i]) * delta_x\n",
    "    \n",
    "axes = fig.add_subplot(2, 1, 2)\n",
    "\n",
    "axes.plot(x, I(x), 'r')\n",
    "# Offset due to indexing above\n",
    "axes.plot(x_star + delta_x / 2.0, I_hat, 'ko')\n",
    "\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"$f(x)$\")\n",
    "axes.set_title(\"Integral and Approximated Integral\")\n",
    "axes.set_xlim((0.0, 2.0 * numpy.pi))\n",
    "axes.set_ylim((-0.1, 2.5))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Error Analysis\n",
    "\n",
    "From before we have a particular quadrature scheme $I_N$ which we can also write as\n",
    "$$\n",
    "    I_N[f] = \\sum^{N-1}_{i=0} w_i f(x_i).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Define the error $E[f]$ such that\n",
    "\n",
    "$$I[f] = I_N[f] + E_N[f]$$\n",
    "\n",
    "The degree of $I_N[f]$ is the integer $N$ such that $E_N[p_i] = 0 \\quad \\forall i \\leq n$ and $\\exists p_{n+1}$ such that $E[p_{n+1}] \\neq 0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Newton-Cotes Quadrature\n",
    "\n",
    "Using $N+1$ equally spaced points, evaluate $f(x)$ at these points and exactly integrate the interpolating polynomial:\n",
    "\n",
    "$$I_N[f] = \\int^b_a P_N(x) dx$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Trapezoidal Rule\n",
    "\n",
    "Use $N = 1$ polynomial to derive the trapezoidal rule."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Trapezoidal rule uses $N = 1$ order polynomials between each point (i.e. piece-wise defined linear polynomials).  The coefficients of the polynomial in each interval are\n",
    "\n",
    "$$p_0 = f(x_i) \\quad \\quad p_1 = \\frac{f(x_{i+1}) - f(x_i)}{x_{i+1} - x_i}$$\n",
    "\n",
    "which gives the interpolating polynomial\n",
    "\n",
    "$$p_1(x) = \\frac{f(x_{i+1}) - f(x_i)}{x_{i+1} - x_i} ( x- x_i) + f(x_i)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Integrating this polynomial we have\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    I_N[f] &= \\int^{x_{i+1}}_{x_i} (p_0 + p_1 (x - x_i)) dx = \\left . p_0 x + p_1 \\left (\\frac{x^2}{2} - x_i x\\right) \\right |^{x_{i+1}}_{x_i} \\\\\n",
    "    &= p_0 \\Delta x + p_1 \\left (\\frac{1}{2} (x_{i+1} + x_i) \\Delta x - x_i \\Delta x\\right) \\\\\n",
    "    &= f(x_i) \\Delta x + (f(x_{i+1}) - f(x_i))\\left (\\frac{1}{2} (x_{i+1} + x_i) - x_i\\right) \\\\\n",
    "    &= f(x_i) \\Delta x + (f(x_{i+1}) - f(x_i)) \\frac{\\Delta x}{2} \\\\\n",
    "    & = \\frac{\\Delta x}{2} (f(x_{i+1}) + f(x_i))\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can also simplify the sum over all the intervals by noting that all but the end points will have total contribution of $\\Delta x$ to the entire sum such that\n",
    "\n",
    "$$\n",
    "    I_N[f] = \\frac{\\Delta x}{2} (f(x_0) + f(x_N) ) + \\sum^{N-1}_{j=1} \\Delta x f(x_j)\n",
    "$$\n",
    "\n",
    "This is known as the composite trapezoidal rule."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Note that this calculates the cummulative integral from 0.0\n",
    "\n",
    "f = lambda x: numpy.sin(x)\n",
    "I = lambda x: 1.0 - numpy.cos(x)\n",
    "x = numpy.linspace(0.0, 2.0 * numpy.pi, 100)\n",
    "\n",
    "num_partitions = 20\n",
    "x_hat = numpy.linspace(0.0, 2.0 * numpy.pi, num_partitions + 1)\n",
    "delta_x = x_hat[1] - x_hat[0]\n",
    "\n",
    "fig = plt.figure()\n",
    "fig.subplots_adjust(hspace=.5)\n",
    "axes = fig.add_subplot(2, 1, 1)\n",
    "\n",
    "axes.plot(x, numpy.zeros(x.shape), 'k--')\n",
    "axes.plot(x, f(x), 'b')\n",
    "\n",
    "for i in range(num_partitions):\n",
    "    axes.plot([x_hat[i], x_hat[i]], [0.0, f(x_hat[i])], 'k--')\n",
    "    axes.plot([x_hat[i + 1], x_hat[i + 1]], [0.0, f(x_hat[i+1])], 'k--')\n",
    "    axes.plot([x_hat[i], x_hat[i + 1]], [f(x_hat[i]), f(x_hat[i+1])], 'k--')\n",
    "    \n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"$f(x)$\")\n",
    "axes.set_title(\"Partition and $f(x)$\")\n",
    "axes.set_xlim((0.0, 2.0 * numpy.pi))\n",
    "axes.set_ylim((-1.1, 1.1))\n",
    "\n",
    "I_hat = numpy.zeros(x_hat.shape)\n",
    "I_hat[0] = (f(x_hat[1]) + f(x_hat[0])) * delta_x / 2.0\n",
    "for i in range(1, num_partitions):\n",
    "    I_hat[i] = I_hat[i - 1] + (f(x_hat[i + 1]) + f(x_hat[i])) * delta_x / 2.0\n",
    "    \n",
    "axes = fig.add_subplot(2, 1, 2)\n",
    "\n",
    "axes.plot(x, I(x), 'r')\n",
    "# Offset due to indexing above\n",
    "axes.plot(x_hat + delta_x, I_hat, 'ko')\n",
    "\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"$f(x)$\")\n",
    "axes.set_title(\"Integral and Approximated Integral\")\n",
    "axes.set_xlim((0.0, 2.0 * numpy.pi))\n",
    "axes.set_ylim((-0.1, 2.5))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Simpson's Rule\n",
    "\n",
    "Simpson's rule uses $N = 2$ order polynomials between each point (i.e. piece-wise defined quadratic polynomials).  \n",
    "\n",
    "The polynomial has the form\n",
    "\n",
    "$$P_2(x) = \\frac{2 f(x_i)}{\\Delta x^2} \\left (x - \\frac{\\Delta x}{2} \\right ) (x - \\Delta x) - \\frac{4 f\\left(x_i + \\frac{\\Delta x}{2}\\right)}{\\Delta x^2}  x (x - \\Delta x) + \\frac{2 f(x_{i+1})}{\\Delta x^2} x \\left (x - \\frac{\\Delta x}{2} \\right )$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Integrating this polynomial we have\n",
    "\n",
    "$$\n",
    "    I_N[f] = \\int^{x_{i+1}}_{x_i} P_2(x) dx = \\frac{\\Delta x}{6} f(x_i) + \\frac{2 \\Delta x}{3} f\\left(x_i + \\frac{\\Delta x}{2} \\right ) + \\frac{\\Delta x}{6} f(x_{i+1})\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can also show this by using the method of undetermined coefficients.  \n",
    "\n",
    "Use the general form of the quadrature rule and determine weights $w_j$ by using functions we know the solution to.  These functions can be any representation of polynomials up to order $N=2$ however the monomials $1$, $x$, $x^2$ are the easiest in this case.\n",
    "\n",
    "$$\n",
    "    I_{\\Delta x}[f] = w_0 f(0) + w_1 f(\\Delta x / 2) + w_2 f(\\Delta x)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$$\\begin{aligned}\n",
    "    &\\text{if}~f = 1:  &I[f] =  \\int^{\\Delta x}_{0} 1 dx = \\Delta x & & I_N[1] &= w_0 + w_1 + w_2 \\\\\n",
    "    &\\text{if}~f = x:  &I[f] =  \\int^{\\Delta x}_{0} x dx = \\frac{\\Delta x^2}{2} & & I_N[x] &= w_1 \\frac{\\Delta x}{2} + w_2\\Delta x\\\\\n",
    "    &\\text{if}~f = x^2:  &I[f] =  \\int^{\\Delta x}_{0} x^2 dx = \\frac{\\Delta x^3}{3} & & I_N[x^2] &= \\frac{\\Delta x^2}{4} w_1 + w_2\\Delta x^2\\\\\n",
    "\\end{aligned}$$\n",
    "\n",
    "We then have the system of equations:\n",
    "$$\\begin{aligned}\n",
    "    w_0 &+& w_1 &+& w_2 &=\\Delta x \\\\\n",
    "        &\\quad& \\frac{\\Delta x}{2} w_1 &+& \\Delta x w_2  &= \\frac{\\Delta x^2}{2} \\\\\n",
    "        &\\quad& \\frac{\\Delta x^2}{4} w_1 &+& \\Delta x^2 w_2 &=\\frac{\\Delta x^3}{6} \\\\\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "or\n",
    "\n",
    "$$\\begin{bmatrix}\n",
    "    1 & 1 & 1 \\\\\n",
    "    0 & \\Delta x / 2 & \\Delta x \\\\\n",
    "    0 & \\Delta x^2 / 4 & \\Delta x^2 \\\\\n",
    "\\end{bmatrix} \\begin{bmatrix}\n",
    "    w_0 \\\\ w_1 \\\\ w_2\n",
    "\\end{bmatrix} = \\begin{bmatrix} \n",
    "    \\Delta x \\\\ \\Delta x^2 / 2 \\\\ \\Delta x^3 / 3\n",
    "\\end{bmatrix} \\Rightarrow \\begin{bmatrix}\n",
    "    1 & 1 & 1 \\\\\n",
    "    0 & 1 / 2 & 1 \\\\\n",
    "    0 & 1 / 4 & 1 \\\\\n",
    "\\end{bmatrix} \\begin{bmatrix}\n",
    "    w_0 \\\\ w_1 \\\\ w_2\n",
    "\\end{bmatrix} = \\begin{bmatrix} \n",
    "    \\Delta x \\\\ \\Delta x / 2 \\\\ \\Delta x / 3\n",
    "\\end{bmatrix} \\Rightarrow \\begin{bmatrix}\n",
    "    1 & 1 & 1 \\\\\n",
    "    0 & 1 / 2 & 1 \\\\\n",
    "    0 & 0 & -1 \\\\\n",
    "\\end{bmatrix} \\begin{bmatrix}\n",
    "    w_0 \\\\ w_1 \\\\ w_2\n",
    "\\end{bmatrix} = \\begin{bmatrix} \n",
    "    \\Delta x \\\\ \\Delta x / 2 \\\\ -\\Delta x / 6\n",
    "\\end{bmatrix}$$\n",
    "\n",
    "Leading to \n",
    "\n",
    "$$ w_2 = \\frac{\\Delta x}{6} \\quad w_1 = \\frac{2}{3} \\Delta x \\quad w_0 = \\frac{\\Delta x}{6}$$\n",
    "\n",
    "Another way to write Simpson's rule is to use intervals of three points (similar to one of the ways we did this last time).  The formulation here effectively has a $\\Delta x$ half of what the intervals show but is easier to program."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Note that this calculates the cummulative integral from 0.0\n",
    "\n",
    "f = lambda x: numpy.sin(x)\n",
    "I = lambda x: 1.0 - numpy.cos(x)\n",
    "x = numpy.linspace(0.0, 2.0 * numpy.pi, 100)\n",
    "\n",
    "num_partitions = 10\n",
    "x_hat = numpy.linspace(0.0, 2.0 * numpy.pi, num_partitions + 1)\n",
    "delta_x = x_hat[1] - x_hat[0]\n",
    "\n",
    "fig = plt.figure()\n",
    "fig.subplots_adjust(hspace=.5)\n",
    "axes = fig.add_subplot(2, 1, 1)\n",
    "\n",
    "axes.plot(x, numpy.zeros(x.shape), 'k--')\n",
    "axes.plot(x, f(x), 'b')\n",
    "\n",
    "for i in range(num_partitions):\n",
    "    axes.plot([x_hat[i], x_hat[i]], [0.0, f(x_hat[i])], 'k--')\n",
    "    axes.plot([x_hat[i + 1], x_hat[i + 1]], [0.0, f(x_hat[i + 1])], 'k--')\n",
    "    coeff = numpy.polyfit((x_hat[i], x_hat[i] + delta_x / 2.0, x_hat[i + 1]), \n",
    "                          (f(x_hat[i]), f(x_hat[i] + delta_x / 2.0), f(x_hat[i+1])), 2)\n",
    "    x_star = numpy.linspace(x_hat[i], x_hat[i+1], 10)\n",
    "    axes.plot(x_star, numpy.polyval(coeff, x_star), 'k--')\n",
    "    \n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"$f(x)$\")\n",
    "axes.set_title(\"Partition and $f(x)$\")\n",
    "axes.set_xlim((0.0, 2.0 * numpy.pi))\n",
    "axes.set_ylim((-1.1, 1.1))\n",
    "\n",
    "I_hat = numpy.zeros(x_hat.shape)\n",
    "I_hat[0] = delta_x * (1.0 / 6.0 * (f(x_hat[0]) + f(x_hat[1])) + 2.0 / 3.0 * f(x_hat[0] + delta_x / 2.0))\n",
    "for i in range(1, num_partitions):\n",
    "    I_hat[i] = I_hat[i - 1] + delta_x * (1.0 / 6.0 * (f(x_hat[i]) + f(x_hat[i+1])) + 2.0 / 3.0 * f(x_hat[i] + delta_x / 2.0))\n",
    "    \n",
    "axes = fig.add_subplot(2, 1, 2)\n",
    "\n",
    "axes.plot(x, I(x), 'r')\n",
    "# Offset due to indexing above\n",
    "axes.plot(x_hat + delta_x, I_hat, 'ko')\n",
    "\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"$f(x)$\")\n",
    "axes.set_title(\"Integral and Approximated Integral\")\n",
    "axes.set_xlim((0.0, 2.0 * numpy.pi))\n",
    "axes.set_ylim((-0.1, 2.5))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Quadrature Accuracy\n",
    "\n",
    "We can also use our polynomial analysis from before to analyze the errors made using both of the aforementioned methods.  From Lagrange's theorem we have the remainder term as before which we can use to look at the error\n",
    "\n",
    "$$R_N(x) = (x - x_0)(x - x_1) \\cdots (x- x_N) \\frac{f^{(N+1)}(c)}{(N+1)!}$$\n",
    "\n",
    "and integrate it to find the form and magnitude of the error on a single interval.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "To find the total error we must sum the error over all the intervals:\n",
    "\n",
    "$$I[f] = \\sum_{i=0}^N \\int^{x_{i+1}}_{x_i} P_N(x) dx + \\sum_{i=0}^N \\int^{x_{i+1}}_{x_i} R_N(x) dx = I_N[f] + E_N[f]$$\n",
    "\n",
    "as we defined before."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Trapezoidal error\n",
    "\n",
    "With $N=1$ we have\n",
    "\n",
    "$$R_1(x) = (x - x_i) (x - x_{i+1}) \\frac{f''(c)}{2}$$\n",
    "\n",
    "Integrating this leads to\n",
    "\n",
    "$$\\int^{x_{i+1}}_{x_i} (x - x_i) (x - x_{i+1}) \\frac{f''(c)}{2} dx = \\frac{\\Delta x^3}{12} f''(c)$$\n",
    "\n",
    "giving us a form for the error.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "If we sum up across all the intervals the total error is\n",
    "\n",
    "$$E_N[f] = -\\frac{\\Delta x^3}{12} \\sum_{i=0}^{N} f''(c_i)$$\n",
    "\n",
    "or more illustrative\n",
    "\n",
    "$$E_N[f] = -\\frac{1}{2} \\Delta x^2 (b - a) \\left [ \\frac{1}{N} \\sum^{N-1}_{i=0} f''(c_i) \\right ]$$\n",
    "\n",
    "where the expression in the brackets is the mean value of the second derivative over the interval $[a,b]$.  This also shows that the trapezoidal rule converges quadratically as $\\Delta x \\rightarrow 0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Simpson's Rule Error\n",
    "\n",
    "Similarly here we have $N = 2$ and \n",
    "\n",
    "$$R_2(x) = (x - x_i) \\left(x - x_i - \\frac{\\Delta x}{2} \\right) (x - x_{i+1}) \\frac{f'''(c)}{3!}$$\n",
    "\n",
    "Integrating and summing the error contributions we find\n",
    "\n",
    "$$E_N[f] = -\\frac{1}{180} (b - a) \\Delta x^4 f^{(4)}(c)$$\n",
    "\n",
    "Interestingly we have gained two orders of accuracy by increasing the polynomial order by only 1!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "##### Example 1:\n",
    "\n",
    "If $f(x) = \\sin \\pi x$ look at the relative accuracy of midpoint, trapezoidal and simpson's rules for a single interval $x\\in[0,1]$.\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    \\text{Exact:}  &I[f] = \\int^1_0 \\sin \\pi x = \\left . \\frac{-\\cos \\pi x}{\\pi} \\right |^1_0 = \\frac{2}{\\pi} \\approx 0.636619772 \\\\\n",
    "    \\text{Midpoint:}  &I_1[f] = \\Delta x f(1/2) = \\sin (\\pi / 2) = 1 \\\\\n",
    "    \\text{Trapezoid:}  &I_1[f] = \\frac{\\Delta x}{2} (\\sin(0) + \\sin(\\pi)) = 0 \\\\\n",
    "    \\text{Simpson's:}  &I_1[f] = \\frac{\\Delta x}{6} \\sin(0) + \\frac{2 \\Delta x}{3} \\sin(\\pi / 2) + \\frac{\\Delta x}{6} \\sin(\\pi) = \\frac{2 \\Delta x}{3} = \\frac{2}{3}\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Compute the error as a function of delta_x for each method\n",
    "f = lambda x: numpy.sin(numpy.pi * x)\n",
    "\n",
    "num_partitions = range(50, 1000, 50)\n",
    "delta_x = numpy.empty(len(num_partitions))\n",
    "error_mid = numpy.empty(len(num_partitions))\n",
    "error_trap = numpy.empty(len(num_partitions))\n",
    "error_simpson = numpy.empty(len(num_partitions))\n",
    "\n",
    "for (j, N) in enumerate(num_partitions):\n",
    "    x_hat = numpy.linspace(0.0, 1.0, N + 1)\n",
    "    delta_x[j] = x_hat[1] - x_hat[0]\n",
    "\n",
    "    # Compute Midpoint\n",
    "    x_star = 0.5 * (x_hat[1:] + x_hat[:-1])\n",
    "    I_hat = 0.0\n",
    "    for i in range(0, N):\n",
    "        I_hat += f(x_star[i]) * delta_x[j]\n",
    "    error_mid[j] = numpy.abs(I_hat - 2.0 / numpy.pi)\n",
    "    \n",
    "    # Compute trapezoid\n",
    "    I_hat = 0.0\n",
    "    for i in range(1, N):\n",
    "        I_hat += (f(x_hat[i + 1]) + f(x_hat[i])) * delta_x[j] / 2.0\n",
    "    error_trap[j] = numpy.abs(I_hat - 2.0 / numpy.pi)\n",
    "    \n",
    "    # Compute simpson's    \n",
    "    I_hat = 0.0\n",
    "    for i in range(0, N):\n",
    "        I_hat += delta_x[j] * (1.0 / 6.0 * (f(x_hat[i]) + f(x_hat[i+1])) + 2.0 / 3.0 * f(x_hat[i] + delta_x[j] / 2.0))\n",
    "    error_simpson[j] = numpy.abs(I_hat - 2.0 / numpy.pi)\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "order_C = lambda delta_x, error, order: numpy.exp(numpy.log(error) - order * numpy.log(delta_x))\n",
    "axes.loglog(delta_x, error_mid, 'ro', label=\"Midpoint\")\n",
    "axes.loglog(delta_x, error_trap, 'bo', label=\"Trapezoid\")\n",
    "axes.loglog(delta_x, error_simpson, 'go', label=\"Simpson's\")\n",
    "axes.loglog(delta_x, order_C(delta_x[0], error_trap[0], 2.0) * delta_x**2.0, 'b--', label=\"2nd Order\")\n",
    "axes.loglog(delta_x, order_C(delta_x[0], error_simpson[0], 4.0) * delta_x**4.0, 'g--', label=\"4th Order\")\n",
    "axes.legend(loc=4)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Recursive Improvement of Accuracy\n",
    "\n",
    "Say we ran the trapezoidal rule with step size $2 \\Delta x$, we then will have\n",
    "\n",
    "$$\\begin{aligned}\n",
    "\t\t\t\\int^{x_2}_{x_0} f(x) dx  &= \\frac{2 \\Delta x}{2} (f_0 + f_2) =  h (f_0 + f_2) \\Rightarrow \\\\\n",
    "\t\t\t\\int^b_a f(x)dx &\\approx I_{2\\Delta x}[f] = \\sum^{N/2-1}_{j=0} \\Delta x (f_{2j} + f_{2j+2}) \\\\\n",
    "\t\t\t&= \\Delta x (f_{0} + f_{2})  + \\Delta x (f_{2} + f_{4})  + \\cdots + \\Delta x (f_{N-2} + f_{N}) \\\\\n",
    "\t\t\t&= \\Delta x\\left ( f_0 + f_N +  2 \\sum^{N/2-1}_{j=1} f_{2j} \\right )\n",
    "  \\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Now compare the two rules for $\\Delta x$ and $2 \\Delta x$:\n",
    "\n",
    "$$I_{\\Delta x}[f] = \\frac{\\Delta x}{2} \\left (f_0 + f_N + 2 \\sum^{N-1}_{j=1} f_j \\right)~~~~~~~~~ I_{2 \\Delta x}[f] = \\Delta x \\left ( f_0 + f_N +  2 \\sum^{N/2-1}_{j=1} f_{2j} \\right )$$\n",
    "\n",
    "$$I_{\\Delta x}[f] = \\frac{1}{2} I_{2\\Delta x} + \\Delta x(f_1 + f_3 + \\cdots + f_{N-1})$$\n",
    "\n",
    "Here we see we can actually reuse the work we did to calculate $Q_{2 \\Delta x}[f]$ to refine the integral."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Arbitrary Intervals (Affine Transforms)\n",
    "\n",
    "Mapping $\\xi \\in [-1,1] \\rightarrow x \\in [a,b]$ can be done through an *affine transform* or *affine map* which is a linear transformation.\n",
    "\n",
    "$$x = \\underbrace{\\frac{b - a}{2}}_{\\text{scaling}} \\xi + \\underbrace{\\frac{a+b}{2}}_{\\text{translation}} ~~~~~ \\text{or} ~~~~~ \\xi = \\left( x - \\frac{a + b}{2}\\right) \\frac{2}{b-a}$$\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    I[f] &= \\int^b_a f(x) dx = \\int^1_{-1} f(x(\\xi)) \\frac{dx}{d\\xi} d\\xi = \\frac{b - a}{2} \\int^1_{-1} f(x(\\xi)) d\\xi\\\\\n",
    "    I_N[f] &= \\sum_i w_i f(x(\\xi_i)) \\left . \\frac{dx}{d\\xi}\\right|_{\\xi_i}\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Optimal Quadrature Methods\n",
    "\n",
    "Can we determine $I_{\\Delta x}[f]$ to maximize degree for a given number of function evaluations (points)?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Generalized Gaussian Quadrature\n",
    "\n",
    "Given $g(x) \\in P_N(x)$ with roots $\\{x_i\\}^N_{i=1}$ we have\n",
    "$$\n",
    "    \\int^1_{-1} w(x) x^i g(x) dx = 0 \\quad \\forall i < N,\n",
    "$$\n",
    "i.e. $g(x)$ is orthogonal to the $x^i$ with respect to the weight function $w(x)$.\n",
    "\n",
    "Recall something similar:\n",
    "$$\n",
    "    <x, y> = \\sum^N_{i=1} x_i y_i = ||x|| \\cdot ||y|| \\cos \\theta.\n",
    "$$\n",
    "If $<x, y> = 0$ then the vectors $x$ and $y$ are orthogonal."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Given the above $g(x)$ there then exists $\\{w_i\\}$ such that\n",
    "\n",
    "$$\\int^1_{-1} w(x) P_j(x) dx = \\sum^N_{i=1} w_i P_j(x_i) \\quad \\forall j \\leq 2 N - 1$$\n",
    "\n",
    "In other words, given a polynomial basis function and weight and orthogonality to all polynomials of order $i < N$ we can exactly integrate polynomials of order $2 N - 1$.  Choosing the correct weighting function and basis leads to a number of useful quadrature approaches:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Gauss-Legendre\n",
    "\n",
    "General Gauss-Legendre quadrature uses $w(x) = 1$ and $g(x) = \\ell_N(x)$ which can be shown to have weights\n",
    "\n",
    "$$w_i = \\frac{2}{(1-x_i^2)(P'_n(x_i))^2}$$\n",
    "\n",
    "and $x_i$ is the $i$th root of $\\ell_N$.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The first few rules yield\n",
    "\n",
    "<table width=\"80%\">\n",
    "    <tr align=\"center\"><th>$$N$$</th> <th align=\"center\">$$x_i$$</th> <th align=\"center\"> $$w_i$$ </th></tr>\n",
    "    <tr align=\"center\"><td>$$1$$</td>           <td> $$0$$ </td> <td> $$2$$ </td> </tr>\n",
    "    <tr align=\"center\"><td>$$2$$</td>           <td> $$\\pm \\sqrt{\\frac{1}{3}}$$ </td> <td> $$1$$ </td> </tr>\n",
    "    <tr align=\"center\"><td rowspan=2>$$3$$</td> <td> $$0$$ </td> <td> $$8/9$$ </td> </tr>\n",
    "    <tr align=\"center\">                     <td> $$\\pm \\sqrt{\\frac{3}{5}}$$ </td> <td> $$5/9$$</td> </tr>\n",
    "    <tr align=\"center\"><td rowspan=2>$$4$$</td> <td> $$\\pm \\sqrt{\\frac{3}{7} - \\frac{2}{7} \\sqrt{\\frac{6}{5}}}$$</td> <td> $$\\frac{18 + \\sqrt{30}}{36}$$ </td> </tr>\n",
    "    <tr align=\"center\">                     <td> $$\\pm \\sqrt{\\frac{3}{7} + \\frac{2}{7} \\sqrt{\\frac{6}{5}}}$$</td> <td>$$\\frac{18 - \\sqrt{30}}{36}$$ </td> </tr>\n",
    "    <tr align=\"center\"><td rowspan=3>$$5$$</td> <td> $$0$$ </td> <td> $$\\frac{128}{225}$$ </td> </tr>\n",
    "    <tr align=\"center\">                     <td> $$\\pm \\frac{1}{3} \\sqrt{5 - 2 \\sqrt{\\frac{10}{7}}}$$</td> <td> $$\\frac{322 + 13\\sqrt{70}}{900}$$</td> </tr>\n",
    "    <tr align=\"center\">                     <td> $$\\pm \\frac{1}{3} \\sqrt{5 + 2 \\sqrt{\\frac{10}{7}}}$$</td> <td> $$\\frac{322 - 13\\sqrt{70}}{900}$$</td> </tr>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "##### Example 2:  2-Point Gauss-Legendre Quadrature\n",
    "\n",
    "Let $N=2$ on $x \\in [-1,1]$\n",
    "\n",
    "$$I_2[f] = w_0 f(x_0) + w_1 f(x_1)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Using undetermined coefficients again we have\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    &\\text{if}~f = 1:  &I[f] =  \\int^{1}_{-1} 1 dx = 2 & & I_2[1] &= w_0 + w_1\\\\\n",
    "    &\\text{if}~f = x:  &I[f] =  \\int^{1}_{-1} x dx = 0 & & I_2[x] &= w_0 x_0 + w_1 x_1\\\\\n",
    "    &\\text{if}~f = x^2:  &I[f] =  \\int^{1}_{-1} x^2 dx = \\frac{2}{3} & & I_2[x^2] &= w_0 x_0^2 + w_1 x_1^2\\\\\n",
    "    &\\text{if}~f = x^3:  &I[f] =  \\int^{1}_{-1} x^3 dx = 0 & & I_2[x^3] &= w_0 x_0^3 + w_1 x_1^3\\\\\n",
    "\\end{aligned}$$\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    &w_0 + w_1 = 2\\\\\n",
    "    &w_0 x_0 + w_1 x_1 = 0\\\\\n",
    "    &w_0 x_0^2 + w_1 x_1^2 = \\frac{2}{3}\\\\\n",
    "    &w_0 x_0^3 + w_1 x_1^3 = 0\\\\\n",
    "\\end{aligned}$$\n",
    "\n",
    "Note that we need to solve for 4 unknowns $x_0$, $x_1$, $w_0$, and $w_1$.  Solving these equations leads to\n",
    "\n",
    "$$x_0 = -\\sqrt{\\frac{1}{3}}, x_1 = \\sqrt{\\frac{1}{3}}  \\quad \\text{and} \\quad w_0 = w_1 = 1 $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Compute Gauss-Legendre based quadrature and affine transforms\n",
    "f = lambda x: numpy.sin(x)\n",
    "I = lambda x: 1.0 - numpy.cos(x)\n",
    "x = numpy.linspace(0.0, 2.0 * numpy.pi, 100)\n",
    "\n",
    "num_partitions = 10\n",
    "x_hat = numpy.linspace(0.0, 2.0 * numpy.pi, num_partitions + 1)\n",
    "delta_x = x_hat[1] - x_hat[0]\n",
    "\n",
    "xi_map = lambda a,b,xi : (b - a) / 2.0 * xi + (a + b) / 2.0\n",
    "xi_0 = -numpy.sqrt(1.0 / 3.0)\n",
    "xi_1 =  numpy.sqrt(1.0 / 3.0)\n",
    "I_hat = numpy.zeros(x_hat.shape)\n",
    "I_hat[0] = (f(xi_map(x_hat[0], x_hat[1], xi_0)) + f(xi_map(x_hat[0], x_hat[1], xi_1))) * delta_x / 2.0\n",
    "for i in range(1, num_partitions):\n",
    "    I_hat[i] = I_hat[i - 1] + (f(xi_map(x_hat[i], x_hat[i+1], xi_0)) + f(xi_map(x_hat[i], x_hat[i+1], xi_1))) * delta_x / 2.0\n",
    "    \n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "axes.plot(x, I(x), 'r')\n",
    "# Offset due to indexing above\n",
    "axes.plot(x_hat + delta_x, I_hat, 'ko')\n",
    "\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"$f(x)$\")\n",
    "axes.set_title(\"Integral and Approximated Integral\")\n",
    "axes.set_xlim((0.0, 2.0 * numpy.pi))\n",
    "axes.set_ylim((-0.1, 2.5))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Other Quadrature Families\n",
    "\n",
    " - Clenshaw-Curtis (Gauss-Chebyshev): If $w(x) = \\frac{1}{\\sqrt{1 - x^2}}$ and $g(x)$ are Chebyshev polynomials then we know the roots of the polynomials to be $x_i = \\cos\\left(\\frac{2i-1}{2N} \\pi \\right)$ (the Chebyshev nodes) and we can derive that $w_i = \\frac{\\pi}{N}$.\n",
    " - Gauss-Hermite:  If $w(x) = e^{-x^2}$ and $g(x)$ are Hermite polynomials $H_i(x)$ then\n",
    "   $$w_i = \\frac{2^{N-1} N! \\sqrt{\\pi}}{N^2 (H_{N-1}(x_i))^2}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "##### Example 3:\n",
    "\n",
    "If $f(x) = e^x$ look at the relative accuracy of midpoint, trapezoidal, Simpson's and 2-point Gauss-Legendre quadrature for a single interval $x \\in [-1,1]$.\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    \\text{Exact:} &I[f] = \\int^1_{-1} e^x = \\left . e^x \\right |^1_{-1} = e - \\frac{1}{e} \\approx 2.350402387 \\\\\n",
    "    \\text{Midpoint:} &I_2[f] = 2 e^0 = 2 \\\\\n",
    "    \\text{Trapezoid:} &I_2[f] = \\frac{2}{2} (e^{-1} + e^1) = e + \\frac{1}{e} = 3.08616127 \\\\\n",
    "    \\text{Simpson's:} &I_2[f] = \\frac{2}{6} e^{-1} + \\frac{4}{3} e^0 + \\frac{2}{6} e^1 = \\frac{4}{3} + \\frac{1}{3} (e^{-1} + e^1) \\approx 2.362053757 \\\\\n",
    "    \\text{Gauss-Legendre:} &I_2[f] = e^{-\\sqrt{\\frac{1}{3}}} + e^{\\sqrt{\\frac{1}{3}}} \\approx 2.342696088\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Compute the error as a function of delta_x for each method\n",
    "f = lambda x: numpy.sin(numpy.pi * x)\n",
    "I = 2.0 / numpy.pi\n",
    "\n",
    "# num_partitions = range(50, 1000, 50)\n",
    "num_partitions = range(5, 50, 5)\n",
    "delta_x = numpy.empty(len(num_partitions))\n",
    "error_trap = numpy.empty(len(num_partitions))\n",
    "error_simpson = numpy.empty(len(num_partitions))\n",
    "error_2 = numpy.empty(len(num_partitions))\n",
    "error_3 = numpy.empty(len(num_partitions))\n",
    "error_4 = numpy.empty(len(num_partitions))\n",
    "\n",
    "for (j, N) in enumerate(num_partitions):\n",
    "    x_hat = numpy.linspace(0.0, 1.0, N)\n",
    "    delta_x[j] = x_hat[1] - x_hat[0]\n",
    "    \n",
    "    # Compute trapezoid\n",
    "    I_hat = 0.0\n",
    "    for i in range(0, N - 1):\n",
    "        I_hat += (f(x_hat[i + 1]) + f(x_hat[i])) * delta_x[j] / 2.0\n",
    "    error_trap[j] = numpy.abs(I_hat - I)\n",
    "    \n",
    "    # Compute simpson's    \n",
    "    I_hat = 0.0\n",
    "    for i in range(0, N - 1):\n",
    "        I_hat += delta_x[j] * (1.0 / 6.0 * (f(x_hat[i]) + f(x_hat[i+1])) + 2.0 / 3.0 * f(x_hat[i] + delta_x[j] / 2.0))\n",
    "    error_simpson[j] = numpy.abs(I_hat - I)\n",
    "    \n",
    "    # Compute Gauss-Legendre 2-point\n",
    "    xi_map = lambda a,b,xi : (b - a) / 2.0 * xi + (a + b) / 2.0\n",
    "    xi = [-numpy.sqrt(1.0 / 3.0), numpy.sqrt(1.0 / 3.0)]\n",
    "    w = [1.0, 1.0]\n",
    "    I_hat = 0.0\n",
    "    for i in range(0, N - 1):\n",
    "        for k in range(len(xi)):\n",
    "            I_hat += f(xi_map(x_hat[i], x_hat[i+1], xi[k])) * w[k]\n",
    "    I_hat *= delta_x[j] / 2.0\n",
    "    error_2[j] = numpy.abs(I_hat - I)\n",
    "    \n",
    "    # Compute Gauss-Legendre 3-point\n",
    "    xi_map = lambda a,b,xi : (b - a) / 2.0 * xi + (a + b) / 2.0\n",
    "    xi = [-numpy.sqrt(3.0 / 5.0), 0.0, numpy.sqrt(3.0 / 5.0)]\n",
    "    w = [5.0 / 9.0, 8.0 / 9.0, 5.0 / 9.0]\n",
    "    I_hat = 0.0\n",
    "    for i in range(0, N - 1):\n",
    "        for k in range(len(xi)):\n",
    "            I_hat += f(xi_map(x_hat[i], x_hat[i+1], xi[k])) * w[k]\n",
    "    I_hat *= delta_x[j] / 2.0\n",
    "    error_3[j] = numpy.abs(I_hat - I)\n",
    "    \n",
    "    # Compute Gauss-Legendre 4-point\n",
    "    xi_map = lambda a,b,xi : (b - a) / 2.0 * xi + (a + b) / 2.0\n",
    "    xi = [-numpy.sqrt(3.0 / 7.0 - 2.0 / 7.0 * numpy.sqrt(6.0 / 5.0)), \n",
    "           numpy.sqrt(3.0 / 7.0 - 2.0 / 7.0 * numpy.sqrt(6.0 / 5.0)),\n",
    "          -numpy.sqrt(3.0 / 7.0 + 2.0 / 7.0 * numpy.sqrt(6.0 / 5.0)),\n",
    "           numpy.sqrt(3.0 / 7.0 + 2.0 / 7.0 * numpy.sqrt(6.0 / 5.0))]\n",
    "    w = [(18.0 + numpy.sqrt(30.0)) / 36.0, (18.0 + numpy.sqrt(30.0)) / 36.0,\n",
    "         (18.0 - numpy.sqrt(30.0)) / 36.0, (18.0 - numpy.sqrt(30.0)) / 36.0]\n",
    "    I_hat = 0.0\n",
    "    for i in range(0, N - 1):\n",
    "        for k in range(len(xi)):\n",
    "            I_hat += f(xi_map(x_hat[i], x_hat[i+1], xi[k])) * w[k]\n",
    "    I_hat *= delta_x[j] / 2.0\n",
    "    error_4[j] = numpy.abs(I_hat - I)\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_trap, 'o', label=\"Trapezoid\")\n",
    "axes.loglog(delta_x, error_simpson, 'o', label=\"Simpson's\")\n",
    "axes.loglog(delta_x, error_2, 'o', label=\"G-L 2-point\")\n",
    "axes.loglog(delta_x, error_3, 'o', label=\"G-L 3-point\")\n",
    "axes.loglog(delta_x, error_4, 'o', label=\"G-L 4-point\")\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_trap[0], 2.0) * delta_x**2.0, 'r--', label=\"2nd Order\")\n",
    "axes.loglog(delta_x, order_C(delta_x[0], error_simpson[0], 4.0) * delta_x**4.0, 'g--', label=\"4th Order\")\n",
    "axes.loglog(delta_x, order_C(delta_x[1], error_3[1], 5) * delta_x**5, 'b--', label=\"5th Order\")\n",
    "axes.loglog(delta_x, order_C(delta_x[1], error_4[1], 7.0) * delta_x**7.0, 'g--', label=\"7th Order\")\n",
    "\n",
    "axes.legend(loc=2)\n",
    "axes.set_xlim((5e-3, delta_x[0]))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## SciPy Integration Routines\n",
    "\n",
    "SciPy has a number of integration routines that we have derived here including general purpose integrators that can control error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "import scipy.integrate as integrate\n",
    "# integrate?"
   ]
  }
 ],
 "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.6.3"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}