{ "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 Wolfram Alpha's answer." ] } ], "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 }