```
using CalculusWithJulia
using Plots
plotly()
using QuadGK
using SymPy
using HCubature
```

# 60 Multi-dimensional integrals

This section uses these add-on packages:

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. 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.

## 60.1 Integration theory

The definition of the multi-dimensional integral is more involved than the one-dimensional case due to the possibly increased complexity of the region. This will require additional steps. 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:

\[ \iint_R f(x) dV, \quad \iint_R fdV, \quad \iint_R f(x_1, \dots, x_n) dx_1 \cdot\cdots\cdot dx_n, \quad\iint_R f(\vec{x}) d\vec{x}. \]

A key fact, requiring proof, is:

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\):

\[ \iint_R (af(x) + bg(x))dV = a\iint_R f(x)dV + b\iint_R g(x) dV. \]

**Disjoint:**

- If \(R\) and \(R'\) are
*disjoint*rectangular regions (possibly sharing a boundary), then the integral over the union is defined by linearity:

\[ \iint_{R \cup R'} f(x) dV = \iint_R f(x)dV + \iint_{R'} f(x) dV. \]

**Monotonicity:**

- 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\).

### 60.1.1 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:

```
f(x,y) = x^2 + 5y^2
f(v) = f(v...) # f accepts a vector
= 0, 1
a0, b0 = 0, 2
a1, b1 hcubature(f, (a0, a1), (b0, b1))
```

`(14.0, 1.7763568394002505e-15)`

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:

```
f(x,y) = 3
f(v) = f(v...)
= 0, 4
a0, b0 = 0, 5 # R is area 20, so V = 60 = 3 ⋅ 20
a1, b1 hcubature(f, (a0, a1), (b0, b1))
```

`(60.0, 7.105427357601002e-15)`

- 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\):

```
f(x,y) = x
f(v) = f(v...)
= 0, 1
a0, b0 = 0, 1
a1, b1 hcubature(f, (a0, a1), (b0, b1))
```

`(0.5, 0.0)`

- 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. 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:

```
d(x, y) = sqrt(x^2 + y^2)
function l(x, y, a)
= atan(y,x)
theta = abs(theta)
atheta if (pi/4 <= atheta < 3pi/4) # this is the y=a or y=-a case
/2)/sin(atheta)
(aelse
/2)/abs(cos(atheta))
(aend
end
```

`l (generic function with 1 method)`

And then

```
f(x,y,a,h) = h * (l(x,y,a) - d(x,y))/l(x,y,a)
= 2, 3
𝒂, 𝒉 f(x,y) = f(x, y, 𝒂, 𝒉) # fix a and h
f(v) = f(v...)
```

`f (generic function with 3 methods)`

We can visualize the volume to be computed, as follows:

```
= ys = range(-1, 1, length=20)
xs surface(xs, ys, f)
```

Trying this, we have:

`hcubature(f, (-𝒂/2, -𝒂/2), (𝒂/2, 𝒂/2))`

`(4.000000009419327, 5.9590510310780554e-8)`

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 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:

```
function f(x, y, r)
if x^2 + y^2 < r^2
sqrt(r^2 - (x^2 + y^2))
else
0.0
end
end
```

`f (generic function with 4 methods)`

**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.

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.

The `Quadrature`

package provides a uniform interface for `QuadGK`

, `HCubature`

, and other numeric integration routines available in `Julia`

.

## 60.2 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 48 cells covering the curve and its interior. So the outer measure is less than \(48\cdot 8/100\), as this is just one possible covering.

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.

## 60.3 Fubini’s theorem

Consider again this figure

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. 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\[ \iint_{R\times S}fdV = \iint_R \left(\iint_S f_x(y) dy\right) dx = \iint_S \left(\iint_R f^y(x) dx\right) dy. \]

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, 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 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.

In Volumes 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 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 \cdot (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 \cdot (1-\frac{2x}{a}) dy) dx = 8 \int_0^{a/2} (h \cdot (1-2x/a) \cdot x) dx = 8 h \cdot (\frac{x^2}{2} \big\lvert_{0}^{a/2} - \frac{2}{a}\frac{x^3}{3}\big\lvert_0^{a/2})= 8 h \cdot (\frac{a^2}{8} - \frac{2}{24}a^2) = \frac{a^2h}{3}. \]

### 60.3.1 `SymPy`

’s `integrate`

The `integrate`

function of `SymPy`

uses various algorithms to symbolically integrate definite (and indefinite) integrals. In the section on integrals 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:

`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 right, 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:

`@syms x::real y::real z::real`

`(x, y, z)`

##### Example

For example, the last integral to compute the volume of a square pyramid, could be computed through

```
@syms a height
8 * integrate(height * (1 - 2x/a), (y, 0, x), (x, 0, a/2))
```

\(\frac{a^{2} height}{3}\)

##### Example

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:

`integrate( y * sin(x^2), (x, y^2, 1), (y, 0, 1))`

\(- \frac{3 \sqrt{2} \sqrt{\pi} \left(\frac{3 \sqrt{2} \cos{\left(1 \right)} \Gamma\left(\frac{3}{4}\right)}{16 \sqrt{\pi} \Gamma\left(\frac{7}{4}\right)} + \frac{3 S\left(\frac{\sqrt{2}}{\sqrt{\pi}}\right) \Gamma\left(\frac{3}{4}\right)}{8 \Gamma\left(\frac{7}{4}\right)}\right) \Gamma\left(\frac{3}{4}\right)}{8 \Gamma\left(\frac{7}{4}\right)} + \frac{3 \sqrt{2} \sqrt{\pi} S\left(\frac{\sqrt{2}}{\sqrt{\pi}}\right) \Gamma\left(\frac{3}{4}\right)}{16 \Gamma\left(\frac{7}{4}\right)} + \frac{9 \Gamma^{2}\left(\frac{3}{4}\right)}{64 \Gamma^{2}\left(\frac{7}{4}\right)}\)

##### 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}\).

```
f(x, y) = x^2
h(x) = x^2
g(x) = 5
integrate(f(x,y), (y, h(x), g(x)), (x, -sqrt(Sym(5)), sqrt(Sym(5))))
```

\(\frac{20 \sqrt{5}}{3}\)

##### Example

Find the volume above the \(x\)-\(y\) plane when a cylinder, \(x^2 + y^2 = 2^2\) is intersected by a plane \(3x + 4y + 5z = 6\).

We solve for \(z = (1/5)\cdot(6 - 3x - 4y)\) and take \(R\) as the disk at the origin of radius \(2\):

```
f(x,y) = 6 - 3x - 4y
g(x) = sqrt(2^2 - x^2)
h(x) = -sqrt(2^2 - x^2)
1//5) * integrate(f(x,y), (y, h(x), g(x)), (x, -2, 2)) (
```

\(\frac{24 \pi}{5}\)

##### Example

Find the volume:

- in the first octant
- 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:

```
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

```
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))
```

\(\frac{500}{27}\)

##### 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?

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 wedge between \(y=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.

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:

```
= 1 # if using r as a symbolic variable specify `positive=true`
r 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)))
```

\(16 - 8 \sqrt{2}\)

##### 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:

```
f(x,y) = x*y
g(y) = sqrt(1 - y^2)
h(y) = y
integrate(f(x,y), (x, h(y), g(y)), (y, 0, sin(PI/4)))
```

\(\frac{1}{16}\)

##### Example: Average value

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.

```
f(x,y) = x^2 + y^2
g(x) = x # solve x - y = 0 for y
h(x) = 0
= integrate(f(x,y), (y, h(x), g(x)), (x, 0, 1))
A = integrate(Sym(1), (y, h(x), g(x)), (x, 0, 1))
B /B A
```

\(\frac{2}{3}\)

(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 = 3x^2\).

So we get a mass of:

```
rho(x,y) = x^2*y^2
g(x) = 2 - x^2
h(x) = -3 + 2x^2
= sqrt(Sym(5//3))
a integrate(rho(x,y), (y, h(x), g(x)), (x, -a, a))
```

\(\frac{460 \sqrt{15}}{729}\)

##### 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:

`integrate(cos(x^2), (y, 0, x), (x, 0, 1))`

\(\frac{\sin{\left(1 \right)}}{2}\)

### 60.3.2 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.

For these notes, we define three operations using Unicode operators entered with `\int[tab]`

, `\iint[tab]`

, `\iiint[tab]`

.

```
# 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
# integrate int_a^b int_h(x)^g(y) f(x,y) dy dx
∬(f, ys, xs) = ∫(x -> ∫(y -> f(x,y), endpoints(ys, x)), xs)
# integrate f(x,y,z) dz dy dx
∭(f, zs, ys, xs) = ∫(
-> ∫(
x -> ∫(
y -> f(x,y,z),
z endpoints(zs, (x,y))),
endpoints(ys,x)),
xs)
```

`∭ (generic function with 1 method)`

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.

##### Example

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.

```
f(x,y) = exp(-x^2 - 2y^2)
f(v) = f(v...)
hcubature(f, (0,0), (3,3)) # (a0, a1), (b0, b1)
```

`(0.5553480979840428, 8.155874399598429e-9)`

```
f(x,y) = exp(-x^2 - 2y^2)
∬(f, (0,3), (0,3)) # (a1, b1), (a0, b0)
```

`0.5553480979874703`

##### Example

Show the area of the unit circle is \(\pi\) using the “Fubini” function.

```
f(x,y) = 1
= ∬(f, (x-> -sqrt(1-x^2), x-> sqrt(1-x^2)), (-1, 1))
a - pi # answer and error a, a
```

`(3.141592655913248, 2.3234547619210844e-9)`

(The error is similar to that returned by `quadgk(x -> 2sqrt(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.

```
f(x,y) = sqrt(1 - x^2 - y^2)
= 2 * ∬(f, (x-> -sqrt(1-x^2), x-> sqrt(1-x^2)), (-1, 1))
a - 4/3*pi a, a
```

`(4.188790207884331, 3.0979405707398655e-9)`

##### Example

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:

`∬((y,x) -> cos(x^2), (y -> y, 1), (0, 1))`

`0.42073549240394825`

Compare to

`sin(1)/2`

`0.42073549240394825`

## 60.4 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.

Here we compute the volumes of these using a triple integral of the form \(\iiint_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:

```
@syms a b c
f(x,y,z) = Sym(1) # need to integrate a symbolic object in integrand or call `sympy.integrate`
integrate(f(x,y,z), (x, 0, a), (y, 0, b), (z, 0, c))
```

\(a b c\)

- Prism. Consider a prism or wedge formed by \(ay + bz = 1\) with \(a,b > 0\) and over the region in the first quadrant \(0 \leq x \leq c\). Find its area.

The function to integrate is \(f(x,y) = (1 - ay)/b\) over the region bounded by \([0,c] \times [0,1/a]\):

```
@syms a b c
f(x,y,z) = Sym(1)
integrate(f(x,y,z), (z, 0, (1 - a*y)/b), (y, 0, 1/a), (x, 0, c))
```

\(\frac{c}{2 a b}\)

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

```
@syms a b c
f(x,y,z) = Sym(1)
integrate(f(x,y,z), (z, 0, (1 - a*x - b*y)/c), (y, 0, (1 - a*x)/b), (x, 0, 1/a))
```

\(\frac{1}{6 a b c}\)

This is \(1/6\)th the volume of the box.

- 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:

\[ \int_{x=-r}^r \int_{y=-\sqrt{r^2 - x^2}}^{\sqrt{r^2-x^2}} \int_0^{a - b(x^2 + y^2)^{1/2}} 1 dz dy dx. \]

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:

\[ \int_{x=-1}^1 \int_{y=-\sqrt{1-x^2}}^{\sqrt{1-x^2}} \int_{z=-\sqrt{1 - x^2 - y^2}}^{\sqrt{1-x^2 - y^2}} 1 dz dy dx. \]

This integral is doable, but `SymPy`

has trouble with it. We will return to this when spherical coordinates are defined.

## 60.5 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 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

\[\begin{align*} dx &= A du + B dv\\ dy &= C du + D dv, \end{align*}\]

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 = -(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 - C(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:

\[\begin{align*} x &= a + mt + \sqrt{1-m^2} v\\ y & = b + \sqrt{1-m^2}t -mv \end{align*}\]

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\).

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\).

```
G(u, v) = u * [cos(v), sin(v)]
G(v) = G(v...)
J(v) = ForwardDiff.jacobian(G, v) # [∇g1', ∇g2']
= 6
n = range(0, 1, length=3n) # radius
us = range(0, 2pi, length=3n) # angle
vs
plot(unzip(G.(us', vs))..., legend = false, aspect_ratio=:equal) # plots constant u lines
plot!(unzip(G.(us, vs'))...) # plots constant v lines
= [us[n],vs[n]]
pt
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.

`showG (generic function with 3 methods)`

The transformation to elliptical coordinates, \(G(u,v) = \langle \cosh(u)\cos(v), \sinh(u)\sin(v)\rangle\), may be viewed similarly:

The transformation \(G(u,v) = v \langle e^u, e^{-u} \rangle\) uses hyperbolic coordinates:

The transformation \(G(u,v) = \langle u^2-v^2, u\cdot v \rangle\) yields a partition of the plane:

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

\[ \| \det\left( \begin{bmatrix} \hat{i} & \hat{j} & \hat{k}\\ u_1 & u_2 & 0\\ v_1 & v_2 & 0 \end{bmatrix} \right) \| = \| \hat{k} \det\left( \begin{bmatrix} u_1 & u_2\\ v_1 & v_2 \end{bmatrix} \right) \| = | \det\left( \begin{bmatrix} u_1 & u_2\\ v_1 & v_2 \end{bmatrix} \right)|. \]

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 Let \(U\) be an open set in \(R^n\), \(G:U \rightarrow R^n\) be an

injectivedifferentiable function withcontinuouspartial 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.

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)\).

### 60.5.1 Two dimensional change of variables

Now we see several examples of two-dimensional transformations.

#### Polar integrals

We have seen 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 1dA\). Let \(G(r, \theta) = r \langle \cos\theta, \sin\theta\rangle\). Then

\[ J_G = \begin{bmatrix} \cos(\theta) & - r\sin(\theta)\\ \sin(\theta) & r\cos(\theta) \end{bmatrix}, \]

with determinant \(r\).

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} -(\rho^2 - 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 from \(2\)D computer graphics are represented in matrix notation:

\[ \begin{bmatrix} x\\ y \end{bmatrix} = \begin{bmatrix} a & b\\ c & d \end{bmatrix} \begin{bmatrix} u\\ v \end{bmatrix} , \]

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\)

**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\).

```
= 2
k G(u, v) = [u + 2v, v]
showG(G)
```