quad method from the integrate submodule of scipy\n",
"Sophisticated numerical integration routines do much of the \"dirty work\" for you. For example, you don't have \n",
"to specify the number of intervals (for simple integrations); and the routines will return an estimate of the error (but you do retain some controls).\n",
"\n",
"Here's how to integrate the function f from above from 0 to 3:\n"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "c8bac3cb",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(48.599999999999994, 5.39568389967826e-13)"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"integrate.quad(f,0,3)"
]
},
{
"cell_type": "markdown",
"id": "75b0be80",
"metadata": {},
"source": [
"The first element of the tuple returned by the function is the approximate value, and the second is \n",
"the estimated error. These are consistent with the exact answer.\n",
"\n",
"If you just want the value of the intergral, just pick of the first value returned by quad:"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "8bcbfb93",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"48.599999999999994"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"integrate.quad(f,0,3)[0]"
]
},
{
"cell_type": "markdown",
"id": "69fb1fef",
"metadata": {},
"source": [
"Here's what happens when quad is used on the more problematic function $g(x)$:"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "becb0487",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/home/mligare/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:1: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.\n",
" If increasing the limit yields no improvement it is advised to analyze \n",
" the integrand in order to determine the difficulties. If the position of a \n",
" local difficulty can be determined (singularity, discontinuity) one will \n",
" probably gain from splitting up the interval and calling the integrator \n",
" on the subranges. Perhaps a special-purpose integrator should be used.\n",
" \"\"\"Entry point for launching an IPython kernel.\n"
]
},
{
"data": {
"text/plain": [
"(1.3311544977190692, 0.0006258410606639012)"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"integrate.quad(g,0,2)"
]
},
{
"cell_type": "markdown",
"id": "917b14a0",
"metadata": {},
"source": [
"The value that is returned is accompanied by a warning: there's something fishy with this function on this \n",
"interval, and quad won't return a value without a warning that the default maximum of the number \n",
"of intervals isn't \"good enough\" to reach the default precision. Therefore, the user must either\n",
"change the default on the number of intervals, or else change the default precision.\n",
"\n",
"After looking at the graph of the function $g(x)$, I choose to increase the allowed number of intervals:"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "8ab7003f",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(1.3311544977210887, 1.870593173189536e-08)"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"integrate.quad(g,0,2,limit=80)"
]
},
{
"cell_type": "markdown",
"id": "f40480a9",
"metadata": {},
"source": [
"I get an output without a pink warning! But note that the number of intervals needed by quad \n",
"is very much smaller than the number needed to achieve a precision 1.9e-08 when I use the simple numerical \n",
"integration I introduced above. What's going on?\n",
"\n",
"Well, quad is more sophisticated than I am. It uses adaptive step sizes, meaning \n",
"that it uses small intervals $\\Delta x$ when necessary, i.e., when the the period of the sin\n",
"function is small, and larger intervals where the period is larger.\n",
"\n",
"How adaptive steps sizes work is not a topic for this course. "
]
},
{
"cell_type": "markdown",
"id": "0e5cd4f8",
"metadata": {},
"source": [
"## Two-dimensional Integration"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "5db6f3fb",
"metadata": {},
"outputs": [],
"source": [
"def f2(y,x):\n",
" return 1\n",
"\n",
"def f3(y,x):\n",
" return x\n",
"\n",
"def ramp(x):\n",
" return 0.5*x\n",
"\n",
"def circle(x):\n",
" return np.sqrt(1 - x**2)\n",
"\n",
"def ncircle(x):\n",
" return -np.sqrt(1 - x**2)"
]
},
{
"cell_type": "markdown",
"id": "2cb40baa",
"metadata": {},
"source": [
"#### Area of a rectangle"
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "086e345a",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(2.0, 2.220446049250313e-14)"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x1 = 0\n",
"x2 = 4\n",
"y1 = 0\n",
"y2 = 0.5\n",
"\n",
"integrate.dblquad(f2, x1, x2, y1, y2)"
]
},
{
"cell_type": "markdown",
"id": "d3fffb19",
"metadata": {},
"source": [
"#### Area of a triangle"
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "e48cb12b",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(4.0, 4.440892098500626e-14)"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x1 = 0\n",
"x2 = 4\n",
"y1 = 0\n",
"y2 = ramp\n",
"\n",
"integrate.dblquad(f2, x1, x2, y1, y2)"
]
},
{
"cell_type": "markdown",
"id": "2a1a955e",
"metadata": {},
"source": [
"$$\n",
"\\int_{\\rm y=10}^{15}\\int_{\\rm x = 0}^4 x\\, dx\\, dy = 40\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "b4e8e39d",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(40.0, 4.440892098500626e-13)"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x1 = 0\n",
"x2 = 4\n",
"y1 = 10\n",
"y2 = 15\n",
"integrate.dblquad(f3, x1, x2, y1, y2)"
]
},
{
"cell_type": "markdown",
"id": "721bd8f9",
"metadata": {},
"source": [
"#### Area of a quadrant of a unit circle\n",
"Answer should be $\\pi/4 = 0.7854$ "
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "bf12afa2",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(0.7853981633974481, 8.833911380179416e-11)"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x1 = 0\n",
"x2 = 1\n",
"y1 = 0\n",
"y2 = circle\n",
"\n",
"integrate.dblquad(f2, x1, x2, y1, y2)"
]
},
{
"cell_type": "markdown",
"id": "720ab776",
"metadata": {},
"source": [
"#### Area of a complete unit circle¶\n",
"Answer should be $\\pi$"
]
},
{
"cell_type": "code",
"execution_count": 20,
"id": "e0422ea2",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(3.1415926535897967, 2.000470900043183e-09)"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x1 = -1\n",
"x2 = 1\n",
"y1 = ncircle\n",
"y2 = circle\n",
"\n",
"integrate.dblquad(f2, x1, x2, y1, y2)"
]
},
{
"cell_type": "markdown",
"id": "ba8a924d",
"metadata": {},
"source": [
"#### Order of limits clarification -- Be careful!\n",
"\n",
"$$\n",
"\\iint_{x=0, y=0}^{x=1, y=2} x y^2 dx\\, dy = \\frac{4}{3}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 21,
"id": "028fdb2c",
"metadata": {},
"outputs": [],
"source": [
"def f4(y,x):\n",
" return x*y**2"
]
},
{
"cell_type": "code",
"execution_count": 22,
"id": "7e1bae9b",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"((1.3333333333333335, 2.954166042199163e-14),\n",
" (1.9999999999999998, 2.2204460492503128e-14))"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"integrate.dblquad(f4, 0, 1, 0, 2), integrate.dblquad(f3, 0, 2, 0, 1)"
]
},
{
"cell_type": "markdown",
"id": "c6b60f80",
"metadata": {},
"source": [
"## Three-dimensional integration"
]
},
{
"cell_type": "markdown",
"id": "0342b213",
"metadata": {},
"source": [
"#### Volume of octant of unit sphere. Resulty should be $\\pi/6 = 0.5236$"
]
},
{
"cell_type": "code",
"execution_count": 23,
"id": "476327b5",
"metadata": {},
"outputs": [],
"source": [
"def f4(x,y,z):\n",
" return 1\n",
"\n",
"def sphere(x,y):\n",
" return np.sqrt(1- x**2 - y**2)"
]
},
{
"cell_type": "code",
"execution_count": 24,
"id": "97d16cae",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(0.5235987755982989, 8.83377815341646e-11)"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x1 = 0\n",
"x2 = 1\n",
"y1 = 0\n",
"y2 = circle\n",
"z1 = 0\n",
"z2 = sphere\n",
"\n",
"integrate.tplquad(f4, x1, x2, y1, y2, z1, z2)"
]
},
{
"cell_type": "markdown",
"id": "2ed6193f",
"metadata": {},
"source": [
"# Version Information \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", "\n", "If `version_information` has been installed system wide (as it has been on Bucknell linux computers with shared file systems), 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": 25, "id": "1762de81", "metadata": {}, "outputs": [], "source": [ "load_ext version_information" ] }, { "cell_type": "code", "execution_count": 26, "id": "5442e12d", "metadata": {}, "outputs": [ { "data": { "application/json": { "Software versions": [ { "module": "Python", "version": "3.7.13 64bit [GCC 7.5.0]" }, { "module": "IPython", "version": "7.31.1" }, { "module": "OS", "version": "Linux 4.9.0 9 amd64 x86_64 with debian 9.13" }, { "module": "numpy", "version": "1.21.5" }, { "module": "scipy", "version": "1.7.3" }, { "module": "matplotlib", "version": "3.5.2" } ] }, "text/html": [ "
| Software | Version |
|---|---|
| Python | 3.7.13 64bit [GCC 7.5.0] |
| IPython | 7.31.1 |
| OS | Linux 4.9.0 9 amd64 x86_64 with debian 9.13 |
| numpy | 1.21.5 |
| scipy | 1.7.3 |
| matplotlib | 3.5.2 |
| Mon Nov 14 14:13:59 2022 EST | |