{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numerical Integration with Simpson's Rule\n",
    "\n",
    "## Formula\n",
    "\n",
    "Unlike the Riemann sum methods, [Stewart](https://www.amazon.ca/Calculus-Early-Transcendentals-James-Stewart/dp/1285741552) does not give a sigma formula for Simpson's Rule. Instead, the textbook defines it as a series:\n",
    "\n",
    "$$ \\int_{a}^b f(x) dx \\approx S_n = \\frac{\\Delta x}{3} [f(x_0) + 4f(x_1) + 2f(x_2) + 4f(x_3) + ... + 2f(x_{n-2}) + 4f(x_{n-1}) + f(x_n)] $$\n",
    "\n",
    "Where $n$ is even and $\\Delta x = \\frac{b - a}{n} $."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The authors say \"note the pattern of the coefficients\" (pg. 520), which feels like a cop-out. What does $x_0$ mean? How do I find the actual values I'm supposed to plug into the function being integrated? Attempting to convert the formula into sigma form made these blindspots obvious.\n",
    "\n",
    "$$ \\frac{\\Delta x}{3} \\sum_{i = 0}^{n} \\Delta x \\ S(i) $$\n",
    "\n",
    "Where $n$ is even and $\\Delta x = \\frac{b - a}{n} $.\n",
    "\n",
    "While it is possible to start the summation at 1 (to match the formulas for the left and right Riemann sums), I think it is simpler to start at 0."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is $S(i)$? This is a special *piecewise* function that determines the coefficient (if any) to multiply $f(x_{i})$ by."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$S(i) =\n",
    "\\begin{cases}\n",
    "f(x_{i}),  & \\text{if $i$ is 0 or $n$} \\\\\n",
    "2 \\ f(x_{i}), & \\text{if $i$ is even} \\\\\n",
    "4 \\ f(x_{i}), & \\text{if $i$ is odd}\n",
    "\\end{cases}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that precedence is important here. I don't know how to represent this in math notation, but programatically the logic is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Will not actually run because n, x_i, and f(x) are not defined\n",
    "def S(i):\n",
    "    if (i == 0) or (i == n):\n",
    "        return f(x_i)\n",
    "    elif (i % 2 == 0):\n",
    "        return 2 * f(x_i)\n",
    "    else:\n",
    "        # odd\n",
    "        return 4 * f(x_i)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, $x_{i}$ is a function that returns the value to plug into $f(x)$ given the current index $i$.\n",
    "\n",
    "$$x_{i} = a + i \\ \\Delta x $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The notation, $x_{i}$, is somewhat confusing because it looks like a variable but is actually a function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Worked Example\n",
    "\n",
    "This is from an old uOttawa Calculus I exam.\n",
    "\n",
    "Approximate the integral $ \\int_{0}^4 \\frac{1}{x+1} dx $ numerically to 4 decimal places with n = 4 subdivisions with Simpson's Rule.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Lower bound\n",
    "a = 0\n",
    "# Upper bound\n",
    "b = 4\n",
    "# Number of intervals\n",
    "n = 4\n",
    "# Difference in x between intervals\n",
    "deltaX = (b - a) / n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Function being integrated\n",
    "def f(x):\n",
    "    return 1.0 / (x + 1.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Determines what value to plug into f(x) for each i\n",
    "def currentX(i):\n",
    "    return a + (i * deltaX)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Function that applies Simpson's coefficients (based on i)\n",
    "def S(i):\n",
    "    # if zeroth and last values\n",
    "    if (i == 0) or (i == n):\n",
    "        # no coefficient\n",
    "        return f(currentX(i))\n",
    "\n",
    "    # if even\n",
    "    elif (i % 2 == 0):\n",
    "        # coefficient of 2\n",
    "        return 2 * f(currentX(i))\n",
    "    \n",
    "    # if odd\n",
    "    else:\n",
    "        # coefficient of 4\n",
    "        return 4 * f(currentX(i))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the return statements are different from the previously defined $S(i)$ function, here $currentX()$ is used to determine the value to submit to f(x). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.622222222222222"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Computation\n",
    "simpsonSum = 0.0\n",
    "\n",
    "# Loop from 0 to n (inclusive)\n",
    "for i in range(0, n+1):\n",
    "    simpsonSum += S(i)\n",
    "\n",
    "# Don't forget to complete the formula by multiplying the sum by (deltaX / 3.0)!\n",
    "(deltaX / 3.0)*simpsonSum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which is roughly equivalent to <a href=\"http://www.wolframalpha.com/input/?i=integrate+1%2F(x+%2B+1)+using+Simpson%27s+rule+with+4+intervals+from+x+%3D+0+to+4\">Wolfram Alpha's answer</a>."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda root]",
   "language": "python",
   "name": "conda-root-py"
  },
  "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}