{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical integration\n",
    "\n",
    "Simple example of numerical integration of a one-dimensional and multi-dimensional integrals using the `quad()` from the `integrate` sub-module of `scipy`.\n",
    "\n",
    "Marty Ligare, \n",
    "August 2020"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import integrate   # importing sub-module of scipy\n",
    "\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Following is an Ipython magic command that puts figures in the  notebook.\n",
    "%matplotlib notebook\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",
    "mpl.style.use('classic')\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": 4,
   "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 np.sin(x-2)/(x-2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3.7500000000000004, 4.1633363423443377e-14)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "value, error = 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": 6,
   "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-6-da1962f33e71>\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[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[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[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/remote/anaconda-3.7-2020-05-28/lib/python3.7/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    350\u001b[0m     \u001b[0;32mif\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    351\u001b[0m         retval = _quad(func, a, b, args, full_output, epsabs, epsrel, limit,\n\u001b[0;32m--> 352\u001b[0;31m                        points)\n\u001b[0m\u001b[1;32m    353\u001b[0m     \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    354\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mpoints\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\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[0;32m/usr/remote/anaconda-3.7-2020-05-28/lib/python3.7/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    461\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[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    462\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[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 463\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[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    464\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    465\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[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: must be real number, not NoneType"
     ]
    }
   ],
   "source": [
    "value, error = integrate.quad(func2,1,3)\n",
    "value, error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Specify troublesome points:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(8.666666666666668, 9.621932880084691e-14)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "value, error = integrate.quad(func2,1,3,points=[2,])\n",
    "value, error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3.210825953605389, 3.5647329017567276e-14)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "value, error = 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": 9,
   "metadata": {},
   "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 np.sqrt(xulim**2-x**2)\n",
    "\n",
    "def q(x,y):\n",
    "    return 0\n",
    "\n",
    "def r(x,y):\n",
    "    return np.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 = integrate.tplquad(func4, xllim, xulim, g, h, q, r)\n",
    "8*value, error, 4*np.pi*xulim**3/3."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With $g$, $h$, $q$, and $r$ defined more concisely:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = lambda x: 0\n",
    "h = lambda x: np.sqrt(xulim**2-x**2)\n",
    "q = lambda x,y: 0\n",
    "r = lambda x,y: np.sqrt(xulim**2-x**2-y**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(33.51032163829113, 3.533511261366584e-10, 33.510321638291124)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xllim = 0\n",
    "xulim = 2\n",
    "value, error = integrate.tplquad(func4, xllim, xulim, g, h, q, r)\n",
    "8*value, error, 4*np.pi*xulim**3/3."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or even more concisely:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(33.51032163829113, 3.533511261366584e-10, 33.510321638291124)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xllim = 0\n",
    "xulim = 2\n",
    "value, error = integrate.tplquad(func4, xllim, xulim, \\\n",
    "                                    lambda x:0, lambda x:np.sqrt(xulim**2-x**2), \\\n",
    "                                    lambda x,y:0, lambda x,y:np.sqrt(xulim**2-x**2-y**2))\n",
    "8*value, error, 4*np.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": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def func5(phi,theta,r):\n",
    "    return r**2*np.sin(theta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(33.510321638291124, 5.556380318779972e-13, 33.510321638291124)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rllim = 0\n",
    "rulim = 2\n",
    "value, error = integrate.tplquad(func5, rllim, rulim, \\\n",
    "                                    lambda theta:0, lambda theta:np.pi, \\\n",
    "                                    lambda theta,phi:0, lambda theta,phi:2.*np.pi)\n",
    "value, error, 4*np.pi*xulim**3/3."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Version information\n",
    "`version_information` is from J.R. Johansson (jrjohansson at gmail.com); see <a href='http://nbviewer.jupyter.org/github/jrjohansson/scientific-python-lectures/blob/master/Lecture-0-Scientific-Computing-with-Python.ipynb'>Introduction to scientific computing with Python</a> for more information and instructions for package installation.\n",
    "\n",
    "`version_information` is installed on the linux network at Bucknell"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext version_information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/json": {
       "Software versions": [
        {
         "module": "Python",
         "version": "3.7.7 64bit [GCC 7.3.0]"
        },
        {
         "module": "IPython",
         "version": "7.16.1"
        },
        {
         "module": "OS",
         "version": "Linux 3.10.0 1062.9.1.el7.x86_64 x86_64 with centos 7.7.1908 Core"
        },
        {
         "module": "numpy",
         "version": "1.18.5"
        },
        {
         "module": "scipy",
         "version": "1.5.2"
        },
        {
         "module": "matplotlib",
         "version": "3.3.0"
        }
       ]
      },
      "text/html": [
       "<table><tr><th>Software</th><th>Version</th></tr><tr><td>Python</td><td>3.7.7 64bit [GCC 7.3.0]</td></tr><tr><td>IPython</td><td>7.16.1</td></tr><tr><td>OS</td><td>Linux 3.10.0 1062.9.1.el7.x86_64 x86_64 with centos 7.7.1908 Core</td></tr><tr><td>numpy</td><td>1.18.5</td></tr><tr><td>scipy</td><td>1.5.2</td></tr><tr><td>matplotlib</td><td>3.3.0</td></tr><tr><td colspan='2'>Fri Aug 07 15:21:07 2020 EDT</td></tr></table>"
      ],
      "text/latex": [
       "\\begin{tabular}{|l|l|}\\hline\n",
       "{\\bf Software} & {\\bf Version} \\\\ \\hline\\hline\n",
       "Python & 3.7.7 64bit [GCC 7.3.0] \\\\ \\hline\n",
       "IPython & 7.16.1 \\\\ \\hline\n",
       "OS & Linux 3.10.0 1062.9.1.el7.x86\\_64 x86\\_64 with centos 7.7.1908 Core \\\\ \\hline\n",
       "numpy & 1.18.5 \\\\ \\hline\n",
       "scipy & 1.5.2 \\\\ \\hline\n",
       "matplotlib & 3.3.0 \\\\ \\hline\n",
       "\\hline \\multicolumn{2}{|l|}{Fri Aug 07 15:21:07 2020 EDT} \\\\ \\hline\n",
       "\\end{tabular}\n"
      ],
      "text/plain": [
       "Software versions\n",
       "Python 3.7.7 64bit [GCC 7.3.0]\n",
       "IPython 7.16.1\n",
       "OS Linux 3.10.0 1062.9.1.el7.x86_64 x86_64 with centos 7.7.1908 Core\n",
       "numpy 1.18.5\n",
       "scipy 1.5.2\n",
       "matplotlib 3.3.0\n",
       "Fri Aug 07 15:21:07 2020 EDT"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "version_information numpy, scipy, matplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}