numerical integration with scipy quad.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A comparison of Scipy's quad integration\n",
    "\n",
    "scipy's quad algorithm is based on the FORTRAN library QUADPACK.\n",
    "\n",
    "You can read more up on QUADPACK if you're feeling ambitious http://people.sc.fsu.edu/~jburkardt/f_src/quadpack/quadpack.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Scipy is part of the Anaconda distribution\n",
    "\n",
    "So make sure you have Anaconda with scipy installed!\n",
    "\n",
    "Refer to the scipy documentation https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html#scipy.integrate.quad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import integrate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    \"\"\"simple function to integrate\"\"\"\n",
    "    return np.sin(x)\n",
    "\n",
    "\n",
    "def trap(f, xmin, xmax, npoints=10):\n",
    "    \"\"\"\n",
    "    computes the integral of f using trapezoid rule\n",
    "    \"\"\"\n",
    "    area = 0\n",
    "    x = np.linspace(xmin, xmax, npoints)\n",
    "    N = len(x)\n",
    "    dx = x[1] - x[0]\n",
    "    \n",
    "    for k in range(1, N):\n",
    "        area += (f(x[k - 1]) + f(x[k])) * dx / 2\n",
    "        \n",
    "    return area"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Integration using trapezoid rule\n",
    "\n",
    "Let's use the trapezoid function with enough points to get something close to the accuracy of scipy's quad algorithm\n",
    "\n",
    "This will return the error and also the execution time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loop, best of 3: 2.9 s per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit trap(f, 0, np.pi, 1500000) - 2   # 2 is the actual integral value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Integration using Scipy's quad algorithm\n",
    "\n",
    "no let's use quad and time how fast and accurate it is."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10000 loops, best of 3: 34 µs per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit integrate.quad(f, 0, np.pi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## WOW! \n",
    "\n",
    "quad is about 80,000 times faster and still an order of magnitude more accurate!!\n",
    "\n",
    "Isn't Scipy amazing!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}