We load some packages and then begin:
","metadata":{}}, {"outputs":[],"cell_type":"code","source":["using Plots, SymPy, LinearAlgebra"],"metadata":{},"execution_count":1}, {"cell_type":"markdown","source":"Let $f: R^2 \\rightarrow R$ be a continuous, real-valued function of two variables. If $f(x,y) \\geq 0$ over a region $D$, then the volume under $f$ above the $x-y$ plan is given by an integral
","metadata":{}}, {"cell_type":"markdown","source":"\n$$\n\\int_D f(x,y) dA\n$$\n","metadata":{}}, {"cell_type":"markdown","source":"More generally, if $f(x,y)$ is possibly negative, then the above is interpreted as the signed volume. As with one-dimensional integrals, this insight can lead to methods to numerically compute the integral above. For now we focus on applying the fundamental theorem of calculus to perform the integration. For this, the double integral must be turned into a sequence of single integrals.
","metadata":{}}, {"cell_type":"markdown","source":"Fubini's theorem gives conditions on when it is possible to turn a double (or higher dimensional integral) into a sequence of iterated integrals. Fubini's theorem will apply if the region $D$ is the rectangle $[a,b] \\times [c,d]$ and the function $f$ is continuous, in which case we have
","metadata":{}}, {"cell_type":"markdown","source":"\n$$\n\\int_D f(x,y) dA = \\int_a^b [\\int_c^d f(x,y) dy] dx = \\int_c^d [\\int_a^b f(x,y) dx] dy.\n$$\n","metadata":{}}, {"cell_type":"markdown","source":"The value on the right can be computed by hand or with SymPy's integrate function within julia.
For example, consider the area under the plane $f(x,y) = 1 - (x + y)/2$ over the region $[0,1] \\times [0,1]$. Fubini's theorem says this is simply
","metadata":{}}, {"cell_type":"markdown","source":"\n$$\n\\int_0^1 [\\int_0^1 (1 - (x+y)/2) dy] dx.\n$$\n","metadata":{}}, {"cell_type":"markdown","source":"This can be computed as:
","metadata":{}}, {"outputs":[{"output_type":"execute_result","data":{"text/latex":["\\begin{equation*}\\frac{1}{2}\\end{equation*}"]},"metadata":{},"execution_count":1}],"cell_type":"code","source":["using SymPy\n@vars x y z real=true\n\nintegrate(1 - (x + y)/2, (y, 0, 1), (x, 0, 1))"],"metadata":{},"execution_count":1}, {"cell_type":"markdown","source":"The limits of integration are specified with the innermost first and outermost last. They are specified with a triple. So the math notation
","metadata":{}}, {"cell_type":"markdown","source":"\n$$\n\\int_a^b (...) dx\n$$\n","metadata":{}}, {"cell_type":"markdown","source":"Becomes integrate( ..., (x, a, b)) for whatever ... means.
More generally, the region may be more complicated than a rectangle. If $D$ is \"simple\", then one can say either $a \\leq x \\leq b$ and $g_1(x) \\leq y \\leq g_2(x)$ (or vice-versa in $x$ and $y$) for some functions $g_1$ and $g_2$, then the double integral over $D$ can be iterated to:
","metadata":{}}, {"cell_type":"markdown","source":"\n$$\n\\int_a^b \\int_{g_1(x)}^{g_2(x)} f(x,y) dy dx.\n$$\n","metadata":{}}, {"cell_type":"markdown","source":"For example, suppose $D$ is the region in the $x-y$ plane described by $0 \\leq x \\leq 1$ and $0 \\leq y \\leq x^2$, then the integral of $(x+y)^2$ over $D$ can be computed by
","metadata":{}}, {"outputs":[{"output_type":"execute_result","data":{"text/latex":["\\begin{equation*}\\frac{29}{70}\\end{equation*}"]},"metadata":{},"execution_count":1}],"cell_type":"code","source":["integrate((x+y)^2, (y, 0, x^2), (x,0,1))"],"metadata":{},"execution_count":1}, {"cell_type":"markdown","source":"If it is more convenient to express the domain $D$ in polar coordinates, then the double integral of $f$ over $D$ can be re-expressed as:
","metadata":{}}, {"cell_type":"markdown","source":"\n$$\n\\int_D f(x,y) dA = \\int_\\theta \\int_r f(r\\cos(\\theta), r\\sin(\\theta)) \\cdot r \\cdot dr \\cdot d\\theta.\n$$\n","metadata":{}}, {"cell_type":"markdown","source":"The $r \\cdot dr \\cdot d\\theta$ matches $dA$ after the transformation, of which more is described later.
","metadata":{}}, {"cell_type":"markdown","source":"For example, consider the problem of integrating $f(x,y) = x^2 + y^2$ over the unit circle. Using rectangular coordinates we have an error:
","metadata":{}}, {"outputs":[],"cell_type":"code","source":["# integrate(x^2 + y^2, (y, -sqrt(1-x^2), sqrt(1-x^2)), (x, -1, 1)) ## uncomment to see error"],"metadata":{},"execution_count":1}, {"cell_type":"markdown","source":"We can see where the error comes from by doing it one at a time:
","metadata":{}}, {"outputs":[{"output_type":"execute_result","data":{"text/latex":["\\begin{equation*}2 x^{2} \\sqrt{- x^{2} + 1} + \\frac{2 \\left(- x^{2} + 1\\right)^{\\frac{3}{2}}}{3}\\end{equation*}"]},"metadata":{},"execution_count":1}],"cell_type":"code","source":["@vars x y real=true\nex = integrate(x^2 + y^2, y)\nex1 = subs(ex, y, sqrt(1-x^2)) - subs(ex, y, -sqrt(1-x^2))"],"metadata":{},"execution_count":1}, {"cell_type":"markdown","source":"Trying integrate(ex1) should work – and does – but for some reason it also fails at times. Here we break this up into two pieces and integrate each:
The second integral gives:
","metadata":{}}, {"outputs":[{"output_type":"execute_result","data":{"text/latex":["\\begin{equation*}- \\frac{x^{3} \\sqrt{- x^{2} + 1}}{6} + \\frac{5 x \\sqrt{- x^{2} + 1}}{12} + \\frac{\\operatorname{asin}{\\left (x \\right )}}{4}\\end{equation*}"]},"metadata":{},"execution_count":1}],"cell_type":"code","source":["outb = integrate(exb, x)"],"metadata":{},"execution_count":1}, {"cell_type":"markdown","source":"(Though there can be errors some times when this is computed.)
","metadata":{}}, {"cell_type":"markdown","source":"The first integral has a condition:
","metadata":{}}, {"outputs":[{"output_type":"execute_result","data":{"text/latex":["\\begin{equation*}2 \\left(\\begin{cases} \\frac{i x^{5}}{4 \\sqrt{x^{2} - 1}} - \\frac{3 i x^{3}}{8 \\sqrt{x^{2} - 1}} - \\frac{i x \\operatorname{acosh}{\\left (\\left|{x}\\right| \\right )}}{8 \\left|{x}\\right|} + \\frac{\\pi x}{16 \\left|{x}\\right|} + \\frac{i x}{8 \\sqrt{x^{2} - 1}} & \\text{for}\\: x^{2} > 1 \\\\- \\frac{x^{5}}{4 \\sqrt{- x^{2} + 1}} + \\frac{3 x^{3}}{8 \\sqrt{- x^{2} + 1}} + \\frac{x \\operatorname{asin}{\\left (\\left|{x}\\right| \\right )}}{8 \\left|{x}\\right|} - \\frac{x}{8 \\sqrt{- x^{2} + 1}} & \\text{otherwise} \\end{cases}\\right)\\end{equation*}"]},"metadata":{},"execution_count":1}],"cell_type":"code","source":["out = integrate(exa, x)"],"metadata":{},"execution_count":1}, {"cell_type":"markdown","source":"We can fish out the value for $|x| < 1$, but it isn't pretty:
","metadata":{}}, {"outputs":[{"output_type":"execute_result","data":{"text/latex":["\\begin{equation*}- \\frac{x^{5}}{2 \\sqrt{- x^{2} + 1}} + \\frac{3 x^{3}}{4 \\sqrt{- x^{2} + 1}} + \\frac{x \\operatorname{asin}{\\left (\\left|{x}\\right| \\right )}}{4 \\left|{x}\\right|} - \\frac{x}{4 \\sqrt{- x^{2} + 1}}\\end{equation*}"]},"metadata":{},"execution_count":1}],"cell_type":"code","source":["outa = 2*collect((out/2).x.args[2].x)[1]"],"metadata":{},"execution_count":1}, {"cell_type":"markdown","source":"With this we can see the issue – evaluating at 1 or -1 gives an NaN. We could take a limit – but this doesn't work – rather we just put in a small value of $h$ and see:
An answer that looks like $\\pi/2$.
","metadata":{}}, {"cell_type":"markdown","source":"Okay, that was too much work and still not that satisfying! Let's see what happens if we change coordinate systems. The unit circle is described by $0 \\leq r \\leq 1$ and $0 \\leq \\theta \\leq 2\\pi$ – a simple region, and we have
","metadata":{}}, {"outputs":[{"output_type":"execute_result","data":{"text/latex":["\\begin{equation*}\\frac{\\pi}{2}\\end{equation*}"]},"metadata":{},"execution_count":1}],"cell_type":"code","source":["@vars r theta real=true\nf(x,y) = x^2 + y^2\nintegrate(f(r*cos(theta), r*sin(theta)) * r, (theta, 0, 2PI), (r, 0, 1)) # note extra \"* r\""],"metadata":{},"execution_count":1}, {"cell_type":"markdown","source":"Triple integrals are conceptually no more difficult than double integrals. Fubini's theorem is used to iterate the integrals, and if a region can be written iteratively, then the process is similar.
","metadata":{}}, {"cell_type":"markdown","source":"For example, compute $\\int\\int\\int_D z dV$ where $D$ is the region between the planes $z=x+y$, $z=3x+5y$ lying over the rectangle $[0,3] \\times [0,2]$. This becomes:
","metadata":{}}, {"outputs":[{"output_type":"execute_result","data":{"text/latex":["\\begin{equation*}294\\end{equation*}"]},"metadata":{},"execution_count":1}],"cell_type":"code","source":["@vars x y z\nintegrate(z, (z, x+y, 3x+5y), (y, 0, 2), (x, 0,3))"],"metadata":{},"execution_count":1}, {"cell_type":"markdown","source":"As polar coordinates make certain double integrals easier, triple integrals can be easier to compute using either cylindrical coordinate systems or spherical. This relates them:
","metadata":{}}, {"cell_type":"markdown","source":"\n$$\n\\int_D f(x,y,z) dV = \\int_\\theta \\int_r \\int_z f(r\\cos(\\theta), r\\sin(\\theta)z) r dz dr d\\theta\n= \\int_\\theta \\int_\\phi \\int_\\rho f(\\rho\\cos\\theta\\sin\\phi, \\rho\\sin\\theta\\sin\\phi, \\rho\\cos\\phi) \\rho^2 \\sin\\phi d\\rho d\\phi d\\theta.\n$$\n","metadata":{}}, {"cell_type":"markdown","source":"For example, computing the integral of $z$ over the region $D$ described by it fitting in the cylinder $x^2 + y^2\\leq 4$ with $0 \\leq z \\leq y$ we have noting that the cylinder is a circle with radius 2:
","metadata":{}}, {"outputs":[{"output_type":"execute_result","data":{"text/latex":["\\begin{equation*}\\pi\\end{equation*}"]},"metadata":{},"execution_count":1}],"cell_type":"code","source":["integrate(z * r, (z, 0, r*sin(theta)), (r, 0, 2), (theta, 0, PI))"],"metadata":{},"execution_count":1}, {"cell_type":"markdown","source":"The subtlety above is the limit of integration on $\\theta$ which is only $[0,\\pi]$, as the condition that $0 \\leq z \\leq y$ only allows values when the plane $z=y$ is positive, or the upper half of the $x-y$ plane.
","metadata":{}}, {"cell_type":"markdown","source":"The integrands for polar, cylindrical, and spherical integrals look like $f \\cdot G$ for some term $G$. (For example, $G=r$ for polar and cylindrical.) The $G$ comes from investigating the change in the \"$dA$\" or \"$dV$\" under the tranformation. What exactly $G$ is follows from the change of variable formula, written in two dimensions here:
","metadata":{}}, {"cell_type":"markdown","source":"\n$$\n\\int_D f(x,y) dx dy = \\int_D f(u,v) G(x,y;u,v) du dv.\n$$\n","metadata":{}}, {"cell_type":"markdown","source":"The $G$ is the determinant of the Jacobian, $J$, that is the determinant of the matrix (in this case) $\\partial x_i/ \\partial u_j$.
","metadata":{}}, {"cell_type":"markdown","source":"For example, suppose we have the transformation $x = r\\sin(\\theta)$ and $y=r\\cos(\\theta)$. Then $J$ is computed with:
\nWe help SymPy out by specifying the type of the output on the comprehension. We see J is a 2 by 2 matrix and its determinant can be taken with det and then simplified through:
To see how spherical coordinates can fit into the change of variable formula, we have:
\nHere is an application where we aren't just verifying an output. Let $D$ be the region between the lines $y=1$, $y=3$ and the parallel lines $x+2y=6$ and $x+2y=10$. Integrate $(x+3y)$.
\nThe change of variable $x = u - 2v$, $y=v$ translates into
","metadata":{}}, {"outputs":[{"output_type":"execute_result","data":{"text/latex":["\\begin{equation*}u\\end{equation*}"]},"metadata":{},"execution_count":1}],"cell_type":"code","source":["@vars u v real=true\nxs = x, y = u - 2v, v\nx+2y"],"metadata":{},"execution_count":1}, {"cell_type":"markdown","source":"So with these coordinate, $D$ becomes a rectangular region $1 \\leq v \\leq 3$ and $6 \\leq u \\leq 10$. We then have:
","metadata":{}}, {"outputs":[{"output_type":"execute_result","data":{"text/latex":["\\begin{equation*}80\\end{equation*}"]},"metadata":{},"execution_count":1}],"cell_type":"code","source":["J = Sym[diff(x, u) for x in xs, u in (u, v)]\nintegrate((x + 3y) * det(J), (v, 1, 3), (u, 6, 10))"],"metadata":{},"execution_count":1}, {"cell_type":"markdown","source":"As an aside, this is an example of a linear transformation called a shear. From the linear algebra viewpoint, the transformation can be written in matrix notation as xs = J * us (if us and xs were vectors and not tuples). That J is the same as the Jacobian is related to the linear nature of the transformation.
Compute the area of the right half of the centroid $r^2 = \\cos(2\\theta)$.
\nThe key is $\\theta$ goes from $-\\pi/4$ to $\\pi/4$ to trace out the region, giving:
","metadata":{}}, {"outputs":[{"output_type":"execute_result","data":{"text/latex":["\\begin{equation*}\\frac{1}{2}\\end{equation*}"]},"metadata":{},"execution_count":1}],"cell_type":"code","source":["integrate(1*r, (r, 0, sqrt(cos(2theta))), (theta, -PI/4, PI/4))"],"metadata":{},"execution_count":1}, {"cell_type":"markdown","source":"Find the center of mass of the domain bounded by $y=1-x^2$ and the $x$ axis with density $\\rho(x,y) = y$.
\nThe center of mass in the $x$ direction is given by $(\\int_D x \\cdot \\rho \\cdot DA)/M$ and similarly with a $y$.
","metadata":{}}, {"cell_type":"markdown","source":"We have:
","metadata":{}}, {"outputs":[{"output_type":"execute_result","data":{"text/plain":["(0, 4/7)"]},"metadata":{},"execution_count":1}],"cell_type":"code","source":["@vars x y real=true\nrho = y\nM = integrate(rho, (y, 0, 1-x^2), (x, -1, 1))\nMx = integrate(x * rho, (y, 0, 1-x^2), (x, -1, 1))\nMy = integrate(y * rho, (y, 0, 1-x^2), (x, -1, 1))\nMx/M, My/M"],"metadata":{},"execution_count":1}, {"cell_type":"markdown","source":"That $M_x$ is $0$ is no surprise, as there is symmetry in both the region and the density to the left and right of the line $x=0$.
","metadata":{}} ], "metadata": { "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "0.6" }, "kernelspec": { "display_name": "Julia 1.0.0", "language": "julia", "name": "julia-1.0" } }, "nbformat": 4, "nbformat_minor": 2 }