{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical integration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy as sp\n",
    "from scipy import integrate   # not included in basic scipy\n",
    "\n",
    "import matplotlib as mpl       # As of July 2017 Bucknell computers use v. 2.x \n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Following is an Ipython magic command that puts figures in the  notebook.\n",
    "# For figures in separate windows, comment out following line and uncomment\n",
    "# the next line\n",
    "# Must come before defaults are changed.\n",
    "%matplotlib notebook\n",
    "#%matplotlib\n",
    "\n",
    "# As of Aug. 2017 reverting to 1.x defaults.\n",
    "# In 2.x text.ustex requires dvipng, texlive-latex-extra, and texlive-fonts-recommended, \n",
    "# which don't seem to be universal\n",
    "# See https://stackoverflow.com/questions/38906356/error-running-matplotlib-in-latex-type1cm?\n",
    "mpl.style.use('classic')\n",
    "        \n",
    "# M.L. modifications of matplotlib defaults using syntax of v.2.0 \n",
    "# More info at http://matplotlib.org/2.0.0/users/deflt_style_changes.html\n",
    "# Changes can also be put in matplotlibrc file, or effected using mpl.rcParams[]\n",
    "plt.rc('figure', figsize = (6, 4.5))            # Reduces overall size of figures\n",
    "plt.rc('axes', labelsize=16, titlesize=14)\n",
    "plt.rc('figure', autolayout = True)             # Adjusts supblot parameters for new size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Single variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def func1(x):      # Continuous function\n",
    "    return x**3\n",
    "\n",
    "def func2(x):      # Discontinuous function\n",
    "    if x< 2:\n",
    "        return x**2\n",
    "    if x>2:\n",
    "        return x**2\n",
    "    \n",
    "def func3(x):    # Function with a singluarity\n",
    "    return sp.sin(x-2)/(x-2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3.7500000000000004, 4.1633363423443377e-14)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "value, error = sp.integrate.quad(func1,1,2)\n",
    "value, error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If there are discontinuities or singularities, quad will fail, eg.,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "must be real number, not NoneType",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-4-228de3762d39>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mvalue\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0merror\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mintegrate\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mquad\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfunc2\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0mvalue\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0merror\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/remote/anaconda-3.6/lib/python3.6/site-packages/scipy/integrate/quadpack.py\u001b[0m in \u001b[0;36mquad\u001b[0;34m(func, a, b, args, full_output, epsabs, epsrel, limit, points, weight, wvar, wopts, maxp1, limlst)\u001b[0m\n\u001b[1;32m    321\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mweight\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    322\u001b[0m         retval = _quad(func, a, b, args, full_output, epsabs, epsrel, limit,\n\u001b[0;32m--> 323\u001b[0;31m                        points)\n\u001b[0m\u001b[1;32m    324\u001b[0m     \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    325\u001b[0m         retval = _quad_weight(func, a, b, args, full_output, epsabs, epsrel,\n",
      "\u001b[0;32m/usr/remote/anaconda-3.6/lib/python3.6/site-packages/scipy/integrate/quadpack.py\u001b[0m in \u001b[0;36m_quad\u001b[0;34m(func, a, b, args, full_output, epsabs, epsrel, limit, points)\u001b[0m\n\u001b[1;32m    386\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mpoints\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    387\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0minfbounds\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 388\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0m_quadpack\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_qagse\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfunc\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mb\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mfull_output\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mepsabs\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mepsrel\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mlimit\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    389\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    390\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0m_quadpack\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_qagie\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfunc\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mbound\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0minfbounds\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mfull_output\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mepsabs\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mepsrel\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mlimit\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: must be real number, not NoneType"
     ]
    }
   ],
   "source": [
    "value, error = sp.integrate.quad(func2,1,3)\n",
    "value, error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Specify troublesome points:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(8.666666666666668, 9.621932880084691e-14)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "value, error = sp.integrate.quad(func2,1,3,points=[2,])\n",
    "value, error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3.210825953605389, 3.5647329017567276e-14)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "value, error = sp.integrate.quad(func3,0,4,points=[2,])\n",
    "value, error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multivariable"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Volume of a sphere I "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ V = \\int_0^R\\int_0^\\sqrt{1-x^2}\\int_0^\\sqrt{1-x^2-y^2} \\, dz dy dx \\equiv \n",
    "  \\int_0^R\\int_{g(x)}^{h(x)}\\int_{q(x,y)}^{r(x,y)}\\, dz dy dx$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With $g$, $h$, $q$, and $r$ defined normally:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def func4(z,y,x):    # ORDER OF ARGUMENTS IMPORTANT\n",
    "    return 1\n",
    "\n",
    "def g(x):\n",
    "    return 0\n",
    "\n",
    "def h(x):\n",
    "    return sp.sqrt(xulim**2-x**2)\n",
    "\n",
    "def q(x,y):\n",
    "    return 0\n",
    "\n",
    "def r(x,y):\n",
    "    return sp.sqrt(xulim**2-x**2-y**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(33.51032163829113, 3.533511261366584e-10, 33.510321638291124)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xllim = 0\n",
    "xulim = 2\n",
    "value, error = sp.integrate.tplquad(func4, xllim, xulim, g, h, q, r)\n",
    "8*value, error, 4*sp.pi*xulim**3/3."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With $g$, $h$, $q$, and $r$ defined more concisely:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = lambda x: 0\n",
    "h = lambda x: sp.sqrt(xulim**2-x**2)\n",
    "q = lambda x,y: 0\n",
    "r = lambda x,y: sp.sqrt(xulim**2-x**2-y**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(33.51032163829113, 3.533511261366584e-10, 33.510321638291124)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xllim = 0\n",
    "xulim = 2\n",
    "value, error = sp.integrate.tplquad(func4, xllim, xulim, g, h, q, r)\n",
    "8*value, error, 4*sp.pi*xulim**3/3."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or even more concisely:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(33.51032163829113, 3.533511261366584e-10, 33.510321638291124)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xllim = 0\n",
    "xulim = 2\n",
    "value, error = sp.integrate.tplquad(func4, xllim, xulim, \\\n",
    "                                    lambda x:0, lambda x:sp.sqrt(xulim**2-x**2), \\\n",
    "                                    lambda x,y:0, lambda x,y:sp.sqrt(xulim**2-x**2-y**2))\n",
    "8*value, error, 4*sp.pi*xulim**3/3."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Volume of a sphere I I\n",
    "\n",
    "$$ V = \\int_0^{2\\pi}\\int_0^\\pi\\int_0^R r^2 \\sin\\theta\\, dr d\\theta d\\phi $$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def func5(phi,theta,r):\n",
    "    return r**2*sp.sin(theta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(33.510321638291124, 5.556380318779972e-13, 33.510321638291124)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rllim = 0\n",
    "rulim = 2\n",
    "value, error = sp.integrate.tplquad(func5, rllim, rulim, \\\n",
    "                                    lambda theta:0, lambda theta:sp.pi, \\\n",
    "                                    lambda theta,phi:0, lambda theta,phi:2.*sp.pi)\n",
    "value, error, 4*sp.pi*xulim**3/3."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Version information\n",
    "\n",
    "version_information is from J.R. Johansson (jrjohansson at gmail.com)\n",
    "See Introduction to scientific computing with Python:\n",
    "http://nbviewer.jupyter.org/github/jrjohansson/scientific-python-lectures/blob/master/Lecture-0-Scientific-Computing-with-Python.ipynb \n",
    "for more information and instructions for package installation.\n",
    "If version_information has been installed system wide (as it has been on linuxremotes), continue with next cell as written. If not, comment out top line in next cell and uncomment the second line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loading extensions from ~/.ipython/extensions is deprecated. We recommend managing extensions like any other Python packages, in site-packages.\n"
     ]
    }
   ],
   "source": [
    "%load_ext version_information\n",
    "\n",
    "#%install_ext http://raw.github.com/jrjohansson/version_information/master/version_information.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/json": {
       "Software versions": [
        {
         "module": "Python",
         "version": "3.6.1 64bit [GCC 4.4.7 20120313 (Red Hat 4.4.7-1)]"
        },
        {
         "module": "IPython",
         "version": "6.1.0"
        },
        {
         "module": "OS",
         "version": "Linux 3.10.0 327.36.3.el7.x86_64 x86_64 with redhat 7.2 Maipo"
        },
        {
         "module": "scipy",
         "version": "0.19.1"
        },
        {
         "module": "matplotlib",
         "version": "2.0.2"
        }
       ]
      },
      "text/html": [
       "<table><tr><th>Software</th><th>Version</th></tr><tr><td>Python</td><td>3.6.1 64bit [GCC 4.4.7 20120313 (Red Hat 4.4.7-1)]</td></tr><tr><td>IPython</td><td>6.1.0</td></tr><tr><td>OS</td><td>Linux 3.10.0 327.36.3.el7.x86_64 x86_64 with redhat 7.2 Maipo</td></tr><tr><td>scipy</td><td>0.19.1</td></tr><tr><td>matplotlib</td><td>2.0.2</td></tr><tr><td colspan='2'>Tue Aug 01 11:18:58 2017 EDT</td></tr></table>"
      ],
      "text/latex": [
       "\\begin{tabular}{|l|l|}\\hline\n",
       "{\\bf Software} & {\\bf Version} \\\\ \\hline\\hline\n",
       "Python & 3.6.1 64bit [GCC 4.4.7 20120313 (Red Hat 4.4.7-1)] \\\\ \\hline\n",
       "IPython & 6.1.0 \\\\ \\hline\n",
       "OS & Linux 3.10.0 327.36.3.el7.x86\\_64 x86\\_64 with redhat 7.2 Maipo \\\\ \\hline\n",
       "scipy & 0.19.1 \\\\ \\hline\n",
       "matplotlib & 2.0.2 \\\\ \\hline\n",
       "\\hline \\multicolumn{2}{|l|}{Tue Aug 01 11:18:58 2017 EDT} \\\\ \\hline\n",
       "\\end{tabular}\n"
      ],
      "text/plain": [
       "Software versions\n",
       "Python 3.6.1 64bit [GCC 4.4.7 20120313 (Red Hat 4.4.7-1)]\n",
       "IPython 6.1.0\n",
       "OS Linux 3.10.0 327.36.3.el7.x86_64 x86_64 with redhat 7.2 Maipo\n",
       "scipy 0.19.1\n",
       "matplotlib 2.0.2\n",
       "Tue Aug 01 11:18:58 2017 EDT"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "version_information scipy, matplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}