use quarto, not Pluto to render pages

quarto/integrals/Project.toml

@@ -0,0 +1,8 @@
+[deps]
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
+Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
+SymPy = "24249f21-da20-56a4-8eb1-6a02cf4ae2e6"
+Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
+UnitfulUS = "7dc9378f-8956-57ef-a780-aa31cc70ff3d"

quarto/integrals/area_between_curves.qmd

@@ -0,0 +1,781 @@
+# Area between two curves
+
+
+```{julia}
+#| echo: false
+
+import Logging
+Logging.disable_logging(Logging.Info) # or e.g. Logging.Info
+Logging.disable_logging(Logging.Warn)
+
+import SymPy
+function Base.show(io::IO, ::MIME"text/html", x::T) where {T <: SymPy.SymbolicObject}
+    println(io, "<span class=\"math-left-align\" style=\"padding-left: 4px; width:0; float:left;\"> ")
+    println(io, "\\[")
+    println(io, sympy.latex(x))
+    println(io, "\\]")
+    println(io, "</span>")
+end
+
+# hack to work around issue
+import Markdown
+import CalculusWithJulia
+function CalculusWithJulia.WeaveSupport.ImageFile(d::Symbol, f::AbstractString, caption; kwargs...)
+    nm = joinpath("..", string(d), f)
+    u = "![$caption]($nm)"
+    Markdown.parse(u)
+end
+
+nothing
+```
+
+This section uses these add-on packages:
+
+
+```{julia}
+using CalculusWithJulia
+using Plots
+using Roots
+using QuadGK
+using SymPy
+```
+
+```{julia}
+#| echo: false
+#| results: "hidden"
+using CalculusWithJulia.WeaveSupport
+
+const frontmatter = (
+        title = "Area between two curves",
+        description = "Calculus with Julia: Area between two curves",
+        tags = ["CalculusWithJulia", "integrals", "area between two curves"],
+);
+nothing
+```
+
+---
+
+
+The definite integral gives the "signed" area between the function $f(x)$ and the $x$-axis over $[a,b]$. Conceptually, this is the area between two curves, $f(x)$ and $g(x)=0$. More generally, this integral:
+
+
+$$
+\int_a^b (f(x) - g(x)) dx
+$$
+
+can be interpreted as the "signed" area between $f(x)$ and $g(x)$ over $[a,b]$. If on this interval $[a,b]$ it is true that $f(x) \geq g(x)$, then this would just be the area, as seen in this figure. The rectangle in the figure has area: $(f(a)-g(a)) \cdot (b-a)$ which could be a term in a left Riemann sum of the integral of $f(x) - g(x)$:
+
+
+```{julia}
+#| hold: true
+#| echo: false
+f1(x) = x^2
+g1(x) = sqrt(x)
+a,b = 1/4, 3/4
+
+xs = range(a, stop=b, length=250)
+ss = vcat(xs, reverse(xs))
+ts = vcat(f1.(xs), g1.(reverse(xs)))
+
+plot(f1, 0, 1, legend=false)
+plot!(g1)
+plot!(ss, ts, fill=(0, :red))
+plot!(xs, f1.(xs), linewidth=5, color=:green)
+plot!(xs, g1.(xs), linewidth=5, color=:green)
+
+
+plot!(xs, f1.(xs), legend=false, linewidth=5, color=:blue)
+plot!(xs, g1.(xs), linewidth=5, color=:blue)
+u,v = .4, .5
+plot!([u,v,v,u,u], [f1(u), f1(u), g1(u), g1(u), f1(u)], color=:black, linewidth=3)
+```
+
+For the figure, we have $f(x) = \sqrt{x}$, $g(x)= x^2$ and $[a,b] = [1/4, 3/4]$. The shaded area is then found by:
+
+
+$$
+\int_{1/4}^{3/4} (x^{1/2} - x^2) dx = (\frac{x^{3/2}}{3/2} - \frac{x^3}{3})\big|_{1/4}^{3/4} = \frac{\sqrt{3}}{4} -\frac{7}{32}.
+$$
+
+#### Examples
+
+
+Find the area bounded by the line $y=2x$ and the curve $y=2 - x^2$.
+
+
+We can plot to see the area in question:
+
+
+```{julia}
+f(x) = 2 - x^2
+g(x) = 2x
+plot(f, -3,3)
+plot!(g)
+```
+
+For this problem we need to identify $a$ and $b$. These are found numerically through:
+
+
+```{julia}
+a,b = find_zeros(x -> f(x) - g(x), -3, 3)
+```
+
+The answer then can be found numerically:
+
+
+```{julia}
+quadgk(x -> f(x) - g(x), a, b)[1]
+```
+
+##### Example
+
+
+Find the integral between $f(x) = \sin(x)$ and $g(x)=\cos(x)$ over $[0,2\pi]$ where $f(x) \geq g(x)$.
+
+
+A plot shows the areas:
+
+
+```{julia}
+𝒇(x) = sin(x)
+𝒈(x) = cos(x)
+plot(𝒇, 0, 2pi)
+plot!(𝒈)
+```
+
+There is a single interval when $f \geq g$ and this can be found algebraically using basic trigonometry, or numerically:
+
+
+```{julia}
+𝒂,𝒃 = find_zeros(x -> 𝒇(x) - 𝒈(x), 0, 2pi)  # pi/4, 5pi/4
+quadgk(x -> 𝒇(x) - 𝒈(x), 𝒂, 𝒃)[1]
+```
+
+##### Example
+
+
+Find the area between $x^n$ and $x^{n+1}$ over $[0,1]$ for $n=1,2,\dots$.
+
+
+We have on this interval $x^n \geq x^{n+1}$, so the integral can be found symbolically through:
+
+
+```{julia}
+@syms x::positive n::positive
+ex = integrate(x^n - x^(n+1), (x, 0, 1))
+together(ex)
+```
+
+Based on this answer, what is the value of this
+
+
+$$
+\frac{1}{2\cdot 3} + \frac{1}{3\cdot 4} + \frac{1}{4\cdot 5} + \cdots?
+$$
+
+This should should be no surprise, given how the areas computed carve up the area under the line $y=x^1$ over $[0,1]$, so the answer should be $1/2$.
+
+
+```{julia}
+p = plot(x, 0, 1, legend=false)
+[plot!(p, x^n, 0, 1) for n in 2:20]
+p
+```
+
+We can check using the `summation` function of `SymPy` which is similar in usage to `integrate`:
+
+
+```{julia}
+summation(1/(n+1)/(n+2), (n, 1, oo))
+```
+
+##### Example
+
+
+Verify [Archimedes'](http://en.wikipedia.org/wiki/The_Quadrature_of_the_Parabola) finding that the area of the parabolic segment is $4/3$rds that of the triangle joining $a$, $(a+b)/2$ and $b$.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+f(x) = 2 - x^2
+a,b = -1, 1/2
+c = (a + b)/2
+xs = range(-sqrt(2), stop=sqrt(2), length=50)
+rxs = range(a, stop=b, length=50)
+rys = map(f, rxs)
+
+
+plot(f, a, b, legend=false, linewidth=3)
+xs = [a,c,b,a]
+plot!(xs, f.(xs), linewidth=3)
+```
+
+For concreteness, let $f(x) = 2-x^2$ and $[a,b] = [-1, 1/2]$, as in the figure. Then the area of the triangle can be computed through:
+
+
+```{julia}
+𝐟(x) = 2 - x^2
+𝐚, 𝐛 = -1, 1/2
+𝐜 = (𝐚 + 𝐛)/2
+
+sac, sab, scb = secant(𝐟, 𝐚, 𝐜), secant(𝐟, 𝐚, 𝐛), secant(𝐟, 𝐜, 𝐛)
+f1(x) = min(sac(x), scb(x))
+f2(x) = sab(x)
+
+A1 = quadgk(x -> f1(x) - f2(x), 𝐚, 𝐛)[1]
+```
+
+As we needed three secant lines, we used the `secant` function from `CalculusWithJulia` to create functions representing each. Once that was done, we used the `max` function to facilitate integrating over the top bounding curve, alternatively, we could break the integral over $[a,c]$ and $[c,b]$.
+
+
+The area of the parabolic segment is more straightforward.
+
+
+```{julia}
+A2 = quadgk(x -> 𝐟(x) - f2(x), 𝐚, 𝐛)[1]
+```
+
+Finally, if Archimedes was right, this relationship should bring about $0$ (or something within round-off error):
+
+
+```{julia}
+A1 * 4/3 - A2
+```
+
+##### Example
+
+
+Find the area bounded by $y=x^4$ and $y=e^x$ when $x^4 \geq e^x$ and $x > 0$.
+
+
+A graph  over $[0,10]$ shows clearly the largest zero, for afterwards the exponential dominates the power.
+
+
+```{julia}
+h1(x) = x^4
+h2(x) = exp(x)
+plot(h1, 0, 10)
+plot!(h2)
+```
+
+There must be another zero, though it is hard to see from the graph over $[0,10]$, as $0^4=0$ and $e^0=1$, so the polynomial must cross below the exponential to the left of $5$. (Otherwise, plotting over $[0,2]$ will clearly reveal the other zero.) We now find these intersection points numerically and then integrate:
+
+
+```{julia}
+#| hold: true
+a,b = find_zeros(x -> h1(x) - h2(x), 0, 10)
+quadgk(x -> h1(x) - h2(x), a, b)[1]
+```
+
+##### Examples
+
+
+The area between $y=\sin(x)$ and $y=m\cdot x$ between $0$ and the first positive intersection depends on $m$ (where $0 \leq m \leq 1$. The extremes are when $m=0$, the area is $2$ and when $m=1$ (the line is tangent at $x=0$), the area is $0$. What is it for other values of $m$? The picture for $m=1/2$ is:
+
+
+```{julia}
+m = 1/2
+plot(sin, 0, pi)
+plot!(x -> m*x)
+```
+
+For a given $m$, the area is found after computing $b$, the intersection point. We express this as a function of $m$ for later reuse:
+
+
+```{julia}
+intersection_point(m) = maximum(find_zeros(x -> sin(x) - m*x, 0, pi))
+a1 = 0
+b1 = intersection_point(m)
+quadgk(x -> sin(x) - m*x, a1, b1)[1]
+```
+
+In general, the area then as a function of `m` is found by substituting `intersection_point(m)` for `b`:
+
+
+```{julia}
+area(m) = quadgk(x -> sin(x) - m*x, 0, intersection_point(m))[1]
+```
+
+A plot shows the relationship:
+
+
+```{julia}
+plot(area, 0, 1)
+```
+
+While here, let's also answer the question of which $m$ gives an area of $1$, or one-half the total? This can be done as follows:
+
+
+```{julia}
+find_zero(m -> area(m) - 1, (0, 1))
+```
+
+(Which is a nice combination of using `find_zeros`, `quadgk` and `find_zero` to answer a problem.)
+
+
+##### Example
+
+
+Find the area bounded by the $x$ axis, the line $x-1$ and the function $\log(x+1)$.
+
+
+A plot shows us the basic area:
+
+
+```{julia}
+j1(x) = log(x+1)
+j2(x) = x - 1
+plot(j1, 0, 3)
+plot!(j2)
+plot!(zero)
+```
+
+The value for "$b$" is found from the intersection point of  $\log(x+1)$ and $x-1$, which is near $2$:
+
+
+```{julia}
+ja = 0
+jb = find_zero(x -> j1(x) - j2(x), 2)
+```
+
+We see that the lower part of the area has a condition: if $x < 1$ then use $0$, otherwise use $g(x)$. We can handle this many different ways:
+
+
+  * break the integral into two pieces and add:
+
+
+```{julia}
+quadgk(x -> j1(x) - zero(x), ja, 1)[1] + quadgk(x -> j1(x) - j2(x), 1, jb)[1]
+```
+
+  * make a new function for the bottom bound:
+
+
+```{julia}
+j3(x) = x < 1 ? 0.0 : j2(x)
+quadgk(x -> j1(x) - j3(x), ja, jb)[1]
+```
+
+  * Turn the picture on its side and integrate in the $y$ variable. To do this, we need to solve for inverse functions:
+
+
+```{julia}
+#| hold: true
+a1=j1(ja)
+b1=j1(jb)
+f1(y)=y+1                # y=x-1, so x=y+1
+g1(y)=exp(y)-1           # y=log(x+1) so e^y = x + 1, x = e^y - 1
+quadgk(y -> f1(y) - g1(y), a1, b1)[1]
+```
+
+:::{.callout-note}
+## Note
+When doing problems by hand this latter style can often reduce the complications, but when approaching the task numerically, the first two styles are generally easier, though computationally more expensive.
+
+:::
+
+#### Integrating in different directions
+
+
+The last example suggested integrating in the $y$ variable. This could have more explanation.
+
+
+It has been noted that different symmetries can aid in computing integrals through their interpretation as areas. For example, if $f(x)$ is odd, then $\int_{-b}^b f(x)dx=0$ and if $f(x)$ is even, $\int_{-b}^b f(x) dx = 2\int_0^b f(x) dx$.
+
+
+Another symmetry of the $x-y$ plane is the reflection through the line $y=x$. This has the effect of taking the graph of $f(x)$ to the graph of $f^{-1}(x)$ and vice versa. Here is an example with $f(x) = x^3$ over $[-1,1]$.
+
+
+```{julia}
+#| hold: true
+f(x) = x^3
+xs = range(-1, stop=1, length=50)
+ys = f.(xs)
+plot(ys, xs)
+```
+
+By switching the order of the `xs` and `ys` we "flip" the graph through the line $x=y$.
+
+
+We can use this symmetry to our advantage. Suppose instead of being given an equation $y=f(x)$, we are given it in "inverse" style: $x = f(y)$, for example suppose we have $x = y^3$. We can plot this as above via:
+
+
+```{julia}
+#| hold: true
+ys = range(-1, stop=1, length=50)
+xs = [y^3 for y in ys]
+plot(xs, ys)
+```
+
+Suppose we wanted the area in the first quadrant between this graph, the $y$ axis and the line $y=1$. What to do? With the problem "flipped" through the $y=x$ line, this would just be $\int_0^1 x^3dx$. Rather than mentally flipping the picture to integrate, instead we can just integrate in the $y$ variable. That is, the area is  $\int_0^1 y^3 dy$. The mental picture for Riemann sums would be have the approximating rectangles laying flat and as a function of $y$, are given a length of $y^3$ and height of "$dy$".
+
+
+---
+
+
+```{julia}
+#| hold: true
+#| echo: false
+f(x) = x^(1/3)
+f⁻¹(x) = x^3
+plot(f, 0, 1, label="f", linewidth=5, color=:blue, aspect_ratio=:equal)
+plot!([0,1,1],[0,0,1], linewidth=1, linestyle=:dash, label="")
+x₀ = 2/3
+Δ = 1/16
+col = RGBA(0,0,1,0.25)
+function box(x,y,Δₓ, Δ, color=col)
+	plot!([x,x+Δₓ, x+Δₓ, x, x], [y,y,y+Δ,y+Δ,y], color=:black, label="")
+	plot!(x:Δₓ:(x+Δₓ), u->y, fillto = u->y+Δ, color=col, label="")
+end
+box(x₀, 0, Δ, f(x₀), col)
+box(x₀+Δ, 0, Δ, f(x₀+Δ), col)
+box(x₀+2Δ, 0, Δ, f(x₀+2Δ), col)
+colᵣ = RGBA(1,0,0,0.25)
+box(f⁻¹(x₀-0Δ), x₀-1Δ, 1 - f⁻¹(x₀-0Δ), Δ, colᵣ)
+box(f⁻¹(x₀-1Δ), x₀-2Δ, 1 - f⁻¹(x₀-1Δ), Δ, colᵣ)
+box(f⁻¹(x₀-2Δ), x₀-3Δ, 1 - f⁻¹(x₀-2Δ), Δ, colᵣ)
+```
+
+The figure above suggests that the area under $f(x)$ over $[a,b]$ could be represented as the area between the curves $f^{-1}(y)$ and $y=b$ from $[f(a), f(b)]$.
+
+
+---
+
+
+For a less trivial problem, consider the area between $x = y^2$ and $x = 2-y$ in the first quadrant.
+
+
+```{julia}
+#| hold: true
+ys = range(0, stop=2, length=50)
+xs = [y^2 for y in ys]
+plot(xs, ys)
+xs = [2-y for y in ys]
+plot!(xs, ys)
+plot!(zero)
+```
+
+We see the bounded area could be described in the "$x$" variable in terms of two integrals, but in the $y$ variable in terms of the difference of two functions with the limits of integration running from $y=0$ to $y=1$. So, this area may be found as follows:
+
+
+```{julia}
+#| hold: true
+f(y) = 2-y
+g(y) = y^2
+a, b = 0, 1
+quadgk(y -> f(y) - g(y), a, b)[1]
+```
+
+## Questions
+
+
+###### Question
+
+
+Find the area enclosed by the curves $y=2-x^2$ and $y=x^2 - 3$.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+f(x) = 2 - x^2
+g(x) = x^2 - 3
+a,b = find_zeros(x -> f(x) - g(x), -10, 10)
+val, _ = quadgk(x -> f(x) - g(x), a, b)
+numericq(val)
+```
+
+###### Question
+
+
+Find the area between $f(x) = \cos(x)$, $g(x) = x$ and the $y$ axis.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+f(x) = cos(x)
+g(x) = x
+a = 0
+b = find_zero(x -> f(x) - g(x), 1)
+val, _ = quadgk(x -> f(x) - g(x), a, b)
+numericq(val)
+```
+
+###### Question
+
+
+Find the area between the line $y=1/2(x+1)$ and half circle $y=\sqrt{1 - x^2}$.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+f(x) = sqrt(1 - x^2)
+g(x) = 1/2 * (x + 1)
+a,b = find_zeros(x -> f(x) - g(x), -1, 1)
+val, _ = quadgk(x -> f(x) - g(x), a, b)
+numericq(val)
+```
+
+###### Question
+
+
+Find the area in the first quadrant between the lines $y=x$, $y=1$, and the curve $y=x^2 + 4$.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+f(x) = x
+g(x) = 1.0
+h(x) = min(f(x), g(x))
+j(x) = x^2 / 4
+a,b = find_zeros(x -> h(x) - j(x), 0, 3)
+val, _ = quadgk(x -> h(x) - j(x), a, b)
+numericq(val)
+```
+
+###### Question
+
+
+Find the area between $y=x^2$ and $y=-x^4$ for $\lvert x \rvert \leq 1$.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+f(x) = x^2
+g(x) = -x^4
+a,b = -1, 1
+val, _ = quadgk(x -> f(x) - g(x), a, b)
+numericq(val)
+```
+
+###### Question
+
+
+Let `f(x) = 1/(sqrt(pi)*gamma(1/2)) * (1 + t^2)^(-1)` and `g(x) = 1/sqrt(2*pi) * exp(-x^2/2)`. These graphs intersect in two points. Find the area bounded by them.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+import SpecialFunctions: gamma
+f(x) = 1/(sqrt(pi)*gamma(1/2)) * (1 + x^2)^(-1)
+g(x) = 1/sqrt(2*pi) * exp(-x^2/2)
+a,b =  find_zeros(x -> f(x) - g(x), -3, 3)
+val, _ = quadgk(x -> f(x) - g(x), a, b)
+numericq(val)
+```
+
+(Where `gamma(1/2)` is a call to the [gamma](http://en.wikipedia.org/wiki/Gamma_function) function.)
+
+
+###### Question
+
+
+Find the area in the first quadrant bounded by the graph of $x = (y-1)^2$, $x=3-y$ and $x=2\sqrt{y}$. (Hint: integrate in the $y$ variable.)
+
+
+```{julia}
+#| hold: true
+#| echo: false
+f(y) = (y-1)^2
+g(y) = 3 - y
+h(y) = 2sqrt(y)
+a = 0
+b = find_zero(y -> f(y) - g(y), 2)
+f1(y) = max(f(y), zero(y))
+g1(y) = min(g(y), h(y))
+val, _ = quadgk(y -> g1(y) - f1(y), a, b)
+numericq(val)
+```
+
+###### Question
+
+
+Find the total area bounded by the lines $x=0$, $x=2$ and the curves $y=x^2$ and $y=x$. This would be $\int_a^b \lvert f(x) - g(x) \rvert dx$.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+f(x) = x^2
+g(x) = x
+a, b = 0, 2
+val, _ = quadgk(x -> abs(f(x) - g(x)), a, b)
+numericq(val)
+```
+
+###### Question
+
+
+Look at the sculpture [Le Tamanoir](https://www.google.com/search?q=Le+Tamanoir+by+Calder.&num=50&tbm=isch&tbo=u&source=univ&sa=X&ved=0ahUKEwiy8eO2tqzVAhVMPz4KHXmgBpgQsAQILQ&biw=1556&bih=878) by Calder. A large scale work. How much does it weigh? Approximately?
+
+
+Let's try to answer that with an educated guess. The right most figure looks to be about 1/5th the total amount. So if we estimate that piece and multiply by 5 we get a good guess. That part looks like an area of metal bounded by two quadratic polynomials. If we compute that area in square inches, then multiply by an assumed thickness of one inch, we have the cubic volume. The density of galvanized steel is 7850 kg/$m^3$ which we convert into pounds/in$^3$ via:
+
+
+```{julia}
+7850 * 2.2 * (1/39.3)^3
+```
+
+The two parabolas, after rotating, might look like the following (with $x$ in inches):
+
+
+$$
+f(x) = x^2/70, \quad g(x) = 35 + x^2/140
+$$
+
+Put this altogether to give an estimated weight in pounds.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+f(x) = x^2/70
+g(x) = 35 + x^2/140
+a,b = find_zeros(x -> f(x) - g(x), -100, 100)
+ar, _ = quadgk(x -> abs(f(x) - g(x)), a, b)
+val = 5 * ar * 7850 * 2.2 * (1/39.3)^3
+numericq(val)
+```
+
+Is the guess that the entire sculpture is more than two tons?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices=["Less than two tons", "More than two tons"]
+answ = 2
+radioq(choices, answ, keep_order=true)
+```
+
+:::{.callout-note}
+## Note
+We used area to estimate weight in this example, but Galileo used weight to estimate area. It is [mentioned](https://www.maa.org/sites/default/files/pdf/cmj_ftp/CMJ/January%202010/3%20Articles/3%20Martin/08-170.pdf) by Martin that in order to estimate the area enclosed by one arch of a cycloid, Galileo cut the arch from from some material and compared the weight to the weight of the generating circle. He concluded the area is close to $3$ times that of the circle, a conjecture proved by Roberval in 1634.
+
+:::
+
+###### Question
+
+
+Formulas from the business world say that revenue is the integral of *marginal revenue* or the additional money from  selling 1 more unit. (This is basically the derivative of profit). Cost is the integral of *marginal cost*, or the cost to produce 1 more. Suppose we have
+
+
+$$
+\text{mr}(x) = 2 - \frac{e^{-x/10}}{1 + e^{-x/10}}, \quad
+\text{mc}(x) = 1 - \frac{1}{2} \cdot \frac{e^{-x/5}}{1 + e^{-x/5}}.
+$$
+
+Find the profit to produce 100 units: $P = \int_0^{100} (\text{mr}(x) - \text{mc}(x)) dx$.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+mr(x) = 2 + exp((-x/10)) / (1 + exp(-x/10))
+mc(x) = 1 + (1/2) * exp(-x/5) / (1 + exp(-x/5))
+a, b = 0, 100
+val, _ = quadgk(x -> mr(x) - mc(x), 0, 100)
+numericq(val)
+```
+
+###### Question
+
+
+Can `SymPy` do what Archimedes did?
+
+
+Consider the following code which sets up the area of an inscribed triangle, `A1`, and the area of a parabolic segment, `A2` for a general parabola:
+
+
+```{julia}
+#| hold: true
+@syms x::real A::real B::real C::real a::real b::real
+c = (a + b) / 2
+f(x) = A*x^2 + B*x + C
+Secant(f, a, b) = f(a) + (f(b)-f(a))/(b-a) * (x - a)
+A1 = integrate(Secant(f, a, c) - Secant(f,a,b), (x,a,c)) + integrate(Secant(f,c,b)-Secant(f,a,b), (x, c, b))
+A2 = integrate(f(x) - Secant(f,a,b), (x, a, b))
+out = 4//3 * A1 - A2
+```
+
+Does `SymPy` get the correct output, $0$, after calling `simplify`?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+yesnoq(true)
+```
+
+###### Question
+
+
+In [Martin](https://www.maa.org/sites/default/files/pdf/cmj_ftp/CMJ/January%202010/3%20Articles/3%20Martin/08-170.pdf) a fascinating history of the cycloid can be read.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+imgfile="figures/cycloid-companion-curve.png"
+caption = """
+Figure from Martin  showing the companion curve to the cycloid. As the generating circle rolls, from ``A`` to ``C``, the original point of contact, ``D``, traces out an arch of the cycloid. The companion curve is that found by congruent line segments. In the figure, when ``D`` was at point ``P`` the line segment ``PQ`` is congruent to ``EF`` (on the original position of the generating circle).
+"""
+ImageFile(:integrals, imgfile, caption)
+```
+
+In particular, it can be read that Roberval proved that the area between the cycloid and its companion curve is half the are of the generating circle. Roberval didn't know integration, so finding the area between two curves required other tricks. One is called "Cavalieri's principle." From the figure above, which of the following would you guess this principle to be:
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices = ["""
+If two regions bounded by parallel lines are such that any parallel between them cuts each region in segments of equal length, then the regions have equal area.
+""",
+           """
+The area of the cycloid is nearly the area of a semi-ellipse with known values, so one can approximate the area of the cycloid with formula for the area of an ellipse
+"""]
+radioq(choices, 1)
+```
+
+Suppose the generating circle has radius $1$, so the area shown is $\pi/2$. The companion curve is then  $1-\cos(\theta)$ (a fact not used by Roberval). The area *under* this curve is then
+
+
+```{julia}
+#| hold: true
+@syms theta
+integrate(1 - cos(theta), (theta, 0, SymPy.PI))
+```
+
+That means the area under **one-half** arch of the cycloid is
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices = ["``\\pi``",
+           "``(3/2)\\cdot \\pi``",
+           "``2\\pi``"
+           ]
+radioq(choices, 2, keep_order=true)
+```
+
+Doubling the answer above gives a value that Galileo had struggled with for many years.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+imgfile="figures/companion-curve-bisects-rectangle.png"
+caption = """
+Roberval, avoiding a trignometric integral, instead used symmetry to show that the area under the companion curve was half the area of the rectangle, which in this figure is ``2\\pi``.
+"""
+ImageFile(:integrals, imgfile, caption)
+```
+

quarto/integrals/center_of_mass.qmd

@@ -0,0 +1,669 @@
+# Center of Mass
+
+
+```{julia}
+#| echo: false
+
+import Logging
+Logging.disable_logging(Logging.Info) # or e.g. Logging.Info
+Logging.disable_logging(Logging.Warn)
+
+import SymPy
+function Base.show(io::IO, ::MIME"text/html", x::T) where {T <: SymPy.SymbolicObject}
+    println(io, "<span class=\"math-left-align\" style=\"padding-left: 4px; width:0; float:left;\"> ")
+    println(io, "\\[")
+    println(io, sympy.latex(x))
+    println(io, "\\]")
+    println(io, "</span>")
+end
+
+# hack to work around issue
+import Markdown
+import CalculusWithJulia
+function CalculusWithJulia.WeaveSupport.ImageFile(d::Symbol, f::AbstractString, caption; kwargs...)
+    nm = joinpath("..", string(d), f)
+    u = "![$caption]($nm)"
+    Markdown.parse(u)
+end
+
+nothing
+```
+
+This section uses these add-on packages:
+
+
+```{julia}
+using CalculusWithJulia
+using Plots
+using Roots
+using QuadGK
+using SymPy
+
+```
+
+```{julia}
+#| echo: false
+#| results: "hidden"
+using CalculusWithJulia.WeaveSupport
+
+const frontmatter = (
+        title = "Center of Mass",
+        description = "Calculus with Julia: Center of Mass",
+        tags = ["CalculusWithJulia", "integrals", "center of mass"],
+);
+nothing
+```
+
+---
+
+
+```{julia}
+#| hold: true
+#| echo: false
+imgfile = "figures/seesaw.png"
+caption = L"""
+
+A silhouette of two children on a seesaw. The seesaw can be balanced
+only if the distance from the central point for each child reflects
+their relative weights, or masses, through the formula $d_1m_1 = d_2
+m_2$. This means if the two children weigh the same the balance will
+tip in favor of the child farther away, and if both are the same
+distance, the balance will tip in favor of the heavier.
+"""
+ImageFile(:integrals, imgfile, caption)
+```
+
+The game of seesaw is one where children earn an early appreciation for the effects of distance and relative weight. For children with equal weights, the seesaw will balance if they sit an equal distance from the center (on opposite sides, of course). However, with unequal weights that isn't the case. If one child weighs twice as much, the other must sit twice as far.
+
+
+The key relationship is that $d_1 m_1 = d_2 m_2$. This come from physics, where the moment about a point is defined by the mass times the distance. This balance relationship says the overall moment balances out.  When this is the case, then the *center of mass* is at the fulcrum point, so there is no impetus to move.
+
+
+The [center](http://en.wikipedia.org/wiki/Center_of_mass) of mass is an old concept that often allows a possibly complicated relationship involving weights to be reduced to a single point. The seesaw is an example: if the center of mass is at the fulcrum the seesaw can balance.
+
+
+In general, we use position of the mass, rather than use distance from some fixed fulcrum. With this, the center of mass  for a finite set of point masses distributed on the real line, is defined by:
+
+
+$$
+\bar{\text{cm}} = \frac{m_1 x_1 + m_2 x_2 + \cdots + m_n x_n}{m_1 + m_2 + \cdots + m_n}.
+$$
+
+Writing $w_i = m_i / (m_1 + m_2 + \cdots + m_n)$, we get the center of mass is just a weighted sum: $w_1 x_1 + \cdots + w_n x_n$, where the $w_i$ are the relative weights.
+
+
+With some rearrangment, we can see that the center of mass satisfies the equation:
+
+
+$$
+w_1 \cdot (x_1 - \bar{\text{cm}}) + w_2 \cdot (x_2 - \bar{\text{cm}}) + \cdots + w_n \cdot (x_n - \bar{\text{cm}}) = 0.
+$$
+
+The center of mass is a balance of the weighted signed distances. This property of the center of mass being a balancing point makes it of intrinsic interest and can be - in the case of sufficient symmetry - easy to find.
+
+
+##### Example
+
+
+A set of weights sits on a dumbbell rack. They are spaced 1 foot apart starting with the 5, then the 10-, 15-, 25-, and 35-pound weights. Where is the center of mass?
+
+
+We begin by letting $m_1=5$, $m_2=10$, $m_3=15$, $m_4=25$ and $m_5=35$. Our positions will be labeled $x_i = i-1$, so the five-pound weight is at position $0$ and the $35$-pound one at $4$. The center of mass is then given by:
+
+
+$$
+\frac{5\cdot 0 + 10\cdot 1 + 15 \cdot 2 + 25 \cdot 3 + 35\cdot 4}{5 + 10 + 15 + 25 + 35} = \frac{255}{90} = 2.833
+$$
+
+If the $25$-pound weight is removed, how does the center of mass shift?
+
+
+We could add the terms again, or just reduce from our sum:
+
+
+$$
+\frac{255 - 25\cdot 3}{90 - 25} = \frac{180}{65} = 2.769...
+$$
+
+The center of mass shifts slightly, but since the removed weight was already close to the center of mass, the movement wasn't much.
+
+
+## Center of mass of figures
+
+
+Consider now a more general problem, the center of mass of a solid figure. We will restrict our attention to figures that can be represented by functions in the $x-y$ plane which are two dimensional. For example, consider the region in the plane bounded by the $x$ axis and the function $1 - \lvert x \rvert$. This is triangle with vertices $(-1,0)$, $(0,1)$, and $(1,0)$.
+
+
+This graph shows that the figure is symmetric:
+
+
+```{julia}
+#| hold: true
+f(x) = 1 - abs(x)
+a, b = -1.5, 1.5
+plot(f, a, b)
+plot!(zero, a, b)
+```
+
+As the center of mass should be a balancing value, we would guess intuitively that the center of mass in the $x$ direction will be $x=0$.
+
+
+But what should the center of mass formula be?
+
+
+As with many formulas that will end up involving a derived integral, we start with a sum approximation. If the region is described as the area under the graph of $f(x)$ between $a$ and $b$, then we can form a Riemann sum approximation, that is a choice of $a = x_0 < x_1 < x_2 \cdots < x_n = b$ and points $c_1$, $\dots$, $c_n$. If all the rectangles are made up of a material of uniform density, say $\rho$, then the mass of each rectangle will be the area times $\rho$, or $\rho f(c_i) \cdot (x_i - x_{i-1})$, for $i = 1, \dots , n$.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+n = 21
+f(x) = 1 - abs(x)
+a, b = -1.5, 1.5
+
+xs = range(-1, stop=1, length=n)
+cs = (xs[1:(end-1)] + xs[2:end]) / 2
+
+p = plot(legend=false);
+
+for i in 1:(n-1)
+   xi, xi1 = xs[i], xs[i+1]
+   plot!(p, [xi, xi1, xi1, xi], [0,0,1,1]*f(cs[i]), linetype=:polygon, color=:red);
+end
+
+for i in 1:(n-1)
+  ci = cs[i]
+  scatter!(p, [ci], [.1], markersize=12*f(cs[i]), color=:orange);
+end
+plot!(p, [-1,1], [0,0])
+
+p
+```
+
+The figure shows the approximating rectangles and circles representing their masses for $n=20$.
+
+
+Generalizing from this figure shows the center of mass for such an approximation will be:
+
+
+$$
+\begin{align*}
+&\frac{\rho f(c_1) (x_1 - x_0) \cdot x_1 + \rho f(c_2) (x_2 - x_1) \cdot x_1 + \cdots + \rho f(c_n) (x_n- x_{n-1}) \cdot x_{n-1}}{\rho f(c_1) (x_1 - x_0)  + \rho f(c_2) (x_2 - x_1)  + \cdots + \rho f(c_n) (x_n- x_{n-1})} \\
+&=\\
+&\quad\frac{f(c_1) (x_1 - x_0) \cdot x_1 + f(c_2) (x_2 - x_1) \cdot x_1 + \cdots + f(c_n) (x_n- x_{n-1}) \cdot x_{n-1}}{f(c_1) (x_1 - x_0)  + f(c_2) (x_2 - x_1)  + \cdots + f(c_n) (x_n- x_{n-1})}.
+\end{align*}
+$$
+
+But the top part is an approximation to the integral $\int_a^b x f(x) dx$ and the bottom part the integral $\int_a^b f(x) dx$. The ratio of these defines the center of mass.
+
+
+> **Center of Mass**: The center of mass (in the $x$ direction) of a region in the $x-y$ plane described by the area under a (positive) function $f(x)$ between $a$ and $b$ is given by
+>
+> $\text{Center of mass} = \text{cm}_x = \frac{\int_a^b xf(x) dx}{\int_a^b f(x) dx}.$
+>
+> For regions described by a more complicated set of equations, the center of mass is found from the same formula where $f(x)$ is the total height in the $x$ direction for a given $x$.
+
+
+
+For the triangular shape, we have by the fact that $f(x) = 1 - \lvert x \rvert$ is an even function that $xf(x)$ will be odd, so the integral around $-1,1$ will be $0$. So the center of mass formula applied to this problem agrees with our expectation.
+
+
+##### Example
+
+
+What about the center of mass of the triangle formed by the line $x=-1$, the $x$ axis and $(1-x)/2$? This too is defined between $a=-1$ and $b=-1$, but the center of mass will be negative, as a graph shows more mass to the left of $0$ than the right:
+
+
+```{julia}
+#| hold: true
+f(x) = (1-x)/2
+plot(f, -1, 1)
+plot!(zero, -1, 1)
+```
+
+The formulas give:
+
+
+$$
+\int_{-1}^1 xf(x) dx = \int_{-1}^1 x\cdot (1-x)/2 = (\frac{x^2}{4} - \frac{x^3}{6})\big|_{-1}^1 = -\frac{1}{3}.
+$$
+
+The bottom integral is just the area (or total mass if the $\rho$ were not canceled) and by geometry is $1/2 (1)(2) = 1$. So $\text{cm}_x = -1/3$.
+
+
+##### Example
+
+
+Find the center of mass formed by the intersection of the parabolas $y=1 - x^2$ and $y=(x-1)^2 - 2$.
+
+
+The center of mass (in the $x$ direction) can be seen to be close to $x=1/2$:
+
+
+```{julia}
+f1(x) = 1 - x^2
+f2(x) = (x-1)^2 -2
+plot(f1, -3, 3)
+plot!(f2, -3, 3)
+```
+
+To find it, we need to find the intersection points, then integrate. We do so numerically.
+
+
+```{julia}
+#| hold: true
+h(x) = f1(x) - f2(x)
+a,b = find_zeros(h, -3, 3)
+top, err = quadgk(x -> x * h(x), a, b)
+bottom, err = quadgk(h, a, b)
+cm = top / bottom
+```
+
+Our guess from the diagram proves correct.
+
+
+:::{.callout-note}
+## Note
+It proves convenient to use the `->` notation for an anonymous function above, as our function `h` is not what is being integrated all the time, but some simple modification. If this isn't palatable, a new function could be defined and passed along to `quadgk`.
+
+:::
+
+##### Example
+
+
+Consider a region bounded by a probability density function. (These functions are non-negative, and integrate to $1$.) The center of mass formula simplifies to $\int xf(x) dx$, as the denominator will be $1$, and the answer is called the *mean*, and often denoted by the Greek letter $\mu$.
+
+
+For the probability density $f(x) = e^{-x}$ for $x\geq 0$ and $0$ otherwise, find the mean.
+
+
+We need to compute $\int_{-\infty}^\infty xf(x) dx$, but in this case since $f$ is $0$ to the left of the origin, we just have:
+
+
+$$
+\mu = \int_0^\infty x e^{-x} dx = -(1+x) \cdot e^{-x} \big|_0^\infty = 1
+$$
+
+For fun, we compare this to the median, which is the value $M$ so that the total area is split in half. That is,  the following formula is satisfied: $\int_0^M f(x) dx = 1/2$. To compute, we have:
+
+
+$$
+\int_0^M e^{-x} dx = -e^{-x} \big|_0^M = 1 - e^{-M}.
+$$
+
+Solving $1/2 = 1 - e^{-M}$ gives $M=\log(2) = 0.69...$, The median is to the left of the mean in this example.
+
+
+:::{.callout-note}
+## Note
+In this example, we used an infinite region, so the idea of "balancing" may be a bit unrealistic, nonetheless, this intuitive interpretation is still a good one to keep this in mind. The point of comparing to the median is that the balancing point is to the right of where the area splits in half. Basically, the center of mass follows in the direction of the area far to the right of the median, as this area is skewed in that direction.
+
+:::
+
+##### Example
+
+
+A figure is formed by transformations of the function $\phi(u) = e^{2(k-1)} - e^{2(k-u)}$, for some fixed $k$, as follows:
+
+
+```{julia}
+k = 3
+phi(u) = exp(2(k-1)) - exp(2(k-u))
+f(u) = max(0, phi(u))
+g(u) = min(f(u+1), f(k))
+
+plot(f, 0, k, legend=false)
+plot!(g, 0, k)
+plot!(zero, 0, k)
+```
+
+(This is basically the graph of $\phi(u)$ and the graph of its shifted value $\phi(u+1)$, only truncated on the top and bottom.)
+
+
+The center of mass of this figure is found with:
+
+
+```{julia}
+h(x) = g(x) - f(x)
+top, _ = quadgk(x -> x*h(x), 0, k)
+bottom, _ = quadgk(h, 0, k)
+top/bottom
+```
+
+This figure has constant slices of length $1$ for fixed values of $y$. If we were to approximate the values with blocks of height $1$, then the center of mass would be to the left of $1$ - for any $k$, but the top most block would have an overhang to the right of $1$ - out to a value of $k$. That is, this figure should balance:
+
+
+```{julia}
+#| echo: false
+u(i) = 1/2*(2k - log(exp(2(k-1)) - i))
+p = plot(legend=false);
+for i in 0:floor(phi(k))
+  x = u(i)
+  plot!(p, [x,x,x-1,x-1], [f(x),f(x)+1, f(x)+1, f(x)], linetype=:polygon, color=:orange);
+end
+xs = range(0, stop=exp(1), length=50)
+plot!(p, f, 0, e, linewidth=5);
+plot!(p, g, 0, 3, linewidth=5)
+p
+```
+
+See this [paper](https://math.dartmouth.edu/~pw/papers/maxover.pdf) and its references for some background on this example and its extensions.
+
+
+### The $y$ direction.
+
+
+We can talk about the center of mass in the $y$ direction too. The approximating picture uses horizontal rectangles - not vertical ones - and if we describe them by $f(y)$, then the corresponding formulas would be
+
+
+> $\text{center of mass} = \text{cm}_y = \frac{\int_a^b y f(y) dy}{\int_a^b f(y) dy}.$
+
+
+
+For example, consider, again, the triangle bounded by the line $x=-1$, the $x$ axis, and the line $y=(1-x)/2$. In terms of describing this in $y$, the function $f(y)=2 -2y$ gives the total length of the horizontal slice (which comes from solving $y=(x-1)/2$for $x$, the general method to find an inverse function, and subtracting $-1$) and the interval is $y=0$ to $y=1$. Thus our center of mass in the $y$ direction will be
+
+
+$$
+\text{cm}_y = \frac{\int_0^1 y (2 - 2y) dy}{\int_0^1 (2 - 2y) dy} = \frac{(2y^2/2 - 2y^3/3)\big|_0^1}{1} = \frac{1}{3}.
+$$
+
+Here the center of mass is below $1/2$ as the bulk of the area is. (The bottom area is just $1$, as known from the area of a triangle.)
+
+
+As seen, the computation of the center of mass in the $y$ direction has an identical formula, though may be more involved if an inverse function must be computed.
+
+
+##### Example
+
+
+More generally, consider a right triangle with vertices $(0,0)$, $(0,a)$, and $(b,0)$. The center of mass of this can be computed with the help of the equation for the line that forms the hypotenuse: $x/b + y/a = 1$. We find the center of mass symbolically in the $y$ variable by solving for $x$ in terms of $y$, then integrating from $0$ to $a$:
+
+
+```{julia}
+@syms a b x y
+eqn = x/b + y/a - 1
+fy = solve(eqn, x)[1]
+integrate(y*fy, (y, 0, a)) / integrate(fy, (y, 0, a))
+```
+
+The answer involves $a$ linearly, but not $b$. If we find the center of mass in $x$, we *could* do something similar:
+
+
+```{julia}
+fx = solve(eqn, y)[1]
+integrate(x*fx, (x, 0, b)) / integrate(fx, (x, 0, b))
+```
+
+But really, we should have just noted that simply by switching the labels $a$ and $b$ in the diagram we could have discovered this formula.
+
+
+:::{.callout-note}
+## Note
+The [centroid](http://en.wikipedia.org/wiki/Centroid) of a region in the plane is just $(\text{cm}_x, \text{cm}_y)$. This last fact says the centroid of the right triangle  is just $(b/3, a/3)$. The centroid can be found by other geometric means. The link shows the plumb line method. For triangles, the centroid is also the intersection point of the medians, the lines that connect a vertex with its opposite midpoint.
+
+:::
+
+##### Example
+
+
+Compute the $x$ and $y$ values of the center of mass of the half circle described by the area below the function $f(x) = \sqrt{1 - x^2}$ and above the $x$-axis.
+
+
+A plot shows the value of cm$_x$ will be $0$ by symmetry:
+
+
+```{julia}
+#| hold: true
+f(x) = sqrt(1 - x^2)
+plot(f, -1, 1)
+```
+
+($f(x)$ is even, so $xf(x)$ will be odd.)
+
+
+However, the value for cm$_y$ will - like the last problem - be around $1/3$. The exact value is compute using slices in the $y$ direction. Solving for $x$ in $y=\sqrt{1-x^2}$, or $x = \pm \sqrt{1-y^2}$, if $f(y) = 2\sqrt{1 - y^2}$. The value is then:
+
+
+$$
+\text{cm}_y = \frac{\int_{0}^1 y 2 \sqrt{1 - y^2}dy}{\int_{0}^1 2\sqrt{1-y^2}} =
+\frac{-(1-x^2)^{3/2}/3\big|_0^1}{\pi/2} = \frac{1}{3}.
+$$
+
+In fact it is exactly $1/3$. The top calculation is done by $u$-substitution, the bottom by using the area formula for a half circle, $\pi r^2/2$.
+
+
+##### Example
+
+
+A disc of radius $2$ is centered at the origin, as a disc of radius $1$ is bored out between $y=0$ and $y=1$. Find the resulting center of mass.
+
+
+A picture shows that this could be complicated, especially for $y > 0$, as we need to describe the length of the red lines below for $-2 < y < 2$:
+
+
+```{julia}
+#| hold: true
+#| echo: false
+a,b = 0, 2pi
+ts = range(a, stop=b, length=50)
+p = plot(t -> 2cos(t), t->2sin(t), a, b, legend=false, aspect_ratio=:equal);
+plot!(p, cos.(ts), 1 .+ sin.(ts), linetype=:polygon, color=:red);
+plot!(p, [-sqrt(3), sqrt(3)], [-1,-1], color=:orange);
+plot!(p, [-sqrt(3), -1],      [1,1],   color=:orange);
+plot!(p, [sqrt(3), 1],        [1,1],   color=:orange);
+p
+
+```
+
+We can see that cm$_x = 0$, by symmetry, but to compute cm$_y$ we need to find $f(y)$, which will depend on the value of $y$ between $-2$ and $2$. The outer circle is $x^2 + y^2 = 4$, the inner circle $x^2 + (y-1)^2 = 1$.  When $y < 0$, $f(y)$ is the distance across the outer circle or, $2\sqrt{4 - y^2}$. When $y \geq 0$, $f(y)$ is *twice* the distance from the bigger circle to the smaller, of $2(\sqrt{4 - y^2} - \sqrt{1 - (y-1)^2})$.
+
+
+We use this to compute:
+
+
+```{julia}
+#| hold: true
+f(y) = y < 0 ? 2*sqrt(4 - y^2) : 2* (sqrt(4 - y^2)- sqrt(1 - (y-1)^2))
+top, _ = quadgk( y -> y * f(y), -2, 2)
+bottom, _ = quadgk( f, -2, 2)
+top/bottom
+```
+
+The nice answer of $-1/3$ makes us think there may be a different way to visualize this. Were we to rearrange the top integral, we could write it as $\int_{-2}^2 y 2 \sqrt{4 -y^2}dy - \int_0^2 2y\sqrt{1 - (y-1)^2}dy$. Call this $A - B$.  The left term, $A$, is part of the center of mass formula for the big circle (which is this value divided by $M=4\pi$), and the right term, $B$, is part of the center of mass formula for the (drilled out) smaller circle (which is this value divided by $m=\pi$. These values are weighted according to $(AM - Bm)/(M-m)$. In this case $A=0$, $B=1$ and $M=4m$, so the answer is $-1/3$.
+
+
+## Questions
+
+
+###### Question
+
+
+Find the center of mass in the $x$ variable for the region bounded by parabola $x=4 - y^2$  and the $y$ axis.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+f(x) = sqrt(4 - x)
+a, b = 0, 4
+top, _ = quadgk(x -> x*f(x), a,b)
+bottom, _ = quadgk(f, a,b)
+val = top/bottom
+numericq(val)
+```
+
+###### Question
+
+
+Find the center of mass in the $x$ variable of the region in the first and fourth quadrants bounded by the ellipse $(x/2)^2 + (y/3)^2 = 1$.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+f(x) = 3 * sqrt(1 - (x/2)^2)
+a, b= 0, 2
+
+top, _ = quadgk(x -> x*f(x), a,b)
+bottom, _ = quadgk(f, a,b)
+val = top/bottom
+numericq(val)
+```
+
+###### Question
+
+
+Find the center of mass of the region in the first quadrant bounded by the function $f(x) = x^3(1-x)^4$.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+f(x) = x^3 * (1-x)^4
+a, b= 0, 1
+
+top, _ = quadgk(x -> x*f(x), a,b)
+bottom, _ = quadgk(f, a,b)
+val = top/bottom
+numericq(val)
+```
+
+###### Question
+
+
+Let $k$ and $\lambda$ be parameters in $(0, \infty)$. The [Weibull](http://en.wikipedia.org/wiki/Weibull_distribution) density is a probability density on $[0, \infty)$ (meaning it is $0$ when $x < 0$ satisfying:
+
+
+$$
+f(x) = \frac{k}{\lambda}\left(\frac{x}{\lambda}\right)^{k-1} \exp(-(\frac{x}{\lambda})^k)
+$$
+
+For $k=2$ and $\lambda = 2$, compute the mean. (The center of mass, assuming the total area is $1$.)
+
+
+```{julia}
+#| hold: true
+#| echo: false
+k, lambda = 2, 2
+f(x) = (k/lambda) * (x/lambda)^(k-1) * exp(-(x/lambda)^k)
+a, b = 0, Inf
+
+top, _ = quadgk(x -> x*f(x), a,b)
+bottom, _ = quadgk(f, a,b)
+val = top/bottom
+numericq(val)
+```
+
+###### Question
+
+
+The [logistic](http://en.wikipedia.org/wiki/Logistic_distribution) density depends on two parameters $m$ and $s$ and is given by:
+
+
+$$
+f(x) = \frac{1}{4s} \text{sech}(\frac{x-\mu}{2s})^2, \quad -\infty < x < \infty.
+$$
+
+(Where $\text{sech}$ is the hyperbolic secant, implemented in `julia` through `sech`.)
+
+
+For $m=2$ and $s=4$ compute the mean, or center of mass, of this density.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+m, s = 2, 4
+f(x) = 1/(4s) * sech((x-m)/2s)^2
+a,b = -Inf, Inf
+
+top, _ = quadgk(x -> x*f(x), a,b)
+bottom, _ = quadgk(f, a,b)
+val = top/bottom
+numericq(val)
+```
+
+###### Question
+
+
+A region is formed by intersecting the area bounded by the circle $x^2 + y^2 = 1$ that lies above the line $y=3/4$. Find the center of mass in the $y$ direction (that of the $x$ direction is $0$ by symmetry).
+
+
+```{julia}
+#| hold: true
+#| echo: false
+f(y) = 2*sqrt(1 - y^2)
+a, b = 3/4, 1
+
+
+top, _ = quadgk(x -> x*f(x), a,b)
+bottom, _ = quadgk(f, a,b)
+val = top/bottom
+numericq(val)
+```
+
+###### Question
+
+
+Find the center of mass in the $y$ direction of the area bounded by the cosine curve and the $x$ axis between $-\pi/2$ and $\pi/2$.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+f(y) = 2*acos(y)
+a, b= 0, 1
+
+top, _ = quadgk(x -> x*f(x), a,b)
+bottom, _ = quadgk(f, a,b)
+val = top/bottom
+numericq(val)
+```
+
+###### Question
+
+
+A penny, nickel, dime and quarter are stacked so that their right most edges align and are centered so that the center of mass in the $y$ direction is $0$. Find the center of mass in the $x$ direction.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+ds = [0.75, 0.835, 0.705, 0.955]
+rs = ds/2
+xs = rs[4] .- rs
+ts = range(0,stop=2pi, length=50)
+p = plot(legend=false, aspect_ratio=:equal);
+for i in 1:4
+  plot!(p, xs[i] .+ rs[i]*cos.(ts), rs[i]*sin.(ts));
+end
+
+p
+```
+
+You will need some specifications, such as these from the [US Mint](http://www.usmint.gov/about_the_mint/?action=coin_specifications)
+
+
+```{eval=false}
+         diameter(in)   weight(gms)
+penny     0.750           2.500
+nickel    0.835           5.000
+dime      0.705           2.268
+quarter   0.955           5.670
+
+```
+
+(Hint: Though this could be done with integration, it is easier to treat each coin as a single point (its centroid) with the given mass and then apply the formula for sums.)
+
+
+```{julia}
+#| hold: true
+#| echo: false
+ds = [0.75, 0.835, 0.705, 0.955]
+ms = [2.5, 5, 2.268, 5.670]
+rs = ds/2
+xs = rs[4] .- rs
+val = sum(ms .* xs) / sum(ms)
+numericq(val)
+```
+

quarto/integrals/improper_integrals.qmd

@@ -0,0 +1,616 @@
+# Improper Integrals
+
+
+```{julia}
+#| echo: false
+
+import Logging
+Logging.disable_logging(Logging.Info) # or e.g. Logging.Info
+Logging.disable_logging(Logging.Warn)
+
+import SymPy
+function Base.show(io::IO, ::MIME"text/html", x::T) where {T <: SymPy.SymbolicObject}
+    println(io, "<span class=\"math-left-align\" style=\"padding-left: 4px; width:0; float:left;\"> ")
+    println(io, "\\[")
+    println(io, sympy.latex(x))
+    println(io, "\\]")
+    println(io, "</span>")
+end
+
+# hack to work around issue
+import Markdown
+import CalculusWithJulia
+function CalculusWithJulia.WeaveSupport.ImageFile(d::Symbol, f::AbstractString, caption; kwargs...)
+    nm = joinpath("..", string(d), f)
+    u = "![$caption]($nm)"
+    Markdown.parse(u)
+end
+
+nothing
+```
+
+This section uses these add-on packages:
+
+
+```{julia}
+using CalculusWithJulia
+using Plots
+using SymPy
+using QuadGK
+```
+
+```{julia}
+#| echo: false
+#| results: "hidden"
+using CalculusWithJulia.WeaveSupport
+
+const frontmatter = (
+        title = "Improper Integrals",
+        description = "Calculus with Julia: Improper Integrals",
+        tags = ["CalculusWithJulia", "integrals", "improper integrals"],
+);
+fig_size=(800, 600)
+
+nothing
+```
+
+---
+
+
+A function $f(x)$ is Riemann integrable over an interval $[a,b]$ if some limit involving Riemann sums exists. This limit will fail to exist if $f(x) = \infty$ in $[a,b]$. As well, the Riemann sum idea is undefined if either $a$ or $b$ (or both) are infinite, so the limit won't exist in this case.
+
+
+To define integrals  with either functions having singularities or infinite domains, the idea of an improper integral is introduced with definitions to handle the two cases above.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+#| cache: true
+### {{{sqrt_graph}}}
+
+function make_sqrt_x_graph(n)
+
+    b = 1
+    a = 1/2^n
+    xs = range(1/2^8, stop=b, length=250)
+    x1s = range(a, stop=b, length=50)
+    f(x) = 1/sqrt(x)
+    val = N(integrate(f, 1/2^n, b))
+    title = "area under f over [1/$(2^n), $b] is $(rpad(round(val, digits=2), 4))"
+
+    plt = plot(f, range(a, stop=b, length=251), xlim=(0,b), ylim=(0, 15), legend=false, size=fig_size, title=title)
+    plot!(plt,  [b, a, x1s...], [0, 0, map(f, x1s)...], linetype=:polygon, color=:orange)
+
+    plt
+
+
+end
+caption = L"""
+
+Area under $1/\sqrt{x}$ over $[a,b]$ increases as $a$ gets closer to $0$. Will it grow unbounded or have a limit?
+
+"""
+n = 10
+anim = @animate for i=1:n
+    make_sqrt_x_graph(i)
+end
+
+imgfile = tempname() * ".gif"
+gif(anim, imgfile, fps = 1)
+
+ImageFile(imgfile, caption)
+```
+
+## Infinite domains
+
+
+Let $f(x)$ be a reasonable function, so reasonable that for any $a < b$ the function is Riemann integrable, meaning $\int_a^b f(x)dx$ exists.
+
+
+What needs to be the case so that we can discuss the integral over the entire real number line?
+
+
+Clearly something. The function $f(x) = 1$ is reasonable by the idea above. Clearly the integral over and $[a,b]$ is just $b-a$, but the limit over an unbounded domain would be $\infty$. Even though limits of infinity can be of interest in some cases, not so here. What will ensure that the area is finite over an infinite region?
+
+
+Or is that even the right question. Now consider $f(x) = \sin(\pi x)$. Over every interval of the type $[-2n, 2n]$ the area is $0$, and over any interval, $[a,b]$ the area never gets bigger than $2$. But still this function does not have a well defined area on an infinite domain.
+
+
+The right question involves a limit. Fix a finite $a$. We define the definite integral over $[a,\infty)$ to be
+
+
+$$
+\int_a^\infty f(x) dx = \lim_{M \rightarrow \infty} \int_a^M f(x) dx,
+$$
+
+when the limit exists. Similarly, we define the definite integral over $(-\infty, a]$ through
+
+
+$$
+\int_{-\infty}^a f(x) dx  = \lim_{M \rightarrow -\infty} \int_M^a f(x) dx.
+$$
+
+For the interval $(-\infty, \infty)$ we have need *both* these limits to exist, and then:
+
+
+$$
+\int_{-\infty}^\infty f(x) dx  = \lim_{M \rightarrow -\infty} \int_M^a f(x) dx + \lim_{M \rightarrow \infty} \int_a^M f(x) dx.
+$$
+
+:::{.callout-note}
+## Note
+When the integral exists, it is said to *converge*. If it doesn't exist, it is said to *diverge*.
+
+:::
+
+##### Examples
+
+
+  * The function $f(x) = 1/x^2$ is integrable over $[1, \infty)$, as this limit exists:
+
+
+$$
+\lim_{M \rightarrow \infty} \int_1^M \frac{1}{x^2}dx = \lim_{M \rightarrow \infty} -\frac{1}{x}\big|_1^M
+= \lim_{M \rightarrow \infty} 1 - \frac{1}{M} = 1.
+$$
+
+  * The function $f(x) = 1/x^{1/2}$ is not integrable over $[1, \infty)$, as this limit fails to exist:
+
+
+$$
+\lim_{M \rightarrow \infty} \int_1^M \frac{1}{x^{1/2}}dx = \lim_{M \rightarrow \infty} \frac{x^{1/2}}{1/2}\big|_1^M
+= \lim_{M \rightarrow \infty} 2\sqrt{M} - 2 = \infty.
+$$
+
+The limit is infinite, so does not exist except in an extended sense.
+
+
+  * The function $x^n e^{-x}$ for $n = 1, 2, \dots$ is integrable over $[0,\infty)$.
+
+
+Before showing this, we recall the fundamental theorem of calculus. The limit existing is the same as saying the limit of $F(M) - F(a)$ exists for an antiderivative of $f(x)$.
+
+
+For this particular problem, it can be shown by integration by parts that for positive, integer values of $n$ that an antiderivative exists of the form $F(x) = p(x)e^{-x}$, where $p(x)$ is a polynomial of degree $n$. But we've seen that for any $n>0$, $\lim_{x \rightarrow \infty} x^n e^{-x} = 0$, so the same is true for any polynomial. So, $\lim_{M \rightarrow \infty} F(M) - F(1) = -F(1)$.
+
+
+  * The function $e^x$ is integrable over $(-\infty, a]$ but not
+
+
+$$
+[a, \infty)
+$$
+
+for any finite $a$. This is because, $F(M) = e^x$ and this has a limit as $x$ goes to $-\infty$, but not $\infty$.
+
+
+  * Let $f(x) = x e^{-x^2}$. This function has an integral over $[0, \infty)$ and more generally $(-\infty, \infty)$. To see, we note that as it is an odd function, the area from $0$ to $M$ is the opposite sign of that from $-M$ to $0$. So $\lim_{M \rightarrow \infty} (F(M) - F(0)) = \lim_{M \rightarrow -\infty} (F(0) - (-F(\lvert M\lvert)))$. We only then need to investigate the one limit. But we can see by substitution with $u=x^2$, that an antiderivative is $F(x) = (-1/2) \cdot e^{-x^2}$. Clearly, $\lim_{M \rightarrow \infty}F(M) = 0$, so the answer is well defined, and the area from $0$ to $\infty$ is just $e/2$. From $-\infty$ to $0$ it is $-e/2$ and the total area is $0$, as the two sides "cancel" out.
+  * Let $f(x) = \sin(x)$. Even though $\lim_{M \rightarrow \infty} (F(M) - F(-M) ) = 0$, this function is not integrable. The fact is we need *both* the limit $F(M)$ and $F(-M)$ to exist as $M$ goes to $\infty$. In this case, even though the area cancels if $\infty$ is approached at the same rate, this isn't sufficient to guarantee the two limits exists independently.
+
+
+  * Will the function $f(x) = 1/(x\cdot(\log(x))^2)$ have an integral over $[e, \infty)$?
+
+
+We first find an antiderivative using the $u$-substitution $u(x) = \log(x)$:
+
+
+$$
+\int_e^M \frac{e}{x \log(x)^{2}} dx
+= \int_{\log(e)}^{\log(M)} \frac{1}{u^{2}} du
+= \frac{-1}{u} \big|_{1}^{\log(M)}
+= \frac{-1}{\log(M)} - \frac{-1}{1}
+= 1 - \frac{1}{M}.
+$$
+
+As $M$ goes to $\infty$, this will converge to $1$.
+
+
+  * The sinc function $f(x) = \sin(\pi x)/(\pi x)$ does not have a nice antiderivative.  Seeing if the limit exists is a bit of a problem. However, this function is important enough that there is a built-in function, `Si`, that computes $\int_0^x \sin(u)/u\cdot du$. This function can be used through `sympy.Si(...)`:
+
+
+```{julia}
+@syms M
+limit(sympy.Si(M), M => oo)
+```
+
+### Numeric integration
+
+
+The `quadgk` function (available through `QuadGK`) is able to accept `Inf` and `-Inf` as endpoints of the interval. For example, this will integrate $e^{-x^2/2}$ over the real line:
+
+
+```{julia}
+f(x) = exp(-x^2/2)
+quadgk(f, -Inf, Inf)
+```
+
+(If may not be obvious, but this is $\sqrt{2\pi}$.)
+
+
+## Singularities
+
+
+Suppose $\lim_{x \rightarrow c}f(x) = \infty$ or $-\infty$. Then a Riemann sum that contains an interval including $c$ will not be finite if the point chosen in the interval is $c$. Though we could choose another point, this is not enough as the definition must hold for any choice of the $c_i$.
+
+
+However, if $c$ is isolated, we can get close to $c$ and see how the area changes.
+
+
+Suppose $a < c$, we define $\int_a^c f(x) dx = \lim_{M \rightarrow c-} \int_a^c f(x) dx$. If this limit exists, the definite integral with $c$ is well defined. Similarly, the integral from $c$ to $b$, where $b > c$, can be defined by a right limit going to $c$. The integral from $a$ to $b$ will exist if both the limits are finite.
+
+
+##### Examples
+
+
+  * Consider the example of the initial illustration, $f(x) = 1/\sqrt{x}$ at $0$. Here $f(0)= \infty$, so the usual notion of a limit won't apply to $\int_0^1 f(x) dx$. However,
+
+
+$$
+\lim_{M \rightarrow 0+} \int_M^1 \frac{1}{\sqrt{x}} dx
+= \lim_{M \rightarrow 0+} \frac{\sqrt{x}}{1/2} \big|_M^1
+= \lim_{M \rightarrow 0+} 2(1) - 2\sqrt{M} = 2.
+$$
+
+:::{.callout-note}
+## Note
+The cases $f(x) = x^{-n}$ for $n > 0$ are tricky to keep straight. For $n > 1$, the functions can be integrated over $[1,\infty)$, but not $(0,1]$. For $0 < n < 1$, the functions can be integrated over $(0,1]$ but not $[1, \infty)$.
+
+:::
+
+  * Now consider $f(x) = 1/x$. Is this integral  $\int_0^1 1/x \cdot dx$ defined? It will be *if* this limit exists:
+
+
+$$
+\lim_{M \rightarrow 0+} \int_M^1 \frac{1}{x} dx
+= \lim_{M \rightarrow 0+} \log(x) \big|_M^1
+= \lim_{M \rightarrow 0+} \log(1) - \log(M) = \infty.
+$$
+
+As the limit does not exist, the function is not integrable around $0$.
+
+
+  * `SymPy` may give answers which do not coincide with our definitions, as it uses complex numbers as a default assumption. In this case it returns the proper answer when integrated from $0$ to $1$ and `NaN` for an integral over $(-1,1)$:
+
+
+```{julia}
+@syms x
+integrate(1/x, (x, 0, 1)), integrate(1/x, (x, -1, 1))
+```
+
+  * Suppose you know $\int_1^\infty x^2 f(x) dx$ exists. Does this imply $\int_0^1 f(1/x) dx$ exists?
+
+
+We need to consider the limit of $\int_M^1 f(1/x) dx$. We try the $u$-substitution $u(x) = 1/x$. This gives $du = -(1/x^2)dx = -u^2 dx$. So, the substitution becomes:
+
+
+$$
+\int_M^1 f(1/x) dx = \int_{1/M}^{1/1} f(u) (-u^2) du = \int_1^{1/M} u^2 f(u) du.
+$$
+
+But the limit as $M \rightarrow 0$ of $1/M$ is the same going to $\infty$, so the right side will converge by the assumption. Thus we get $f(1/x)$ is integrable over $(0,1]$.
+
+
+### Numeric integration
+
+
+So far our use of the `quadgk` function specified the region to integrate via `a`, `b`, as in `quadgk(f, a, b)`. In fact, it can specify values in between for which the function should not be sampled. For example, were we to integrate $1/\sqrt{\lvert x\rvert}$ over $[-1,1]$, we would want to avoid $0$ as a point to sample. Here is how:
+
+
+```{julia}
+#| hold: true
+f(x) = 1 / sqrt(abs(x))
+quadgk(f, -1, 0, 1)
+```
+
+Just trying `quadgk(f, -1, 1)` leads to a `DomainError`, as `0` will be one of the points sampled. The general call is like `quadgk(f, a, b, c, d,...)` which integrates over $(a,b)$ and $(b,c)$ and $(c,d)$, $\dots$. The algorithm is not supposed to evaluate the function at the endpoints of the intervals.
+
+
+## Probability applications
+
+
+A probability density is a function $f(x) \geq 0$ which is integrable on $(-\infty, \infty)$ and for which $\int_{-\infty}^\infty f(x) dx =1$. The cumulative distribution function is defined by $F(x)=\int_{-\infty}^x f(u) du$.
+
+
+Probability densities are good example of using improper integrals.
+
+
+  * Show that $f(x) = (1/\pi) (1/(1 + x^2))$ is a probability density function.
+
+
+We need to show that the integral exists and is $1$. For this, we use the fact that $(1/\pi) \cdot \tan^{-1}(x)$ is an antiderivative. Then we have:
+
+
+$$
+\lim_{M \rightarrow \infty} F(M) = (1/\pi) \cdot \pi/2
+$$
+
+and as $\tan^{-1}(x)$ is odd, we must have $F(-\infty) = \lim_{M \rightarrow -\infty} f(M) = -(1/\pi) \cdot \pi/2$. All told, $F(\infty) - F(-\infty) = 1/2 - (-1/2) = 1$.
+
+
+  * Show that $f(x) = 1/(b-a)$ for $a \leq x \leq b$ and $0$ otherwise is a probability density.
+
+
+The integral for $-\infty$ to $a$ of $f(x)$ is just an integral of the constant $0$, so will be $0$. (This is the only constant with finite area over an infinite domain.) Similarly, the integral from $b$ to $\infty$ will be $0$. This means:
+
+
+$$
+\int_{-\infty}^\infty f(x) dx = \int_a^b \frac{1}{b-a} dx = 1.
+$$
+
+(One might also comment that $f$ is Riemann integrable on any $[0,M]$ despite being discontinuous at $a$ and $b$.)
+
+
+  * Show that if $f(x)$ is a probability density then so is $f(x-c)$ for any $c$.
+
+
+We have by the $u$-substitution
+
+
+$$
+\int_{-\infty}^\infty f(x-c)dx = \int_{u(-\infty)}^{u(\infty)} f(u) du = \int_{-\infty}^\infty f(u) du = 1.
+$$
+
+The key is that we can use the regular $u$-substitution formula provided $\lim_{M \rightarrow \infty} u(M) = u(\infty)$ is defined. (The *informal* notation $u(\infty)$ is defined by that limit.)
+
+
+  * If $f(x)$ is a probability density, then so is $(1/h) f((x-c)/h)$ for any $c, h > 0$.
+
+
+Again, by a $u$ substitution with, now, $u(x) = (x-c)/h$, we have $du = (1/h) \cdot dx$ and the result follows just as before:
+
+
+$$
+\int_{-\infty}^\infty \frac{1}{h}f(\frac{x-c}{h})dx = \int_{u(-\infty)}^{u(\infty)} f(u) du = \int_{-\infty}^\infty f(u) du = 1.
+$$
+
+  * If $F(x) = 1 - e^{-x}$, for $x \geq 0$, and $0$ otherwise,  find $f(x)$.
+
+
+We want to just say $F'(x)= e^{-x}$ so $f(x) = e^{-x}$. But some care is needed. First, that isn't right. The derivative for $x<0$ of $F(x)$ is $0$, so $f(x) = 0$ if $x < 0$. What about for $x>0$? The derivative is $e^{-x}$, but is that the right answer? $F(x) = \int_{-\infty}^x f(u) du$, so we have to at least discuss if the $-\infty$ affects things. In this case, and in general the answer is *no*. For any $x$ we can find $M < x$ so that we have $F(x) = \int_{-\infty}^M f(u) du + \int_M^x f(u) du$. The first part is a constant, so will have derivative $0$, the second will have derivative $f(x)$, if the derivative exists (and it will exist at $x$ if the derivative is continuous in a neighborhood of $x$).
+
+
+Finally, at $x=0$ we have an issue, as $F'(0)$ does not exist. The left limit of the secant line approximation is $0$, the right limit of the secant line approximation is $1$. So, we can take $f(x) = e^{-x}$ for $x > 0$ and $0$ otherwise, noting that redefining $f(x)$ at a point will not effect the integral as long as the point is finite.
+
+
+## Questions
+
+
+###### Question
+
+
+Is $f(x) = 1/x^{100}$ integrable around $0$?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+yesnoq("no")
+```
+
+###### Question
+
+
+Is $f(x) = 1/x^{1/3}$ integrable around $0$?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+yesnoq("yes")
+```
+
+###### Question
+
+
+Is $f(x) = x\cdot\log(x)$ integrable on $[1,\infty)$?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+yesnoq("no")
+```
+
+###### Question
+
+
+Is $f(x) = \log(x)/ x$ integrable on $[1,\infty)$?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+yesnoq("no")
+```
+
+###### Question
+
+
+Is $f(x) = \log(x)$ integrable on $[1,\infty)$?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+yesnoq("no")
+```
+
+###### Question
+
+
+Compute the integral $\int_0^\infty 1/(1+x^2) dx$.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+f(x) = 1/(1+x^2)
+a, b= 0, Inf
+val, _ = quadgk(f, a, b)
+numericq(val)
+```
+
+###### Question
+
+
+Compute the the integral $\int_1^\infty \log(x)/x^2 dx$.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+f(x) =log(x)/x^2
+a, b= 1, Inf
+val, _ = quadgk(f, a, b)
+numericq(val)
+```
+
+###### Question
+
+
+Compute the integral $\int_0^2 (x-1)^{2/3} dx$.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+f(x) = cbrt((x-1)^2)
+val, _ = quadgk(f , 0, 1, 2)
+numericq(val)
+```
+
+###### Question
+
+
+From the relationship that if $0 \leq f(x) \leq g(x)$ then $\int_a^b f(x) dx \leq \int_a^b g(x) dx$ it can be deduced that
+
+
+  * if $\int_a^\infty f(x) dx$ diverges, then so does $\int_a^\infty g(x) dx$.
+  * if $\int_a^\infty g(x) dx$ converges, then so does $\int_a^\infty f(x) dx$.
+
+
+Let $f(x) = \lvert \sin(x)/x^2 \rvert$.
+
+
+What can you say about $\int_1^\infty f(x) dx$, as $f(x) \leq 1/x^2$ on $[1, \infty)$?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices =[
+"It is convergent",
+"It is divergent",
+"Can't say"]
+answ = 1
+radioq(choices, answ, keep_order=true)
+```
+
+---
+
+
+Let $f(x) = \lvert \sin(x) \rvert / x$.
+
+
+What can you say about $\int_1^\infty f(x) dx$, as $f(x) \leq 1/x$ on $[1, \infty)$?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices =[
+"It is convergent",
+"It is divergent",
+"Can't say"]
+answ = 3
+radioq(choices, answ, keep_order=true)
+```
+
+---
+
+
+Let $f(x) = 1/\sqrt{x^2 - 1}$. What can you say about $\int_1^\infty f(x) dx$, as $f(x) \geq 1/x$ on $[1, \infty)$?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices =[
+"It is convergent",
+"It is divergent",
+"Can't say"]
+answ = 2
+radioq(choices, answ, keep_order=true)
+```
+
+---
+
+
+Let $f(x) = 1 + 4x^2$. What can you say about $\int_1^\infty f(x) dx$, as $f(x) \leq 1/x^2$ on $[1, \infty)$?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices =[
+"It is convergent",
+"It is divergent",
+"Can't say"]
+answ = 2
+radioq(choices, answ, keep_order=true)
+```
+
+---
+
+
+Let $f(x) = \lvert \sin(x)^{10}\rvert/e^x$. What can you say about $\int_1^\infty f(x) dx$, as $f(x) \leq e^{-x}$ on $[1, \infty)$?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices =[
+"It is convergent",
+"It is divergent",
+"Can't say"]
+answ = 1
+radioq(choices, answ, keep_order=true)
+```
+
+###### Question
+
+
+The difference between "blowing up" at $0$ versus being integrable at $\infty$ can be seen to be related through the $u$-substitution $u=1/x$. With this $u$-substitution, what becomes of $\int_0^1 x^{-2/3} dx$?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices = [
+"``\\int_1^\\infty u^{2/3}/u^2 \\cdot du``",
+"``\\int_0^1 u^{2/3} \\cdot du``",
+"``\\int_0^\\infty 1/u \\cdot du``"
+]
+answ = 1
+radioq(choices, answ)
+```
+
+###### Question
+
+
+The antiderivative of $f(x) = 1/\pi \cdot 1/\sqrt{x(1-x)}$ is $F(x)=(2/\pi)\cdot \sin^{-1}(\sqrt{x})$.
+
+
+Find $\int_0^1 f(x) dx$.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+f(x) = 1/pi * 1/sqrt(x*(1-x))
+a, b= 0, 1
+val, _ = quadgk(f, a, b)
+numericq(val)
+```
+

quarto/integrals/integration_by_parts.qmd

@@ -0,0 +1,732 @@
+# Integration By Parts
+
+
+```{julia}
+#| echo: false
+
+import Logging
+Logging.disable_logging(Logging.Info) # or e.g. Logging.Info
+Logging.disable_logging(Logging.Warn)
+
+import SymPy
+function Base.show(io::IO, ::MIME"text/html", x::T) where {T <: SymPy.SymbolicObject}
+    println(io, "<span class=\"math-left-align\" style=\"padding-left: 4px; width:0; float:left;\"> ")
+    println(io, "\\[")
+    println(io, sympy.latex(x))
+    println(io, "\\]")
+    println(io, "</span>")
+end
+
+# hack to work around issue
+import Markdown
+import CalculusWithJulia
+function CalculusWithJulia.WeaveSupport.ImageFile(d::Symbol, f::AbstractString, caption; kwargs...)
+    nm = joinpath("..", string(d), f)
+    u = "![$caption]($nm)"
+    Markdown.parse(u)
+end
+
+nothing
+```
+
+This section uses these add-on packages:
+
+
+```{julia}
+using CalculusWithJulia
+using Plots
+using SymPy
+```
+
+```{julia}
+#| echo: false
+#| results: "hidden"
+using CalculusWithJulia.WeaveSupport
+using QuadGK
+
+frontmatter = (
+    title = "Integration By Parts",
+    description = "Calculus with Julia: Integration By Parts",
+    tags = ["CalculusWithJulia", "integrals", "integration by parts"],
+);
+
+nothing
+```
+
+---
+
+
+So far we have seen that the *derivative* rules lead to *integration rules*. In particular:
+
+
+  * The sum rule $[au(x) + bv(x)]' = au'(x) + bv'(x)$ gives rise to an integration rule: $\int (au(x) + bv(x))dx = a\int u(x)dx + b\int v(x))dx$. (That is, the linearity of the derivative means the integral has linearity.)
+  * The chain rule $[f(g(x))]' = f'(g(x)) g'(x)$ gives $\int_a^b f(g(x))g'(x)dx=\int_{g(a)}^{g(b)}f(x)dx$. That is, substitution reverses the chain rule.
+
+
+Now we turn our attention to the implications of the *product rule*: $[uv]' = u'v + uv'$. The resulting technique is called integration by parts.
+
+
+The following illustrates integration by parts of the integral $(uv)'$ over $[a,b]$ [original](http://en.wikipedia.org/wiki/Integration_by_parts#Visualization).
+
+
+```{julia}
+#| echo: false
+let
+    ## parts picture
+    u(x) = sin(x*pi/2)
+    v(x) = x
+    xs = range(0, stop=1, length=50)
+    a,b = 1/4, 3/4
+    p = plot(u, v, 0, 1, legend=false)
+    plot!(p, zero, 0, 1)
+    scatter!(p, [u(a), u(b)], [v(a), v(b)], color=:orange, markersize=5)
+
+    plot!(p, [u(a),u(a),0, 0, u(b),u(b),u(a)],
+          [0, v(a), v(a), v(b), v(b), 0, 0],
+          linetype=:polygon, fillcolor=:orange, alpha=0.25)
+    annotate!(p, [(0.65, .25, "A"), (0.4, .55, "B")])
+    annotate!(p, [(u(a),v(a) + .08, "(u(a),v(a))"), (u(b),v(b)+.08, "(u(b),v(b))")])
+end
+```
+
+The figure is a parametric plot of $(u,v)$ with the points $(u(a), v(a))$ and $(u(b), v(b))$ marked. The difference $u(b)v(b) - u(a)v(a) = u(x)v(x) \mid_a^b$ is shaded.  This area breaks into two pieces, $A$ and $B$, partitioned by the curve. If $u$ is increasing and the curve is parameterized by $t \rightarrow u^{-1}(t)$, then $A=\int_{u^{-1}(a)}^{u^{-1}(b)} v(u^{-1}(t))dt$. A $u$-substitution with $t = u(x)$ changes this into the integral $\int_a^b v(x) u'(x) dx$. Similarly, for increasing $v$, it can be seen that $B=\int_a^b u(x) v'(x) dx$. This suggests a relationship between the integral of $u v'$, the integral of $u' v$ and the value $u(b)v(b) - u(a)v(a)$.
+
+
+In terms of formulas, by the fundamental theorem of calculus:
+
+
+$$
+u(x)\cdot v(x)\big|_a^b = \int_a^b [u(x) v(x)]' dx = \int_a^b u'(x) \cdot v(x) dx + \int_a^b u(x) \cdot v'(x) dx.
+$$
+
+This is re-expressed as
+
+
+$$
+\int_a^b u(x) \cdot v'(x) dx = u(x) \cdot v(x)\big|_a^b -  \int_a^b v(x) \cdot u'(x) dx,
+$$
+
+Or, more informally, as $\int udv = uv - \int v du$.
+
+
+This can sometimes be confusingly written as:
+
+
+$$
+\int f(x) g'(x) dx = f(x)g(x) - \int f'(x) g(x) dx.
+$$
+
+(The confusion coming from the fact that the indefinite integrals are only defined up to a constant.)
+
+
+How does this help? It allows us to differentiate parts of an integral in hopes it makes the result easier to integrate.
+
+
+An illustration can clarify.
+
+
+Consider the integral $\int_0^\pi x\sin(x) dx$. If we let $u=x$ and $dv=\sin(x) dx$, then $du = 1dx$ and $v=-\cos(x)$. The above then says:
+
+
+$$
+\begin{align*}
+\int_0^\pi x\sin(x) dx &= \int_0^\pi u dv\\
+&= uv\big|_0^\pi - \int_0^\pi v du\\
+&= x \cdot (-\cos(x)) \big|_0^\pi - \int_0^\pi (-\cos(x)) dx\\
+&= \pi (-\cos(\pi)) - 0(-\cos(0)) + \int_0^\pi \cos(x) dx\\
+&= \pi + \sin(x)\big|_0^\pi\\
+&= \pi.
+\end{align*}
+$$
+
+The technique means one part is differentiated and one part integrated. The art is to break the integrand up into a piece that gets easier through differentiation and a piece that doesn't get much harder through integration.
+
+
+#### Examples
+
+
+Consider $\int_1^2 x \log(x) dx$. We might try differentiating the $\log(x)$ term, so we set
+
+
+$$
+u=\log(x) \text{ and  } dv=xdx
+$$
+
+Then we get
+
+
+$$
+du = \frac{1}{x} dx \text{ and } v = \frac{x^2}{2}.
+$$
+
+Putting together gives:
+
+
+$$
+\begin{align*}
+\int_1^2 x \log(x) dx
+&= (\log(x) \cdot \frac{x^2}{2}) \big|_1^2 - \int_1^2 \frac{x^2}{2} \frac{1}{x} dx\\
+&= (2\log(2) - 0) - (\frac{x^2}{4})\big|_1^2\\
+&= 2\log(2) - (1 - \frac{1}{4}) \\
+&= 2\log(2) - \frac{3}{4}.
+\end{align*}
+$$
+
+##### Example
+
+
+This related problem, $\int \log(x) dx$, uses the same idea, though perhaps harder to see at first glance, as setting `dv=dx` is almost too simple to try:
+
+
+$$
+\begin{align*}
+u &= \log(x) & dv &= dx\\
+du &= \frac{1}{x}dx & v &= x
+\end{align*}
+$$
+
+$$
+\begin{align*}
+\int \log(x) dx
+&= \int u dv\\
+&= uv - \int v du\\
+&= (\log(x) \cdot x) - \int x \cdot \frac{1}{x} dx\\
+&= x \log(x) - \int dx\\
+&= x \log(x) - x
+\end{align*}
+$$
+
+Were this a definite integral problem, we would have written:
+
+
+$$
+\int_a^b \log(x) dx = (x\log(x))\big|_a^b - \int_a^b dx = (x\log(x) - x)\big|_a^b.
+$$
+
+##### Example
+
+
+Sometimes integration by parts is used two or more times. Here we let $u=x^2$ and $dv = e^x dx$:
+
+
+$$
+\int_a^b x^2 e^x dx = (x^2 \cdot e^x)\big|_a^b - \int_a^b 2x e^x dx.
+$$
+
+But we can do $\int_a^b x e^xdx$ the same way:
+
+
+$$
+\int_a^b x e^x = (x\cdot e^x)\big|_a^b - \int_a^b 1 \cdot e^xdx = (xe^x - e^x)\big|_a^b.
+$$
+
+Combining gives the answer:
+
+
+$$
+\int_a^b x^2 e^x dx
+= (x^2 \cdot e^x)\big|_a^b - 2( (xe^x - e^x)\big|_a^b ) =
+e^x(x^2  - 2x - 1) \big|_a^b.
+$$
+
+In fact, it isn't hard to see that an integral of $x^m e^x$, $m$ a positive integer, can be handled in this manner. For example, when $m=10$, `SymPy` gives:
+
+
+```{julia}
+@syms 𝒙
+integrate(𝒙^10 * exp(𝒙), 𝒙)
+```
+
+The general answer is $\int x^n e^xdx = p(x) e^x$, where $p(x)$ is a polynomial of degree $n$.
+
+
+##### Example
+
+
+The same technique is attempted for  this integral, but ends differently. First in the following we let $u=\sin(x)$ and $dv=e^x dx$:
+
+
+$$
+\int e^x \sin(x)dx = \sin(x) e^x - \int \cos(x) e^x dx.
+$$
+
+Now we let $u = \cos(x)$ and again $dv=e^x dx$:
+
+
+$$
+\int e^x \sin(x)dx = \sin(x) e^x - \int \cos(x) e^x dx = \sin(x)e^x - \cos(x)e^x - \int (-\sin(x))e^x dx.
+$$
+
+But simplifying this gives:
+
+
+$$
+\int e^x \sin(x)dx = - \int e^x \sin(x)dx +  e^x(\sin(x) - \cos(x)).
+$$
+
+Solving for the "unknown" $\int e^x \sin(x) dx$ gives:
+
+
+$$
+\int e^x \sin(x) dx = \frac{1}{2} e^x (\sin(x) - \cos(x)).
+$$
+
+##### Example
+
+
+Positive integer powers of trigonometric functions can be addressed by this technique. Consider $\int \cos(x)^n dx$. We let $u=\cos(x)^{n-1}$ and $dv=\cos(x) dx$. Then $du = (n-1)\cos(x)^{n-2}(-\sin(x))dx$ and $v=\sin(x)$. So,
+
+
+$$
+\begin{align*}
+\int \cos(x)^n dx &= \cos(x)^{n-1} \cdot (\sin(x)) - \int (\sin(x)) ((n-1)\sin(x) \cos(x)^{n-2}) dx \\
+&= \sin(x) \cos(x)^{n-1} + (n-1)\int \sin^2(x) \cos(x)^{n-1} dx\\
+&= \sin(x) \cos(x)^{n-1} + (n-1)\int (1 - \cos(x)^2) \cos(x)^{n-2} dx\\
+&= \sin(x) \cos(x)^{n-1} + (n-1)\int \cos(x)^{n-2}dx - (n-1)\int \cos(x)^n dx.
+\end{align*}
+$$
+
+We can then solve for the unknown ($\int \cos(x)^{n}dx$) to get this *reduction formula*:
+
+
+$$
+\int \cos(x)^n dx = \frac{1}{n}\sin(x) \cos(x)^{n-1} + \frac{n-1}{n}\int \cos(x)^{n-2}dx.
+$$
+
+This is called a reduction formula as it reduces the problem from an integral with a power of $n$ to one with a power of $n - 2$, so could be repeated until the remaining indefinite integral required knowing either $\int \cos(x) dx$ (which is $-\sin(x)$) or $\int \cos(x)^2 dx$, which by a double angle formula application, is $x/2 - \sin(2x)/4$.
+
+
+`SymPy` is quite able to do this repeated bookkeeping. For example with $n=10$:
+
+
+```{julia}
+integrate(cos(𝒙)^10, 𝒙)
+```
+
+##### Example
+
+
+The visual interpretation of integration by parts breaks area into two pieces, the one labeled "B" looks like it would be labeled "A" for an inverse function for $f$. Indeed, integration by parts gives a means to possibly find antiderivatives for inverse functions.
+
+
+Let $uv = x f^{-1}(x)$. Then we have $[uv]' = u'v + uv' = f^{-1}(x) + x [f^{-1}(x)]'$. So, up to a constant $uv = \int [uv]'dx = \int f^{-1}(x) + \int  x [f^{-1}(x)]'$. Re-expressing gives:
+
+
+$$
+\begin{align*}
+\int f^{-1}(x) dx
+&= xf^{-1}(x) - \int x [f^{-1}(x)]' dx\\
+&= xf^{-1}(x) - \int f(u) du.\\
+\end{align*}
+$$
+
+The last line follows from the $u$-substitution: $u=f^{-1}(x)$ for then $du = [f^{-1}(x)]' dx$ and $x=f(u)$.
+
+
+We use this to find an antiderivative for $\sin^{-1}(x)$:
+
+
+$$
+\begin{align*}
+\int \sin^{-1}(x) dx &= x \sin^{-1}(x) - \int \sin(u) du \\
+&=  x \sin^{-1}(x) + \cos(u) \\
+&= x \sin^{-1}(x) + \cos(\sin^{-1}(x)).
+\end{align*}
+$$
+
+Using right triangles to simplify, the last value $\cos(\sin^{-1}(x))$ can otherwise be written as $\sqrt{1 - x^2}$.
+
+
+##### Example
+
+
+The [trapezoid](http://en.wikipedia.org/wiki/Trapezoidal_rule) rule is an approximation to the definite integral like a Riemann sum, only instead of approximating the area above $[x_i, x_i + h]$ by a rectangle with height $f(c_i)$ (for some $c_i$), it uses a trapezoid formed by the left and right endpoints. That is, this area is used in the estimation: $(1/2)\cdot (f(x_i) + f(x_i+h)) \cdot h$.
+
+
+Even though we suggest just using `quadgk` for numeric integration, estimating the error in this approximation is still of some theoretical interest.
+
+
+Recall, just using *either* $x_i$ or $x_{i-1}$ for $c_i$ gives an error that is "like" $1/n$, as $n$ gets large, though the exact rate depends on the function and the length of the interval.
+
+
+This [proof](http://www.math.ucsd.edu/~ebender/20B/77_Trap.pdf) for the error estimate is involved, but is reproduced here, as it nicely integrates many of the theoretical concepts of integration discussed so far.
+
+
+First, for convenience, we consider the interval $x_i$ to $x_i+h$. The actual answer over this is just $\int_{x_i}^{x_i+h}f(x) dx$. By a $u$-substitution with $u=x-x_i$ this becomes $\int_0^h f(t + x_i) dt$. For analyzing this we integrate once by parts using $u=f(t+x_i)$ and $dv=dt$. But instead of letting $v=t$, we choose to add - as is our prerogative - a constant of integration $A$, so $v=t+A$:
+
+
+$$
+\begin{align*}
+\int_0^h f(t + x_i) dt &= uv \big|_0^h - \int_0^h v du\\
+&= f(t+x_i)(t+A)\big|_0^h - \int_0^h (t + A) f'(t + x_i) dt.
+\end{align*}
+$$
+
+We choose $A$ to be $-h/2$, any constant is possible, for then the term $f(t+x_i)(t+A)\big|_0^h$ becomes $(1/2)(f(x_i+h) + f(x_i)) \cdot h$, or the trapezoid approximation. This means, the error over this interval - actual minus estimate - satisfies:
+
+
+$$
+\text{error}_i = \int_{x_i}^{x_i+h}f(x) dx - \frac{f(x_i+h) -f(x_i)}{2} \cdot h  = - \int_0^h (t + A) f'(t + x_i) dt.
+$$
+
+For this, we *again* integrate by parts with
+
+
+$$
+\begin{align*}
+u  &= f'(t + x_i)  & dv &= (t + A)dt\\
+du &= f''(t + x_i) & v  &= \frac{(t + A)^2}{2} + B
+\end{align*}
+$$
+
+Again we added a constant of integration, $B$, to  $v$. The  error becomes:
+
+
+$$
+\text{error}_i = -(\frac{(t+A)^2}{2} + B)f'(t+x_i)\big|_0^h + \int_0^h (\frac{(t+A)^2}{2} + B) \cdot f''(t+x_i) dt.
+$$
+
+With $A=-h/2$, $B$ is chosen so $(t+A)^2/2 + B = 0$, or  $B=-h^2/8$. The error becomes
+
+
+$$
+\text{error}_i = \int_0^h \left(\frac{(t-h/2)^2}{2} - \frac{h^2}{8}\right) \cdot f''(t + x_i) dt.
+$$
+
+Now, we assume the $\lvert f''(t)\rvert$ is bounded by $K$ for any $a \leq t \leq b$. This will be true, for example, if the second derivative is assumed to exist and be continuous.  Using this fact about definite integrals $\lvert \int_a^b g dx\rvert \leq \int_a^b \lvert g \rvert dx$ we have:
+
+
+$$
+\lvert \text{error}_i  \rvert \leq K \int_0^h \lvert (\frac{(t-h/2)^2}{2} - \frac{h^2}{8}) \rvert dt.
+$$
+
+But what is the function in the integrand? Clearly it is a quadratic in $t$. Expanding gives $1/2 \cdot (t^2 - ht)$.  This is negative over $[0,h]$ (and $0$ at these endpoints, so the integral above is just:
+
+
+$$
+\frac{1}{2}\int_0^h (ht - t^2)dt = \frac{1}{2} (\frac{ht^2}{2} - \frac{t^3}{3})\big|_0^h = \frac{h^3}{12}
+$$
+
+This gives the bound: $\vert \text{error}_i \rvert \leq K h^3/12$. The *total* error may be less, but is not more than the value found by adding up the error over each of the $n$ intervals. As our bound does not depend on the $i$, we have this sum satisfies:
+
+
+$$
+\lvert \text{error}\rvert \leq n \cdot \frac{Kh^3}{12} = \frac{K(b-a)^3}{12}\frac{1}{n^2}.
+$$
+
+So the error is like $1/n^2$, in contrast to the $1/n$ error of the Riemann sums. One way to see this, for the Riemann sum it takes twice as many terms to half an error estimate, but for the trapezoid rule only $\sqrt{2}$ as many, and for Simpson's rule, only $2^{1/4}$ as many.
+
+
+## Area related to parameterized curves
+
+
+The figure introduced to motivate the integration by parts formula also suggests that areas described parametrically (by a pair of functions $x=u(t), y=v(t)$ for $a \le t \le b$) can have their area computed.
+
+
+When $u(t)$ is strictly *increasing*, and hence having an inverse function, then re-parameterizing by $\phi(t) = u^{-1}(t)$ gives a $x=u(u^{-1}(t))=t, y=v(u^{-1}(t))$ and integrating this gives the area by $A=\int_a^b v(t) u'(t) dt$
+
+
+However, the correct answer requires understanding a minus sign.  Consider the area enclosed by $x(t) = \cos(t), y(t) = \sin(t)$:
+
+
+```{julia}
+#| echo: false
+let
+    r(t) = [cos(t), sin(t)]
+    p=plot_parametric(0..2pi, r, aspect_ratio=:equal, legend=false)
+    for t ∈ (pi/4, 3pi/4, 5pi/4, 7pi/4)
+	    quiver!(unzip([r(t)])..., quiver=Tuple(unzip([0.1*r'(t)])))
+    end
+    ti, tj = pi/3, pi/3+0.1
+    plot!([cos(tj), cos(ti), cos(ti), cos(tj), cos(tj)], [0,0,sin(tj), sin(tj),0])
+    quiver!([0],[0], quiver=Tuple(unzip([r(ti)])))
+    quiver!([0],[0], quiver=Tuple(unzip([r(tj)])))
+    p
+end
+```
+
+We added a rectangle for a Riemann sum for $t_i = \pi/3$ and $t_{i+1} = \pi/3 + \pi/8$. The height of this rectangle if $y(t_i)$, the base is of length $x(t_i) - x(t_{i+1})$ *given* the orientation of how the circular curve is parameterized (counter clockwise here).
+
+
+Taking this Riemann sum approach, we can approximate the area under the curve parameterized by $(u(t), v(t))$ over the time range $[t_i, t_{i+1}]$ as a rectangle with height $y(t_i)$ and base $x(t_{i}) - x(t_{i+1})$. Then we get, as expected:
+
+
+$$
+\begin{align*}
+A &\approx \sum_i y(t_i) \cdot (x(t_{i}) - x(t_{i+1}))\\
+  &= - \sum_i y(t_i) \cdot (x(t_{i+1}) - x(t_{i}))\\
+  &= - \sum_i y(t_i) \cdot \frac{x(t_{i+1}) - x(t_i)}{t_{i+1}-t_i} \cdot (t_{i+1}-t_i)\\
+  &\approx -\int_a^b y(t) x'(t) dt.
+\end{align*}
+$$
+
+So with a counterclockwise rotation, the actual answer for the area includes a minus sign. If the area is traced out in a *clockwise* manner, there is no minus sign.
+
+
+This is a case of [Green's Theorem](https://en.wikipedia.org/wiki/Green%27s_theorem#Area_calculation) to be taken up in [Green's Theorem, Stokes' Theorem, and the Divergence Theorem](file:///Users/verzani/julia/CalculusWithJulia/html/integral_vector_calculus/stokes_theorem.html).
+
+
+##### Example
+
+
+Apply the formula to a parameterized circle to ensure, the signed area is properly computed. If we use $x(t) = r\cos(t)$ and $y(t) = r\sin(t)$ then we have the motion is counterclockwise:
+
+
+```{julia}
+#| hold: true
+@syms 𝒓 t
+𝒙 = 𝒓 * cos(t)
+𝒚 = 𝒓 * sin(t)
+-integrate(𝒚 * diff(𝒙, t), (t, 0, 2PI))
+```
+
+We see the expected answer for the area of a circle.
+
+
+##### Example
+
+
+Apply the formula to find the area under one arch of a cycloid, parameterized by $x(t) = t - \sin(t), y(t) = 1 - \cos(t)$.
+
+
+Working symbolically, we have one arch given by the following described in a *clockwise* manner, so we use $\int y(t) x'(t) dt$:
+
+
+```{julia}
+#| hold: true
+@syms t
+𝒙 = t - sin(t)
+𝒚 = 1 - cos(t)
+integrate(𝒚 * diff(𝒙, t), (t, 0, 2PI))
+```
+
+([Galileo](https://mathshistory.st-andrews.ac.uk/Curves/Cycloid/) was thwarted in finding this answer exactly and resorted to constructing one from metal to *estimate* the value.)
+
+
+##### Example
+
+
+Consider the example $x(t) = \cos(t) + t\sin(t), y(t) = \sin(t) - t\cos(t)$ for $0 \leq t \leq 2\pi$.
+
+
+```{julia}
+#| echo: false
+let
+    x(t) = cos(t) + t*sin(t)
+    y(t) = sin(t) - t*cos(t)
+    ts = range(0,2pi, length=100)
+    plot(x.(ts), y.(ts))
+end
+```
+
+How much area is enclosed by this curve and the $x$ axis? The area is described in a counterclockwise manner, so we have:
+
+
+```{julia}
+#| hold: true
+let
+    x(t) = cos(t) + t*sin(t)
+    y(t) = sin(t) - t*cos(t)
+    yx′(t) = -y(t) * x'(t)  # yx\prime[tab]
+    quadgk(yx′, 0, 2pi)
+end
+```
+
+This particular problem could also have been done symbolically, but many curves will need to have a numeric approximation used.
+
+
+## Questions
+
+
+###### Question
+
+
+In the integral of $\int \log(x) dx$ we let $u=\log(x)$ and $dv=dx$. What are $du$ and $v$?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices = [
+"``du=1/x dx \\quad v = x``",
+"``du=x\\log(x) dx\\quad v = 1``",
+"``du=1/x dx\\quad v = x^2/2``"]
+answ = 1
+radioq(choices, answ)
+```
+
+###### Question
+
+
+In the integral $\int \sec(x)^3 dx$ we let $u=\sec(x)$ and $dv = \sec(x)^2 dx$. What are $du$ and $v$?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices = [
+"``du=\\sec(x)\\tan(x)dx \\quad v=\\tan(x)``",
+"``du=\\csc(x) dx \\quad v=\\sec(x)^3 / 3``",
+"``du=\\tan(x)  dx \\quad v=\\sec(x)\\tan(x)``"
+]
+answ = 1
+radioq(choices, answ)
+```
+
+###### Question
+
+
+In the integral $\int e^{-x} \cos(x)dx$ we let $u=e^{-x}$ and $dv=\cos(x) dx$. What are $du$ and $v$?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices = [
+"``du=-e^{-x} dx \\quad v=\\sin(x)``",
+"``du=-e^{-x} dx \\quad v=-\\sin(x)``",
+"``du=\\sin(x)dx \\quad v=-e^{-x}``"
+]
+answ = 1
+radioq(choices, answ)
+```
+
+###### Question
+
+
+Find the value of $\int_1^4 x \log(x) dx$. You can integrate by parts.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+f(x) = x*log(x)
+a,b = 1,4
+val,err = quadgk(f, a, b)
+numericq(val)
+```
+
+###### Question
+
+
+Find the value of $\int_0^{\pi/2} x\cos(2x) dx$. You can integrate by parts.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+f(x) = x*cos(2x)
+a,b = 0, pi/2
+val,err = quadgk(f, a, b)
+numericq(val)
+```
+
+###### Question
+
+
+Find the value of $\int_1^e (\log(x))^2 dx$. You can integrate by parts.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+f(x) = log(x)^2
+a,b = 1,exp(1)
+val,err = quadgk(f, a, b)
+numericq(val)
+```
+
+###### Question
+
+
+Integration by parts can be used to provide "reduction" formulas, where an antiderivative is written in terms of another antiderivative with a lower power. Which is the proper reduction formula for $\int (\log(x))^n dx$?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices = [
+"``x(\\log(x))^n - n \\int (\\log(x))^{n-1} dx``",
+"``\\int (\\log(x))^{n+1}/(n+1) dx``",
+"``x(\\log(x))^n - \\int (\\log(x))^{n-1} dx``"
+]
+answ = 1
+radioq(choices, answ)
+```
+
+###### Question
+
+
+The [Wikipedia](http://en.wikipedia.org/wiki/Integration_by_parts) page has a rule of thumb with an acronym LIATE to indicate what is a good candidate to be "$u$": **L**og function, **I**nverse functions, **A**lgebraic functions ($x^n$), **T**rigonmetric functions, and **E**xponential functions.
+
+
+Consider the integral $\int x \cos(x) dx$. Which letter should be tried first?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices = ["L", "I", "A", "T", "E"]
+answ = 3
+radioq(choices, answ, keep_order=true)
+```
+
+---
+
+
+Consider the integral $\int x^2\log(x) dx$. Which letter should be tried first?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices = ["L", "I", "A", "T", "E"]
+answ = 1
+radioq(choices, answ, keep_order=true)
+```
+
+---
+
+
+Consider the integral $\int x^2 \sin^{-1}(x) dx$. Which letter should be tried first?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices = ["L", "I", "A", "T", "E"]
+answ = 2
+radioq(choices, answ, keep_order=true)
+```
+
+---
+
+
+Consider the integral $\int e^x \sin(x) dx$. Which letter should be tried first?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices = ["L", "I", "A", "T", "E"]
+answ = 4
+radioq(choices, answ, keep_order=true)
+```
+
+###### Question
+
+
+Find an antiderivative for $\cos^{-1}(x)$ using the integration by parts formula.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices = [
+"``x\\cos^{-1}(x)-\\sqrt{1 - x^2}``",
+"``x^2/2 \\cos^{-1}(x) - x\\sqrt{1-x^2}/4 - \\cos^{-1}(x)/4``",
+"``-\\sin^{-1}(x)``"]
+answ = 1
+radioq(choices, answ)
+```
+

quarto/integrals/mean_value_theorem.qmd

@@ -0,0 +1,420 @@
+# Mean value theorem for integrals
+
+
+```{julia}
+#| echo: false
+
+import Logging
+Logging.disable_logging(Logging.Info) # or e.g. Logging.Info
+Logging.disable_logging(Logging.Warn)
+
+import SymPy
+function Base.show(io::IO, ::MIME"text/html", x::T) where {T <: SymPy.SymbolicObject}
+    println(io, "<span class=\"math-left-align\" style=\"padding-left: 4px; width:0; float:left;\"> ")
+    println(io, "\\[")
+    println(io, sympy.latex(x))
+    println(io, "\\]")
+    println(io, "</span>")
+end
+
+# hack to work around issue
+import Markdown
+import CalculusWithJulia
+function CalculusWithJulia.WeaveSupport.ImageFile(d::Symbol, f::AbstractString, caption; kwargs...)
+    nm = joinpath("..", string(d), f)
+    u = "![$caption]($nm)"
+    Markdown.parse(u)
+end
+
+nothing
+```
+
+This section uses these add-on packages:
+
+
+```{julia}
+using CalculusWithJulia
+using Plots
+using QuadGK
+```
+
+```{julia}
+#| echo: false
+#| results: "hidden"
+using CalculusWithJulia.WeaveSupport
+
+const frontmatter = (
+        title = "Mean value theorem for integrals",
+        description = "Calculus with Julia: Mean value theorem for integrals",
+        tags = ["CalculusWithJulia", "integrals", "mean value theorem for integrals"],
+);
+nothing
+```
+
+---
+
+
+## Average value of a function
+
+
+Let $f(x)$ be a continuous function over the interval $[a,b]$ with $a < b$.
+
+
+The average value of $f$ over $[a,b]$ is defined by:
+
+
+$$
+\frac{1}{b-a} \int_a^b f(x) dx.
+$$
+
+If $f$ is a constant, this is just the contant value, as would be expected. If $f$ is *piecewise* linear, then this is the weighted average of these constants.
+
+
+#### Examples
+
+
+##### Example: average velocity
+
+
+The average velocity between times $a < b$, is simply the change in position during the time interval divided by the change in time. In notation, this would be $(x(b) - x(a)) / (b-a)$. If $v(t) = x'(t)$ is the velocity, then by the second part of the fundamental theorem of calculus, we have, in agreement with the definition above, that:
+
+
+$$
+\text{average velocity} = \frac{x(b) - x(a)}{b-a} = \frac{1}{b-a} \int_a^b v(t) dt.
+$$
+
+The average speed is the change in *total* distance over time, which is given by
+
+
+$$
+\text{average speed} =  \frac{1}{b-a} \int_a^b \lvert v(t)\rvert dt.
+$$
+
+Let $\bar{v}$ be the average velocity. Then we have $\bar{v} \cdot(b-a) = x(b) - x(a)$, or the change in position can be written as a constant ($\bar{v}$) times the time, as though we had a constant velocity.  This is an old intuition. [Bressoud](http://www.math.harvard.edu/~knill/teaching/math1a_2011/exhibits/bressoud/) comments on the special case known to scholars at Merton College around $1350$ that the distance traveled by an object under uniformly increasing velocity starting at $v_0$ and ending at $v_t$ is equal to the distance traveled by an object with constant velocity of $(v_0 + v_t)/2$.
+
+
+##### Example
+
+
+What is the average value of $f(x)=\sin(x)$ over $[0, \pi]$?
+
+
+$$
+\text{average} = \frac{1}{\pi-0} \int_0^\pi \sin(x) dx = \frac{1}{\pi} (-\cos(x)) \big|_0^\pi = \frac{2}{\pi}
+$$
+
+Visually, we have:
+
+
+```{julia}
+plot(sin, 0, pi)
+plot!(x -> 2/pi)
+```
+
+##### Example
+
+
+What is the average value of the function $f$ which is $1$ between $[0,3]$, $2$ between $(3,5]$ and $1$ between $(5,6]$?
+
+
+Though not continuous, $f(x)$ is integrable as it contains only jumps. The integral from $[0,6]$ can be computed with geometry: $3\cdot 3 + 2 \cdot 2 + 1 \cdot 1 = 14$. The average then is $14/(6-0) = 7/3$.
+
+
+##### Example
+
+
+What is the average value of the function $e^{-x}$ between $0$ and $\log(2)$?
+
+
+$$
+\begin{align*}
+\text{average} = \frac{1}{\log(2) - 0} \int_0^{\log(2)} e^{-x} dx\\
+&= \frac{1}{\log(2)} (-e^{-x}) \big|_0^{\log(2)}\\
+&= -\frac{1}{\log(2)} (\frac{1}{2} - 1)\\
+&= \frac{1}{2\log(2)}.
+\end{align*}
+$$
+
+Visualizing, we have
+
+
+```{julia}
+plot(x -> exp(-x), 0, log(2))
+plot!(x -> 1/(2*log(2)))
+```
+
+## The mean value theorem for integrals
+
+
+If $f(x)$ is assumed integrable, the average value of $f(x)$ is defined, as above. Re-expressing gives that there exists a $K$ with
+
+
+$$
+K \cdot (b-a) = \int_a^b f(x) dx.
+$$
+
+When we assume that $f(x)$ is continuous, we can describe $K$ as a value in the range of $f$:
+
+
+> **The mean value theorem for integrals**: Let $f(x)$ be a continuous function on $[a,b]$ with $a < b$. Then there exists $c$ with $a \leq c \leq b$ with
+>
+> $f(c) \cdot (b-a) = \int_a^b f(x) dx.$`
+
+
+
+The proof comes from the intermediate value theorem and the extreme value theorem. Since $f$ is continuous on a closed interval, there exists values $m$ and $M$ with $f(c_m) = m \leq f(x) \leq M=f(c_M)$, for some $c_m$ and $c_M$ in the interval $[a,b]$. Since $m \leq f(x) \leq M$, we must have:
+
+
+$$
+m \cdot (b-a) \leq K\cdot(b-a) \leq M\cdot(b-a).
+$$
+
+So in particular $K$ is in $[m, M]$. But $m$ and $M$ correspond to values of $f(x)$, so by the intermediate value theorem, $K=f(c)$ for some $c$ that must lie in between $c_m$ and $c_M$, which means as well that it must be in $[a,b]$.
+
+
+##### Proof of second part of Fundamental Theorem of Calculus
+
+
+The mean value theorem is exactly what is needed to prove formally the second part of the Fundamental Theorem of Calculus. Again, suppose $f(x)$ is continuous on $[a,b]$ with $a < b$. For any $a < x < b$, we define $F(x) = \int_a^x f(u) du$. Then the derivative of $F$ exists and is $f$.
+
+
+Let $h>0$. Then consider the forward difference $(F(x+h) - F(x))/h$. Rewriting gives:
+
+
+$$
+\frac{\int_a^{x+h} f(u) du - \int_a^x f(u) du}{h}  =\frac{\int_x^{x+h} f(u) du}{h} = f(\xi(h)).
+$$
+
+The value $\xi(h)$ is just the $c$ corresponding to a given value in $[x, x+h]$ guaranteed by the mean value theorem.  We only know that $x \leq \xi(h) \leq x+h$. But this is plenty - it says that $\lim_{h \rightarrow 0+} \xi(h) = x$. Using the fact that $f$ is continuous and the known properties of limits of compositions of functions this gives $\lim_{h \rightarrow 0+} f(\xi(h)) = f(x)$. But this means that the (right) limit of the secant line expression exists and is equal to $f(x)$, which is what we want to prove. Repeating a similar argument when $h < 0$, finishes the proof.
+
+
+The basic notion used is simply that for small $h$, this expression is well approximated by the left Riemann sum taken over $[x, x+h]$:
+
+
+$$
+f(\xi(h)) \cdot h = \int_x^{x+h} f(u) du.
+$$
+
+## Questions
+
+
+###### Question
+
+
+Between $0$ and $1$ a function is constantly $1$. Between $1$ and $2$ the function is constantly $2$. What is the average value of the function over the interval $[0,2]$?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+f(x) = x < 1 ? 1.0 : 2.0
+a,b = 0, 2
+val, _ = quadgk(f, a, b)
+numericq(val/(b-a))
+```
+
+###### Question
+
+
+Between $0$ and $2$ a function is constantly $1$. Between $2$ and $3$ the function is constantly $2$. What is the average value of the function over the interval $[0,3]$?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+f(x) = x < 2 ? 1.0 : 2.0
+a, b= 0, 2
+val, _ = quadgk(f, a, b)
+numericq(val/(b-a))
+```
+
+###### Question
+
+
+What integral will show the intuition of the Merton College scholars that the distance traveled by an object under uniformly increasing velocity starting at $v_0$ and ending at $v_t$ is equal to the distance traveled by an object with constant velocity of $(v_0 + v_t)/2$.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices = [
+"``\\int_0^t v(u) du = v^2/2 \\big|_0^t``",
+"``\\int_0^t (v(0) + v(u))/2 du = v(0)/2\\cdot t + x(u)/2\\ \\big|_0^t``",
+"``(v(0) + v(t))/2 \\cdot \\int_0^t du = (v(0) + v(t))/2 \\cdot t``"
+]
+answ = 1
+radioq(choices, answ)
+```
+
+###### Question
+
+
+Find the average value of $\cos(x)$ over the interval $[-\pi/2, \pi/2]$.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+f(x) = cos(x)
+a,b = -pi/2,pi/2
+val, _ = quadgk(f, a, b)
+val = val/(b-a)
+numericq(val)
+```
+
+###### Question
+
+
+Find the average value of $\cos(x)$ over the interval $[0, \pi]$.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+f(x) = cos(x)
+a,b = 0, pi
+val, _ = quadgk(f, a, b)
+val = val/(b-a)
+numericq(val)
+```
+
+###### Question
+
+
+Find the average value of $f(x) = e^{-2x}$ between $0$ and $2$.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+f(x) = exp(-2x)
+a, b = 0, 2
+val, _ = quadgk(f, a, b)
+val = val/(b-a)
+numericq(val)
+```
+
+###### Question
+
+
+Find the average value of $f(x) = \sin(x)^2$ over the $0$, $\pi$.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+f(x) = sin(x)^2
+a, b = 0, pi
+val, _ = quadgk(f, a, b)
+val = val/(b-a)
+numericq(val)
+```
+
+###### Question
+
+
+Which is bigger? The average value of $f(x) = x^{10}$ or the average value of $g(x) = \lvert x \rvert$ over the interval $[0,1]$?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices = [
+L"That of $f(x) = x^{10}$.",
+L"That of $g(x) = \lvert x \rvert$."]
+answ = 2
+radioq(choices, answ)
+```
+
+###### Question
+
+
+Define a family of functions over the interval $[0,1]$ by $f(x; a,b) = x^a \cdot (1-x)^b$. Which has a greater average, $f(x; 2,3)$ or $f(x; 3,4)$?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices = [
+"``f(x; 2,3)``",
+"``f(x; 3,4)``"
+]
+n1, _ = quadgk(x -> x^2 *(1-x)^3, 0, 1)
+n2, _ = quadgk(x -> x^3 *(1-x)^4, 0, 1)
+answ = 1 + (n1 < n2)
+radioq(choices, answ)
+```
+
+###### Question
+
+
+Suppose the average value of $f(x)$ over $[a,b]$ is $100$. What is the average value of $100 f(x)$ over $[a,b]$?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+numericq(100 * 100)
+```
+
+###### Question
+
+
+Suppose $f(x)$ is continuous and positive on $[a,b]$.
+
+
+  * Explain why for any $x > a$ it must be that:
+
+
+$$
+F(x) = \int_a^x f(x) dx > 0
+$$
+
+```{julia}
+#| hold: true
+#| echo: false
+choices = [
+L"Because the mean value theorem says this is $f(c) (x-a)$ for some $c$ and both terms are positive by the assumptions",
+"Because the definite integral is only defined for positive area, so it is always positive"
+]
+answ = 1
+radioq(choices, answ)
+```
+
+  * Explain why $F(x)$ is increasing.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices = [
+L"By the extreme value theorem, $F(x)$ must reach its maximum, hence it must increase.",
+L"By the intermediate value theorem, as $F(x) > 0$, it must be true that $F(x)$ is increasing",
+L"By the fundamental theorem of calculus, part I, $F'(x) = f(x) > 0$, hence $F(x)$ is increasing"
+]
+answ = 3
+radioq(choices, answ)
+```
+
+###### Question
+
+
+For $f(x) = x^2$, which is bigger: the average of the function $f(x)$ over $[0,1]$ or the geometric mean which is the exponential of the average of the logarithm of $f$ over the same interval?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+f(x) = x^2
+a,b = 0, 1
+val1 = quadgk(f, a, b)[1] / (b-a)
+val2 = exp(quadgk(x -> log(f(x)), a, b)[1] / (b - a))
+
+choices = [
+L"The average of $f$",
+L"The exponential of the average of $\log(f)$"
+]
+answ = val1 > val2 ? 1 : 2
+radioq(choices, answ)
+```
+

quarto/integrals/partial_fractions.qmd

@@ -0,0 +1,566 @@
+# Partial Fractions
+
+
+```{julia}
+#| echo: false
+
+import Logging
+Logging.disable_logging(Logging.Info) # or e.g. Logging.Info
+Logging.disable_logging(Logging.Warn)
+
+import SymPy
+function Base.show(io::IO, ::MIME"text/html", x::T) where {T <: SymPy.SymbolicObject}
+    println(io, "<span class=\"math-left-align\" style=\"padding-left: 4px; width:0; float:left;\"> ")
+    println(io, "\\[")
+    println(io, sympy.latex(x))
+    println(io, "\\]")
+    println(io, "</span>")
+end
+
+# hack to work around issue
+import Markdown
+import CalculusWithJulia
+function CalculusWithJulia.WeaveSupport.ImageFile(d::Symbol, f::AbstractString, caption; kwargs...)
+    nm = joinpath("..", string(d), f)
+    u = "![$caption]($nm)"
+    Markdown.parse(u)
+end
+
+nothing
+```
+
+```{julia}
+using CalculusWithJulia
+using SymPy
+```
+
+```{julia}
+#| echo: false
+#| results: "hidden"
+using CalculusWithJulia.WeaveSupport
+
+const frontmatter = (
+        title = "Partial Fractions",
+        description = "Calculus with Julia: Partial Fractions",
+        tags = ["CalculusWithJulia", "integrals", "partial fractions"],
+);
+nothing
+```
+
+---
+
+
+Integration is facilitated when an antiderivative for $f$ can be found, as then definite integrals can be evaluated through the fundamental theorem of calculus.
+
+
+However, despite integration being an algorithmic procedure, integration is not. There are "tricks" to try, such as substitution and integration by parts. These work in some cases. However, there are classes of functions for which algorithms exist. For example, the `SymPy` `integrate` function mostly implements an algorithm that decides if an elementary function has an antiderivative. The [elementary](http://en.wikipedia.org/wiki/Elementary_function) functions include exponentials, their inverses (logarithms), trigonometric functions, their inverses, and powers, including $n$th roots. Not every elementary function will have an antiderivative comprised of (finite) combinations of elementary functions. The typical example is $e^{x^2}$, which has no simple antiderivative, despite its ubiquitousness.
+
+
+There are classes of functions where an (elementary) antiderivative can always be found. Polynomials provide a case. More surprisingly, so do their ratios, *rational functions*.
+
+
+## Partial fraction decomposition
+
+
+Let $f(x) = p(x)/q(x)$, where  $p$ and $q$ are polynomial functions with real coefficients. Further, we assume without comment that $p$ and $q$ have no common factors. (If they did, we can divide them out, an act which has no effect on the integrability of $f(x)$.
+
+
+The function $q(x)$ will factor over the real numbers. The fundamental theorem of algebra can be applied to say that $q(x)=q_1(x)^{n_1} \cdots q_k(x)^{n_k}$ where $q_i(x)$ is a linear or quadratic polynomial and $n_k$ a positive integer.
+
+
+> **Partial Fraction Decomposition**: There are unique polynomials $a_{ij}$ with degree $a_{ij} <$ degree $q_i$ such that
+>
+> $$
+> \frac{p(x)}{q(x)} = a(x) + \sum_{i=1}^k \sum_{j=1}^{n_i} \frac{a_{ij}(x)}{q_i(x)^j}.
+> $$
+
+
+
+The method is attributed to John Bernoulli, one of the prolific Bernoulli brothers who put a stamp on several areas of math. This Bernoulli was a mentor to Euler.
+
+
+This basically says that each factor $q_i(x)^{n_i}$ contributes a term like:
+
+
+$$
+\frac{a_{i1}(x)}{q_i(x)^1} + \frac{a_{i2}(x)}{q_i(x)^2} + \cdots + \frac{a_{in_i}(x)}{q_i(x)^{n_i}},
+$$
+
+where each $a_{ij}(x)$ has degree less than the degree of $q_i(x)$.
+
+
+The value of this decomposition is that the terms $a_{ij}(x)/q_i(x)^j$ each have an antiderivative, and so the sum of them will also have an antiderivative.
+
+
+:::{.callout-note}
+## Note
+Many calculus texts will give some examples for finding a partial fraction decomposition. We push that work off to `SymPy`, as for all but the easiest cases - a few are in the problems - it can be a bit tedious.
+
+:::
+
+In `SymPy`, the `apart` function will find the partial fraction decomposition when a factorization is available. For example, here we see $n_i$ terms for each power of $q_i$
+
+
+```{julia}
+@syms a::real b::real c::real A::real B::real x::real
+```
+
+```{julia}
+apart((x-2)*(x-3) / (x*(x-1)^2*(x^2 + 2)^3))
+```
+
+### Sketch of proof
+
+
+A standard proof uses two facts of number systems: the division algorithm and a representation of the greatest common divisor in terms of sums, extended to polynomials. Our sketch shows how these are used.
+
+
+Take one of the factors of the denominators, and consider this representation of the rational function $P(x)/(q(x)^k Q(x))$ where there are no common factors to any of the three polynomials.
+
+
+Since $q(x)$ and $Q(x)$ share no factors, [Bezout's](http://tinyurl.com/kd6prns) identity says there exists polynomials $a(x)$ and $b(x)$  with:
+
+
+$$
+a(x) Q(x) + b(x) q(x) = 1.
+$$
+
+Then dividing by $q^k(x)Q(x)$ gives the decomposition
+
+
+$$
+\frac{1}{q(x)^k Q(x)} = \frac{a(x)}{q(x)^k} + \frac{b(x)}{q(x)^{k-1}Q(x)}.
+$$
+
+So we get by multiplying the $P(x)$:
+
+
+$$
+\frac{P(x)}{q(x)^k Q(x)} = \frac{A(x)}{q(x)^k} + \frac{B(x)}{q(x)^{k-1}Q(x)}.
+$$
+
+This may look more complicated, but what it does is peel off one term (The first) and leave something which is smaller, in this case by a factor of $q(x)$. This process can be repeated pulling off a power of a factor at a time until nothing is left to do.
+
+
+What remains is to establish that we can take $A(x) = a(x)\cdot P(x)$ with a degree less than that of $q(x)$.
+
+
+In Proposition 3.8 of [Bradley](http://www.m-hikari.com/imf/imf-2012/29-32-2012/cookIMF29-32-2012.pdf) and Cook we can see how. Recall the division algorithm, for example, says there are $q_k$ and $r_k$ with $A=q\cdot q_k + r_k$ where the degree of $r_k$ is less than that of $q$, which is linear or quadratic. This is repeatedly applied below:
+
+
+$$
+\begin{align*}
+\frac{A}{q^k} &= \frac{q\cdot q_k + r_k}{q^k}\\
+&= \frac{r_k}{q^k} + \frac{q_k}{q^{k-1}}\\
+&= \frac{r_k}{q^k} + \frac{q \cdot q_{k-1} + r_{k-1}}{q^{k-1}}\\
+&= \frac{r_k}{q^k} + \frac{r_{k-1}}{q^{k-1}} + \frac{q_{k-1}}{q^{k-2}}\\
+&= \frac{r_k}{q^k} + \frac{r_{k-1}}{q^{k-1}} + \frac{q\cdot q_{k-2} + r_{k-2}}{q^{k-2}}\\
+&= \cdots\\
+&= \frac{r_k}{q^k} + \frac{r_{k-1}}{q^{k-1}} + \cdots + q_1.
+\end{align*}
+$$
+
+So the term $A(x)/q(x)^k$ can be expressed in terms of a sum where the numerators or each term have degree less than $q(x)$, as expected by the statement of the theorem.
+
+
+## Integrating the terms in a partial fraction decomposition
+
+
+We discuss, by example, how each type of possible term in a partial fraction decomposition has an antiderivative. Hence, rational functions will *always* have an antiderivative that can be computed.
+
+
+### Linear factors
+
+
+For $j=1$, if $q_i$ is linear, then $a_{ij}/q_i^j$ must look like a constant over a linear term, or something like:
+
+
+```{julia}
+p = a/(x-c)
+```
+
+This has a logarithmic antiderivative:
+
+
+```{julia}
+integrate(p, x)
+```
+
+For $j > 1$, we have powers.
+
+
+```{julia}
+@syms j::positive
+integrate(a/(x-c)^j, x)
+```
+
+### Quadratic factors
+
+
+When $q_i$ is quadratic, it looks like $ax^2 + bx + c$. Then $a_{ij}$ can be a constant or a linear polynomial. The latter can be written as $Ax + B$.
+
+
+The integral of the following general form is presented below:
+
+
+$$
+\frac{Ax +B }{(ax^2  + bx + c)^j},
+$$
+
+With `SymPy`, we consider a few cases of the following form, which results from a shift of `x`
+
+
+$$
+\frac{Ax + B}{((ax)^2 \pm 1)^j}
+$$
+
+This can be done by finding a $d$ so that $a(x-d)^2 + b(x-d) + c = dx^2 + e = e((\sqrt{d/e}x^2 \pm 1)$.
+
+
+The integrals of the type $Ax/((ax)^2 \pm 1)$ can completed by $u$-substitution, with $u=(ax)^2 \pm 1$.
+
+
+For example,
+
+
+```{julia}
+integrate(A*x/((a*x)^2 + 1)^4, x)
+```
+
+The integrals of the type $B/((ax)^2\pm 1)$ are completed by trigonometric substitution and various reduction formulas. They can get involved, but are tractable. For example:
+
+
+```{julia}
+integrate(B/((a*x)^2 + 1)^4, x)
+```
+
+and
+
+
+```{julia}
+integrate(B/((a*x)^2 - 1)^4, x)
+```
+
+---
+
+
+In [Bronstein](http://www-sop.inria.fr/cafe/Manuel.Bronstein/publications/issac98.pdf) this characterization can be found - "This method, which dates back to Newton, Leibniz and Bernoulli, should not be used in practice, yet it remains the method found in most calculus texts and is often taught. Its major drawback is the factorization of the denominator of the integrand over the real or complex numbers." We can also find the following formulas which formalize the above exploratory calculations ($j>1$ and $b^2 - 4c < 0$ below):
+
+
+$$
+\begin{align*}
+\int \frac{A}{(x-a)^j} &= \frac{A}{1-j}\frac{1}{(x-a)^{1-j}}\\
+\int \frac{A}{x-a}     &= A\log(x-a)\\
+\int \frac{Bx+C}{x^2 + bx + c}     &= \frac{B}{2} \log(x^2 + bx + c) + \frac{2C-bB}{\sqrt{4c-b^2}}\cdot \arctan\left(\frac{2x+b}{\sqrt{4c-b^2}}\right)\\
+\int \frac{Bx+C}{(x^2 + bx + c)^j} &= \frac{B' x + C'}{(x^2 + bx + c)^{j-1}} + \int \frac{C''}{(x^2 + bx + c)^{j-1}}
+\end{align*}
+$$
+
+The first returns a rational function; the second yields a logarithm term; the third yields a logarithm and an arctangent term; while the last, which has explicit constants available, provides a reduction that can be recursively applied;
+
+
+That is integrating $f(x)/g(x)$, a rational function, will yield an output that looks like the following, where the functions are polynomials:
+
+
+$$
+\int f(x)/g(x) = P(x) + \frac{C(x)}{D{x}} + \sum v_i \log(V_i(x)) + \sum w_j \arctan(W_j(x))
+$$
+
+(Bronstein also sketches the modern method which is to use a Hermite reduction to express $\int (f/g) dx = p/q + \int (g/h) dx$, where $h$ is square free (the "`j`" are all $1$). The latter can be written over the complex numbers as logarithmic terms of the form $\log(x-a)$, the "`a`s"found following a method due to Trager and Lazard, and Rioboo, which is mentioned in the SymPy documentation as the method used.)
+
+
+#### Examples
+
+
+Find an antiderivative for $1/(x\cdot(x^2+1)^2)$.
+
+
+We have a partial fraction decomposition is:
+
+
+```{julia}
+q = (x * (x^2 + 1)^2)
+apart(1/q)
+```
+
+We see three terms. The first and second will be done by $u$-substitution, the third by a logarithm:
+
+
+```{julia}
+integrate(1/q, x)
+```
+
+---
+
+
+Find an antiderivative of $1/(x^2 - 2x-3)$.
+
+
+We again just let `SymPy` do the work. A partial fraction decomposition is given by:
+
+
+```{julia}
+𝒒 =  (x^2 - 2x - 3)
+apart(1/𝒒)
+```
+
+We see what should yield two logarithmic terms:
+
+
+```{julia}
+integrate(1/𝒒, x)
+```
+
+:::{.callout-note}
+## Note
+`SymPy` will find $\log(x)$ as an antiderivative for $1/x$, but more generally, $\log(\lvert x\rvert)$ is one.
+
+:::
+
+##### Example
+
+
+The answers found can become quite involved. [Corless](https://arxiv.org/pdf/1712.01752.pdf), Moir, Maza, and Xie use this example which at first glance seems tame enough:
+
+
+```{julia}
+ex = (x^2 - 1) / (x^4 + 5x^2 + 7)
+```
+
+But the integral is something best suited to a computer algebra system:
+
+
+```{julia}
+integrate(ex, x)
+```
+
+## Questions
+
+
+###### Question
+
+
+The partial fraction decomposition of $1/(x(x-1))$ must be of the form $A/x + B/(x-1)$.
+
+
+What is $A$? (Use `SymPy` or just put the sum over a common denominator and solve for $A$ and $B$.)
+
+
+```{julia}
+#| hold: true
+#| echo: false
+val = -1
+numericq(val)
+```
+
+What is $B$?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+val = 1
+numericq(val)
+```
+
+###### Question
+
+
+The following gives the partial fraction decomposition for a rational expression:
+
+
+$$
+\frac{3x+5}{(1-2x)^2} = \frac{A}{1-2x} + \frac{B}{(1-2x)^2}.
+$$
+
+Find $A$ (being careful with the sign):
+
+
+```{julia}
+#| hold: true
+#| echo: false
+numericq(-3/2)
+```
+
+Find $B$:
+
+
+```{julia}
+#| hold: true
+#| echo: false
+numericq(13/2)
+```
+
+###### Question
+
+
+The following specifies the general partial fraction decomposition for a rational expression:
+
+
+$$
+\frac{1}{(x+1)(x-1)^2} = \frac{A}{x+1} + \frac{B}{x-1} + \frac{C}{(x-1)^2}.
+$$
+
+Find $A$:
+
+
+```{julia}
+#| hold: true
+#| echo: false
+numericq(1/4)
+```
+
+Find $B$:
+
+
+```{julia}
+#| hold: true
+#| echo: false
+numericq(-1/4)
+```
+
+Find $C$:
+
+
+```{julia}
+#| hold: true
+#| echo: false
+numericq(1/2)
+```
+
+###### Question
+
+
+Compute the following exactly:
+
+
+$$
+\int_0^1 \frac{(x-2)(x-3)}{(x-4)^2\cdot(x-5)} dx
+$$
+
+Is $-6\log(5) - 5\log(3) - 1/6 + 11\log(4)$ the answer?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+yesnoq("yes")
+```
+
+###### Question
+
+
+In the assumptions for the partial fraction decomposition is the fact that $p(x)$ and $q(x)$ share no common factors. Suppose, this isn't the case and in fact we have:
+
+
+$$
+\frac{p(x)}{q(x)} = \frac{(x-c)^m s(x)}{(x-c)^n t(x)}.
+$$
+
+Here $s$ and $t$ are polynomials such that $s(c)$ and $t(c)$ are non-zero.
+
+
+If $m > n$, then why can we cancel out the $(x-c)^n$ and not have a concern?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices = [
+"`SymPy` allows it.",
+L"The value $c$ is a removable singularity, so the integral will be identical.",
+L"The resulting function has an identical domain and is equivalent for all $x$."
+]
+answ = 2
+radioq(choices, answ, keep_order=true)
+```
+
+If $m = n$, then why can we cancel out the $(x-c)^n$ and not have a concern?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices = [
+"`SymPy` allows it.",
+L"The value $c$ is a removable singularity, so the integral will be identical.",
+L"The resulting function has an identical domain and is equivalent for all $x$."
+]
+answ = 2
+radioq(choices, answ, keep_order=true)
+```
+
+If $m < n$, then why can we cancel out the $(x-c)^n$ and not have a concern?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices = [
+"`SymPy` allows it.",
+L"The value $c$ is a removable singularity, so the integral will be identical.",
+L"The resulting function has an identical domain and is equivalent for all $x$."
+]
+answ = 3
+radioq(choices, answ, keep_order=true)
+```
+
+##### Question
+
+
+The partial fraction decomposition, as presented, factors the denominator polynomial into linear and quadratic factors over the real numbers. Alternatively, factoring over the complex numbers is possible, resulting in terms like:
+
+
+$$
+\frac{a + ib}{x - (\alpha + i \beta)} + \frac{a - ib}{x - (\alpha - i \beta)}
+$$
+
+How to see that these give rise to real answers on integration is the point of this question.
+
+
+Breaking the terms up over $a$ and $b$ we have:
+
+
+$$
+\begin{align*}
+I &= \frac{a}{x - (\alpha + i \beta)} + \frac{a}{x - (\alpha - i \beta)} \\
+II &= i\frac{b}{x - (\alpha + i \beta)} - i\frac{b}{x - (\alpha - i \beta)}
+\end{align*}
+$$
+
+Integrating $I$ leads to two logarithmic terms, which are combined to give:
+
+
+$$
+\int I dx = a\cdot \log((x-(\alpha+i\beta)) \cdot (x - (\alpha-i\beta)))
+$$
+
+This involves no complex numbers, as:
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices = ["The complex numbers are complex conjugates, so the term in the logarithm will simply be ``x - 2\\alpha x + \\alpha^2 + \\beta^2``",
+           "The ``\\beta`` are ``0``, as the polynomials in question are real"]
+radioq(choices, 1)
+```
+
+The term $II$ benefits from this computation (attributed to Rioboo by [Corless et. al](https://arxiv.org/pdf/1712.01752.pdf))
+
+
+$$
+\frac{d}{dx} i \log(\frac{X+iY}{X-iY}) = 2\frac{d}{dx}\arctan(\frac{X}{Y})
+$$
+
+Applying this with $X=x - \alpha$ and $Y=-\beta$ shows that $\int II dx$ will be
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices = ["``-2b\\arctan((x - \\alpha)/(\\beta))``",
+           "``2b\\sec^2(-(x-\\alpha)/(-\\beta))``"]
+radioq(choices, 1)
+```
+

quarto/integrals/process.jl

@@ -0,0 +1,41 @@
+
+fnames = [
+          "area",
+          "ftc",
+
+          "substitution",
+          "integration_by_parts",
+          "partial_fractions", # XX add in trig integrals (cos()sin() stuff? mx or ^m... XXX
+          "improper_integrals", ##
+
+          "mean_value_theorem",
+          "area_between_curves",
+          "center_of_mass",
+          "volumes_slice",
+          #"volumes_shell", ## XXX add this in if needed, but not really that excited to now XXX
+          "arc_length",
+          "surface_area"
+]
+
+
+
+function process_file(nm, twice=false)
+    include("$nm.jl")
+    mmd_to_md("$nm.mmd")
+    markdownToHTML("$nm.md")
+    twice && markdownToHTML("$nm.md")
+end
+
+process_files(twice=false) = [process_file(nm, twice) for nm in fnames]
+
+
+
+
+"""
+## TODO integrals
+
+* add in volumes shell???
+* mean value theorem is light?
+* could add surface area problems
+
+"""

quarto/integrals/riemann.js

@@ -0,0 +1,31 @@
+const b = JXG.JSXGraph.initBoard('jsxgraph', {
+    boundingbox: [-0.5,0.3,1.5,-1/4], axis:true
+});
+
+var g = function(x) { return x*x*x*x +  10*x*x - 60* x + 100}
+var f = function(x) {return 1/Math.sqrt(g(x))};
+
+var type = "right";
+var l = 0;
+var r = 1;
+var rsum = function() {
+    return JXG.Math.Numerics.riemannsum(f,n.Value(), type, l, r);
+};
+var n =   b.create('slider', [[0.1, -0.05],[0.75,-0.05], [2,1,50]],{name:'n',snapWidth:1});
+
+
+
+
+var graph = b.create('functiongraph', [f, l, r]);
+var os = b.create('riemannsum',
+                  [f,
+                   function(){ return n.Value();},
+                   type, l, r
+                  ],
+                  {fillColor:'#ffff00', fillOpacity:0.3});
+
+
+
+b.create('text', [0.1,0.25, function(){
+    return 'Riemann sum='+(rsum().toFixed(4));
+}]);

quarto/integrals/substitution.qmd

@@ -0,0 +1,893 @@
+# Substitution
+
+
+```{julia}
+#| echo: false
+
+import Logging
+Logging.disable_logging(Logging.Info) # or e.g. Logging.Info
+Logging.disable_logging(Logging.Warn)
+
+import SymPy
+function Base.show(io::IO, ::MIME"text/html", x::T) where {T <: SymPy.SymbolicObject}
+    println(io, "<span class=\"math-left-align\" style=\"padding-left: 4px; width:0; float:left;\"> ")
+    println(io, "\\[")
+    println(io, sympy.latex(x))
+    println(io, "\\]")
+    println(io, "</span>")
+end
+
+# hack to work around issue
+import Markdown
+import CalculusWithJulia
+function CalculusWithJulia.WeaveSupport.ImageFile(d::Symbol, f::AbstractString, caption; kwargs...)
+    nm = joinpath("..", string(d), f)
+    u = "![$caption]($nm)"
+    Markdown.parse(u)
+end
+
+nothing
+```
+
+This section uses these add-on packages:
+
+
+```{julia}
+using CalculusWithJulia
+using Plots
+using SymPy
+
+```
+
+```{julia}
+#| echo: false
+#| results: "hidden"
+using CalculusWithJulia.WeaveSupport
+
+const frontmatter = (
+        title = "Substitution",
+        description = "Calculus with Julia: Substitution",
+        tags = ["CalculusWithJulia", "integrals", "substitution"],
+);
+nothing
+```
+
+---
+
+
+The technique of $u$-[substitution](https://en.wikipedia.org/wiki/Integration_by_substitution) is derived from reversing the chain rule: $[f(g(x))]' = f'(g(x)) g'(x)$.
+
+
+Suppose that $g$ is continuous and $u(x)$ is differentiable with $u'(x)$ being Riemann integrable. Then both these integrals are defined:
+
+
+$$
+\int_a^b g(u(t)) \cdot u'(t) dt, \quad \text{and}\quad \int_{u(a)}^{u(b)} g(x) dx.
+$$
+
+We wish to show they are equal.
+
+
+Let $G$ be an antiderivative of $g$, which exists as $g$ is assumed to be continuous. (By the Fundamental Theorem part I.) Consider the composition $G \circ u$. The chain rule gives:
+
+
+$$
+[G \circ u]'(t) = G'(u(t)) \cdot u'(t) = g(u(t)) \cdot u'(t).
+$$
+
+So,
+
+
+$$
+\begin{align*}
+\int_a^b g(u(t)) \cdot u'(t) dt &= \int_a^b (G \circ u)'(t) dt\\
+&= (G\circ u)(b) - (G\circ u)(a) \quad\text{(the FTC, part II)}\\
+&= G(u(b)) - G(u(a)) \\
+&= \int_{u(a)}^{u(b)} g(x) dx. \quad\text{(the FTC part II)}
+\end{align*}
+$$
+
+That is, this substitution formula applies:
+
+
+> $\int_a^b g(u(x)) u'(x) dx = \int_{u(a)}^{u(b)} g(x) dx.$
+
+
+
+Further, for indefinite integrals,
+
+
+> $\int f(g(x)) g'(x) dx = \int f(u) du.$
+
+
+
+We have seen a special case of substitution where $u(x) = x-c$ in the formula $\int_{a-c}^{b-c} g(x) dx= \int_a^b g(x-c)dx$.
+
+
+The main use of this is to take complicated things inside of the function $g$ out of the function (the $u(x)$) by renaming them, then accounting for the change of name.
+
+
+Some examples are in order.
+
+
+Consider:
+
+
+$$
+\int_0^{\pi/2} \cos(x) e^{\sin(x)} dx.
+$$
+
+Clearly the $\sin(x)$ inside the exponential is an issue. If we let $u(x) = \sin(x)$, then $u'(x) = \cos(x)$, and this becomes
+
+
+$$
+\int_0^2 u\prime(x) e^{u(x)} dx =
+\int_{u(0)}^{u(\pi/2)} e^x dx = e^x \big|_{\sin(0)}^{\sin(\pi/2)} = e^1 - e^0.
+$$
+
+This all worked, as the problem was such that it was more or less obvious what to choose for $u$ and $G$.
+
+
+### Integration by substitution
+
+
+The process of identifying the result of the chain rule in the function to integrate is not automatic, but rather a bit of an art. The basic step is to try some values and hope one works. Typically, this is taught by "substituting" in some value for part of the expression (basically the $u(x)$) and seeing what happens.
+
+
+In the above problem, $\int_0^{\pi/2} \cos(x) e^{\sin(x)} dx$, we might just rename $\sin(x)$ to be $u$ (suppressing the "of $x$ part). Then we need to rewrite the "$dx$" part of the integral. We know in this case that $du/dx = \cos(x)$. In terms of differentials, this gives $du = \cos(x) dx$. But this allows us to substitute in with $u$ and $du$ as is possible:
+
+
+$$
+\int_0^{\pi/2} \cos(x) e^{\sin(x)} dx = \int_0^{\pi/2}  e^{\sin(x)} \cdot \cos(x) dx = \int_{u(0)}^{u(\pi)} e^u du.
+$$
+
+---
+
+
+Let's illustrate with a new problem: $\int_0^2 4x e^{x^2} dx$.
+
+
+Again, we see that the $x^2$ inside the exponential is a complication. Letting $u = x^2$ we have $du = 2x dx$. We  have $4xdx$ in the original problem, so we will end up with $2du$:
+
+
+$$
+\int_0^2 4x e^{x^2} dx = 2\int_0^2 e^{x^2} \cdot 2x dx = 2\int_{u(0)}^{u(2)} e^u du = 2 \int_0^4 e^u du =
+2 e^u\big|_{u=0}^4 = 2(e^4 - 1).
+$$
+
+---
+
+
+Consider now $\int_0^1 2x^2 \sqrt{1 + x^3} dx$. Here we see that the $1 + x^3$ makes the square root term complicated. If we call this $u$, then what is $du$? Clearly, $du = 3x^2 dx$, or $(1/3)du = x^2 dx$, so we can rewrite this as:
+
+
+$$
+\int_0^1 2x^2 \sqrt{1 + x^3} dx = \int_{u(0)}^{u(1)} 2 \sqrt{u} (1/3) du = 2/3 \cdot \frac{u^{3/2}}{3/2} \big|_1^2 =
+\frac{4}{9} \cdot(2^{3/2} - 1).
+$$
+
+---
+
+
+Consider $\int_0^{\pi} \cos(x)^3 \sin(x) dx$. The $\cos(x)$ function inside the $x^3$ function is complicated. We let $u(x) = \cos(x)$ and see what that implies: $du = \sin(x) dx$, which we see is part of the question. So the above becomes:
+
+
+$$
+\int_0^{\pi} \cos(x)^3 \sin(x) dx = \int_{u(0)}^{u(\pi)} u^3 du= \frac{u^4}{4}\big|_0^0 = 0.
+$$
+
+Changing limits leaves the two endpoints the same, which means the total area after substitution is $0$. A graph of this function shows that about $\pi/2$ the function has odd-like symmetry, so the answer of $0$ is supported by the plot:
+
+
+```{julia}
+#| hold: true
+f(x) = cos(x)^3 * sin(x)
+plot(f, 0, 1pi)
+```
+
+---
+
+
+Consider $\int_1^e \log(x)/x dx$. There isn't really an "inside" function here, but instead just a tricky $\log(x)$. If we let $u=\log(x)$, what happens? We get $du = 1/x \cdot dx$, which we see present in the original. So with this, we have:
+
+
+$$
+\int_1^e \frac{\log(x)}{x} dx = \int_{u(1)}^{u(e)} u du = \frac{u^2}{2}\big|_0^1 = \frac{1}{2}.
+$$
+
+##### Example: Transformations
+
+
+We say that the area intrinsically discussed in the definite integral $A=\int_a^b f(x-c) dx$ is unaffected by shifts, in that $A = \int_{a-c}^{b-c} f(x) dx$. What about more general transformations? For example: if $g(x) = (1/h) \cdot f((x-c)/h)$ for values $c$ and $h$ what is the integral over $a$ to $b$ in terms of the function $f(x)$?
+
+
+If $A = \int_a^b (1/h) \cdot f((x-c)/h) dx$ then we let $u = (x-c)/h$. With this, $du = 1/h \cdot dx$. This allows a straight substitution:
+
+
+$$
+A = \int_a^b \frac{1}{h} f(\frac{x-c}{h}) dx = \int_{(a-c)/h}^{(b-c)/h} f(u) du.
+$$
+
+So the answer is: the area under the transformed function over $a$ to $b$ is the area of the function over the transformed region.
+
+
+For example, consider the "hat" function $f(x) = 1 - \lvert x \rvert $ when $-1 \leq x \leq 1$ and $0$ otherwise. The area under $f$ is just $1$ - the graph forms a triangle with base of length $2$ and height $1$. If we take any values of $c$ and $h$, what do we find for the area under the curve of the transformed function?
+
+
+Let $u(x) = (x-c)/h$ and $g(x) = h f(u(x))$. Then, as $du = 1/h dx$
+
+
+$$
+\begin{align}
+\int_{c-h}^{c+h} g(x) dx
+&= \int_{c-h}^{c+h} h f(u(x)) dx\\
+&= \int_{u(c-h)}^{u(c+h)} f(u) du\\
+&= \int_{-1}^1 f(u) du\\
+&= 1.
+\end{align}
+$$
+
+So the area of this transformed function is still $1$. The shifting by $c$ we know doesn't effect the area, the scaling by $h$ inside of $f$ does, but is balanced out by the multiplication by $h$ outside of $f$.
+
+
+##### Example: Speed versus velocity
+
+
+The "velocity" of an object includes a sense of direction in addition to the sense of magnitude. The "speed" just includes the sense of magnitude. Speed is always non-negative, whereas velocity is a signed quantity.
+
+
+As mentioned previously, position is the integral of velocity, as expressed precisely through this equation:
+
+
+$$
+x(t) = \int_0^t v(u) du - x(0).
+$$
+
+What is the integral of speed?
+
+
+If $v(t)$ is the velocity, the $s(t) = \lvert v(t) \rvert$ is the speed. If integrating either $s(t)$ or $v(t)$, the  integrals would agree when $v(t) \geq 0$. However,  when $v(t) \leq 0$, the position back tracks so $x(t)$ decreases, where the integral of $s(t)$ would only increase.
+
+
+This integral
+
+
+$$
+td(t) = \int_0^t s(u) du = \int_0^t \lvert v(u) \rvert du,
+$$
+
+Gives the *total distance* traveled.
+
+
+To illustrate with a simple example, if a car drives East for one hour at 60 miles per hour, then heads back West for an hour at 60 miles per hour, the car's position after one hour is $x(2) = x(0)$, with a change in position $x(2) - x(0) = 0$. Whereas, the total distance traveled is $120$ miles. (Gas is paid on total distance, not change in position!). What are the formulas for speed and velocity? Clearly $s(t) = 60$, a constant, whereas here $v(t) = 60$ for $0 \leq t \leq 1$ and $-60$ for $1 < t \leq 2$.
+
+
+Suppose $v(t)$ is given by $v(t) = (t-2)^3/3 - 4(t-2)/3$. If $x(0)=0$ Find the position after 3 time units and the total distance traveled.
+
+
+We let $u(t) = t - 2$ so $du=dt$. The position is given by
+
+
+$$
+\int_0^3  ((t-2)^3/3 - 4(t-2)/3) dt = \int_{u(0)}^{u(3)} (u^3/3 - 4/3 u) du =
+(\frac{u^4}{12} - \frac{4}{3}\frac{u^2}{2}) \big|_{-2}^1 = \frac{3}{4}.
+$$
+
+The speed is similar, but we have to work harder:
+
+
+$$
+\int_0^3 \lvert v(t) \rvert dt = \int_0^3  \lvert ((t-2)^3/3 - 4(t-2)/3) \rvert dt  =
+\int_{-2}^1 \lvert u^3/3 - 4u/3 \rvert du.
+$$
+
+But $u^3/3 - 4u/3 = (1/3) \cdot u(u-1)(u+2)$, so between $-2$ and $0$ it is positive and between $0$ and $1$ negative, so this integral is:
+
+
+$$
+\begin{align*}
+\int_{-2}^0 (u^3/3 - 4u/3 ) du + \int_{0}^1 -(u^3/3 - 4u/3) du
+&= (\frac{u^4}{12} - \frac{4}{3}\frac{u^2}{2}) \big|_{-2}^0 - (\frac{u^4}{12} - \frac{4}{3}\frac{u^2}{2}) \big|_{0}^1\\
+&= \frac{4}{3} - -\frac{7}{12}\\
+&= \frac{23}{12}.
+\end{align*}
+$$
+
+##### Example
+
+
+In probability, the normal distribution plays an outsized role. This distribution is characterized by a family of *density* functions:
+
+
+$$
+f(x; \mu, \sigma) = \frac{1}{\sqrt{2\pi}}\frac{1}{\sigma} \exp(-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^2).
+$$
+
+Integrals involving this function are typically transformed by substitution. For example:
+
+
+$$
+\begin{align*}
+\int_a^b f(x; \mu, \sigma) dx
+&= \int_a^b  \frac{1}{\sqrt{2\pi}}\frac{1}{\sigma} \exp(-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^2) dx \\
+&= \int_{u(a)}^{u(b)} \frac{1}{\sqrt{2\pi}} \exp(-\frac{1}{2}u^2) du \\
+&= \int_{u(a)}^{u(b)} f(u; 0, 1) du,
+\end{align*}
+$$
+
+where $u = (x-\mu)/\sigma$, so $du = (1/\sigma) dx$.
+
+
+This shows that integrals involving a normal density with parameters $\mu$ and $\sigma$ can be computed using the *standard* normal density with $\mu=0$ and $\sigma=1$. Unfortunately, there is no elementary antiderivative for $\exp(-u^2/2)$, so integrals for the standard normal must be numerically approximated.
+
+
+There is a function `erf` in the `SpecialFunctions` package (which is loaded by `CalculusWithJulia`) that computes:
+
+
+$$
+\int_0^x \frac{2}{\sqrt{\pi}} \exp(-t^2) dt
+$$
+
+A further change of variables by $t = u/\sqrt{2}$ (with $\sqrt{2}dt = du$) gives:
+
+
+$$
+\begin{align*}
+\int_a^b f(x; \mu, \sigma) dx &=
+\int_{t(u(a))}^{t(u(b))} \frac{\sqrt{2}}{\sqrt{2\pi}} \exp(-t^2) dt\\
+&= \frac{1}{2} \int_{t(u(a))}^{t(u(b))} \frac{2}{\sqrt{\pi}} \exp(-t^2) dt
+\end{align*}
+$$
+
+Up to a factor of $1/2$ this is `erf`.
+
+
+So we would have, for example, with $\mu=1$,$\sigma=2$ and $a=1$ and $b=3$ that:
+
+
+$$
+\begin{align*}
+t(u(a)) &= (1 - 1)/2/\sqrt{2} = 0\\
+t(u(b)) &= (3 - 1)/2/\sqrt{2} = \frac{1}{\sqrt{2}}\\
+\int_1^3 f(x; 1, 2)
+&= \frac{1}{2} \int_0^{1/\sqrt{2}} \frac{2}{\sqrt{\pi}} \exp(-t^2) dt.
+\end{align*}
+$$
+
+Or
+
+
+```{julia}
+1/2 * erf(1/sqrt(2))
+```
+
+:::{.callout-note}
+## The `Distributions` package
+The above calculation is for illustration purposes. The add-on package `Distributions` makes much quicker work of such a task for the normal distribution and many other distributions from probability and statistics.
+
+:::
+
+## SymPy and substitution
+
+
+The `integrate` function in `SymPy` can handle most problems which involve substitution. Here are a few examples:
+
+
+  * This integral, $\int_0^2 4x/\sqrt{x^2 +1}dx$, involves a substitution for $x^2 + 1$:
+
+
+```{julia}
+@syms x::real t::real
+integrate(4x / sqrt(x^2 + 1), (x, 0, 2))
+```
+
+  * This integral, $\int_e^{e^2} 1/(x\log(x)) dx$ involves a substitution of $u=\log(x)$. Here we see the answer:
+
+
+```{julia}
+#| hold: true
+f(x) = 1/(x*log(x))
+integrate(f(x), (x, sympy.E, sympy.E^2))
+```
+
+(We used `sympy.E)` - and not `e` - to avoid any conversion to floating point, which could yield an inexact answer.)
+
+
+The antiderivative is interesting here; it being an *iterated* logarithm.
+
+
+```{julia}
+integrate(1/(x*log(x)), x)
+```
+
+### Failures...
+
+
+Not every integral problem lends itself to solution by substitution. For example, we can use substitution to evaluate the integral of $xe^{-x^2}$, but for $e^{-x^2}$ or $x^2e^{-x^2}$. The first has no familiar antiderivative, the second is done by a different technique.
+
+
+Even when substitution can be used, `SymPy` may not be able to algorithmically identify it. The main algorithm used can determine if expressions involving rational functions, radicals, logarithms, and exponential functions is integrable. Missing from this list are absolute values.
+
+
+For some such problems, we can help `SymPy` out - by breaking the integral into pieces where we know the sign of the expression.
+
+
+For substitution problems, we can also help out. For example, to find an antiderivative for
+
+
+$$
+\int(1 + \log(x)) \sqrt{1 + (x\log(x))^2} dx
+$$
+
+A quick attempt with `SymPy` turns up nothing:
+
+
+```{julia}
+𝒇(x) = (1 + log(x)) * sqrt(1 + (x*log(x))^2 )
+integrate(𝒇(x), x)
+```
+
+But were we to try $u=x\log(x)$, we'd see that this simplifies to $\int \sqrt{1 + u^2} du$, which has some hope of having an antiderivative.
+
+
+We can help `SymPy` out by substitution:
+
+
+```{julia}
+u(x) = x * log(x)
+@syms w dw
+ex = 𝒇(x)
+ex₁ = ex(u(x) => w, diff(u(x),x) => dw)
+```
+
+This verifies the above. Can it be integrated in `w`? The "`dw`" is only for familiarity, `SymPy` doesn't use this, so we set it to 1 then integrate:
+
+
+```{julia}
+ex₂ = ex₁(dw => 1)
+ex₃ = integrate(ex₂, w)
+```
+
+Finally, we put back in the `u(x)` to get an antiderivative.
+
+
+```{julia}
+ex₃(w => u(x))
+```
+
+:::{.callout-note}
+## Note
+Lest it be thought this is an issue with `SymPy`, but not other systems, this example was [borrowed](http://faculty.uml.edu/jpropp/142/Integration.pdf) from an illustration for helping Mathematica.
+
+:::
+
+## Trigonometric substitution
+
+
+Wait, in the last example an antiderivative for $\sqrt{1 + u^2}$ was found. But how? We haven't discussed this yet.
+
+
+This can be found using *trigonometric* substitution. In this example, we know that $1 + \tan(\theta)^2$ simplifies to $\sec(\theta)^2$, so we might *try* a substitution of $\tan(u)=x$. This would simplify $\sqrt{1 + x^2}$ to $\sqrt{1 + \tan(u)^2} = \sqrt{\sec(u)^2}$ which is $\lvert \sec(u) \rvert$. What of $du$? The chain rule gives $\sec(u)^2du = dx$. In short we get:
+
+
+$$
+\int \sqrt{1 + x^2} dx = \int \sec(u)^2 \lvert \sec(u) \rvert du = \int \sec(u)^3 du,
+$$
+
+if we know $\sec(u) \geq 0$.
+
+
+This leaves still the question of integrating $\sec(u)^3$, which we aren't (yet) prepared to discuss, but we see that this type of substitution can re-express an integral in a new way that may pay off.
+
+
+#### Examples
+
+
+Let's see some examples where a trigonometric substitution is all that is needed.
+
+
+##### Example
+
+
+Consider $\int 1/(1+x^2) dx$. This is an antiderivative of some function, but if that isn't observed, we might notice the $1+x^2$ and try to simplify that. First, an attempt at a $u$-substitution:
+
+
+Letting $u = 1+x^2$ we get $du = 2xdx$ which gives $\int (1/u) (2x) du$. We aren't able to address the "$2x$" part successfully, so this attempt is for naught.
+
+
+Now we try a trigonometric substitution, taking advantage of the identity $1+\tan(x)^2 = \sec(x)^2$. Letting $\tan(u) = x$ yields $\sec(u)^2 du = dx$ and we get:
+
+
+$$
+\int \frac{1}{1+x^2} dx = \int \frac{1}{1 + \tan(u)^2} \sec(u)^2 du = \int 1 du = u.
+$$
+
+But $\tan(u) = x$, so in terms of $x$, an antiderivative is just $\tan^{-1}(x)$, or the arctangent. Here we verify with `SymPy`:
+
+
+```{julia}
+integrate(1/(1+x^2), x)
+```
+
+The general form allows $a^2 + (bx)^2$ in the denominator (squared so both are positive and the answer is nicer):
+
+
+```{julia}
+#| hold: true
+@syms a::real, b::real, x::real
+integrate(1 / (a^2  + (b*x)^2), x)
+```
+
+##### Example
+
+
+The expression $1-x^2$ can be attacked by the substitution $\sin(u) =x$ as then $1-x^2 = 1-\cos(u)^2 = \sin(u)^2$. Here we see this substitution being used successfully:
+
+
+$$
+\begin{align*}
+\int \frac{1}{\sqrt{9 - x^2}} dx &= \int \frac{1}{\sqrt{9 - (3\sin(u))^2}} \cdot 3\cos(u) du\\
+&=\int \frac{1}{3\sqrt{1 - \sin(u)^2}}\cdot3\cos(u) du \\
+&= \int du \\
+&= u \\
+&= \sin^{-1}(x/3).
+\end{align*}
+$$
+
+Further substitution allows the following integral to be solved for an antiderivative:
+
+
+```{julia}
+#| hold: true
+@syms a::real, b::real
+integrate(1 / sqrt(a^2 - b^2*x^2), x)
+```
+
+##### Example
+
+
+The expression $x^2 - 1$ is a bit different, this lends itself to $\sec(u) = x$ for a substitution, for $\sec(u)^2 - 1 = \tan(u)^2$. For example, we try $\sec(u) = x$ to integrate:
+
+
+$$
+\begin{align*}
+\int \frac{1}{\sqrt{x^2 - 1}} dx &= \int \frac{1}{\sqrt{\sec(u)^2 - 1}} \cdot \sec(u)\tan(u) du\\
+&=\int \frac{1}{\tan(u)}\sec(u)\tan(u) du\\
+&= \int \sec(u) du.
+\end{align*}
+$$
+
+This doesn't seem that helpful, but the antiderivative to $\sec(u)$ is $\log\lvert (\sec(u) + \tan(u))\rvert$, so we can proceed to get:
+
+
+$$
+\begin{align*}
+\int \frac{1}{\sqrt{x^2 - 1}} dx &= \int \sec(u) du\\
+&= \log\lvert (\sec(u) + \tan(u))\rvert\\
+&= \log\lvert x + \sqrt{x^2-1} \rvert.
+\end{align*}
+$$
+
+SymPy gives a different representation using the arccosine:
+
+
+```{julia}
+#| hold: true
+@syms a::positive, b::positive, x::real
+integrate(1 / sqrt(a^2*x^2 - b^2), x)
+```
+
+##### Example
+
+
+The equation of an ellipse is $x^2/a^2 + y^2/b^2 = 1$. Suppose $a,b>0$. The area under the function $b \sqrt{1 - x^2/a^2}$ between $-a$ and $a$ will then be half the area of the ellipse. Find the area enclosed by the ellipse.
+
+
+We need to compute:
+
+
+$$
+2\int_{-a}^a b \sqrt{1 - x^2/a^2} dx =
+4 b \int_0^a\sqrt{1 - x^2/a^2} dx.
+$$
+
+Letting $\sin(u) = x/a$ gives $a\cos(u)du = dx$ and an antiderivative is found with:
+
+
+$$
+4 b \int_0^a \sqrt{1 - x^2/a^2} dx = 4b \int_0^{\pi/2} \sqrt{1-u^2} a \cos(u) du
+= 4ab \int_0^{\pi/2} \cos(u)^2 du
+$$
+
+The identify $\cos(u)^2 = (1 +  \cos(2u))/2$ makes this tractable:
+
+
+$$
+\begin{align*}
+4ab \int \cos(u)^2 du
+&= 4ab\int_0^{\pi/2}(\frac{1}{2} + \frac{\cos(2u)}{2}) du\\
+&= 4ab(\frac{1}{2}u + \frac{\sin(2u)}{4})\big|_0^{\pi/2}\\
+&= 4ab (\pi/4 + 0) = \pi ab.
+\end{align*}
+$$
+
+Keeping in mind that that a circle with radius $a$ is an ellipse with $b=a$, we see that this gives the correct answer for a circle.
+
+
+## Questions
+
+
+###### Question
+
+
+For $\int \sin(x) \cos(x) dx$, let $u=\sin(x)$. What is the resulting substitution?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices = [
+"``\\int u du``",
+"``\\int u (1 - u^2) du``",
+"``\\int u \\cos(x) du``"
+]
+answ = 1
+radioq(choices, answ)
+```
+
+###### Question
+
+
+For $\int \tan(x)^4 \sec(x)2 dx$ what $u$-substitution makes this easy?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices = [
+"``u=\\tan(x)``",
+"``u=\\tan(x)^4``",
+"``u=\\sec(x)``",
+"``u=\\sec(x)^2``"
+]
+answ = 1
+radioq(choices, answ)
+```
+
+###### Question
+
+
+For $\int x \sqrt{x^2 - 1} dx$ what $u$ substitution makes this easy?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices = [
+"``u=x^2 - 1``",
+"``u=x^2``",
+"``u=\\sqrt{x^2 - 1}``",
+"``u=x``"
+]
+answ = 1
+radioq(choices, answ)
+```
+
+###### Question
+
+
+For $\int x^2(1-x)^2 dx$ will the substitution $u=1-x$ prove effective?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+yesnoq("no")
+```
+
+What about expanding the factored polynomial to get a fourth degree polynomial, will this prove effective?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+yesnoq("yes")
+```
+
+###### Question
+
+
+For  $\int (\log(x))^3/x dx$ the substitution $u=\log(x)$ reduces this to what?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices = [
+"``\\int u^3 du``",
+"``\\int u du``",
+"``\\int u^3/x du``"
+]
+answ = 1
+radioq(choices, answ)
+```
+
+###### Question
+
+
+For $\int \tan(x) dx$ what substitution will prove effective?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices = [
+"``u=\\cos(x)``",
+"``u=\\sin(x)``",
+"``u=\\tan(x)``"
+]
+answ = 1
+radioq(choices, answ)
+```
+
+###### Question
+
+
+Integrating $\int_0^1 x \sqrt{1 - x^2} dx$ can be done by using the $u$-substitution $u=1-x^2$. This yields an integral
+
+
+$$
+\int_a^b \frac{-\sqrt{u}}{2} du.
+$$
+
+What are $a$ and $b$?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices = [
+"``a=0,~ b=1``",
+"``a=1,~ b=0``",
+"``a=0,~ b=0``",
+"``a=1,~ b=1``"
+]
+answ = 1
+radioq(choices, answ)
+```
+
+###### Question
+
+
+The integral $\int \sqrt{1 - x^2} dx$ lends itself to what substitution?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices = [
+"``\\sin(u) = x``",
+"``\\tan(u) = x``",
+"``\\sec(u) = x``",
+"``u = 1 - x^2``"
+]
+answ = 1
+radioq(choices, answ)
+```
+
+###### Question
+
+
+The integral $\int x/(1+x^2) dx$ lends itself to what substitution?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices = [
+"``u = 1 + x^2``",
+"``\\sin(u) = x``",
+"``\\tan(u) = x``",
+"``\\sec(u) = x``"
+]
+answ = 1
+radioq(choices, answ)
+```
+
+###### Question
+
+
+The integral $\int dx / \sqrt{1 - x^2}$ lends itself to what substitution?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices = [
+"``\\sin(u) = x``",
+"``\\tan(u) = x``",
+"``\\sec(u) = x``",
+"``u = 1 - x^2``"
+]
+answ = 1
+radioq(choices, answ)
+```
+
+###### Question
+
+
+The integral $\int dx / \sqrt{x^2 - 16}$ lends itself to what substitution?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices = [
+"``4\\sec(u) = x``",
+"``\\sec(u) = x``",
+"``4\\sin(u) = x``",
+"``\\sin(u) = x``"]
+answ = 1
+radioq(choices, answ)
+```
+
+###### Question
+
+
+The integral $\int dx / (a^2 + x^2)$ lends itself to what substitution?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices = [
+"``\\tan(u) = x``",
+"``\\tan(u) = x``",
+"``a\\sec(u) = x``",
+"``\\sec(u) = x``"]
+answ = 1
+radioq(choices, answ)
+```
+
+###### Question
+
+
+The integral $\int_{1/2}^1 \sqrt{1 - x^2}dx$ can be approached with the substitution $\sin(u) = x$ giving:
+
+
+$$
+\int_a^b \cos(u)^2 du.
+$$
+
+What are $a$ and $b$?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices =[
+"``a=\\pi/6,~ b=\\pi/2``",
+"``a=\\pi/4,~ b=\\pi/2``",
+"``a=\\pi/3,~ b=\\pi/2``",
+"``a=1/2,~ b= 1``"
+]
+answ =1
+radioq(choices, answ)
+```
+
+###### Question
+
+
+How would we verify that $\log\lvert (\sec(u) + \tan(u))\rvert$ is an antiderivative for $\sec(u)$?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices = [
+L"We could differentiate $\sec(u)$.",
+L"We could differentiate $\log\lvert (\sec(u) + \tan(u))\rvert$ "]
+answ = 2
+radioq(choices, answ)
+```
+

quarto/integrals/surface_area.qmd

@@ -0,0 +1,614 @@
+# Surface Area
+
+
+```{julia}
+#| echo: false
+
+import Logging
+Logging.disable_logging(Logging.Info) # or e.g. Logging.Info
+Logging.disable_logging(Logging.Warn)
+
+import SymPy
+function Base.show(io::IO, ::MIME"text/html", x::T) where {T <: SymPy.SymbolicObject}
+    println(io, "<span class=\"math-left-align\" style=\"padding-left: 4px; width:0; float:left;\"> ")
+    println(io, "\\[")
+    println(io, sympy.latex(x))
+    println(io, "\\]")
+    println(io, "</span>")
+end
+
+# hack to work around issue
+import Markdown
+import CalculusWithJulia
+function CalculusWithJulia.WeaveSupport.ImageFile(d::Symbol, f::AbstractString, caption; kwargs...)
+    nm = joinpath("..", string(d), f)
+    u = "![$caption]($nm)"
+    Markdown.parse(u)
+end
+
+nothing
+```
+
+This section uses these add-on packages:
+
+
+```{julia}
+using CalculusWithJulia
+using Plots
+using SymPy
+using QuadGK
+```
+
+```{julia}
+#| echo: false
+#| results: "hidden"
+using CalculusWithJulia.WeaveSupport
+
+const frontmatter = (
+        title = "Surface Area",
+        description = "Calculus with Julia: Surface Area",
+        tags = ["CalculusWithJulia", "integrals", "surface area"],
+);
+fig_size=(800, 600)
+nothing
+```
+
+---
+
+
+## Surfaces of revolution
+
+
+```{julia}
+#| hold: true
+#| echo: false
+imgfile = "figures/gehry-hendrix.jpg"
+caption = """
+
+The exterior of the Jimi Hendrix Museum in Seattle has the signature
+style of its architect Frank Gehry. The surface is comprised of
+patches. A general method to find the amount of material to cover the
+surface - the surface area - might be to add up the area of *each* of the
+patches. However, in this section we will see for surfaces of
+revolution, there is an easier way. (Photo credit to
+[firepanjewellery](http://firepanjewellery.com/).)
+"""
+
+ImageFile(:integrals, imgfile, caption)
+```
+
+> The surface area generated by rotating the graph of $f(x)$ between $a$ and $b$  about the $x$-axis is given by the integral
+>
+> $$
+> \int_a^b 2\pi f(x) \cdot \sqrt{1 + f'(x)^2} dx.
+> $$
+>
+> If the curve is parameterized by $(g(t), f(t))$ between $a$ and $b$ then the surface area is
+>
+> $$
+> \int_a^b 2\pi f(t) \cdot \sqrt{g'(t)^2 + f'(t)^2} dx.
+> $$
+>
+> These formulas do not add in the surface area of either of the ends.
+
+
+
+```{julia}
+#| hold: true
+#| echo: false
+F₀(u,v) = [u, u*cos(v), u*sin(v)] # a cone
+us = range(0, 1, length=25)
+vs = range(0, 2pi, length=25)
+ws = unzip(F₀.(us, vs')) # make square
+surface(ws..., legend=false)
+plot!([-0.5,1.5], [0,0],[0,0])
+```
+
+The above figure shows a cone (the line $y=x$) presented as a surface of revolution about the $x$-axis.
+
+
+To see why this formula is as it is, we look at the parameterized case, the first one being a special instance with $g(t) =t$.
+
+
+Let a partition of $[a,b]$ be given by $a = t_0 < t_1 < t_2 < \cdots < t_n =b$. This breaks the curve into a collection of line segments. Consider the line segment connecting $(g(t_{i-1}), f(t_{i-1}))$ to $(g(t_i), f(t_i))$. Rotating this around the $x$ axis will generate something approximating a disc, but in reality will be the frustum of a cone. What will be the surface area?
+
+
+Consider a right-circular cone parameterized by an angle $\theta$ and the largest radius $r$ (so that the height satisfies $r/h=\tan(\theta)$). If this cone were made of paper, cut up a side, and layed out flat, it would form a sector of a circle, whose area would be $R\gamma$ where $R$ is the radius of the circle (also the side length of our cone), and $\gamma$ an angle that we can figure out from $r$ and $\theta$. To do this, we note that the arc length of the circle's edge is $R\gamma$ and also the circumference of the bottom of the cone so $R\gamma = 2\pi r$. With all this, we can solve to get $A = \pi r^2/\sin(\theta)$. But we have a frustum of a cone with radii $r_0$ and $r_1$, so the surface area is a difference: $A = \pi (r_1^2 - r_0^2) /\sin(\theta)$.
+
+
+Relating this to our values in terms of $f$ and $g$, we have $r_1=f(t_i)$, $r_0 = f(t_{i-1})$, and $\sin(\theta) = \Delta f / \sqrt{(\Delta g)^2 + (\Delta f)^2}$, where $\Delta f = f(t_i) - f(t_{i-1})$ and similarly for $\Delta g$.
+
+
+Putting this altogether we get that the surface area generarated by rotating the line segment around the $x$ axis is
+
+
+$$
+\text{sa}_i = \pi (f(t_i)^2 - f(t_{i-1})^2) \cdot \sqrt{(\Delta g)^2 + (\Delta f)^2} / \Delta f =
+\pi (f(t_i) + f(t_{i-1})) \cdot \sqrt{(\Delta g)^2 + (\Delta f)^2}.
+$$
+
+(This is $2 \pi$ times the average radius times the slant height.)
+
+
+As was done in the derivation of the formula for arc length, these pieces are multiplied both top and bottom by $\Delta t = t_{i} - t_{i-1}$. Carrying the bottom inside the square root and noting that by the mean value theorem $\Delta g/\Delta t = g(\xi)$ and $\Delta f/\Delta t = f(\psi)$ for some $\xi$ and $\psi$ in $[t_{i-1}, t_i]$, this becomes:
+
+
+$$
+\text{sa}_i = \pi (f(t_i) + f(t_{i-1})) \cdot \sqrt{(g'(\xi))^2 + (f'(\psi))^2} \cdot (t_i - t_{i-1}).
+$$
+
+Adding these up, $\text{sa}_1 + \text{sa}_2 + \cdots + \text{sa}_n$, we get a Riemann sum approximation to the integral
+
+
+$$
+\text{SA} = \int_a^b 2\pi f(t) \sqrt{g'(t)^2 + f'(t)^2} dt.
+$$
+
+If we assume integrability of the integrand, then as our partition size goes to zero, this approximate surface area converges to the value given by the limit.  (As with arc length, this needs a technical adjustment to the Riemann integral theorem as here we are evaluating the integrand function at four points ($t_i$, $t_{i-1}$, $\xi$ and $\psi$) and not just at some $c_i$.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+#| cache: true
+## {{{approximate_surface_area}}}
+
+xs,ys = range(-1, stop=1, length=50), range(-1, stop=1, length=50)
+f(x,y)= 2 - (x^2 + y^2)
+
+dr = [1/2, 3/4]
+df = [f(dr[1],0), f(dr[2],0)]
+
+function sa_approx_graph(i)
+    p = plot(xs, ys, f, st=[:surface], legend=false)
+    for theta in range(0, stop=i/10*2pi, length=10*i )
+        path3d!(p,sin(theta)*dr, cos(theta)*dr, df)
+    end
+    p
+end
+n = 10
+
+anim = @animate for i=1:n
+    sa_approx_graph(i)
+end
+
+imgfile = tempname() * ".gif"
+gif(anim, imgfile, fps = 1)
+
+
+caption = L"""
+
+Surface of revolution of $f(x) = 2 - x^2$ about the $y$ axis. The lines segments are the images of rotating the secant line connecting $(1/2, f(1/2))$ and $(3/4, f(3/4))$. These trace out the frustum of a cone which approximates the corresponding surface area of the surface of revolution. In the limit, this approximation becomes exact and a formula for the surface area of surfaces of revolution can be used to compute the value.
+
+"""
+
+ImageFile(imgfile, caption)
+```
+
+#### Examples
+
+
+Lets see that the surface area of an open cone follows from this formula, even though we just saw how to get this value.
+
+
+A cone be be envisioned as rotating the function $f(x) = x\tan(\theta)$ between $0$ and $h$ around the $x$ axis. This integral yields the surface area:
+
+
+$$
+\begin{align*}
+\int_0^h 2\pi f(x) \sqrt{1 + f'(x)^2}dx
+&= \int_0^h 2\pi x \tan(\theta) \sqrt{1 + \tan(\theta)^2}dx   \\
+&= (2\pi\tan(\theta)\sqrt{1 + \tan(\theta)^2} x^2/2 \big|_0^h \\
+&= \pi \tan(\theta) \sec(\theta) h^2                          \\
+&= \pi r^2 / \sin(\theta).
+\end{align*}
+$$
+
+(There are many ways to express this, we used $r$ and $\theta$ to match the work above. If the cone is parameterized by a height $h$ and radius $r$, then the surface area of the sides is $\pi r\sqrt{h^2 + r^2}$. If the base is included, there is an additional $\pi r^2$ term.)
+
+
+##### Example
+
+
+Let the graph of $f(x) = x^2$ from $x=0$ to $x=1$ be rotated around the $x$ axis. What is the resulting surface area generated?
+
+
+$$
+\text{SA} = \int_a^b 2\pi f(x) \sqrt{1 + f'(x)^2}dx = \int_0^1 2\pi x^2 \sqrt{1 + (2x)^2} dx.
+$$
+
+This integral is done by a trig substitution, but gets involved. We let `SymPy` do it:
+
+
+```{julia}
+@syms x
+F = integrate(2 * PI * x^2 * sqrt(1 + (2x)^2), x)
+```
+
+We show `F`, only to demonstrate that indeed the integral is a bit involved. The actual surface area follows from a *definite* integral, which we get through the fundamental theorem of calculus:
+
+
+```{julia}
+F(1) - F(0)
+```
+
+### Plotting surfaces of revolution
+
+
+The commands to plot a surface of revolution will be described more clearly later; for now we present them as simply a pattern to be followed in case plots are desired. Suppose the curve in the $x-y$ plane is given parametrically by $(g(u), f(u))$ for $a \leq u \leq b$.
+
+
+To be concrete, we parameterize the circle centered at $(6,0)$ with radius $2$ by:
+
+
+```{julia}
+g(u) = 6 + 2sin(u)
+f(u) = 2cos(u)
+a, b = 0, 2pi
+```
+
+The plot of this curve is:
+
+
+```{julia}
+#| hold: true
+us = range(a, b, length=100)
+plot(g.(us), f.(us), xlims=(-0.5, 9), aspect_ratio=:equal, legend=false)
+plot!([0,0],[-3,3], color=:red, linewidth=5)    # y axis emphasis
+plot!([3,9], [0,0], color=:green, linewidth=5)  # x axis emphasis
+```
+
+Though parametric plots have a convenience constructor, `plot(g, f, a, b)`, we constructed the points with `Julia`'s broadcasting notation, as we will need to do  for a surface of revolution. The `xlims` are adjusted to show the $y$ axis, which is emphasized with a layered line. The line is drawn by specifying two points, $(x_0, y_0)$ and $(x_1, y_1)$ in the form `[x0,x1]` and `[y0,y1]`.
+
+
+Now, to rotate this about the $y$ axis, creating a surface plot, we have the following pattern:
+
+
+```{julia}
+S(u,v) = [g(u)*cos(v), g(u)*sin(v), f(u)]
+us = range(a, b, length=100)
+vs = range(0, 2pi, length=100)
+ws = unzip(S.(us, vs'))   # reorganize data
+surface(ws..., zlims=(-6,6), legend=false)
+
+plot!([0,0], [0,0], [-3,3], color=:red, linewidth=5) # y axis emphasis
+```
+
+The `unzip` function is not part of base `Julia`, rather part of `CalculusWithJulia`. This function rearranges data into a form consumable by the plotting methods like `surface`. In this case, the result of `S.(us,vs')` is a grid (matrix) of points, the result of `unzip` is three grids of values, one for the $x$ values, one for the $y$ values, and one for the $z$ values.  A manual adjustment to the `zlims` is used, as `aspect_ratio` does not have an effect with the `plotly()` backend and errors on 3d graphics with `pyplot()`.
+
+
+To rotate this about the $x$ axis, we have this pattern:
+
+
+```{julia}
+#| hold: true
+S(u,v) = [g(u), f(u)*cos(v), f(u)*sin(v)]
+us = range(a, b, length=100)
+vs = range(0, 2pi, length=100)
+ws = unzip(S.(us,vs'))
+surface(ws..., legend=false)
+
+plot!([3,9], [0,0],[0,0], color=:green, linewidth=5) # x axis emphasis
+```
+
+The above pattern covers the case of  rotating the graph of a function $f(x)$ of $a,b$ by taking $g(t)=t$.
+
+
+##### Example
+
+
+Rotate the graph of $x^x$ from $0$ to $3/2$ around the $x$ axis. What is the surface area generated?
+
+
+We work numerically for this one, as no antiderivative is forthcoming. Recall, the accompanying `CalculusWithJulia` package defines `f'` to return the automatic derivative through the `ForwardDiff` package.
+
+
+```{julia}
+#| hold: true
+f(x) = x^x
+a, b = 0, 3/2
+val, _ = quadgk(x -> 2pi * f(x) * sqrt(1 + f'(x)^2), a, b)
+val
+```
+
+(The function is not defined at $x=0$ mathematically, but is on the computer to be $1$, the limiting value. Even were this not the case, the `quadgk` function doesn't evaluate the function at the points `a` and `b` that are specified.)
+
+
+```{julia}
+#| hold: true
+g(u) = u
+f(u) = u^u
+S(u,v) = [g(u)*cos(v), g(u)*sin(v), f(u)]
+us = range(0, 3/2, length=100)
+vs = range(0, pi, length=100)  # not 2pi (to see inside)
+ws = unzip(S.(us,vs'))
+surface(ws..., alpha=0.75)
+```
+
+We compare this answer to that of the frustum of a cone with radii $1$ and $(3/2)^2$, formed by rotating the line segment connecting $(0,f(0))$ with $(3/2,f(3/2))$. From looking at the graph of the surface, these values should be comparable. The surface area of the cone part is $\pi (r_1^2 + r_0^2) / \sin(\theta) = \pi (r_1 + r_0) \cdot \sqrt{(\Delta h)^2 + (r_1-r_0)^2}$.
+
+
+```{julia}
+#| hold: true
+f(x) = x^x
+r0, r1 = f(0), f(3/2)
+pi * (r1 + r0) * sqrt((3/2)^2 + (r1-r0)^2)
+```
+
+##### Example
+
+
+What is the surface area generated by Gabriel's Horn, the solid formed by rotating $1/x$ for $x \geq 1$ around the $x$ axis?
+
+
+$$
+\text{SA} = \int_a^b 2\pi f(x) \sqrt{1 + f'(x)^2}dx =
+\lim_{M \rightarrow \infty} \int_1^M 2\pi \frac{1}{x} \sqrt{1 + (-1/x^2)^2} dx.
+$$
+
+We do this with `SymPy`:
+
+
+```{julia}
+@syms M
+ex = integrate(2PI * (1/x) * sqrt(1 + (-1/x)^2), (x, 1, M))
+```
+
+The limit as $M$ gets large is of interest. The only term that might get out of hand is `asinh(M)`.  We check its limit:
+
+
+```{julia}
+limit(asinh(M), M => oo)
+```
+
+So indeed it does. There is nothing to balance this out, so the integral will be infinite, as this shows:
+
+
+```{julia}
+limit(ex, M => oo)
+```
+
+This figure would have infinite surface, were it possible to actually construct an infinitely long solid. (But it has been shown to have *finite* volume.)
+
+
+##### Example
+
+
+The curve described parametrically by $g(t) = 2(1 + \cos(t))\cos(t)$ and $f(t) = 2(1 + \cos(t))\sin(t)$ from $0$ to $\pi$ is rotated about the $x$ axis. Find the resulting surface area.
+
+
+The graph shows half a heart, the resulting area will resemble an apple.
+
+
+```{julia}
+#| hold: true
+g(t) = 2(1 + cos(t)) * cos(t)
+f(t) = 2(1 + cos(t)) * sin(t)
+plot(g, f, 0, 1pi)
+```
+
+The integrand simplifies to $8\sqrt{2}\pi \sin(t) (1 + \cos(t))^{3/2}$. This lends itself to $u$-substitution with $u=\cos(t)$.
+
+
+$$
+\begin{align*}
+\int_0^\pi 8\sqrt{2}\pi \sin(t) (1 + \cos(t))^{3/2}
+&= 8\sqrt{2}\pi \int_1^{-1} (1 + u)^{3/2} (-1) du\\
+&= 8\sqrt{2}\pi (2/5) (1+u)^{5/2} \big|_{-1}^1\\
+&= 8\sqrt{2}\pi (2/5) 2^{5/2} = \frac{2^7 \pi}{5}.
+\end{align*}
+$$
+
+## The first Theorem of Pappus
+
+
+The [first](http://tinyurl.com/le3lvb9) theorem of Pappus provides a simpler means to compute the surface area if the distance the centroid is from the axis ($\rho$) and the arc length of the curve ($L$) are both known. In that case, the surface area satisfies:
+
+
+$$
+\text{SA} = 2 \pi \rho L
+$$
+
+That is, the surface area is simply the circumference of the circle traced out by the centroid of the curve times the length of the curve - the distances rotated are collapsed to that of just the centroid.
+
+
+##### Example
+
+
+The surface area of of an open cone can be computed, as the arc length is $\sqrt{h^2 + r^2}$ and the centroid of the line is a distance $r/2$ from the axis. This gives SA$=2\pi (r/2) \sqrt{h^2 + r^2} = \pi r \sqrt{h^2 + r^2}$.
+
+
+##### Example
+
+
+We can get the surface area of a torus from this formula.
+
+
+The torus is found by rotating the curve $(x-b)^2 + y^2 = a^2$ about the $y$ axis. The centroid is $b$, the arc length $2\pi a$, so the surface area is $2\pi (b) (2\pi a) = 4\pi^2 a b$.
+
+
+A torus with $a=2$ and $b=6$
+
+
+```{julia}
+#| hold: true
+#| echo: false
+a,b = 2, 6
+F₀(u,v) = [a*(cos(u) + b)*cos(v), a*(cos(u) + b)*sin(v), a*sin(u)]
+us = vs = range(0, 2pi, length=35)
+ws = unzip(F₀.(us, vs'))
+surface(ws..., legend=false, zlims=(-12,12))
+```
+
+##### Example
+
+
+The surface area of sphere will be SA$=2\pi \rho (\pi r) = 2 \pi^2 r \cdot \rho$. What is $\rho$? The centroid of an arc formula can be derived in a manner similar to that of the centroid of a region. The formulas are:
+
+
+$$
+\begin{align}
+\text{cm}_x &= \frac{1}{L} \int_a^b g(t) \sqrt{g'(t)^2 + f'(t)^2} dt\\
+\text{cm}_y &= \frac{1}{L} \int_a^b f(t) \sqrt{g'(t)^2 + f'(t)^2} dt.
+\end{align}
+$$
+
+Here, $L$ is the arc length of the curve.
+
+
+For the sphere parameterized by $g(t) = r \cos(t)$, $f(t) = r\sin(t)$, we get that these become
+
+
+$$
+\text{cm}_x = \frac{1}{L}\int_0^\pi r\cos(t) \sqrt{r^2(\sin(t)^2 + \cos(t)^2)} dt =  \frac{1}{L}r^2 \int_0^\pi \cos(t) = 0.
+$$
+
+$$
+\text{cm}_y =  \frac{1}{L}\int_0^\pi r\sin(t) \sqrt{r^2(\sin(t)^2 + \cos(t)^2)} dt =  \frac{1}{L}r^2 \int_0^\pi \sin(t) = \frac{1}{\pi r} r^2 \cdot 2 = \frac{2r}{\pi}.
+$$
+
+Combining this, we see that the surface area of a sphere is $2 \pi^2 r (2r/\pi) = 4\pi r^2$, by Pappus' Theorem.
+
+
+## Questions
+
+
+##### Questions
+
+
+The graph of $f(x) = \sin(x)$ from $0$ to $\pi$ is rotated around the $x$ axis. After a $u$-substitution, what integral would  give the surface area generated?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices = [
+"``-\\int_1^{-1} 2\\pi \\sqrt{1 + u^2} du``",
+"``-\\int_1^{_1} 2\\pi u \\sqrt{1 + u^2} du``",
+"``-\\int_1^{_1} 2\\pi u^2 \\sqrt{1 + u} du``"
+]
+answ = 1
+radioq(choices, answ)
+```
+
+Though the integral can be computed by hand, give a numeric value.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+f(x) = sin(x)
+a, b = 0, pi
+val, _ = quadgk(x -> 2pi* f(x) * sqrt(1 + f'(x)^2), a, b)
+numericq(val)
+```
+
+##### Questions
+
+
+The graph of $f(x) = \sqrt{x}$ from $0$ to $4$ is rotated around the $x$ axis. Numerically find the surface area generated?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+f(x) = sqrt(x)
+a, b = 0, 4
+val, _ = quadgk(x -> 2pi* f(x) * sqrt(1 + f'(x)^2), a, b)
+numericq(val)
+```
+
+##### Questions
+
+
+Find the surface area generated by revolving the graph of the function $f(x) = x^3/9$ from $x=0$ to $x=2$ around the $x$ axis. This can be done by hand or numerically.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+f(x) = x^3/9
+a, b = 0, 2
+val, _ = quadgk(x -> 2pi* f(x) * sqrt(1 + f'(x)^2), a, b)
+numericq(val)
+```
+
+##### Questions
+
+
+(From Stewart.) If a loaf of bread is in the form of a sphere of radius $1$, the amount of crust for a slice depends on the width, but not where in the loaf it is sliced.
+
+
+That is this integral with $f(x) = \sqrt{1 - x^2}$ and $u, u+h$ in $[-1,1]$ does not depend on $u$:
+
+
+$$
+A = \int_u^{u+h} 2\pi f(x) \sqrt{1 + f'(x)^2} dx.
+$$
+
+If we let $f(x) = y$ then $f'(x) = x/y$. With this, what does the integral above come down to after cancellations:
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices = [
+"``\\int_u^{u+h} 2\\pi dx``",
+"``\\int_u^{u_h} 2\\pi y dx``",
+"``\\int_u^{u_h} 2\\pi x dx``"
+]
+answ = 1
+radioq(choices, answ)
+```
+
+##### Questions
+
+
+Find the surface area of the dome of sphere generated by rotating the the curve generated by $g(t) = \cos(t)$ and $f(t) = \sin(t)$ for $t$ in $0$ to $\pi/6$.
+
+
+Numerically find the value.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+g(t) = cos(t)
+f(t) = sin(t)
+a, b = 0, pi/4
+val, _ = quadgk(t -> 2pi* f(t) * sqrt(g'(t)^2 + f'(t)^2), a, b)
+numericq(val)
+```
+
+##### Questions
+
+
+The [astroid](http://www-history.mcs.st-and.ac.uk/Curves/Astroid.html) is parameterized by $g(t) = a\cos(t)^3$ and $f(t) = a \sin(t)^3$. Let $a=1$ and rotate the curve from $t=0$ to $t=\pi$ around the $x$ axis. What is the surface area generated?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+g(t) = cos(t^3)
+f(t) = sin(t^3)
+a, b = 0, pi
+val, _ = quadgk(t -> 2pi* f(t) * sqrt(g'(t)^2 + f'(t)^2), a, b)
+numericq(val)
+```
+
+##### Questions
+
+
+For the curve   parameterized by $g(t) = a\cos(t)^5$ and $f(t) = a \sin(t)^5$. Let $a=1$ and rotate the curve from $t=0$ to $t=\pi$ around the $x$ axis. Numerically find  the surface area generated?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+g(t) = cos(t^5)
+f(t) = sin(t^5)
+a, b = 0, pi
+val, _ = quadgk(t -> 2pi* f(t) * sqrt(g'(t)^2 + f'(t)^2), a, b)
+numericq(val)
+```
+

quarto/integrals/volumes_slice.qmd

@@ -0,0 +1,877 @@
+# Volumes by slicing
+
+
+```{julia}
+#| echo: false
+
+import Logging
+Logging.disable_logging(Logging.Info) # or e.g. Logging.Info
+Logging.disable_logging(Logging.Warn)
+
+import SymPy
+function Base.show(io::IO, ::MIME"text/html", x::T) where {T <: SymPy.SymbolicObject}
+    println(io, "<span class=\"math-left-align\" style=\"padding-left: 4px; width:0; float:left;\"> ")
+    println(io, "\\[")
+    println(io, sympy.latex(x))
+    println(io, "\\]")
+    println(io, "</span>")
+end
+
+# hack to work around issue
+import Markdown
+import CalculusWithJulia
+function CalculusWithJulia.WeaveSupport.ImageFile(d::Symbol, f::AbstractString, caption; kwargs...)
+    nm = joinpath("..", string(d), f)
+    u = "![$caption]($nm)"
+    Markdown.parse(u)
+end
+
+nothing
+```
+
+This section uses these add-on packages:
+
+
+```{julia}
+using CalculusWithJulia
+using Plots
+using QuadGK
+using Unitful, UnitfulUS
+using Roots
+using SymPy
+```
+
+```{julia}
+#| echo: false
+#| results: "hidden"
+using CalculusWithJulia.WeaveSupport
+
+import LinearAlgebra: norm
+
+const frontmatter = (
+        title = "Volumes by slicing",
+        description = "Calculus with Julia: Volumes by slicing",
+        tags = ["CalculusWithJulia", "integrals", "volumes by slicing"],
+);
+nothing
+```
+
+---
+
+
+```{julia}
+#| hold: true
+#| echo: false
+imgfile = "figures/michelin-man.jpg"
+caption = """
+
+Hey Michelin Man, how much does that costume weigh?
+
+"""
+ImageFile(:integrals, imgfile, caption)
+```
+
+An ad for a summer job says work as the Michelin Man! Sounds promising, but how much will that costume weigh? A very hot summer may make walking around in a heavy costume quite uncomfortable.
+
+
+A back-of-the envelope calculation would start by
+
+
+  * Mentally separating out each "tire" and lining them up one by one.
+  * Counting the number of "tires" (or rings), say $n$.
+  * Estimating the radius for each tire, say $r_i$ for $1 \leq i \leq n$.
+  * Estimating the height for each tire, say $h_i$ for $1 \leq i \leq n$
+
+
+Then the volume would be found by adding:
+
+
+$$
+V = \pi \cdot r_1^2 \cdot h_1 + \pi \cdot r_2^2 \cdot h_2 + \cdot + \pi \cdot r_n^2 \cdot h_n.
+$$
+
+The weight would come by multiplying the volume by some appropriate density.
+
+
+Looking at the sum though, we see the makings of an approximate integral. If the heights were to get infinitely small, we might expect this to approach something like $V=\int_a^b \pi r(h)^2  dh$.
+
+
+In fact, we have in general:
+
+
+> **Volume of a figure with a known cross section**: The volume of a solid with known cross-sectional area $A_{xc}(x)$ from $x=a$ to $x=b$ is given by
+>
+> $V = \int_a^b A_{xc}(x) dx.$
+>
+> This assumes $A_{xc}(x)$ is integrable.
+
+
+
+This formula is derived by approximating the volume by "slabs" with volume $A_{xc}(x) \Delta x$ and using the Riemann integral's definition to pass to the limit. The discs of the Michelin man are an example, where the cross-sectional area is just that of a circle, or $\pi r^2$.
+
+
+## Solids of revolution
+
+
+We begin with some examples of a special class of solids - solids of revolution. These have an axis of symmetry from which the slabs are then just circular disks.
+
+
+Consider the volume contained in this glass, it will depend on the radius at different values of $x$:
+
+
+```{julia}
+#| hold: true
+#| echo: false
+imgfile = "figures/integration-glass.jpg"
+caption = L"""
+
+A wine glass oriented so that it is seen as generated by revolving a
+curve about the $x$ axis. The radius of revolution varies as a function of $x$
+between about $0$ and $6.2$cm.
+
+"""
+ImageFile(:integrals, imgfile, caption)
+```
+
+If $r(x)$ is the radius as a function of $x$, then the cross sectional area is $\pi r(x)^2$ so the volume is given by:
+
+
+$$
+V = \int_a^b \pi r(x)^2 dx.
+$$
+
+:::{.callout-note}
+## Note
+The formula is for a rotation around the $x$-axis, but can easily be generalized to rotating around any line (say the $y$-axis or $y=x$, ...) just by adjusting what $r(x)$ is taken to be.
+
+:::
+
+For a numeric example, we consider the original Red [Solo](http://en.wikipedia.org/wiki/Red_Solo_Cup) Cup. The dimensions of the cup were basically: a top diameter of $d_1 = 3~ \frac{3}{4}$ inches, a bottom diameter of $d_0 = 2~ \frac{1}{2}$ inches and a height of $h = 4~ \frac{3}{4}$ inches.
+
+
+The central axis is straight down. If we rotate the cup so this is the $x$-axis, then we can get
+
+
+$$
+r(x) = \frac{d_0}{2} + \frac{d_1/2 - d_0/2}{h}x = \frac{5}{4} + \frac{5}{38}x
+$$
+
+The volume in cubic inches will be:
+
+
+$$
+V = \int_0^h \pi r(x)^2 dx
+$$
+
+This is
+
+
+```{julia}
+d0, d1, h = 2.5, 3.75, 4.75
+rad(x) = d0/2 + (d1/2 - d0/2)/h * x
+vol, _ = quadgk(x -> pi * rad(x)^2, 0, h)
+```
+
+So $36.9 \text{in}^3$. How many ounces is that? It is useful to know that 1 [gallon](http://en.wikipedia.org/wiki/Gallon) of water is defined as $231$ cubic inches, contains $128$ ounces, and weighs $8.34$ pounds.
+
+
+So our cup holds this many ounces:
+
+
+```{julia}
+ozs = vol / 231 * 128
+```
+
+Full it is about $20$ ounces, though this doesn't really account for the volume taken up by the bottom of the cup, etc.
+
+
+If you are poor with units, `Julia` can provide some help through the `Unitful` package. Here the additional `UnitfulUS` package must also be included, as was done above, to access fluid ounces:
+
+
+```{julia}
+vol * u"inch"^3 |> us"floz"
+```
+
+Before Solo "squared" the cup, the Solo cup had markings that - [some thought](http://www.snopes.com/food/prepare/solocups.asp) - indicated certain volume amounts.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+imgfile = "figures/red-solo-cup.jpg"
+caption = "Markings on the red Solo cup indicated various volumes"
+ImageFile(:integrals, imgfile, caption)
+```
+
+What is the height for $5$ ounces (for a glass of wine)? $12$ ounces (for a beer unit)?
+
+
+Here the volume is fixed, but the height is not. For $v$ ounces, we need to convert to cubic inches. The conversion is $1$ ounce is $231/128 \text{in}^3$.
+
+
+So we need to solve $v \cdot (231/128) = \int_0^h\pi r(x)^2 dx$ for $h$ when $v=5$ and $v=12$.
+
+
+Let's express volume as a function of $h$:
+
+
+```{julia}
+Vol(h) = quadgk(x -> pi * rad(x)^2, 0, h)[1]
+```
+
+Then to solve we have:
+
+
+```{julia}
+v₅ = 5
+h5 = find_zero(h -> Vol(h) - v₅ * 231 / 128, 4)
+```
+
+and
+
+
+```{julia}
+v₁₂ = 12
+h12 = find_zero(h -> Vol(h) - v₁₂ * 231 / 128, 4)
+```
+
+As a percentage of the total height, these are:
+
+
+```{julia}
+h5/h, h12/h
+```
+
+:::{.callout-note}
+## Note
+Were performance at issue, Newton's method might also have been considered here, as the derivative is easily computed by the fundamental theorem of calculus.
+
+:::
+
+##### Example
+
+
+By rotating the line segment $x/r + y/h=1$ that sits in the first quadrant around the $y$ axis, we will generate a right-circular cone. The volume of which can be expressed through the above formula by noting the radius, as a function of $y$, will be $R = r(1 - y/h)$. This gives the well-known volume of a cone:
+
+
+```{julia}
+#| hold: true
+@syms r h x y
+R = r*(1 - y/h)
+integrate(pi*R^2, (y, 0, h))
+```
+
+The frustum of a cone is simply viewed as a cone with its top cut off. If the original height would have been $h_0$ and the actual height $h_1$, then the volume remaining is just $\int_{h_0}^h \pi r(y)^2 dy = \pi h_1 r^2/3 - \pi h_0 r^2/3 = \pi r^2 (h_1-h_0)/3$.
+
+
+It is not unusual to parameterize a cone by the angle $\theta$ it makes and the height. Since $r/h=\tan\theta$, this gives the formula $V = \pi/3\cdot h^3\tan(\theta)^2$.
+
+
+##### Example
+
+
+[Gabriel's](http://tinyurl.com/8a6ygv) horn is a geometric figure of mathematics - but not the real world - which has infinite height, but not volume! The figure is found by rotating the curve $y=1/x$ around the $x$ axis from $1$ to $\infty$. If the volume formula holds, what is the volume of this "horn?"
+
+
+```{julia}
+radius(x) = 1/x
+quadgk(x -> pi*radius(x)^2, 1, Inf)[1]
+```
+
+That is a value very reminiscent of $\pi$, which it is as $\int_1^\infty 1/x^2 dx = -1/x\big|_1^\infty=1$.
+
+
+:::{.callout-note}
+## Note
+The interest in this figure is that soon we will be able to show that it has **infinite** surface area, leading to the [paradox](http://tinyurl.com/osawwqm) that it seems possible to fill it with paint, but not paint the outside.
+
+:::
+
+##### Example
+
+
+A movie studio hand is asked to find a prop vase to be used as a [Ming vase](http://en.wikipedia.org/wiki/Chinese_ceramics) in an upcoming scene. The dimensions specified are for the outside diameter in centimeters and are given by
+
+
+$$
+d(h) = \begin{cases}
+2 \sqrt{26^2 - (h-20)^2} & 0 \leq h \leq 44\\
+20 \cdot e^{-(h - 44)/10} & 44 < h \leq 50.
+\end{cases}
+$$
+
+If the vase were solid, what would be the volume?
+
+
+We define `d` using a ternary operator to handle the two cases:
+
+
+```{julia}
+d(h) = h <= 44 ? 2*sqrt(26^2 - (h-20)^2) : 20 * exp(-(h-44)/10)
+rad(h) = d(h)/2
+```
+
+The volume in cm$^3$ is then:
+
+
+```{julia}
+Vₜ, _ = quadgk(h -> pi * rad(h)^2, 0, 50)
+```
+
+For the actual shoot, the vase is to be filled with ash, to simulate a funeral urn. (It will then be knocked over in a humorous manner, of course.) How much ash is needed if the vase has walls that are 1/2 centimeter thick
+
+
+We need to subtract $0.5$ from the radius and `a` then recompute:
+
+
+```{julia}
+V_int, _ = quadgk(h -> pi * (rad(h) - 1/2)^2, 1/2, 50)
+```
+
+A liter of volume is $1000 \text{cm}^3$. So this is about $68$ liters, or more than 15 gallons. Perhaps the dimensions given were bit off.
+
+
+While we are here, to compute the actual volume of the material in the vase could be done by subtraction.
+
+
+```{julia}
+Vₜ - V_int
+```
+
+### The washer method
+
+
+Returning to the Michelin Man, in our initial back-of-the-envelope calculation we didn't account for the fact that a tire isn't a disc, as it has its center cut out. Returning, suppose $R_i$ is the outer radius and $r_i$ the inner radius. Then each tire has volume
+
+
+$$
+\pi R_i^2 h_i - \pi r_i^2 h_i = \pi (R_i^2 - r_i^2) h_i.
+$$
+
+Rather than use $\pi r(x)^2$ for a cross section, we would use $\pi (R(x)^2 - r(x)^2)$.
+
+
+In general we call a shape like the tire a "washer" and use this formula for a washer's cross section $A_{xc}(x) = \pi(R(x)^2 - r(x)^2)$.
+
+
+Then the volume for the solid of revolution whose cross sections are washers would be:
+
+
+$$
+V = \int_a^b \pi \cdot (R(x)^2 - r(x)^2) dx.
+$$
+
+##### Example
+
+
+An artist is working with a half-sphere of material, and wishes to bore out a conical shape. What would be the resulting volume, if the two figures are modeled by
+
+
+$$
+R(x) = \sqrt{1^2 - (x-1)^2}, \quad r(x) = x,
+$$
+
+with $x$ ranging from $x=0$ to $1$?
+
+
+The answer comes by integrating:
+
+
+```{julia}
+#| hold: true
+Rad(x) = sqrt(1 - (x-1)^2)
+rad(x) = x
+V, _ = quadgk(x -> pi*(Rad(x)^2 - rad(x)^2), 0, 1)
+```
+
+## Solids with known cross section
+
+
+The Dart cup company now produces the red solo cup with a [square](http://www.solocup.com/products/squared-plastic-cup/) cross section. Suppose the dimensions are the same: a top diameter of $d_1 = 3 3/4$ inches, a bottom diameter of $d_0 = 2 1/2$ inches and a height of $h = 4 3/4$ inches. What is the volume now?
+
+
+The difference, of course, is that cross sections now have area $d^2$, as opposed to $\pi r^2$. This leads to some difference, which we quantify, as follows:
+
+
+```{julia}
+#| hold: true
+d0, d1, h = 2.5, 3.75, 4.75
+d(x) = d0 + (d1 - d0)/h * x
+vol, _ = quadgk(x -> d(x)^2, 0, h)
+vol / 231 * 128
+```
+
+This shape would have more volume - the cross sections are bigger. Presumably the dimensions have changed. Without going out and buying a cup, let's assume the cross-sectional diameter remained the same, not the diameter. This means the largest dimension is the same. The cross section diameter is $\sqrt{2}$ larger. What would this do to the area?
+
+
+We could do this two ways: divide $d_0$ and $d_1$ by $\sqrt{2}$ and recompute. However, each cross section of this narrower cup would simply be $\sqrt{2}^2$ smaller, so the total volume would change by $2$, or be 13 ounces. We have $26.04$ is too big, and $13.02$ is too small, so some other overall dimensions are used.
+
+
+##### Example
+
+
+For a general cone, we use this [definition](http://en.wikipedia.org/wiki/Cone):
+
+
+> A cone is the solid figure bounded by a base in a plane and by a surface (called the lateral surface) formed by the locus of all straight line segments joining the apex to the perimeter of the base.
+
+
+
+Let $h$ be the distance from the apex to the base. Consider cones with the property that all planes parallel to the base intersect the cone with the same shape, though perhaps a different scale. This figure shows an example, with the rays coming from the apex defining the volume.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+h = 5
+R, r, rho = 1, 1/4, 1/4
+f(t) = (R-r) * cos(t) + rho * cos((R-r)/r * t)
+g(t) = (R-r) * sin(t) - rho * sin((R-r)/r * t)
+ts = range(0, 2pi, length=100)
+
+p = plot(f.(ts), g.(ts), zero.(ts), legend=false)
+for t ∈ range(0, 2pi, length=25)
+    plot!(p, [0,f(t)], [0,g(t)], [h, 0], linecolor=:red)
+end
+p
+```
+
+A right circular cone is one where this shape is a circle. This definition can be more general, as a square-based right pyramid is also such a cone. After possibly reorienting the cone in space so the base is at $u=0$ and the apex at $u=h$ the volume of the cone can be found from:
+
+
+$$
+V = \int_0^h A_{xc}(u) du.
+$$
+
+The cross sectional area $A_{xc}(u)$ satisfies a formula in terms of $A_{xc}(0)$, the area of the base:
+
+
+$$
+A_{xc}(u)  = A_{xc}(0) \cdot (1 - \frac{u}{h})^2
+$$
+
+So the integral becomes:
+
+
+$$
+V = \int_0^h A_{xc}(u) du = A_{xc}(0) \int_0^h (1 - \frac{u}{h})^2 du = A_{xc}(0) \int_0^1 v^2 \frac{1}{h} dv = A_{xc}(0) \frac{h}{3}.
+$$
+
+This gives a general formula for the volume of such cones.
+
+
+### Cavalieri's method
+
+
+[Cavalieri's](http://tinyurl.com/oda9xd9) Principle is "Suppose two regions in three-space (solids) are included between two parallel planes. If every plane parallel to these two planes intersects both regions in cross-sections of equal area, then the two regions have equal volumes." (Wikipedia).
+
+
+With the formula for the volume of solids based on cross sections, this is a trivial observation, as the functions giving the cross-sectional area are identical. Still, it can be surprising. Consider a sphere with an interior cylinder bored out of it.  (The [Napkin](http://tinyurl.com/o237v83) ring problem.) The bore has height $h$ - for larger radius spheres this means very wide bores.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+#The following illustrates $R=5$ and $h=8$.
+
+R =5; h1 = 2*4
+
+theta = asin(h1/2/R)
+thetas = range(-theta, stop=theta, length=100)
+ts = range(-pi, stop=pi, length=100)
+y = h1/4
+
+p = plot(legend=false, aspect_ratio=:equal);
+plot!(p, R*cos.(ts), R*sin.(ts));
+plot!(p, R*cos.(thetas), R*sin.(thetas), color=:orange);
+
+plot!(p, [R*cos.(theta), R*cos.(theta)], [h1/2, -h1/2], color=:orange);
+plot!(p, [R*cos.(theta), sqrt(R^2 - y^2)], [y, y], color=:orange)
+
+plot!(p, [0, R*cos.(theta)], [0,0], color=:red);
+plot!(p,[ 0, R*cos.(theta)], [0,h1/2], color=:red);
+
+annotate!(p, [(.5, -2/3, "sqrt(R²- (h/2)²)"),
+              (R*cos.(theta)-.6, h1/4, "h/2"),
+              (1.5, 1.75*tan.(theta), "R")])
+
+p
+
+```
+
+The small orange line is rotated, so using the washer method we get the cross sections given by $\pi(r_0^2 - r_i^2)$, the outer and inner radii, as a function of $y$.
+
+
+The outer radii has points $(x,y)$ satisfying $x^2 + y^2 = R^2$, so is $\sqrt{R^2 - y^2}$. The inner radii has a constant value, and as indicated in the figure, is $\sqrt{R^2 - (h/2)^2}$, by the Pythagorean theorem.
+
+
+Thus the cross sectional area is
+
+
+$$
+\pi( (\sqrt{R^2 - y^2})^2 - (\sqrt{R^2 - (h/2)^2})^2 )
+= \pi ((R^2 - y^2) - (R^2 - (h/2)^2))
+= \pi ((\frac{h}{2})^2 - y^2)
+$$
+
+As this does not depend on $R$, and the limits of integration would always be $-h/2$ to $h/2$ by Cavalieri's principle, the volume of the solid will be independent of $R$ too.
+
+
+To actually compute this volume, we take $R=h/2$, so that the bore hole is just a line of no volume, the resulting volume is then that of a sphere with radius $h/2$, or $4/3\pi(h/2)^3 = \pi h^3/6$.
+
+
+## The second theorem of Pappus
+
+
+The second theorem of [Pappus](http://tinyurl.com/l43vw4) says that if a plane figure $F$ is rotated around an axis to form a solid of revolution, the total volume can be written as $2\pi r A(F)$, where $r$ is the distance the centroid is from the axis of revolution, and $A(F)$ is the area of the plane figure. In short, the distance traveled by the centroid times the area.
+
+
+This can make some computations trivial. For example, we can make a torus (or donut) by rotating the circle $(x-2)^2 + y^2 = 1$ about the $y$ axis. As the centroid is clearly $(2, 0)$, with $r=2$ in the above formula, and the area of the circle is $\pi 1^2$, the volume of the donut is $2\pi(2)(\pi) = 4\pi^2$.
+
+
+##### Example
+
+
+Above, we found the volume of a cone, as it is a solid of revolution, through the general formula. However, parameterizing the cone as the revolution of a triangle with vertices $(0,0)$, $(r, 0)$, and $(0,h)$ and using the formula for the center of mass in the $x$ direction of such a triangle, $r/3$, we get that the volume of a cone with height $h$ and radius $r$ is  $2\pi (r/3)\cdot (rh/2) = \pi r^2 h/3$, in agreement with the calculus based computation.
+
+
+## Questions
+
+
+###### Question
+
+
+Consider this big Solo cup:
+
+
+```{julia}
+#| hold: true
+#| echo: false
+imgfile ="figures/big-solo-cup.jpg"
+caption = " Big solo cup. "
+ImageFile(:integrals, imgfile, caption)
+```
+
+It has approximate dimensions: smaller radius 5 feet, upper radius 8 feet and height 15 feet. How many gallons is it? At $8$ pounds a gallon this would be pretty heavy!
+
+
+Two facts are useful:
+
+
+  * a cubic foot is 7.48052 gallons
+  * the radius as a function of height is $r(h) = 5  + (3/15)\cdot h$
+
+
+```{julia}
+#| hold: true
+#| echo: false
+gft = 7.48052
+rad(h) = 5 + (3/15)*h
+a,err = quadgk(h -> pi*rad(h)^2, 0, 15)
+val = a*gft
+numericq(val, 1e1)
+```
+
+###### Question
+
+
+In *Glass Shape Influences Consumption Rate* for Alcoholic [Beverages](http://www.plosone.org/article/info%3Adoi%2F10.1371%2Fjournal.pone.0043007) the authors demonstrate that the shape of the glass can have an effect on the rate of consumption, presumably people drink faster when they aren't sure how much they have left. In particular, they comment that people have difficulty judging the half-finished-by-volume mark.
+
+
+This figure shows some of the wide variety of beer-serving glasses:
+
+
+```{julia}
+#| hold: true
+#| echo: false
+imgfile ="figures/beer_glasses.jpg"
+caption = "A variety of different serving glasses for beer."
+ImageFile(:integrals, imgfile, caption)
+```
+
+We work with metric units, as there is a natural relation between volume in cm$^3$ and liquid measure ($1$ liter = $1000$ cm$^3$, so a $16$-oz pint glass is roughly $450$ cm$^3$.)
+
+
+Let two glasses be given as follows. A typical pint glass with linearly increasing radius:
+
+
+$$
+r(h) = 3 + \frac{1}{5}h, \quad 0 \leq h \leq b;
+$$
+
+and a curved-edge one:
+
+
+$$
+s(h) = 3 + \log(1 + h), \quad 0 \leq h \leq b
+$$
+
+The following functions find the volume as a function of height, $h$:
+
+
+```{julia}
+r1(h) = 3 + h/5
+s1(h) = 2 + log(1 + h)
+r_vol(h) = quadgk(x -> pi*r1(x)^2, 0, h)[1]
+s_vol(h) = quadgk(x -> pi*s1(x)^2, 0, h)[1]
+```
+
+  * For the straight-sided glass find $h$ so that the volume is $450$.
+
+
+```{julia}
+#| echo: false
+h450 = find_zero(h -> r_vol(h) - 450, 10)
+numericq(h450)
+```
+
+  * For the straight-sided glass find $h$ so that the volume is $225$ (half full).
+
+
+```{julia}
+#| echo: false
+h225 = find_zero(h -> r_vol(h) - 225, 10)
+numericq(h225)
+```
+
+  * For the straight-sided glass, what is the percentage of the total height when the glass is half full. (For a cylinder it would just be 50.
+
+
+```{julia}
+#| echo: false
+numericq(h225/450 * 100, 2, units="percent")
+```
+
+  * People often confuse the half-way by height amount for the half way by volume, as it is for the cylinder. Take the height for the straight-sided glass filled with $450$ mm, divide it by $2$, then compute the percentage of volume at the half way height to the original.
+
+
+```{julia}
+#| echo: false
+numericq(r_vol(h450/2)/450*100, 2, units="percent")
+```
+
+---
+
+
+  * For the curved-sided glass find $h$ so that the volume is $450$.
+
+
+```{julia}
+#| echo: false
+h_450 = find_zero(h -> s_vol(h) - 450, 10)
+numericq(h_450)
+```
+
+  * For the curved-sided glass find $h$ so that the volume is $225$ (half full).
+
+
+```{julia}
+#| echo: false
+h_225 = find_zero(h -> s_vol(h) - 225, 10)
+numericq(h_225)
+```
+
+  * For the curved-sided glass, what is the percentage of the total height when the glass is half full. (For a cylinder it would just be 50.
+
+
+```{julia}
+#| echo: false
+numericq(h_225/450 * 100, 2, units="percent")
+```
+
+  * People often confuse the half-way by height amount for the half way by volume, as it is for the cylinder. Take the height for the curved-sided glass filled with $450$ mm, divide it by $2$, then compute the percentage of volume at the half way height to the original.
+
+
+```{julia}
+#| echo: false
+numericq(s_vol(h_450/2)/450*100, 2, units="percent")
+```
+
+###### Question
+
+
+A right [pyramid](http://en.wikipedia.org/wiki/Pyramid_%28geometry%29) has its apex (top point) above the centroid of its base, and for our purposes, each of its cross sections. Suppose a pyramid has square base of dimension $w$ and height of dimension $h$.
+
+
+Will this integral give the volume:
+
+
+$$
+V = \int_0^h w^2 (1 - \frac{y}{h})^2 dy?
+$$
+
+```{julia}
+#| hold: true
+#| echo: false
+yesnoq("yes")
+```
+
+What is the volume?
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices = [
+"``1/3 \\cdot b\\cdot h``",
+"``1/3 \\cdot w^2\\cdot h``",
+"``l\\cdot w \\cdot h/ 3``"
+]
+answ = 2
+radioq(choices, answ)
+```
+
+###### Question
+
+
+An ellipsoid is formed by rotating the region in the first and second quadrants bounded by the ellipse $(x/2)^2 + (y/3)^2=1$ and the $x$ axis around the $x$ axis. What is the volume of this ellipsoid? Find it numerically.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+f(x) = 3*sqrt( 1 - (x/2)^2 )
+val, _ = quadgk(x -> pi * f(x)^2, -2, 2)
+numericq(val)
+```
+
+###### Question
+
+
+An ellipsoid is formed by rotating the region in the first and second quadrants bounded by the ellipse $(x/a)^2 + (y/b)^2=1$ and the $x$ axis around the $x$ axis. What is the volume of this ellipsoid? Find it symbolically.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+choices = [
+"``4/3 \\cdot \\pi a b^2``",
+"``4/3 \\cdot \\pi a^2 b``",
+"``\\pi/3 \\cdot a b^2``"
+]
+answ = 1
+radioq(choices, answ)
+```
+
+###### Question
+
+
+A solid is generated by rotating the region enclosed by the graph $y=\sqrt{x}$, the lines $x=1$, $x=2$, and $y=1$ about the $x$ axis. Find the volume of the solid.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+Ra(x) = sqrt(x)
+ra(x) = 1
+a,b=1,2
+val, _ = quadgk(x -> pi * (Ra(x)^2 - ra(x)^2), a,b)
+numericq(val)
+```
+
+###### Question
+
+
+The region enclosed by the graphs of $y=x^3 - 1$ and $y=x-1$ are rotated around the $y$ axis. What is the volume of the solid?
+
+
+```{julia}
+#| hold: true
+@syms x
+plot(x^3 - 1, 0, 1, legend=false)
+plot!(x-1)
+```
+
+```{julia}
+#| hold: true
+#| echo: false
+Ra(y) = cbrt(y+1)
+ra(y) = y + 1
+a,b = 0, 1
+val, _ = quadgk(x -> pi * (Ra(x)^2 - ra(x)^2), a,b)
+numericq(val)
+```
+
+###### Question
+
+
+Rotate the region bounded by $y=e^x$, the line $x=\log(2)$ and the first quadrant about the line $x=\log(2)$.
+
+
+(Be careful, the radius in the formula $V=\int_a^b \pi r(u)^2 du$ is from the line $x=\log(2)$.)
+
+
+```{julia}
+#| hold: true
+#| echo: false
+a, b = 0, exp(log(2))
+ra(y) = log(2) - log(y)
+val, _ = quadgk(y -> pi * ra(y)^2, a, b)
+numericq(val)
+```
+
+###### Question
+
+
+Find the volume of rotating the region bounded by the line $y=x$, $x=1$ and the $x$-axis around the line $y=x$. (The Theorem of Pappus is convenient and the fact that the centroid of the triangular region lies at $(2/3, 1/3)$.)
+
+
+```{julia}
+#| hold: true
+#| echo: false
+cm=[2/3, 1/3]
+c = [1/2, 1/2]
+rr = norm(cm - c)
+A = 1/2 * 1 * 1
+val = 2pi * rr * A
+numericq(val)
+```
+
+###### Question
+
+
+Rotate the region bounded by the line $y=x$ and the function $f(x) = x^2$ about the line $y=x$. What is the resulting volume?
+
+
+You can integrate in the length along the line $y=x$ ($u$ from $0$ to $\sqrt{2}$). The radius then can be found by intersecting the line perpendicular line to $y=x$ at $u$ to the curve $f(x)$. This will do so:
+
+
+```{julia}
+theta = pi/4  ## we write y=x as y = x * tan(pi/4) for more generality, as this allows other slants.
+
+f(x) = x^2
+𝒙(u) = find_zero(x -> u*sin(theta) - 1/tan(theta) * (x - u*cos(theta)) - f(x), (u*cos(theta), 1))
+𝒓(u) = sqrt((u*cos(theta) - 𝒙(u))^2 + (u*sin(theta) - f(𝒙(u)))^2)
+```
+
+(Though in this case you can also find `r(u)` using the quadratic formula.)
+
+
+With this, find the volume.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+a, b = 0, sqrt(2)
+val, _ = quadgk(u -> pi*𝒓(u)^2, a, b)
+numericq(val)
+```
+
+---
+
+
+Repeat (find the volume) only this time with the function $f(x) = x^{20}$.
+
+
+```{julia}
+#| hold: true
+#| echo: false
+a, b = 0, sqrt(2)
+f(x) = x^20
+xval(u) = find_zero(x -> u*sin(theta) - 1/tan(theta) * (x - u*cos(theta)) - f(x), (0,sqrt(2)))
+rad(u) = sqrt((u*cos(theta) - xval(u))^2 + (u*sin(theta) - f(xval(u)))^2)
+val, _ = quadgk(u -> pi*rad(u)^2, a, b)
+numericq(val)
+```
+