Multiple Integral

A multiple integral is a generalization of the usual integral in one dimension to functions of multiple variables in higher-dimensional spaces, e.g.

\[\int \int f(x,y) \,dx \, dy,\]

which is an integral of a function over a two-dimensional region. The most common multiple integrals are double and triple integrals, involving two or three variables, respectively. Since the world has three spatial dimensions, many of the fundamental equations of physics involve multiple integration (e.g. with respect to each spatial variable).

The double integral of the graphed function corresponds to the volume contained underneath the surface corresponding to the function.

Often, multiple integrals are written with a single integral symbol but the notation still implies the correct number of integrals, as in some of the following:

\[\int f\:dA , \quad \int f \: d^3 x , \quad \int f(x_1, \ldots, x_n) \:dx_1\ldots dx_n, \quad \int_{S^2} f, \]

referring to integration over an area \(A\), all of \(\mathbb{R}^3\), all of \(\mathbb{R}^n\), and the surface of the sphere \(S^2\), respectively. Multiple integrals are powerful tools because they allow us to do all the things that can be done in one dimension by integrating, like finding average values or work done, in multiple dimensions.

As in partial differentiation where all but one variable is treated as constant when a partial derivative is taken, in multiple integration all but one variable is treated as constant and integrals are taken iteratively. That is, given an integral

\[\int \int f(x,y) \,dy\,dx,\]

one should treat \(x\) as a constant and perform the integration with respect to \(y,\) and then perform the integration with respect to \(x.\) Note that the integral limits for \(y\) may be functions of \(x\) in this case; these correspond to performing the integration over some region that is not a rectangle in the \(xy\)-plane.

According to Fubini's theorem, in most cases one can interchange the order of integration between \(x\) and \(y\) as desired. This is often very useful, as some integrals can only be evaluated easily in one order. Fubini's theorem fails, however, when \(\int \int \big|f(x,y)\big| dx\, dy = \infty\).

Integration over the same region in 2D can be "sliced" differently by performing \(dx\) or \(dy\) integration first.

Geometrically, a double integral corresponds to the volume under some surface in \(\mathbb{R}^3\). The double integral of the function \(f(x,y) = 1\) therefore corresponds to the area of the region of integration.

Compute the double integral

\[\int_0^1 \int_0^1 \frac{x^2 + y^2}{2} \, dx\,dy.\]

---

First perform the integration with respect to \(x\), and then with respect to \(y\):

\[\int_0^1 \int_0^1 \frac{x^2 + y^2}{2} \: dx\:dy = \left. \int_0^1 \left(\frac{x^3}{6} + \frac{y^2 x}{2} \right) \right|_0^1 dy = \int_0^1 \frac16 + \frac{y^2}{2} \:dy = \left. \frac{y + y^3}{6} \right|_0^1 = \frac13. \]

This integral geometrically corresponds to a sum of the given function over the square with side lengths \(1\). \(_\square\)

Evaluate the triple integral \[\int _{ 0 }^{ 2 }{ \int _{ -3 }^{ 0 }{ \int _{ -1 }^{ 1 }{ \left( { x }^{ 2 }+yz \right) dz\, dy\, dx } } }. \]

---

Since \(z\) is the innermost differential, we integrate with respect to that first. We need to treat all other variables as constants, and integrate with respect to \(z\). This will go as follows: \[\int _{ 0 }^{ 2 }{ \int _{ -3 }^{ 0 }{ { \left[ { x }^{ 2 }z+\frac { 1 }{ 2 } y{ z }^{ 2 } \right] }_{ z=-1 }^{ z=1 }\, dy } \, dx }. \] Just like any other integral, we will plug in the top bound plugged in for \(z\) and then subtract the bottom bound plugged in for \(z\), and get \[\int _{ 0 }^{ 2 }{ \int _{ -3 }^{ 0 }{ \left[ { x }^{ 2 }(1)+\frac { 1 }{ 2 } y({ 1 }^{ 2 })-\left( { { x }^{ 2 }(-1)+\frac { 1 }{ 2 } y{ (-1) }^{ 2 } } \right) \right] dy }\, dx}. \] This simplifies to \[\int _{ 0 }^{ 2 }{ \int _{ -3 }^{ 0 }{ 2{ x }^{ 2 }dy }\, dx }=\int_0^2 6x^2 dx=16.\ _\square \]

Compute the double integral

\[\int_0^1 \int_y^1 e^{-x^2} dx\, dy.\]

---

Since \(e^{-x^2}\) is always positive and finite, and the integration region is a finite area, the integral of the absolute value is finite, so Fubini's theorem applies: the order of integration may be exchanged. This is ideal because the Gaussian integral cannot be computed directly. Note that the integration is over the region:

Region of integration is the triangle given in blue.

Changing the order of integration, the integral becomes \[\int_0^1 \int_0^x e^{-x^2} dy\, dx = \int_0^1 xe^{-x^2} dx = \frac12 \big(1 - e^{-1}\big).\] Note the change in the integral limits coming from the change of variables; after the change of variables, the integration is first performed up from the \(x\)-axis \(y=0\) to the line \(y=x\), and then across from \(x=0\) to \(x=1\). \(_\square\)

Compute the double integrals

\[\int_0^1 \int_0^1 \frac{x^2-y^2}{(x^2+y^2)^2}\, dx\, dy,\quad \int_0^1 \int_0^1 \frac{x^2-y^2}{(x^2+y^2)^2}\, dy\, dx.\]

Comment on the results.

---

Conveniently, the integrand can be written as a mixed second partial derivative:

\[\frac{x^2-y^2}{(x^2+y^2)^2} = -\frac{\partial^2}{\partial x\partial y} \arctan \left(\frac yx\right).\]

A careful analysis of the values of the partial derivatives at the points \((0,1)\) and the origin gives the first double integral

\[ \begin{align} \int_0^1 \int_0^1 \frac{x^2-y^2}{(x^2+y^2)^2} dx\, dy &=\int_0^1 \left(\biggl.-\frac{\partial}{\partial y} \arctan \left(\frac yx\right)\biggr|_{x=1} +\biggl.\frac{\partial}{\partial y} \arctan \left(\frac yx\right)\biggr|_{x=0} \right)dy \\ &= \int_0^1 \left(-\frac{\partial}{\partial y}\arctan (y)+\frac{\partial}{\partial y} \arctan (\infty) \right)dy \\ &= -\arctan(1) \\ &= -\frac{\pi}{4}. \end{align} \]

Performing the other integral yields by symmetry (note the negative sign in the numerator of the integrand)

\[\int_0^1 \int_0^1 \frac{x^2-y^2}{(x^2+y^2)^2}\, dy\, dx = \frac{\pi}{4}.\]

There is a factor of \(-1\) discrepancy between the two different orderings of the integral. This is because the conditions for Fubini's theorem is violated. Near the origin, \(\big(x^2+y^2\big)^2\) diverges rapidly and thus the integral \(\int_0^1 \int_0^1\left| \frac{x^2-y^2}{(x^2+y^2)^2} \right| dx\, dy\) diverges. \(_\square\)

The following problem can be solved with Riemann sums and double integrals. Note that the boundaries of the inner integral can be variables.

For certain functions, it may be advantageous to change the coordinate system used in evaluating integrals. For instance, polar coordinates \(x=r\cos \theta, y=r\sin \theta\) may be used. In this case, however, one must be careful, because the area element \(dA\) is different in Cartesian and polar coordinates. Specifically,

\[\int \int f(x,y)\, dx\, dy = \int \int f(r \cos \theta ,r \sin \theta) \left|\frac{\partial(x,y)}{\partial(r,\theta)}\right| dr\, d\theta = \int \int f(r\cos \theta, r\sin \theta) r\, dr\, d\theta.\]

The expression in the second term is the Jacobian determinant

\[\def\arraystretch{1.5} \left|\frac{\partial(x,y)}{\partial(r,\theta)}\right| = \text{det} \begin{pmatrix} \frac{\partial x}{\partial r} & \frac{\partial x}{\partial \theta} \\ \frac{\partial y}{\partial r} & \frac{\partial y}{\partial \theta} \end{pmatrix} = \text{det} \begin{pmatrix} \cos \theta & -r\sin \theta \\ \sin \theta & r \cos \theta \end{pmatrix} = r\left(\cos^2 \theta + \sin^2 \theta\right) = r.\]

The shaded region, an area element \(dA\), decomposed into its side lengths \(dr\) and \(r\,d\theta\). For infinitesimal \(dr\) and \(d\theta\), \(dA = r\,dr\,d\theta\) is an approximate rectangle.

Integrate the function \(f(r,\theta) = \frac{1}{r}\) over the annulus bounded by the circles \(x^2+y^2 = 4\) and \(x^2+y^2 = 1\).

Annular region described in the problem

---

In polar coordinates, the boundaries of the region are \(r=1\) and \(r=2\), with \(\theta\) ranging from \(0\) to \(2\pi\) over the annulus. The integral is therefore \[\int_0^{2\pi} \int_1^2 \frac{1}{r} r\, dr\, d\theta = 2\pi \int_1^2 dr = 2\pi.\ _\square\]

In a higher number of dimensions, the procedure is no different; the integrals along each variable are done iteratively, and the integral limits may depend on all variables which have not yet been integrated, corresponding to some particular region of integration. In three dimensions, the triple integral \(\int \int \int dx\, dy\, dz\) defined by integrating \(1\) gives the volume of a particular region.

(2012 AIME I, Problem 8)
Cube \(ABCDEFGH\), labeled as shown below, has edge length \(1\) and is cut by a plane passing through vertex \(D\) and the midpoints \(M\) and \(N\) of \(\overline{AB}\) and \(\overline{CG},\) respectively. The plane divides the cube into two solids. The volume of the larger of the two solids can be written in the form \(\tfrac{p}{q}\), where \(p\) and \(q\) are coprime positive integers. Find \(p+q\).

Cube referenced in the problem, from 2012 AIME I

---

Define the point \(M\) to be the origin in Cartesian coordinates. The plane is defined by the two vectors \(\overline{MD}\) and \(\overline{MN}\), which can be written in these coordinates:

\[\overline{MD} = \left(-\frac12 ,1,0\right), \quad \overline{MN} = \left(\frac12, 1, \frac12 \right).\]

The normal vector to the plane is the cross product of the two:

\[\overline{MD} \times \overline{MN} = \left(\frac12 , \frac14, -1\right),\]

so the equation of the plane is

\[\frac12 x + \frac14 y - z = 0.\]

Looking at the smaller solid bounded by this plane and the cube, the \(z\)-coordinate is bounded below by \(0\) and above by the plane \(\frac12 x +\frac14 y\). The bounds on the \(x\)-coordinate are found by intersecting the plane with the \(z=0\) plane: \(x = -\frac12 y\) is the lower bound, while \(x=1\) at the edge of the cube is the upper bound. Within these constraints there are no further constraints on \(y\). Thus the volume of the smaller solid is

\[ \begin{align} \int_0^1 \int_{-\frac12 y}^1 \int_0^{\frac12 x + \frac14 y} \:dz\:dx\:dy &= \int_0^1 \int_{-\frac12 y}^1 \left(\frac12 x + \frac14 y\right) \:dx\:dy \\ &= \int_0^1 \Biggl. \left(\frac{x^2}{4} + \frac14 xy\right)\Bigg|^1_{-\frac12 y} dy\\ &= \int_0^1 \left(\frac14 - \frac14 y - \frac{1}{16} y^2 + \frac{1}{8}y^2 \right) dy \\ &= \frac14 - \frac18 - \frac{1}{48} + \frac{1}{24} = \frac{7}{48}. \end{align} \]

Subtracting from the volume of the unit cube, one finds that the volume of the larger solid is \(\frac{41}{48},\) so \(p+q = 89\). \(_\square\)

As in two dimensions, certain problems may be much easier using a coordinate system adapted to the problem at hand. The natural generalization of polar coordinates in 2D are cylindrical coordinates \((r,\theta, z)\), where \(r\) and \(\theta\) are defined the same way as in 2D and the \(z\)-coordinate is left untouched.

Compute the volume bounded by the cylinders \(x^2+y^2 = 1\) and \(x^2+y^2=4\) and the planes \(z=0\), \(z=1\).

---

The relevant integral can be computed immediately:

\[\int_0^1 \int_0^{2\pi} \int_1^2 r \: dr \: d\theta \:dz = 3\pi.\ _\square\]

Coordinates that are frequently useful in problems with spherical geometry are spherical coordinates \((\rho, \theta, \phi)\), defined by \(x = \rho \sin \phi \cos \theta, y = \rho \sin \phi \sin \theta, z = \rho \cos \phi\). Note that in physics the definitions of \(\theta\) and \(\phi\) are usually interchanged. Geometrically, \(\rho\) is the radius of a point from the origin, \(\theta\) is the usual polar angle, and \(\phi\) is the angle between a vector and the \(z\)-axis.

In spherical coordinates, the Jacobian determinant method used above shows that

\[dx\:dy\:dz = r^2 \sin \phi \: d\rho \: d\theta \: d\phi.\]

Towards the north pole of the sphere, pieces of volume must get smaller by the \(\sin \phi\) factor because at the north pole they shrink to a single point.

Compute the integral of the function \(r \sin \frac{\theta}{4}\) over the upper hemisphere of the unit sphere.

---

Computing from the definition,

\[ \begin{align} \int_0^{\pi/2} \int_0^{2\pi} \int_0^1 r^3 \sin \phi \sin \frac{\theta}{4} \:d\rho\: d\theta\: d\phi &= \int_0^{\pi/2} \int_0^{2\pi} \frac14 \sin \phi \sin \frac{\theta}{4} \: d\theta\: d\phi \\ &= \int_0^{\pi/2} \sin \phi \: d\phi \\ &= 1. \ _\square \end{align}\]

Nothing changes when these methods are generalized to more than three dimensions:

Evaluate the quadruple integral in four dimensions:

\[\int_0^1 \int_0^1 \int_0^1 \int_0^1 (w+x+y+z)\, dw\, dx\, dy\, dz.\]

---

The procedure for evaluating each integral iteratively is the same as previously:

\[ \begin{align} \int_0^1 \int_0^1 \int_0^1 \int_0^1 (w+x+y+z)\, dw\, dx\, dy\, dz &= \int_0^1 \int_0^1 \int_0^1 \left(\frac12 + x+y+z\right) dx\, dy\, dz \\ &= \int_0^1 \int_0^1 \left(\frac12 + \frac12 + y + z\right) dy\, dz \\ &= \int_0^1 \left(\frac12 + \frac12 + \frac12 + z\right) dz \\ &= 2. \ _\square \end{align} \]

The coordinates of the center of mass of a two-dimensional object (or three-dimensional object with translation symmetry in one dimension) are given by

\[\vec{x}_{\text{COM}} = \frac{1}{M} \int \sigma(\vec{x}) \vec{x}\, dA,\]

with \(M\) the total mass of the object and \(\sigma(\vec{x})\) the mass area density.

Find the \((x,y)\) coordinates of the center of mass of the upper half of a disk of radius \(1\) with mass area density \(\sigma = x^2+y^2\).

---

First, compute the total mass of the disk using polar coordinates:

\[M =\int \sigma\, dA = \frac12 \int_0^{2\pi} \int_0^1 r^3 \:dr\:d\theta = (\pi) \frac14 = \frac{\pi}{4}.\]

Now each coordinate can be computed separately. By symmetry, the \(x\)-coordinate is zero. The \(y\)-coordinate is

\[y_{_\text{COM}} = \frac{4}{\pi} \int_0^{\pi} \int_0^1 r^4 \sin \theta \:dr\:d\theta = \frac{4}{5\pi} \int_0^{\pi} \sin \theta \,d\theta = \frac{8}{5\pi}.\ _\square\]

Gauss's law in electricity and magnetism reads

\[\iint \vec{E} \cdot \vec{dA} = \frac{Q}{\epsilon_0},\]

where \(E\) is a vector field in three dimensions, \(A\) is some closed surface, \(Q\) is the total charge enclosed by that surface, and \(\epsilon_0\) is a constant. If \(Q\) is a spherically symmetric charge distribution, derive the electric field outside the charge distribution.

---

On an imaginary sphere at radius \(r\) away from the center of the charge distribution, the electric field must be the same everywhere by spherical symmetry. So the left-hand side turns into

\[\iint \vec{E} \cdot \vec{dA} = E \int dA = E \left(4\pi r^2\right).\]

The surface integral turns into a double integral over the surface once the \(E\) is removed; as discussed previously, this yields the area of the surface. Thus, Gauss's law says

\[E (4\pi r^2) = \frac{Q}{\epsilon_0} \implies E = \frac{Q}{4\pi \epsilon_0 r^2}.\]

This is Coulomb's law for the electric field outside a spherically symmetric ball of charge. \(_\square\)

A particle in quantum mechanics confined to two dimensions has the wavefunction

\[\psi(x,y) = Ce^{-x^2-y^2}.\]

Find the constant \(C\) given the normalization condition \(\int |\psi|^2 dx dy = 1\).

---

Computing the double integral in polar coordinates,

\[\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} C^2 e^{-2x^2-2y^2} dx\, dy = \int_0^{2\pi} \int_0^{\infty} C^2 re^{-2r^2} dr\, d\theta = C^2 2\pi \left(\frac14 \right) = C^2 \frac{\pi}{2}.\]

Therefore, the constant \(C\) is

\[C = \sqrt{\frac{2}{\pi}}.\ _\square\]