tests/inputs/Trapezoid_Rule.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Basic Numerical Integration: the Trapezoid Rule"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A simple illustration of the trapezoid rule for definite integration:\n",
    "\n",
    "$$\n",
    "\\int_{a}^{b} f(x)\\, dx \\approx \\frac{1}{2} \\sum_{k=1}^{N} \\left( x_{k} - x_{k-1} \\right) \\left( f(x_{k}) + f(x_{k-1}) \\right).\n",
    "$$\n",
    "<br>\n",
    "First, we define a simple function and sample it between 0 and 10 at 200 points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    return (x-3)*(x-5)*(x-7)+85\n",
    "\n",
    "x = np.linspace(0, 10, 200)\n",
    "y = f(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Choose a region to integrate over and take only a few points in that region"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "a, b = 1, 8 # the left and right boundaries\n",
    "N = 5 # the number of points\n",
    "xint = np.linspace(a, b, N)\n",
    "yint = f(xint)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot both the function and the area below it in the trapezoid approximation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.plot(x, y, lw=2)\n",
    "plt.axis([0, 9, 0, 140])\n",
    "plt.fill_between(xint, 0, yint, facecolor='gray', alpha=0.4)\n",
    "plt.text(0.5 * (a + b), 30,r\"$\\int_a^b f(x)dx$\", horizontalalignment='center', fontsize=20);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the integral both at high accuracy and with the trapezoid approximation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The integral is: 565.2499999999999 +/- 6.275535646693696e-12\n",
      "The trapezoid approximation with 5 points is: 559.890625\n"
     ]
    }
   ],
   "source": [
    "from __future__ import print_function\n",
    "from scipy.integrate import quad\n",
    "integral, error = quad(f, a, b)\n",
    "integral_trapezoid = sum( (xint[1:] - xint[:-1]) * (yint[1:] + yint[:-1]) ) / 2\n",
    "print(\"The integral is:\", integral, \"+/-\", error)\n",
    "print(\"The trapezoid approximation with\", len(xint), \"points is:\", integral_trapezoid)"
   ]
  }
 ],
 "metadata": {
  "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.4.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}