Multiple Integration

Multiple integration extends the concept of the definite integral to functions of several variables defined on regions in ${\mathbb{R}}^n$. This lesson develops the theory with analytical rigor, progressing from Riemann sums to the Change of Variables Theorem.

1. Riemann Integration in ${\mathbb{R}}^n$

Let $B\subset {\mathbb{R}}^n$ be an $n$-dimensional closed box defined by the Cartesian product of intervals: $B=[a_1,b_1]\times [a_2,b_2]\times \cdots \times [a_n,b_n].$ The $n$-dimensional volume (measure) of $B$ is $\text{vol}(B)={\prod }_{i=1}^n(b_i-a_i)$.

Partitions and Riemann Sums

A partition $P$ of $B$ is a collection $\{B_J\}$ of sub-boxes obtained by partitioning each interval $[a_i,b_i]$ into sub-intervals. The mesh size of $P$, denoted $\Vert P\Vert$, is the maximum diameter of any sub-box $B_J$.

For a bounded function $f:B\to \mathbb{R}$, we define the Riemann sum relative to a partition $P$ and a set of sample points $\Xi =\{{\xi }_J\mid {\xi }_J\in B_J\}$ as: $S(f,P,\Xi )={\sum }_{J}f({\xi }_J)\text{vol}(B_J).$

The Definite Integral

The function $f$ is Riemann integrable on $B$ if there exists a number $I\in \mathbb{R}$ such that for every $ϵ>0$, there exists $\delta >0$ such that for any partition $P$ with $\Vert P\Vert <\delta$ and any choice of sample points $\Xi$: $|S(f,P,\Xi )-I|<ϵ$ We denote this value as ${\int }_{B}f(\mathbf{x})dV$ or $\int \cdots {\int }_{B}f(x_1,\dots ,x_n)d{x}_1\dots d{x}_n$.

2. Iterated Integrals and Fubini’s Theorem

Fubini’s Theorem provides the theoretical justification for evaluating multiple integrals via successive univariate integrations.

Theorem (Fubini): Let $A\subset {\mathbb{R}}^n$ and $B\subset {\mathbb{R}}^m$ be closed boxes, and let $f:A\times B\to \mathbb{R}$ be continuous. Then: ${\int }_{A\times B}f(x,y)d(x,y)={\int }_A\left({\int }_{B}f(x,y)dy\right)dx={\int }_B\left({\int }_{A}f(x,y)dx\right)dy.$

Rigorous Caveat: Continuity is a sufficient but not necessary condition. More generally, if $f$ is Lebesgue integrable on $A\times B$, the iterated integrals exist and match the double integral almost everywhere. However, if $f$ is not in $L^1$ (absolutely integrable), the order of integration can lead to different results, notably in cases like $f(x,y)=\frac{x^2-y^2}{(x^2+y^2{)}^2}$, where the integral depends on the path to the origin.

3. The Change of Variables Theorem

The most powerful tool in multi-variable integration is the general Change of Variables Theorem, which generalizes the $u$-substitution of single-variable calculus.

Theorem: Let $U\subset {\mathbb{R}}^n$ be an open set and $\Phi :U\to {\mathbb{R}}^n$ be a $C^1$ diffeomorphism (injective with a continuous inverse). If $f$ is integrable on the region $\Omega =\Phi (U)$, then: ${\int }_{\Omega }f(\mathbf{x})d\mathbf{x}={\int }_{U}f(\Phi (\mathbf{u}))|\det (D\Phi (\mathbf{u}))|d\mathbf{u},$ where $D\Phi (\mathbf{u})$ is the Jacobian matrix of $\Phi$ at $\mathbf{u}$.

The Jacobian Determinant

The term $|\det (D\Phi (\mathbf{u}))|$ represents the local volume expansion factor.

Proof Sketch: At any point ${\mathbf{u}}_0\in U$, the map $\Phi$ can be approximated by its linearization: $\Phi (\mathbf{u})\approx \Phi ({\mathbf{u}}_0)+D\Phi ({\mathbf{u}}_0)(\mathbf{u}-{\mathbf{u}}_0)$ Under a linear transformation represented by matrix $A$, the volume of the image of a set $S$ is $\text{vol}(A(S))=|\det A|\text{vol}(S)$. Thus, a small $n$-box in the $u$-space with volume $d{V}_u$ is mapped to a region in the $x$-space with volume approximately $|\det (D\Phi )|d{V}_u$. Summing these infinitesimal volumes yields the integral identity.

4. Standard Coordinate Systems

Applying the Change of Variables theorem to common geometries leads to specific Jacobian factors:

Polar Coordinates (${\mathbb{R}}^2$)

Mapping: $\Phi (r,\theta )=(r\cos \theta ,r\sin \theta )$ Jacobian: $\det \left(\begin{array}{cc}\cos \theta & -r\sin \theta \\ \sin \theta & r\cos \theta \end{array}\right)=r$ $dA=rdrd\theta .$

Cylindrical Coordinates (${\mathbb{R}}^3$)

Mapping: $x=r\cos \theta ,y=r\sin \theta ,z=z$ Jacobian: $r$ $dV=rdrd\theta dz.$

Spherical Coordinates (${\mathbb{R}}^3$)

Mapping: $x=\rho \sin \varphi \cos \theta ,y=\rho \sin \varphi \sin \theta ,z=\rho \cos \varphi$ Where $\rho$ is radius, $\varphi$ is the polar angle (from $z$-axis), and $\theta$ is the azimuthal angle. Jacobian determinant: ${\rho }^2\sin \varphi$ $dV={\rho }^2\sin \varphi d\rho d\varphi d\theta .$

5. Improper Multiple Integrals

An integral over an unbounded region $D$ is defined via a sequence of bounded regions $D_k$ such that $D_k\subset D_{k+1}$ and $\cup D_k=D$: ${\int }_{D}fdV={\lim }_{k\to \infty }{\int }_{D_k}fdV.$ For the integral to be well-defined (independent of the sequence $D_k$), the function must be absolutely integrable ($|f|\in L^1$). A classic example is the Gaussian integral: ${\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{e}^{-(x^2+y^2)}dxdy={\int }_0^{2\pi }{\int }_0^{\infty }{e}^{-r^2}rdrd\theta =\pi .$

6. Numerical Methods: Monte Carlo Integration

In high dimensions ($n>10$), traditional grid-based methods (like Simpson’s rule) suffer from the “curse of dimensionality,” where the number of points required grows exponentially. Monte Carlo integration provides an alternative.

By the Law of Large Numbers, the integral of $f$ over a region $\Omega$ of volume $V$ is approximated by: ${\int }_{\Omega }fdV\approx \frac{V}{N}{\sum }_{i=1}^{N}f({\mathbf{x}}_i)$ where ${\mathbf{x}}_i$ are points sampled uniformly from $\Omega$. The error scales as $O(1/\sqrt{N})$, regardless of the dimension $n$.

Implementation in Python

The following example uses scipy.integrate.tplquad to compute the mass of a solid with a variable density function $\rho (x,y,z)=x^2+y^2+z^2$ over the unit cube.

import numpy as np
from scipy.integrate import tplquad

# Density function
def density(z, y, x):
    return x**2 + y**2 + z**2

# Bounds for the unit cube [0, 1]x[0, 1]x[0, 1]
# Note: tplquad takes arguments in order (func, q, r, gfun, hfun, qfun, rfun)
# where q, r are outer bounds, g, h are middle, q, r are inner.
mass, error = tplquad(density, 0, 1, lambda x: 0, lambda x: 1, 
                      lambda x, y: 0, lambda x, y: 1)

print(f"Computed Mass: {mass:.6f}")
print(f"Estimated Error: {error:.6e}")

# Theoretical value: Integral of x^2 + y^2 + z^2 from 0 to 1
# = [x^3/3] + [y^3/3] + [z^3/3] = 1/3 + 1/3 + 1/3 = 1.0