{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Quick background\n", "\n", "# Numeric integration with Julia\n", "\n", "To get started, we load the `MTH229` package:" ], "id": "4cdb7cb5-d742-4e41-ac6e-462a4eaa14e7" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "using MTH229\n", "using Plots\n", "plotly()" ], "id": "41042258" }, { "cell_type": "markdown", "metadata": {}, "source": [ "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.\n", "However, we don’t pursue that here, opting instead for *numeric*\n", "approximations to definite integrals.\n", "\n", "To compute numeric approximations to definite integrals, we can appeal\n", "to the *definition* of the definite integral. For continuous\n", "non-negative $f(x)$, the definite integral is the area under the graph\n", "of $f$ over the interval $[a,b]$. For possibly negative functions, the\n", "indefinite integral is found by the *signed* area under $f$. This area\n", "can be directly *approximated* using Riemann sums, or some other\n", "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": "7737daad-a2b2-421c-ac14-1e63cad05f3d" }, { "cell_type": "code", "execution_count": 3, "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": "380f017d" }, { "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": "f3fa76b0-f8f0-43af-b6b5-270def760529" }, { "cell_type": "code", "execution_count": 4, "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": "9161e50a" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Bigger $n$s mean better approximations. This figure overlays the picture\n", "with $n=10$." ], "id": "a369cb70-25f1-40c6-bb52-7610b407ed08" }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/html": [ "" ] } } ], "source": [ "riemann_plot!(f, a, b, 10; fill=(\"blue\", 0.25, 0))" ], "id": "336e0c73" }, { "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": "64f25cea-3818-43db-8515-e77ec383d7fa" }, { "cell_type": "code", "execution_count": 6, "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": "7e52709d" }, { "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 to compute an approximate sum is similar, it is fairly\n", "easy to wrap the computations in a function for convenience. This is\n", "done for you in the `MTH229` package with the `riemann` function. The\n", "above would be done through:" ], "id": "30b9df85-86c6-4abf-8141-834be278673f" }, { "cell_type": "code", "execution_count": 7, "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": "4c10e2a3" }, { "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": "b6c5af86-8ca2-4f12-979f-e15f9135c8d3" }, { "cell_type": "code", "execution_count": 8, "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": "d87c0ee2" }, { "cell_type": "markdown", "metadata": {}, "source": [ "And for Simpson’s method:" ], "id": "351c7945-7ec3-4b39-b892-bc66608d7dcb" }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/plain": [ "0.3333333333333334" ] } } ], "source": [ "riemann(f, 0, 1, 1000, method=\"simpsons\")" ], "id": "9f657bc4" }, { "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": "e93403b0-ee5d-4de8-b5c9-c102d10e2229" }, { "cell_type": "code", "execution_count": 10, "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": "3e913976" }, { "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", "- `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": "bd0a2f6e-2842-4141-a0d7-05e9cbe3b7b3" }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/plain": [ "(4.6467837188706085, 6.719391001417838e-8)" ] } } ], "source": [ "f(x) = sqrt(x)\n", "a, b = 0, 4\n", "dL(x) = sqrt(1 + f'(x)^2)\n", "answer, err = quadgk(dL, a, b)" ], "id": "1337623c" }, { "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": "cb158b25-f438-4813-97cf-14f20856387c" }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/plain": [ "(212174.79024304505, 0.0)" ] } } ], "source": [ "r(h) = (37^2 - h^2)^(1/2)\n", "a, b = -37, 37\n", "dv(h) = pi * r(h)^2\n", "answer, err = quadgk(dv, a, b)" ], "id": "b7fbbfba" }, { "cell_type": "markdown", "metadata": {}, "source": [ "The volume in cubic inches, then is:" ], "id": "44260424-b342-4bec-a43c-97269b2a2216" }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/plain": [ "12.947700103145083" ] } } ], "source": [ "answer / (2.54 * 10)^3" ], "id": "424d77ca" }, { "cell_type": "markdown", "metadata": {}, "source": [ "------------------------------------------------------------------------" ], "id": "6336c1d0-e4a9-46aa-930c-2769b982e3d2" }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "# Your commands go here" ], "id": "2f3ee123" } ], "nbformat": 4, "nbformat_minor": 5, "metadata": { "kernelspec": { "name": "julia-1.11", "display_name": "Julia 1.11.2", "language": "julia", "path": "/Users/jverzani/Library/Jupyter/kernels/julia-1.11" }, "language_info": { "name": "julia", "file_extension": ".jl", "mimetype": "application/julia", "version": "1.11.2" } } }