The definition of the definite integral, $\int_a^b f(x)dx$, is based on Riemann sums.
We review, using a more general form than [previously](../integrals/area.html). Consider a bounded function $f$ over $[a,b]$. A partition, $P$, is based on $a = x_0 < x_1 < \cdots < x_n = b$. For each subinterval $[x_{i-1}, x_{i}]$ take $m_i(f) = \inf_{u \text{ in } [x_{i-1},x_i]} f(u)$ and $M_i(f) = \sup_{u \text{ in } [x_{i-1},x_i]} f(u)$. (When $f$ is continuous, $m_i$ and $M_i$ are realized at points of $[x_{i-1},x_i]$, though that isn't assumed here. The use of "$\sup$" and "$\inf$" is a mathematically formal means to replace this in general.) Let $\Delta x_i = x_i - x_{i-1}$. Form the sums $m(f, P) = \sum_i m_i(f) \Delta x_i$ and $M(f, P) = \sum_i M_i(f) \Delta x_i$. These are the *lower* and *upper* Riemann sums for a partition. A general Riemann sum would be formed by selecting $c_i$ from $[x_{i-1}, x_i]$ and forming $S(f,P) = \sum f(c_i) \Delta x_i$. It will be the case that $m(f,P) \leq S(f,P) \leq M(f,P)$, as this is true for *each* sub-interval of the partition.
If, as the largest diameter ($\Delta x_i$) of the partition $P$ goes to $0$, the upper and lower sums converge to the same limit, then $f$ is called Riemann integrable over $[a,b]$. If $f$ is Riemann integrable, any Riemann sum will converge to the definite integral as the partitioning shrinks.
Continuous functions are known to be Riemann integrable, as are functions with only finitely many discontinuities, though this isn't the most general case of integrable functions, which will be stated below.
In practice, we don't typically compute integrals using a limit of a partition, though the approach may provide direction to numeric answers, as the Fundamental Theorem of Calculus relates the definite integral with an antiderivative of the integrand.
The multidimensional case will prove to be similar where a Riemann sum is used to define the value being discussed, but a theorem of Fubini will allow the computation of integrals using the Fundamental Theorem of Calculus.
caption = """How to estimate the volume contained within the Chrysler Building? One way might be to break the building up into tall vertical blocks based on its skyline; compute the volume of each block using the formula of volume as area of the base times the height; and, finally, adding up the computed volumes This is the basic idea of finding volumes under surfaces using Riemann integration."""
Computing the volume of a nano-block construction of the Chrysler building is easier than trying to find an actual tree at the Chrysler building, as we can easily compute the volume of columns of equal-sized blocks. Riemann sums are similar.
Computing the volume of a nano-block construction of the Chrysler building is easier than trying to find an actual tree at the Chrysler building, as we can easily compute the volume of columns of equal-sized blocks. Riemann sums are similar.
The definition of the multi-dimensional integral is more involved then the one-dimensional case due to the possibly increased complexity of the region. This will require additional [steps](https://math.okstate.edu/people/lebl/osu4153-s16/chapter10-ver1.pdf). The basic approach is as follows.
First, let $R = [a_1, b_1] \times [a_2, b_2] \times \cdots \times [a_n, b_n]$ be a closed rectangular region. If $n=2$, this is a rectangle, and if $n=3$, a box. We begin by defining integration over closed rectangular regions. For each side, a partition $P_i$ is chosen based on $a_i = x_{i0} < x_{i1} < \cdots < x_{ik} = b_i$. Then a sub-rectangular region would be of the form $R' = P_{1j_1} \times P_{2j_2} \times \cdots \times P_{nj_n}$, where $P_{ij_i}$ is one of the partitioning sub intervals of $[a_i, b_i]$. Set $\Delta R' = \Delta P_{1j_1} \cdot \Delta P_{2j_2} \cdot\cdots\cdot\Delta P_{nj_n}$ to be the $n$-dimensional volume of the sub-rectangular region.
For each sub-rectangular region, we can define $m(f,R')$ to be $\inf_{u \text{ in } R'} f(u)$ and $M(f, R') = \sup_{u \text{ in } R'} f(u)$. If we enumerate all the sub-rectangular regions, we can define $m(f, P) = \sum_i m(f, R_i) \Delta R_i$ and $M(f,P) = \sum_i M(f, R_i)\Delta R_i$, as in the one-dimensional case. These are upper and lower sums, and, as before, would bound the Riemann sum formed by choosing any $c_i$ in $R_i$ and computing $S(f,P) = \sum_i f(c_i) \Delta R_i$.
As with the one-dimensional case, $f$ is Riemann integrable over $R$ if the limits of $m(f,P)$ and $M(f,P)$ exist and are identical as the diameter of the partition (defined as the largest diameter of each side) goes to $0$. If the limits are equal, then so is the limit of any Riemann sum.
When $f$ is Riemann integrable over a rectangular region $R$, we denote the limit by any of:
> Any continuous function, $f$, is Riemann integrable over a closed, bounded rectangular region.
---
As with one-dimensional integrals, from the Riemann sum definition, several familiar properties for integrals follow. Let $V(R)$ be the volume of $R$ found by multiplying the side-lengths together.
**Constants:**
* A constant is Riemann integrable and: $\iint_R c dV = c V(R)$.
**Linearity:**
* For integrable $f$ and $g$ and constants $a$ and $b$:
* As $f$ is bounded, let $m \leq f(x) \leq M$ for all $x$ in $R$. Then
$$
m V(R) \leq \iint_R f(x) dV \leq MV(R).
$$
* If $f$ and $g$ are integrable *and* $f(x) \leq g(x)$, then the integrals have the same property, namely $\iint_R f dV \leq \iint_R gdV$.
* If $S \subset R$, both closed rectangles, then if $f$ is integrable over $R$ it will be also over $S$ and, when $f\geq 0$, $\iint_S f dV \leq \iint_R fdV$.
**Triangle inequality:**
* If $f$ is bounded and integrable, then $|\iint_R fdV| \leq \iint_R |f| dV$.
### HCubature
To numerically compute multidimensional integrals over rectangular regions in `Julia` is efficiently done with the `HCubature` package. The `hcubature` function is defined for $n$-dimensional integrals, so the integrand is specified through a function which takes a vector as an input. The region to integrate over is of rectangular form. It is specified by a tuple of left endpoints and a tuple of right endpoints. The order is in terms of the order of the vector.
To elaborate, if we think of $f(\vec{x}) = f(x_1, x_2, \dots, x_n)$ and we are integrating over $[a_1, b_1] \times \cdots \times [a_n, b_n]$, then the region would be specified through two tuples: `(a1, a2, ..., an)` and `(b1, b2, ..., bn)`.
To illustrate, to integrate the function $f(x,y) = x^2 + 5y^2$ over the region $[0,1] \times [0,2]$ using `HCubature`'s `hcubature` function, we would proceed as follows:
```{julia}
#| hold: true
f(x,y) = x^2 + 5y^2
f(v) = f(v...) # f accepts a vector
a0, b0 = 0, 1
a1, b1 = 0, 2
hcubature(f, (a0, a1), (b0, b1))
```
The computed value and a worst case estimate for the error is returned, in a manner similar to the `quadgk` function (from the `QuadGK` package) used previously for one-dimensional numeric integrals.
The order above is `x` then `y`, which is clear from the first definition of `f` and as belabored in the tuples passed to `hcubature`. A more convenient use is to just put the constants into the function call, as in `hcubature(f, (0,0), (1,2))`.
##### Example
Let's verify the numeric approach works for figures where an answer is known from the geometry of the problem.
* A constant function $c=f(x,y)$. In this case, the volume is simply a box, so the volume will come from multiplying the three dimensions. Here is an example:
```{julia}
#| hold: true
f(x,y) = 3
f(v) = f(v...)
a0, b0 = 0, 4
a1, b1 = 0, 5 # R is area 20, so V = 60 = 3 ⋅ 20
hcubature(f, (a0, a1), (b0, b1))
```
* A wedge. Let $f(x,y) = x$ and $R= [0,1] \times [0,1]$. The the volume is a wedge, and should be half the value of the unit cube, or simply $1/2$:
```{julia}
#| hold: true
f(x,y) = x
f(v) = f(v...)
a0, b0 = 0, 1
a1, b1 = 0, 1
hcubature(f, (a0, a1), (b0, b1))
```
* The volume of a right square pyramid is $V=(1/3)a^2 h$, or a third of an enclosing box. We computed this area previously using the method of [slices](../integrals/volumes_slice.html). Here we do it thinking of the pyramid as the volume formed by the surface over the region $[-a,a] \times [-a,a]$ generated by $f(x,y) = h \cdot (l(x,y) - d(x,y))/l(x,y)$ where $d(x,y)$ is the distance to the origin, or $\sqrt{x^2 + y^2}$ and $l(x,y)$ is the length of the line segment from the origin to the boundary of $R$ that goes through $(x,y)$.
Identifying a formula for this is a bit tricky. Here we use a brute force approach; later we will simplify this. Using polar coordinates, we know $r\cos(\theta) = a$ describes the line $x=a$ and $r\sin(\theta)=a$ describes the line $y=a$. Using the square, we have to alternate between these depending on where $\theta$ is (e.g., between $-\pi/4$ and $\pi/4$ it would be $r\cos(\theta)=a$ or $a/\cos(\theta)$ is $l(x,y)$. We write a function for this:
```{julia}
𝒅(x, y) = sqrt(x^2 + y^2)
function 𝒍(x, y, a)
theta = atan(y,x)
atheta = abs(theta)
if (pi/4 <= atheta < 3pi/4) # this is the y=a or y=-a case
(a/2)/sin(atheta)
else
(a/2)/abs(cos(atheta))
end
end
```
And then
```{julia}
𝒇(x,y,a,h) = h * (𝒍(x,y,a) - 𝒅(x,y))/𝒍(x,y,a)
𝒂, 𝒉 = 2, 3
𝒇(x,y) = 𝒇(x, y, 𝒂, 𝒉) # fix a and h
𝒇(v) = 𝒇(v...)
```
We can visualize the volume to be computed, as follows:
```{julia}
#| hold: true
xs = ys = range(-1, 1, length=20)
surface(xs, ys, 𝒇)
```
Trying this, we have:
```{julia}
hcubature(𝒇, (-𝒂/2, -𝒂/2), (𝒂/2, 𝒂/2))
```
The answer agrees with that known from the formula, $4 = (1/3)a^2 h$, but the answer takes a long time to be produce. The `hcubature` function is slow with functions defined in terms of conditions. For this problem, volumes by [slicing](../integrals/volumes_slice.html) is more direct. But also symmetry can be used, were we able to compute the volume above the triangular region formed by the $x$-axis, the line $x=a/2$ and the line $y=x$, which would be $1/8$th the total volume. (As then $l(x,y,a) = (a/2)/\sin(\tan^{-1}(y,x))$.).
* The volume of a sphere is $4/3 \pi r^3$. We could verify this by integrating $z = f(x,y) = \sqrt{r^2 - (x^2 + y^2)}$ over $R = \{(x,y): x^2 + y^2 \leq r^2\}$. *However*, this is not a *rectangular* region, so we couldn't directly proceed.
We might try integrating a function with a condition:
```{julia}
#| hold: true
function f(x, y, r)
if x^2 + y^2 < r
sqrt(z - x^2 + y^2)
else
0.0
end
end
```
**But** `hcubature` is **very** slow to integrate such functions. We will see our instincts are good – this is the approach taken to discuss integrals over general regions – but this is not practical here. There are two alternative approaches to be discussed: approach the integral *iteratively* or *transform* the circular region into a rectangular region and integrate. Before doing so, we discuss how the integral is developed for more general regions.
:::{.callout-note}
## Note
The approach above takes a nice smooth function and makes it non smooth at the boundary. In general this is not a good idea for numeric solutions, as many algorithms work better with assumptions of smoothness.
:::
:::{.callout-note}
## Note
The `Quadrature` package provides a uniform interface for `QuadGK`, `HCubature`, and other numeric integration routines available in `Julia`.
:::
## Integrals over more general regions
To proceed further, it is necessary to discuss certain types of sets that will be used to describe the boundaries of regions that can be integrated over, though we don't dig into the details.
Let the *measure* of a rectangular region be its volume and for any subset of $S \subset R^n$, define the *outer* measure of $S$ by $m^*(S) = \inf\sum_{j=1}^\infty V(R_j)$ where the infimum is taken over all closed, countable, rectangles with $S \subset \cup_{j=1}^\infty R_j$.
In two dimensions, if $S$ is viewed on a grid, then this would be *area* of the smallest collection of cells that contain any part of $S$. This is the smallest this value takes as the grid becomes infinite.
For the following graph, there are $100$ cells each of area $8/100$. Their are 58 cells covering the curve and its interior. So the outer measure is less than $58\cdot 8/100$, as this is just one possible covering.
```{julia}
#| hold: true
#| echo: false
function cassini(theta)
a, b = .75, 1
A = 1; B = -2a^2*cos(2theta)
C = a^4 - b^4
(-B - sqrt(B^2 - 4*A*C))/(2A)
end
polar_plot(r, a, b) = plot(t -> r(t)*cos(t), t->r(t)*sin(t), a, b, legend=false, linewidth=3)
p = polar_plot(cassini, 0, 8pi)
n=10
a1,b1 = -1, 1
a2, b2 = -2, 2
for a in range(a1, b1, length=n+1)
for b in range(a2, b2, length=n+1)
plot!(p, [a,a],[a2, b2], alpha=0.75)
plot!(p, [a1,b1],[b,b], alpha=0.75)
end
end
p
```
A set has measure $0$ if the outer measure is $0$. An alternate definition, among other characterizations, is a set has measure $0$ if for any $\epsilon > 0$ there exists rectangular regions $R_1, R_2, \dots, R_n$ (for some $n$) with $\sum V(R_i) < \epsilon$. Measure zero sets have many properties not discussed here.
For now, let's see that graph of $y=f(x)$ over $[a,b]$, as a two dimensional set, has measure zero when $f(x)$ has a bounded derivative ($|f'|$ bounded by $M$). Fix some $\epsilon>0$. Take $n$ with $2M(b-a)^2/n < \epsilon$, then divide $[a,b]$ into $n$ equal length intervals (of length $\delta = (b-a)/n)$. For each interval, we consider the box $[a_i, b_i] \times [f(a_i)-\delta M, f(a_i) + \delta M]$. By the mean value theorem, we have $|f(x) - f(a_i)| \leq |b_i-a_i|M$ so $f(a_i) - \delta M \leq f(x) \leq f(a_i) + \delta M$, so the curve will stay in the boxes. These boxes have total area $n \cdot \delta \cdot 2\delta M = 2M(b-a)^2/n$, an area less than $\epsilon$.
The above can be extended to any graph of a continuous function over $[a,b]$.
For a function $f$ the set of discontinuities in $R$ is all points where $f$ is not continuous. A formal definition is often given in terms of oscillation. Let $o(f, \vec{x}, \delta) = \sup_{\{\vec{y} : \| \vec{y}-\vec{x}\| < \delta\}}f(\vec{y}) - \inf_{\{\vec{y}: \|\vec{y}-\vec{x}\|<\delta\}}f(\vec{y})$. A function is discontinuous at $\vec{x}$ if the limit as $\delta \rightarrow 0+$ (which must exist) is not $0$.
With this, we can state the Riemann-Lebesgue theorem on integrable functions:
> Let $R$ be a closed, rectangular region, and $f:R^n \rightarrow R$ a bounded function. Then $f$ is Riemann integrable over $R$ if and only if the set of discontinuities is a set of measure $0$.
It was said at the outset we would generalize the regions we can integrate over, but this theorem generalizes the functions. We can tie the two together as follows. Define the integral over any *bounded* set $S$ with boundary of measure $0$. Bounded means $S$ is contained in some bounded rectangle $R$. Let $f$ be defined on $S$ and extend it to be $0$ on points in $R$ that are not in $S$. If this extended function is integrable over $R$, then we can define the integral over $S$ in terms of that. This is why the *boundary* of $S$ must have measure zero, as in general it is among the set of discontinuities of the extend function $f$. Such regions are also called Jordan regions.
## Fubini's theorem
Consider again this figure
```{julia}
#| hold: true
#| echo: false
function cassini(theta)
a, b = .75, 1
A = 1; B = -2a^2*cos(2theta)
C = a^4 - b^4
(-B - sqrt(B^2 - 4*A*C))/(2A)
end
polar_plot(r, a, b) = plot(t -> r(t)*cos(t), t->r(t)*sin(t), a, b, legend=false, linewidth=3)
p = polar_plot(cassini, 0, 8pi)
n=10
a1,b1 = -1, 1
a2, b2 = -2, 2
for a in range(a1, b1, length=n+1)
for b in range(a2, b2, length=n+1)
plot!(p, [a,a],[a2, b2], alpha=0.75)
plot!(p, [a1,b1],[b,b], alpha=0.75)
end
end
p
```
Let $C_i$ enumerate all the cells shown, assume $f$ is extended to be $0$ outside the region, and let $c_i$ be a point in the cell. Then the Riemann sum $\sum_i f(c_i) V(C_i)$ can be visualized three identical ways:
* as a linear sum over the indices $i$, as written, leading to $\iint_R f(x) dV$.
* by indexing the cells by row ($i$) and column ($j$) and summing as $\sum_i (\sum_j f(x_{ij}, y_{ij}) \Delta y_j) \Delta x_i$.
* by indexing the cells by row ($i$) and column ($j$) and summing as $\sum_j (\sum_i f(x_{ij}, y_{ij}) \Delta x_i) \Delta y_j$.
The last two suggest that their limit will be *iterated* integrals of the form $\int_{-1}^1 (\int_{-2}^2 f(x,y) dy) dx$ and $\int_{-2}^2 (\int_{-1}^1 f(x,y) dx) dy$.
By "iterated" we mean performing two different definite integrals. For example, to compute $\int_{-1}^1 (\int_{-2}^2 f(x,y) dy) dx$ the first task would be to compute $I(x) = \int_{-2}^2 f(x,y) dy$. Like partial derivatives, this integrates in $y$ while treating $x$ as a constant. Once the interior integral is computed, then the integral $\int_{-1}^1 I(x) dx$ would be computed to find the answer.
The question then: under what conditions will the three integrals be equal?
> [Fubini](https://math.okstate.edu/people/lebl/osu4153-s16/chapter10-ver1.pdf). Let $R \times S$ be a closed rectangular region in $R^n \times R^m$. Suppose $f$ is bounded. Define $f_x(y) = f(x,y)$ and $f^y(x) = f(x,y)$ where $x$ is in $R^n$ and $y$ in $R^m$. *If* $f_x$ and $f^y$ are integrable then
Similarly, if $f^y$ is integrable for all $y$, then $\iint_{R\times S}fdV =\iint_S \iint_R f(x,y) dx dy$.
An immediate corollary is that the above holds for continuous functions when $R$ and $S$ are bounded, the case described here.
The case of continuous functions was known to [Euler](https://en.wikipedia.org/wiki/Fubini%27s_theorem#History), Lebesgue (1904) discussed bounded functions, as in our statement, and Fubini and Tonnelli (1907 and 1909) generalized the statement to more general functions than continuous functions, thereby earning naming rights.
In [Ferzola](https://doi.org/10.2307/2687130) we can read a summary of Euler's thinking of 1769 when trying to understand the integral of a function $f(x,y)$ over a bounded domain $R$ enclosed by arcs in the $x$-$y$ plane. (That is, the area below $g(x)$ and above $h(x)$ over the interval $[a,b]$.) Euler wrote the answer as $\int_a^b dx (\int_{g(x)}^{h(x)} f(x,y)dy)$. Ferzola writes that Euler saw this integral yielding a *volume* as the integral $\int_{g(x)}^{h(x)} f(x,y)dy$ gives the area of a slice (parallel to the $y$ axis) and integrating in $x$ adds these slices to give a volume. This is the typical usage of Fubini's theorem today.
```{julia}
#| hold: true
#| echo: false
imgfile ="figures/strang-slicing.png"
caption = L"""Figure 14.2 of Strang illustrating the slice when either $x$ is fixed or $y$ is fixed. The inner integral computes the shared area, the outer integral adds the areas up to compute volume."""
In [Volumes](../integrals/volumes_slice.html) the formula for a volume with a known cross-sectional area is given by $V = \int_a^b CA(x) dx$. The inner integral, $\int_{R_x} f(x,y) dy$ is a function depending on $x$ that yields the area of the slice (where $R_x$ is the region sliced by the line of constant $x$ value). This is consistent with Euler's view of the iterated integral.
A domain, as described above, is known as a [normal](https://en.wikipedia.org/wiki/Multiple_integral#Normal_domains_on_R2) domain. Using Fubini's theorem to integrate iteratively, employing the fundamental theorem of calculus at each step, is the standard approach.
For example, we return to the problem of a square pyramid, only now using symmetry, we integrate only over the triangular region between $0 \leq x \leq a/2$ and $0 \leq y \leq x$. The answer is then (the $8$ by symmetry)
$$
V = 8 \int_0^{a/2} \int_0^x h(l(x,y) - d(x,y))/l(x,y) dy dx.
$$
But, using similar triangles, we have $d/x = l/(a/2)$ so $(l-d)/l = 1 - 2x/a$. Continuing, our answer becomes
$$
V = 8 \int_0^{a/2} (\int_0^x h(1-\frac{2x}{a}) dy) dx =
The `integrate` function of `SymPy` uses various algorithms to symbolically integrate definite (and indefinite) integrals. In the section on [integrals](../integrals/ftc.html) its use for one-dimensional integrals was shown. For multi-dimensional integrals the usage is similar, the syntax following, somewhat, the Fubini-like notation.
For example, to perform the integral
$$
\int_a^b \int_{h(x)}^{g(x)} f(x,y) dy dx
$$
the call would look like:
```{julia}
#| eval: false
integrate(f(x,y), (y, h(x), g(x)), (x, a, b))
```
That is, the variable to integrate and the endpoints are passed as tuples. (Unlike `hcubature` which always uses two tuples to specify the bounds, `integrate` uses $n$ tuples to specify an $n$-dimensional integral.) The iteration happens from left to write, so in the above the `y` integral is done (and, as seen, may depend on the variable `x`) and then the `x` integral is performed. The above uses `f(x,y)`, `h(x)` and `g(x)`, but these may be simple symbolic expressions and not function calls using symbolic variables.
We define `x` and `y` below for use throughout:
```{julia}
@syms x::real y::real z::real
```
##### Example
For example, the last integral to compute the volume of a square pyramid, could be computed through
Find the integral $\int_0^1\int_{y^2}^1 y \sin(x^2) dx dy$.
Without concerning ourselves with what or why, we just translate:
```{julia}
#| hold: true
integrate( y * sin(x^2), (x, y^2, 1), (y, 0, 1))
```
##### Example
Find the volume enclosed by $y = x^2$, $y = 5$, $z = x^2$, and $z = 0$.
The limits on $z$ say this is the volume under the surface $f(x,y) = x^2$, over the region defined by $y=5$ and $y = x^2$. The region is a parabola with $y$ running from $x^2$ to $5$, while $x$ ranges from $-\sqrt{5}$ to $\sqrt{5}$.
* bounded by $x+y+z = 10$, $2x + 3y = 20$, and $x + 3y = 10$
The first plane can be expressed as $z = f(x,y) = 10 - x - y$ and the volume is that below the surface of $f$ over the region $R$ formed by the two lines and the $x$ and $y$ axes. Plotting that we have:
```{julia}
#| hold: true
g1(x) = (20 - 2x)/3
g2(x) = (10 - x)/3
plot(g1, 0, 20)
plot!(g2, 0, 20)
```
We see the intersection is when $x=10$, so this becomes
```{julia}
#| hold: true
f(x,y) = 10 - x - y
h(x) = (10 - x)/3
g(x) = (20 - 3x)/3
integrate(f(x,y), (y, h(x), g(x)), (x, 0, 10))
```
##### Example
Let $r=1$ and define three cylinders along the $x$, $y$, and $z$ axes by: $y^2+z^2 = r^2$, $x^2 + z^2 = r^2$, and $x^2 + y^2 = r^2$. What is the enclosed [volume](http://mathworld.wolfram.com/SteinmetzSolid.html)?
Using the cylinder along the $z$ axis, we have the volume sits above and below the disk $R = x^2 + y^2 \leq r^2$. By symmetry, we can double the volume that sits above the disk to answer the question.
Using symmetry, we can tell that the the wedge between $x=0$, $y=x$, and $x^2 + y^2 \leq 1$ (corresponding to a polar angle in $[0,\pi/4]$ in $R$ contains $1/8$ the volume of the top, so $1/16$ of the total.
```{julia}
#| hold: true
#| echo: false
rad(theta) = 1
plot(t -> rad(t)*cos(t), t -> rad(t)*sin(t), 0, pi/4, legend=false, linewidth=3)
Over this wedge the height is given by the cylinder along the $y$ axis, $x^2 + z^2 = r^2$. We *could* break this wedge into a triangle and a semicircle to integrate piece by piece. However, from the figure we can integrate in the $y$ direction on the outside, and use only one intergral:
```{julia}
#| hold: true
r = 1 # if using r as a symbolic variable specify `positive=true`
f(x,y) = sqrt(r^2 - x^2)
16 * integrate(f(x,y), (x, y, sqrt(r^2-y^2)), (y, 0, r*cos(PI/4)))
```
##### Example
Find the volume under $f(x,y) = xy$ in the cone swept out by $r(\theta) = 1$ as $\theta$ goes between $[0, \pi/4]$.
The region $R$, the same as the last one. As seen, it can be described in two pieces as a function of $x$, but needs only $1$ as a function of $y$, so we use that below:
The average value of a function, $f(x,y)$, over a region $R$ is the integral of $f$ over $R$ divided by the area of $R$. It can be computed through two integrals, as below.
let $R$ be the region in the first quadrant bounded by $x - y = 0$ and $f(x,y) = x^2 + y^2$. Find the average value.
```{julia}
#| hold: true
f(x,y) = x^2 + y^2
g(x) = x # solve x - y = 0 for y
h(x) = 0
A = integrate(f(x,y), (y, h(x), g(x)), (x, 0, 1))
B = integrate(Sym(1), (y, h(x), g(x)), (x, 0, 1))
A/B
```
(We integrate `Sym(1)` and not just `1`, as we either need to have a symbolic value for the first argument or use the `sympy.integrate` method directly.)
##### Example: Density
The area of a region $R$ can be computed by $\iint_R 1 dA$. If the region is physical, say a disc, then its mass can be of interest. If the mass is uniform with density $\rho$, then the mass would be $\iint_R \rho dA$. If the mass is non uniform, say it is a function $\rho(x,y)$, then the integral to find the mass becomes $\iint_R \rho(x,y) dA$. (In a Riemann sum, the term $\rho(c_{ij}) \Delta x_i\Delta y_j$ would be the mass of a constant-density solid, the integral just adds these up to find total mass.)
Find the mass of a disc bounded by the two parabolas $y=2 - x^2$ and $y = -3 + 2x^2$ with density function given by $\rho(x,y) = x^2y^2$.
First we need the intersection points of the two parabolas. Solving $2-x^2 = -3 + 2x^2$ for $x$ yields: $5 = x^2$.
So we get a mass of:
```{julia}
#| hold: true
rho(x,y) = x^2*y^2
g(x) = 2 - x^2
h(x) = -3 + 2x^2
a = sqrt(Sym(5))
integrate(rho(x,y), (y, h(x), g(x)), (x, -a, a))
```
##### Example (Strang)
Integrate $\int_0^1 \int_y^1 \cos(x^2) dx dy$ avoiding the *impossible* integral of $\cos(x^2)$. As the integrand is continuous, Fubini's Theorem allows the interchange of the variable of integraton. The region, $R$, is a triangle in the first quadrant below the line $y=x$ and left of the line $x=1$. So we have:
$$
\int_0^1 \int_0^x \cos(x^2) dy dx
$$
We can integrate this, as the interior integral leaves $x \cos(x^2)$ to integrate:
```{julia}
integrate(cos(x^2), (y, 0, x), (x, 0, 1))
```
### A "Fubini" function
The computationally efficient way to perform multiple integrals numerically would be to use `hcubature`. However, this function is defined only for *rectangular* regions. In the event of non-rectangular regions, the suggested performant way would be to find a suitable transformation (below).
However, for simple problems, where ease of expressing a region is preferred to computational efficiency, something can be implemented using repeated uses of `quadgk`. Again, this isn't recommended, save for its relationship to how iteration is approached algebraically.
In the `CalculusWithJulia` package, the `fubini` function is provided. For these notes, we define three operations using Unicode operators entered with `\int[tab]`, `\iint[tab]`, `\iiint[tab]`. (Using this, better shows the mechanics involved.)
```{julia}
# adjust endpoints when expressed as a functions of outer variables
callf(f::Number, x) = f
callf(f, x) = f(x...)
endpoints(ys, x) = callf.(ys, Ref(x))
# integrate f(x) dx
∫(@nospecialize(f), xs) = quadgk(f, xs...)[1] # @nospecialize is not necessary, but offers a speed boost
The one issue with the above is the tolerances should increase for each outer integration. Otherwise, a certain precision is expected for a noisier integrand. This isn't an issue for many integrands, but can become one.
Compare the integral of $f(x,y) = \exp(-x^2 -2y^2)$ over the region $R=[0,3]\times[0,3]$ using `hcubature` and the above.
```{julia}
#| hold: true
f(x,y) = exp(-x^2 - 2y^2)
f(v) = f(v...)
hcubature(f, (0,0), (3,3)) # (a0, a1), (b0, b1)
```
```{julia}
#| hold: true
f(x,y) = exp(-x^2 - 2y^2)
∬(f, (0,3), (0,3)) # (a1, b1), (a0, b0)
```
##### Example
Show the area of the unit circle is $\pi$ using the "Fubini" function.
```{julia}
#| hold: true
f(x,y) = 1
a = ∬(f, (x-> -sqrt(1-x^2), x-> sqrt(1-x^2)), (-1, 1))
a, a - pi # answer and error
```
(The error is similar to that returned by `quadgk(x -> sqrt(1-x^2), -1, 1)`.)
###### Example
Show the volume of a sphere of radius $1$ is $4/3\pi = 4/3\pi\cdot 1^3$ by doubling the integral of $f(x,y) = \sqrt{1-x^2-y^2}$ over $R$, the unit disk.
Numeric integrals don't need to worry about integrands without antiderivatives. Their concerns are highly oscillatory integrands. Here we compute $\int_0^1 \int_y^1 \cos(x^2) dx dy$ directly. The limits are in a different order than the "Fubini" function expects, so we switch the variables:
```{julia}
∬((y,x) -> cos(x^2), (y -> y, 1), (0, 1))
```
Compare to
```{julia}
sin(1)/2
```
## Triple integrals
Triple integrals are identical in theory to double integrals, though the computations can be more involved and the regions more complicated to describe. The main regions (emphasized by Strang) to understand are: box, prism, cylinder, cone, tetrahedron, and sphere.
```{julia}
#| echo: false
let
ts = range(0, pi/2, length=50)
O = [0,0,0]
bx, by,bz = [1, 0, 0], [0,2,0], [0,0,3]
p = plot(unzip([O])..., legend=false, title="box",
axis=nothing)
arrow!(p, O, bx), arrow!(p, O, by), arrow!(p, O, bz)
p = plot(unzip([O])..., legend=false, title="tetrahedron",
axis=nothing)
arrow!(p, O, bx), arrow!(p, O, by), arrow!(p, O, bz)
arrow!(p, bx, by-bx)
arrow!(p, bx, bz-bx), arrow!(p, by, bz-by)
push!(ps, p)
p = plot(unzip([O])..., legend=false, title="cone",
camera=(70,20), axis=nothing)
arrow!(p, O, bx), arrow!(p, O, by/2), arrow!(p, O, bz)
arrow!(p, bx, bz-bx), arrow!(p, by/2, bz-by/2)
for h in range(0.1, 0.9, length=5)
z = 3*h
r = 1 - h
plot!(p, r*cos.(ts), r*sin.(ts), z .+ 0*ts)
end
push!(ps, p)
p = plot(unzip([O])..., legend=false, title="sphere",
camera=(70,20), axis=nothing)
Os = 0 * ts
plot!(p, cos.(ts), sin.(ts), Os)
plot!(p, cos.(ts), Os, sin.(ts))
plot!(p, Os, cos.(ts), sin.(ts))
for h in .2:.2:.8
r = sqrt(1 - h^2)
plot!(p, r*cos.(ts), r*sin.(ts), h .+ Os)
end
push!(ps, p)
l = @layout [a b; c d; e]
plot(ps..., layout=l)
end
```
Here we compute the volumes of these using a triple integral of the form $\iint_R 1 dV$.
* Box. Consider the box-like, or "rectangular," region $[0,a]\times [0,b] \times [0,c]$. This has volume $abc$ which we see here using Fubini's theorem:
```{julia}
#| hold: true
@syms a b c
f(x,y,z) = Sym(1) # need to integrate a symbolic object in integrand or call `sympy.integrate`
Which, as expected, is half the volume of the box $[0,c] \times [0, 1/a] \times [0, 1/b]$.
* Tetrahedron. Consider the volume formed by $x,y,z \geq 0$ and bounded by $ax+by+cz = 1$ where $a,b,c \geq 0$. The volume is a tetrahedron. The base in the $x$-$y$ plane is a triangle with vertices $(1/a, 0, 0)$ and $(0, 1/b, 0)$.
(The third easy-to-find point is $(0, 0, 1/c)$). The line connecting the points in the $x$-$y$ plane is $ax + by = 1$. With this, the integral to compute the volume is
* Cone. Consider a cone formed by the function $z = f(x,y) = a - b(x^2+y^2)^{1/2}$ ($a,b > 0$) and the $x$-$y$ plane. This will have radius $r = a/b$ and height $a$. The volume is given by this integral:
This integral is doable, but `SymPy` has trouble with it. We will return to this when cylindrical coordinates are defined.
* Sphere. The sphere $x^2 + y^2 + z^2 \leq 1$ has a known volume. Can we compute it using integration? In Cartesian coordinates, we can describe the region $x^2 + y^2 \leq 1$ and then the $z$-limits will follow:
This integral is doable, but `SymPy` has trouble with it. We will return to this when spherical coordinates are defined.
## Change of variables
The change of variables, or substitution, formula from first-semester calculus is expressed, under assumptions, by:
$$
\int_{g(R)} f(x) dx = \int_R (f\circ g)(u)g'(u) du.
$$
The derivation comes from reversing the chain rule. When using it, we start on the right hand side and typically write $x = g(u)$ and from here derive an expression involving differentials: $dx = g'(u) du$ and the rest follows. In practice, this is used to simplify the integrand in the search for an antiderivative, as $(f\circ g)$ is generally more complicated than $f$ alone.
In higher dimensions, we will see that change of variables can not only simplify the integrand, but is also of great use to simplify the region to integrate over. We mentioned, for example, that to use `hcubature` efficiently over a non-rectangular region, a transformation–-or change of variables–-is needed. The key to the multi-dimensional formula is understanding what should replace $dx = g'(u) du$. We take a bit of a circuitous route to get there.
In [Katz](http://www.jstor.org/stable/2689856) a review of the history of "change of variables" from Euler to Cartan is given. We follow Lagrange's formal analysis to derive the change of variable formula in two dimensions.
We view $R$ in two coordinate systems $(x,y)$ and $(u,v)$. We have that
where $A = \partial{x}/\partial{u}$, $B = \partial{x}/\partial{v}$, $C= \partial{y}/\partial{u}$, and $D = \partial{y}/\partial{v}$. Lagrange, following Euler, first sets $x$ to be constant (as is done in iterated integration). Hence, $dx = 0$ and so $du = -C(B/A) dv$ and, after substitution, $dy = (D-C(B/A))dv$. Then Lagrange set $y$ to be a constant, so $dy = 0$ and hence $dv=0$ so $dx = Adu$. The area "element" $dx dy = A du \cdot (D - (B/A)) dv = (AD - BC) du dv$. Since areas and volumes are non-negative, the absolute value is used. With this, we have "$dxdy = |AD-BC|du dv$" as the analog of $dx = g'(u) du$.
The expression $AD - BC$ was also derived by Euler, by related means. Lagrange extended the analysis to 3 dimensions. Before doing so, it is helpful to understand the problem from a geometric perspective. Euler was attempting to understand the effects of the following change of variable:
Euler knew this to be a clockwise *rotation* by an angle $\theta$ with $\cos(\theta) = m$, a *reflection* through the $x$ axis, and a translation by $\langle a, b\rangle$. All these *should* preserve the area represented by $dx dy$, so he was *expecting* $dx dy = dt dv$.
```{julia}
#| hold: true
#| echo: false
imgfile ="figures/euler-rotation.png"
caption = "Figure from Katz showing rotation of Euler."
The figure, taken from Katz, shows the translation, and rotation that should preserve area on a differential scale.
However Euler knew $dx = (\partial{g_1}/\partial{t}) dt + (\partial{g_1}/\partial{v}) dv$ and $dy = (\partial{g_2}/{\partial{t}}) dt + (\partial{g_2}/\partial{v}) dv$. Just multiplying gives $dx dy = m\sqrt{1-m^2} dt dt + (1-m^2) dv dt -m^2 dt dv -m\sqrt{1-m^2} dv dv$, a result that didn't make sense physically as $dt dt$ and $dv dv$ have no meaning in integration and $1 - m^2 - m^2$ is not $1$ as expected. Euler, like Lagrange, used a formal trick to proceed, but the geometric insight that the incremental areas for a change of variable should be related and for this change of variable identical is correct.
The following illustrates the polar-coordinate transformation $\langle x,y\rangle = G(r, \theta) = r \langle \cos\theta, \sin\theta\rangle$.
plot!(unzip(G.(us, vs'))...) # plots constant v lines
pt = [us[n],vs[n]]
arrow!(G(pt), J(pt)*[1,0], color=:blue)
arrow!(G(pt), J(pt)*[0,1], color=:blue)
```
This graphic shows the image of the box $[0,1] \times [0, 2\pi]$ under the transformation. The `plot` commands draw lines for values of constant `u` or constant `v`. If $G(u,v) = \langle g_1(u,v), g_2(u,v)\rangle$, then the Taylor expansion for $g_i$ is $g_i(u+du, v+dv) \approx g_i(u,v) + (\nabla{g_i})^T \cdot \langle du, dv \rangle$ and combining $G(u+du, v+dv) \approx G(u,v) + J_G(u,v) \langle du, dv \rangle$. The vectors added above represent the images when $u$ is constant (so $du=0$) and when $v$ is constant (so $dv=0$). The two arrows define a parallelogram whose area gives the change of area undergone by the unit square under the transformation. The area is $|\det(J_G)|$, the absolute value of the determinant of the Jacobian.
```{julia}
#| echo: false
function showG(G, a=1, b=1;a0=0, b0=0, an = 3, bn=3, n=5, lambda=1/2, k1=1, k2=1)
p = plot(unzip(G.(us', vs))..., legend = false, aspect_ratio=:equal) # plots constant u lines
plot!(p,unzip(G.(us, vs'))...) # plots constant v lines
pt = [us[k1 * n],vs[k2*n]]
P, U, V = G(pt...), lambda * J(pt)*[1,0], lambda * J(pt)*[0,1]
arrow!(P, U, color=:blue, linewidth=2)
arrow!(P+V, U, color=:red, linewidth=1)
arrow!(P, V, color=:blue, linewidth=2)
arrow!(P+U, V, color=:red, linewidth=1)
p
end
```
The tranformation to elliptical coordinates, $G(u,v) = \langle \cosh(u)\cos(v), \sinh(u)\sin(v)\rangle$, may be viewed similarly:
```{julia}
#| hold: true
#| echo: false
G(u,v) = [cosh(u)*cos(v), sinh(u)*sin(v)]
showG(G, 1, 2pi)
```
The transformation $G(u,v) = v \langle e^u, e^{-u} \rangle$ uses hyperbolic coordinates:
```{julia}
#| hold: true
#| echo: false
G(u,v) = v * [exp(u), exp(-u)]
showG(G, 1, 2pi, bn = 6, k2=4)
```
The transformation $G(u,v) = \langle u^2-v^2, u\cdot v \rangle$ yields a partition of the plane:
```{julia}
#| hold: true
#| echo: false
G(u,v) = [u^2 - v^2, u*v]
showG(G, 1, 1)
```
The arrows are the images of the standard unit vectors. We see some transformations leave these *orthogonal* and some change the respective lengths. The area of the associated parallelogram can be found using the determinant of an accompanying matrix. For two dimensions, using the cross product formulation on the embedded vectors, the area is
Using the fact that the two vectors involved are columns in the Jacobian of the transformation, this is just $|\det(J_G)|$. For $3$ dimensions, the determinant gives the volume of the 3-dimensional parallelepiped in the same manner. This holds for higher dimensions.
The absolute value of the determinant of the Jacobian is the multiplying factor that is seen in the change of variable formula for all dimensions:
> [Change of variable](https://en.wikipedia.org/wiki/Integration_by_substitution#Substitution_for_multiple_variables) Let $U$ be an open set in $R^n$, $G:U \rightarrow R^n$ be an *injective* differentiable function with *continuous* partial derivatives. If $f$ is continuous and compactly supported, then
>
> $$
> \iint_{G(S)} f(\vec{x}) dV = \iint_S (f \circ G)(\vec{u}) |\det(J_G)(\vec{u})| dU.
> $$
For the one-dimensional case, there is no absolute value, but there the interval is reversed, producing "negative" area. This is not the case here, where $S$ is parameterized to give positive volume.
:::{.callout-note}
## Note
The term "functional determinant" is found for the value $\det(J_G)$, as is the notation $\partial(x_1, x_2, \dots x_n)/\partial(u_1, u_2, \dots, u_n)$.
:::
### Two dimensional change of variables
Now we see several examples of two-dimensional transformations.
#### Polar integrals
We have [seen](../differentiable_vector_calculus/polar_coordinates.html) how to compute area in polar coordinates through the formula $A = \int (1/2) r^2(\theta) d\theta$. This formula can be derived as follows. Consider a region $R$ parameterized in polar coordinates by $r(\theta)$ for $a \leq \theta \leq b$. The area of this region would be $\iint_R fdA$. Let $G(r, \theta) = r \langle \cos\theta, \sin\theta\rangle$. Then
That is, for *polar coordinates* $dx dy = r dr d\theta$ ($r \geq 0$).
So by the change of variable formula, we have:
$$
A = \iint_R 1 dx dy = \int_a^b \int_0^{r(\theta)} 1 r dr d\theta = \int_a^b \frac{r^2(\theta)}{2} d\theta.
$$
The key is noting that the region, $S$, described by $\theta$ running from $a$ to $b$ and $r$ running from $0$ to $r(\theta)$, maps onto $R$ through the change of variables. As polar coordinates is just a renaming, this is clear to see.
---
Now consider finding the volume of a sphere using polar coordinates. We have, with $\rho$ being the radius:
$$
V = 2 \iint_R \sqrt{\rho^2 - x^2 - y^2} dy dx,
$$
where $R$ is the disc of radius $\rho$. Using polar coordinates, we have $x^2 + y^2 = r^2$ and the expression becomes:
$$
V = 2 \int_0^{2\pi} \int_0^\rho \sqrt{\rho^2 - r^2} r dr d\theta = 2 \int_0^{2\pi} -(1 - r^2)^{3/2}\frac{1}{3} \mid_0^\rho d\theta = 2\int_0^{2\pi} \frac{\rho^3}{3}d\theta = \frac{4\pi\rho^3}{3}.
$$
##### Linear transformations
Some [transformations](https://en.wikipedia.org/wiki/Transformation_matrix#Examples_in_2D_computer_graphics) from $2$D computer graphics are represented in matrix notation:
or $G(u,v) = \langle au+bv, cu+dv\rangle$. The Jacobian of this *linear* transformation is the matrix itself.
Some common transformations are:
* **Stretching** or $G(u,v) = \langle ku, v \rangle$ or $G(u,v) = \langle u, kv\rangle$ for some $k >0$. The former stretching the $x$ axis, the latter the $y$. These have Jacobian determinant $k$
```{julia}
#| hold: true
#| echo: false
k = 2
G(u,v) = [k*u, v]
showG(G, 1, 1)
```
* **Rotation**. Let $\theta$ be a clockwise rotation parameter, then $G(u,v) = \langle\cos\theta u + \sin\theta v, -\sin\theta u + \cos\theta v\rangle$ will be the transform. The Jacobian is $1$. This figure rotates by $\pi/6$:
* **Shearing**. Let $k > 0$ and $G(u,v) = \langle u + kv, v \rangle$. This transformation is shear parallel to the $x$ axis. (Use $G(u,v) = \langle u, ku+v\rangle$ for the $y$ axis). A shear has Jacobian $1$.
```{julia}
#| hold: true
k = 2
G(u, v) = [u + 2v, v]
showG(G)
```
* **Reflection** If $\vec{l} = \langle l_x, l_y \rangle$ with norm $\|\vec{l}\|$. The reflection through the line in the direction of $\vec{l}$ through the origin is defined, using a matrix, by:
For some simple cases: $\langle l_x, l_y \rangle = \langle 1, 1\rangle$, the diagonal, this is $G(u,v) = (1/2) \langle 2v, 2u \rangle$; $\langle l_x, l_y \rangle = \langle 0, 1\rangle$ (the $y$-axis) this is $G(u,v) = \langle -u, v\rangle$.
* A translation by $\langle a ,b \rangle$ would be given by $G(u,v) = \langle u+a, y+b \rangle$ and would have Jacobian determinant $1$.
As an example, consider the transformation of reflecting through the line $x = 1/2$. Let $\vec{ab} = \langle 1/2, 0\rangle$. This would be found by translating by $-\vec{ab}$ then reflecting through the $y$ axis, then translating by $\vec{ab}$:
```{julia}
#| hold: true
T(u, v, a, b) = [u+a, v+b]
G(u, v) = [-u, v]
@syms u v
a,b = 1//2, 0
x1, y1 = T(u,v, -a, -b)
x2, y2 = G(x1, y1)
x, y = T(x2, y2, a, b)
```
##### Triangle
Consider the problem of integrating $f(x,y)$ over the triangular region bounded by $y=x$, $y=0$, and $x=1$. Such an integral may be computed through Fubini's theorem through $\int_0^1 \int_0^x f(x,y) dy dx$ or $\int_0^1 \int_y^1 f(x,y) dx dy$, but *if* these can not be computed, and a numeric option is preferred, a transformation so that the integral is over a rectangle is preferred.
For this, the transformation $x = u$, $y=uv$ for $(u,v)$ in $[0,1] \times [0,1]$ is possible:
So, $\iint_R f(x,y) dA = \int_0^1\int_0^1 f(u, uv) u du dv$. Here we illustrate with a generic monomial:
```{julia}
@syms x y n::positive m::positive
monomial(x,y) = x^n*y^m
integrate(monomial(x,y), (y, 0, x), (x, 0, 1))
```
And compare with:
```{julia}
#| hold: true
@syms u v
integrate(monomial(u, u*v)*u, (u,0,1), (v,0,1))
```
###### Composition of transformations
What about other triangles, say the triangle bounded by $x=0$, $y=0$ and $y-x=1$?
This can be seen as a reflection through the line $x=1/2$ of the triangle above. If $G_1$ represents the mapping from $U [0,1]\times[0,1]$ into the triangle of the last problem, and $G_2$ represents the reflection through the line $x=1/2$, then the transformation $G_2 \circ G_1$ will map the box $U$ into the desired region. By the chain rule, we have:
Consider the region, $R$, bounded by $y=0$, $x=e^{-n}$, $x=e^n$, and $y=1/x$. An integral over this region may be computed with the help of the transform $G(u,v) = v \langle e^u, e^{-u}\rangle$ which takes the box $[-n, n] \times [0,1]$ onto $R$.
With this, we compute $\iint_R x^2 y^3 dA$ using `SymPy` to compute the Jacobian:
The center of mass is a balancing point of a region with density $\rho(x,y)$. In two dimensions it is a point $\langle \bar{x}, \bar{y}\rangle$. These are found by the following formulas:
$$
A = \iint_R \rho(x,y) dA, \quad \bar{x} = \frac{1}{A} \iint_R x \rho(x,y) dA, \quad
\bar{y} = \frac{1}{A} \iint_R y \rho(x,y) dA.
$$
The $x$ value can be seen in terms of Fubini by integrating in $y$ first:
$$
\iint_R x \rho(x,y) dA = \int_{x=a}^b (\int_{y=h(x)}^{g(x)} \rho(x,y) dy) dx.
$$
The inner integral is the mass of a slice at a value along the $x$ axis. The center of mass is formed then by the mass times the distance from the origin. The center of mass is a "balance" point, in the sense that $\iint_R (x - \bar{x}) dA = 0$ and $\iint_R (y-\bar{y})dA = 0$.
For example, the center of mass of the upper half *unit* disc will have a centroid with $\bar{x} = 0$, by symmetry. We can see this by integrating in *Cartesian* coordinates, as follows
$$
\iint_R x dA = \int_{y=0}^1 \int_{x=-\sqrt{1-y^2}}^{\sqrt{1 - y^2}} x dx dy.
$$
The inner integral is $0$ as it an integral of an *odd* function over an interval symmetric about $0$.
The value of $\bar{y}$ is found using polar coordinate transformation from:
$$
\iint_R y dA = \int_{r=0}^1 \int_{\theta=0}^{\pi} (r\sin(\theta))r d\theta dr =
\int_{r=0}^1 r^2 dr \int_{\theta=0}^{\pi}\sin(\theta) = \frac{1}{3} \cdot 2.
The moment of [inertia](https://en.wikipedia.org/wiki/Moment_of_inertia) of a point mass about an axis is $I = mr^2$ where $m$ is the mass and $r$ the distance to the axis. The moment of inertia of a body is the sum of the moment of inertia of each piece. If $R$ is a region in the $x$-$y$ plane with density $\rho(x,y)$ and the axis is the $y$ axis, then an approximate moment of inertia would be $\sum (x_i)^2\rho(x_i, y_i)\Delta x_i \Delta y_i$ which would lead to $I = \iint_R x^2\rho(x,y) dA$.
Let $R$ be the half disc contained by $x^2 + y^2 = 1$ and $y \geq 0$. Let $\rho(x,y) = xy^2$. Find the moment of inertia.
$$
R
$$
is best described in polar coordinates, so we try to compute
$$
\int_0^1 \int_{-\pi/2}^{\pi/2} (r\cos(\theta))^2 (r\cos(\theta))(r\sin(\theta)) r d\theta dr.
$$
That requires integrating $\sin^2(\theta)\cos^3(\theta)$, a doable task, but best left to SymPy:
(Strang) Find the moment of inertia about the $y$ axis of the unit square tilted *counter*-clockwise an angle $0 \leq \alpha \leq \pi/2$.
The counterclockwise rotation of the unit square is $G(u,v) = \langle \cos(\alpha)u-\sin(\alpha)v, \sin(\alpha)u + \cos(\alpha) v\rangle$. This comes from the above formula for clockwise rotation using $-\alpha$. This transformation has Jacobian determinant $1$, as the area is not deformed. With this, we have
$$
\iint_R x^2 dA = \iint_{G(U)} (f\circ G)(u) |\det(J_G(u))| dU,
Let $R$ be a ring with inner radius $4$ and outer radius $5$. Find its moment of inertia about the $y$ axis.
The integral to compute is:
$$
\iint_R x^2 dA,
$$
with domain that is easy to describe in polar coordinates:
```{julia}
#| hold: true
@syms r theta
x = r*cos(theta)
integrate(x^2 * r, (r, 4, 5), (theta, 0, 2PI))
```
### Three dimensional change of variables
The change of variables formula is no different between dimensions $2$ and $3$ (or higher), but the question of suitable transformation is more involved as the dimensions increase. We stick here to a few widely used ones.
#### Cylindrical coordinates
Polar coordinates describe the $x$-$y$ plane in terms of a radius $r$ and angle $\theta$. *Cylindrical* coordinates describe the $x-y-z$ plane in terms of $r, \theta$, and $z$. A transformation is:
where $r=a/b$. This is *much* simpler in Cylindrical coordinates, as the region is described by the rectangle in $(r, \theta)$: $[0, \sqrt{b/a}] \times [0, 2\pi]$ and the $z$ range is from $0$ to $a - b r$.
The volume then is:
$$
\int_{theta=0}^{2\pi} \int_{r=0}^{a/b} \int_{z=0}^{a - br} 1 r dz dr d\theta =
2\pi \int_{r=0}^{a/b} (a-br)r dr = \frac{\pi a^3}{3b^2}.
$$
This is in agreement with $\pi r^2 h/3$.
---
Find the centroid for the cone. First in the $x$ direction, $\iint_R x dV$ is found by:
That this is $0$ is no surprise. The same will be true for the $y$ direction, as the figure is symmetric about the plane $y=0$ and $x=0$. However, the $z$ direction is different:
```{julia}
#| hold: true
@syms r theta z a b
f(x,y,z) = z
x = r*cos(theta)
y = r*sin(theta)
Jac = r
A = integrate(f(x,y,z) * Jac, (z, 0, a - b*r), (r, 0, a/b), (theta, 0, 2PI))
B = integrate(1 * Jac, (z, 0, a - b*r), (r, 0, a/b), (theta, 0, 2PI))
A, B, A/B
```
The answer depends on the height through $a$, but *not* the size of the base, parameterized by $b$. To finish, the centroid is $\langle 0, 0, a/4\rangle$.
##### Example
A sphere of radius $2$ is intersected by a cylinder of radius $1$ along the $z$ axis. Find the volume of the intersection.
We have $x^2 + y^2 + z^2 = 4$ or $z^2 = 4 - r^2$ in cylindrical coordinates. The integral then is:
Spherical coordinates describe a point in space by a radius from the origin, $r$ or $\rho$; a azimuthal angle $\theta$ in $[0, 2\pi]$ and an *inclination* angle $\phi$ (also called polar angle) in $[0, \pi]$. The $z$ axis is the direction of the zenith and gives a reference line to define the inclination angle. The $x$-$y$ plane is the reference plane, with the $x$ axis giving a reference direction for the azimuth measurement.
The exact formula to relate $(\rho, \theta, \phi)$ to $(x,y,z)$ is given by
$$
G(\rho, \theta, \phi) = \rho \langle
\sin(\phi)\cos(\theta),
\sin(\phi)\sin(\theta),
\cos(\phi)
\rangle.
$$
```{julia}
#| hold: true
#| echo: false
imgfile = "figures/spherical-coordinates.png"
caption = "Figure showing the parameterization by spherical coordinates. (Wikipedia)"
Computing the volume of a sphere is a challenge (for SymPy) in Cartesian coordinates, but a breeze in spherical coordinates. Using $r^2\sin(\phi)$ as the multiplying factor, the volume is simply:
Compute the volume of the ellipsoid, $R$, described by $(x/a)^2 + (y/v)^2 + (z/c)^2 \leq 1$.
We first change variables via $G(u,v,w) = \langle ua, vb, wc \rangle$. This maps the unit sphere, $S$, given by $u^2 + v^2 + w^2 \leq 1$ into the ellipsoid. Then
$$
\iint_R 1 dV = \iint_S 1 |\det(J_G)| dU
$$
But the Jacobian is a constant:
```{julia}
#| hold: true
@syms u v w a b c
G(u,v,w) = [u*a, v*b, w*c]
det(G(u,v,w).jacobian([u,v,w]))
```
So the answer is $abc V(S) = 4\pi abc/3$
## Questions
###### Question
Suppose $f(x,y) = f_1(x)f_2(y)$ and $R = [a_1, b_1] \times [a_2,b_2]$ is a rectangular region. Is this true?
$$
\iint_R f dA = (\int_{a_1}^{b_1} f_1(x) dx) \cdot (\int_{a_2}^{b_2} f_2(y) dy).
$$
```{julia}
#| hold: true
#| echo: false
choices = [
L"Yes. As an inner integral $\int_{a^2}^{b_2} f(x,y) dy = f_1(x) \int_{a_2}^{b_2} f_2(y) dy$.",
"No."
]
answ = 1
radioq(choices, answ)
```
###### Question
Which integrals of the following are $0$ by symmetry? Let $R$ be the unit disc.
$$
a = \iint_R x dA, \quad b = \iint_R (x^2 + y^2) dA, \quad c = \iint_R xy dA
$$
```{julia}
#| hold: true
#| echo: false
choices = [
L"Both $a$ and $b$",
L"Both $a$ and $c$",
L"Both $b$ and $c$"
]
answ = 2
radioq(choices, answ)
```
###### Question
Let $R$ be the unit disc. Which integrals can be found from common geometric formulas (e.g., known formulas for the sphere, cone, pyramid, ellipse, ...)
$$
a = \iint_R (1 - (x^2+y2)) dA, \quad
b = \iint_R (1 - \sqrt{x^2 + y^2}) dA, \quad
c = \iint_R (1 - (x^2 + y^2)^2 dA
$$
```{julia}
#| hold: true
#| echo: false
choices = [
L"Both $a$ and $b$",
L"Both $a$ and $c$",
L"Both $b$ and $c$"
]
answ = 1
radioq(choices, answ)
```
###### Question
Let the region $R$ be described by: in the first quadrant and bounded by $x^3 + y^3 = 1$. What integral below will **not** find the area of $R$?
```{julia}
#| hold: true
#| echo: false
choices = [
raw" ``\int_0^1 \int_0^{(1-x^3)^{1/3}} 1\cdot dy dx``",
raw" ``\int_0^1 \int_0^{(1-y^3)^{1/3}} 1\cdot dy dx``"
]
answ = 3
radioq(choices, answ, keep_order=true)
```
###### Question
Let $R$ be a triangular region with vertices $(0,0), (2,0), (1, b)$ where $b \geq 0$. What integral below computes the area of :R?
```{julia}
#| hold: true
#| echo: false
choices = [
raw" ``\int_0^b\int_{y/b}^{2-y/b} dx dy``",
raw" ``\int_0^2\int_0^{bx} dy dx``",
raw" ``\int_0^2 \int_0^{2b - bx} dy dx``"
]
answ = 1
radioq(choices, answ)
```
###### Question
Let $f(x) \geq 0$ be an integrable function. The area under $f(x)$ over $[a,b]$, $\int_a^b f(x) dx$, is equivalent to?
```{julia}
#| hold: true
#| echo: false
choices = [
raw" ``\int_a^b \int_0^{f(x)} dy dx``",
raw" ``\int_a^b \int_0^{f(x)} dx dy``",
raw" ``\int_0^{f(x)} \int_a^b dx dy``"
]
answ = 1
radioq(choices, answ)
```
###### Question
The region $R$ contained within $|x| + |y| = 1$ is square, but not rectangular (in the sense of integration). What transformation of $S = [-1/2,1/2] \times [-1/2,1/2]$ will have $G(S) = R$?
Let $G(u,v) = \langle \cosh(u)\cos(v), \sinh(u)\sin(v) \rangle$. Using `ForwardDiff` find the determinant of the Jacobian at $[1,2]$.
```{julia}
#| hold: true
#| echo: false
G(u,v) = [cosh(u)*cos(v), sinh(u)*sin(v)]
pt = [1,2]
val = det(ForwardDiff.jacobian(v -> G(v...), [1,2]))
numericq(val)
```
###### Question
Let $G(u, v) = \langle \cosh(u)\cos(v), \sinh(u)\sin(v) \rangle$. Compute the determinant of the Jacobian symbolically:
```{julia}
#| hold: true
#| echo: false
choices = [
raw" ``\sin^{2}{\left (v \right )} \cosh^{2}{\left (u \right )} + \cos^{2}{\left (v \right )} \sinh^{2}{\left (u \right )}``",
raw" ``1``",
raw" ``\sinh(u)\cosh(v)``"
]
answ = 1
radioq(choices, answ)
```
###### Question
Compute the determinant of the Jacobian of the composition of a clockwise rotation by $\theta$, a reflection through the $x$ axis, and then a translation by $\langle a,b\rangle$, using the fact that the Jacobian determinant of *compositions* can be written as product of determinants of the individual Jacobians.
```{julia}
#| hold: true
#| echo: false
choices = [
L"It is $1$, as each is area preserving",
L"It is $r$, as the rotation uses polar coordinates",
L"It is $r^2 \sin(\phi)$, as the rotations use spherical coordinates"
]
answ = 1
radioq(choices, answ)
```
###### Question
A wedge, $R$, is specified by $0 \leq r \leq a$, $0 \leq \theta \leq b$.
```{julia}
#| hold: true
#| echo: false
@syms r theta a b
x = r*cos(theta)
y = r*sin(theta)
A = integrate(r, (r, 0, a), (theta, 0, b))
B = integrate(x * r, (r, 0, a), (theta, 0, b))
C = integrate(y * r, (r, 0, a), (theta, 0, b))
```
What does `A` compute?
```{julia}
#| hold: true
#| echo: false
choices = [
L"The area of $R$",
L"The value $\bar{x}$ of the centroid",
L"The value $\bar{y}$ of the centroid",
L"The moment of inertia of $R$ about the $x$ axis"
]
answ = 1
radioq(choices, answ, keep_order=true)
```
What does $B/A$ compute?
```{julia}
#| hold: true
#| echo: false
choices = [
L"The area of $R$",
L"The value $\bar{x}$ of the centroid",
L"The value $\bar{y}$ of the centroid",
L"The moment of inertia of $R$ about the $x$ axis"
]
answ = 2
radioq(choices, answ, keep_order=true)
```
###### Question
According to [Katz](http://www.jstor.org/stable/2689856) in 1899 Cartan formalized the subject of differential forms (elements such as $dx$ or $du$). Using the rules $dtdt = 0 = dv=dv$ and $dv dt = - dt dv$, what is the product of $dx=mdt + dv\sqrt{1-m^2}$ and $dy=dt\sqrt{1-m^2}-mdv$?
