{
"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"
}
}
}