{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# \n", "\n", "# Numeric integration with Julia\n", "\n", "To get started, we load the `MTH229` package:" ], "id": "1c1d3a1e-ba4f-4a51-8cf1-71d2b8e9f646" }, { "cell_type": "code", "execution_count": 0, "metadata": {}, "outputs": [], "source": [ "using MTH229\n", "using Plots\n", "plotly()" ], "id": "2" }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Quick background\n", "\n", "Read more about this material here:\n", "[integration](http://mth229.github.io/integration.html).\n", "\n", "In many cases the task of evaluating a definite integral is made easy by\n", "the fundamental theorem of calculus, one part of which states that for a\n", "continuous function $f$ with (any) antiderivative, $F$, the definite\n", "integral can be computed through:\n", "\n", "$$\n", "\\int_a^b f(x) dx = F(b) - F(a).\n", "$$\n", "\n", "That is, the definite integral is found by evaluating a related function\n", "at the endpoints, $a$ and $b$.\n", "\n", "The `SymPy` package can compute many antiderivatives using a version of\n", "the [Risch algorithm](http://en.wikipedia.org/wiki/Risch_algorithm) that\n", "works for [elementary\n", "functions](http://en.wikipedia.org/wiki/Elementary_function). `SymPy`’s\n", "`integrate` function can be used to find an indefinite integral:" ], "id": "320269e7-8090-4e96-9e46-731005273bf6" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/plain": [ " 3\n", "x \n", "──\n", "3 " ] } } ], "source": [ "f(x) = x^2\n", "@syms x\n", "integrate(f(x), x)" ], "id": "6" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Or a definite integral by passing in values `a` and `b` through the\n", "grouping `(x,a,b)`:" ], "id": "195b47c0-a91f-4ee1-b478-fe0a02ffb433" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/plain": [ "1/3" ] } } ], "source": [ "integrate(f(x), (x, 0, 1)) # returns a \"symbolic\" number" ], "id": "8" }, { "cell_type": "markdown", "metadata": {}, "source": [ "However, this only works *if* there is a known antiderivative\n", "$F(x)$—which is not always the case. If not, what to do?\n", "\n", "In this case, we can appeal to the *definition* of the definite\n", "integral. For continuous non-negative $f(x)$, the definite integral is\n", "the area under the graph of $f$ over the interval $[a,b]$. For possibly\n", "negative functions, the indefinite integral is found by the *signed*\n", "area under $f$. This area can be directly *approximated* using Riemann\n", "sums, or some other approximation scheme.\n", "\n", "The Riemann approximation for a definite integral uses approximating\n", "rectangles. The following pattern will compute a Riemann sum with\n", "equal-sized partitions using right-hand endpoints:" ], "id": "4268acd5-6b7e-4132-9bfb-fc5f473dba25" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/plain": [ "0.44000000000000006" ] } } ], "source": [ "f(x) = x^2\n", "a, b, n = 0, 1, 5 # 5 partitions of [0,1] requested\n", "delta = (b - a)/n # size of partition\n", "xs = a .+ (1:n) * delta # a < x₁ < ⋯ < xₙ₋₁ < xₙ = b\n", "sum(f.(xs)) * delta # a new function `sum` to add up values in a container" ], "id": "10" }, { "cell_type": "markdown", "metadata": {}, "source": [ "That value isn’t very close to $1/3$. But we only took $n=5$ rectangles\n", "$-$ clearly there will be some error." ], "id": "9cce8d62-50cc-40bf-8e36-9bcc889ad624" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/html": [ "" ] } } ], "source": [ "plot(f, a, b; legend=false)\n", "riemann_plot!(f, a, b, n; fill=(\"red\", 0.25, 0))" ], "id": "12" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Bigger $n$s mean better approximations. This figure overlays the picture\n", "with $n=10$." ], "id": "01964be9-491e-4845-b490-9f65a7c1d877" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/html": [ "" ] } } ], "source": [ "riemann_plot!(f, a, b, 10; fill=(\"blue\", 0.25, 0))" ], "id": "14" }, { "cell_type": "markdown", "metadata": {}, "source": [ "With large `n` the figures become too crowded to illustrate, but we can\n", "easily use a large value of `n` in our computations. With $n=50,000$ we\n", "have the first $4$ decimal points are accurate:" ], "id": "eccc5c25-a9fd-4f14-a655-24cc216f111d" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/plain": [ "0.33334333340000005" ] } } ], "source": [ "f(x) = x^2\n", "a, b, n = 0, 1, 50_000 # 50,000 partitions of [0,1] requested\n", "delta = (b - a)/n\n", "xs = a .+ (1:n) * delta\n", "sum(f.(xs)) * delta" ], "id": "16" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that only the first two lines of the calculation needed changing to\n", "adjust to a new problem.\n", "\n", "------------------------------------------------------------------------\n", "\n", "As the pattern is similar, it is fairly easy to wrap the computations in\n", "a function for convenience. This is done for you in the `MTH229` package\n", "with the `riemann` function. The above would be done through:" ], "id": "03602194-68aa-4823-a3bb-545ac2898b19" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/plain": [ "0.33334333339999955" ] } } ], "source": [ "f(x) = x^2\n", "a, b, n = 0, 1, 50_000\n", "riemann(f, a, b, n)" ], "id": "18" }, { "cell_type": "markdown", "metadata": {}, "source": [ "------------------------------------------------------------------------\n", "\n", "The Riemann sum is slow to converge here. There are faster algorithms\n", "both mathematically and computationally. We will briefly discuss two:\n", "the [trapezoid](https://en.wikipedia.org/wiki/Trapezoidal_rule) rule,\n", "which replaces rectangles with trapezoids; and\n", "[Simpson’s](https://en.wikipedia.org/wiki/Simpson%27s_rule)rule which is\n", "a quadratic approximation. Each is invoked by passing a value to the\n", "`method` argument or `riemann`:" ], "id": "56e98710-088b-4ac7-9d7c-021e4898f5e0" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/plain": [ "0.3333334999999997" ] } } ], "source": [ "f(x) = x^2\n", "riemann(f, 0, 1, 1000, method=\"trapezoid\")" ], "id": "20" }, { "cell_type": "markdown", "metadata": {}, "source": [ "And for Simpson’s method:" ], "id": "19f66e91-4ec5-44bb-b687-987c3b474196" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/plain": [ "0.3333333333333334" ] } } ], "source": [ "riemann(f, 0, 1, 1000, method=\"simpsons\")" ], "id": "22" }, { "cell_type": "markdown", "metadata": {}, "source": [ "------------------------------------------------------------------------\n", "\n", "For real-world use, `Julia` has the `QuadGK` package and its function\n", "`quadgk`, which is loaded with `MTH229`. By using a different approach\n", "altogether, this function is much more efficient and estimates the\n", "potential error in the approximation. It is used quite easily—no $n$ is\n", "needed, as the algorithm is adaptive:" ], "id": "19832a51-3e78-475b-9f65-629c3ae90dbe" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/plain": [ "(0.3333333333333333, 0.0)" ] } } ], "source": [ "f(x) = x^2\n", "answer, err = quadgk(f, 0, 1)" ], "id": "24" }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `quadgk` function returns two values, an answer and an estimated\n", "maximum possible error. The answer is the first number, clearly it is\n", "$1/3$, and the estimated maximum error is the second. In this case it is\n", "small ($10^{-17}$) $-$ basically `0`.\n", "\n", "------------------------------------------------------------------------\n", "\n", "In summary we consider here these functions to find definite integrals:\n", "\n", "- `integrate`—symbolically find a definite integral using the\n", " fundamental theorem of calculus, **if possible**.\n", "\n", "- `riemann`—approximate a definite integral using either Riemann sums\n", " or a related method.\n", "\n", "- `quadgk`—use Gauss quadrature approach to efficiently find\n", " approximations to definite integrals.\n", "\n", "## Different interpretations of other integrals\n", "\n", "The integral can represent other quantities besides the area under a\n", "curve. We give two examples: arc-length and certain volumes of\n", "revolution.\n", "\n", "A formula to compute the length of the graph of a differentiable\n", "function $f(x)$ from $a$ to $b$ is given by the formula:\n", "\n", "$$\n", "\\int_a^b \\sqrt{1 + f'(x)^2} dx\n", "$$\n", "\n", "That is, a function *related* to $f$ is integrated and this has a\n", "different interpretation than the area under $f$.\n", "\n", "For example, the arc-length of the square root function between $[0, 4]$\n", "is given by:" ], "id": "a14cdbef-d749-4902-b102-846d6a8ea245" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/plain": [ "(4.6467837188706085, 6.719391001417838e-8)" ] } } ], "source": [ "f(x) = sqrt(x)\n", "dL(x) = sqrt(1 + f'(x)^2)\n", "answer, err = quadgk(dL, 0, 4)" ], "id": "26" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Application: Volume of glasses.\n", "\n", "We discuss an application of the integral to finding volumes—not areas.\n", "\n", "A *solid of revolution* is a figure with rotational symmetry around some\n", "axis, such as a soda can, a snow cone, a red solo cup, and other common\n", "objects. A formula for the volume of an object with rotational symmetry\n", "can be written in terms of an integral based on a function, $r(h)$,\n", "which specifies the radius for various values of $h$.\n", "\n", "> If the radius as a function of height is given by $r(h)$, the the\n", "> volume is $\\int_a^b \\pi r(h)^2 dh$.\n", "\n", "So for example, a baseball has a overall diameter of $2\\cdot 37$mm, but\n", "if we place the center at the origin, its rotational radius is given by\n", "$r(h) = (37^2 - h^2)^{1/2}$ for $-37 \\leq h \\leq 37$. The volume in\n", "mm$^3$ is given by:" ], "id": "297b6e8c-5750-421d-9017-494c20e1e28e" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/plain": [ "(212174.79024304505, 0.0)" ] } } ], "source": [ "r(h) = (37^2 - h^2)^(1/2)\n", "dv(h) = pi * r(h)^2\n", "answer, err = quadgk(dv, -37, 37)" ], "id": "28" }, { "cell_type": "markdown", "metadata": {}, "source": [ "The volume in cubic inches, then is:" ], "id": "e6578b5d-4836-4b0a-ae53-4d187facce1b" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/plain": [ "12.947700103145083" ] } } ], "source": [ "answer / (2.54 * 10)^3" ], "id": "30" }, { "cell_type": "markdown", "metadata": {}, "source": [ "------------------------------------------------------------------------" ], "id": "1c5bd6f0-10c5-48f7-ae30-bd32051598b1" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# Your commands go here\n" ], "id": "32" } ], "nbformat": 4, "nbformat_minor": 5, "metadata": { "kernel_info": { "name": "julia" }, "kernelspec": { "name": "julia", "display_name": "Julia", "language": "julia" }, "language_info": { "name": "julia", "codemirror_mode": "julia", "version": "1.10.0" } } }