Multiple Integrals — Double & Triple

Fubini's theorem says the two orders agree whenever $f$ is continuous on a closed bounded region, or more generally when $\iint_R \abs{f}\,\dd{A} < \infty$ (Tonelli–Fubini for measurable functions). Geometrically, you're computing the same volume by slicing it two different ways — vertical strips vs horizontal strips. They must agree because volume doesn't depend on how you slice.

It can fail when $f$ is unbounded and $\iint\abs{f}=\infty$. Classic counterexample: $f(x,y) = (x^2-y^2)/(x^2+y^2)^2$ on $[0,1]^2$. The two iterated integrals give $\pi/4$ and $-\pi/4$ — different! The double integral doesn't exist in the Riemann/Lebesgue sense because $\iint\abs{f}=\infty$.

In polar coordinates, a "cell" with $\Delta r$ and $\Delta\theta$ is an annular sector. Its area is approximately $r\,\Delta r\,\Delta\theta$ — not $\Delta r\,\Delta\theta$. The extra factor of $r$ is exactly $\abs{\det J_\Phi}$ where $\Phi(r,\theta)=(r\cos\theta,r\sin\theta)$. Sectors at large $r$ are much bigger than sectors near the origin, even though both have the same $\Delta r,\Delta\theta$. The simulation paints each sector with its actual Cartesian area, which is why outer rings appear wider.

Operationally: when you switch to polar mode, the readout shows $S_n = \sum f(r_i,\theta_j)\cdot r_i\,\Delta r\,\Delta\theta$. Drop the $r$ and you get a wrong answer that systematically under-weights large-$r$ contributions.

Physics: mass $m=\iiint_E \rho\dd{V}$ and center of mass $\bar x = \frac{1}{m}\iiint_E x\rho\dd{V}$; moments of inertia $I_z = \iiint_E (x^2+y^2)\rho\dd{V}$; electric/gravitational potentials are 3D integrals of source densities; flux through a surface is a double integral; total charge from charge density. Probability: for jointly continuous $(X,Y)$ with density $f$, $P((X,Y)\in B) = \iint_B f\dd{A}$ and $E[g(X,Y)] = \iint g f\dd{A}$. The normalizing constant of a 2D Gaussian is exactly the polar trick above. Engineering: finite-element methods reduce PDE problems to assembling many local 2D/3D integrals over mesh elements; computing volumes/areas of CAD parts; structural moments. Statistics / ML: Bayesian posteriors are integrals over high-dimensional parameter spaces — Monte Carlo and MCMC exist because deterministic quadrature collapses past $\sim 6$ dimensions.

For a smooth map $\Phi:S\to R$, near any point $\Phi$ is well-approximated by its derivative — a linear map represented by the Jacobian matrix $J_\Phi$. A linear map sends a unit square to a parallelogram whose area equals $\abs{\det J}$. So $\abs{\det J_\Phi(u,v)}$ is the local area amplification factor at $(u,v)$. The change-of-variables formula is "weight each tiny patch in $S$ by how much it stretches when mapped to $R$, then integrate".

In 3D the same idea works with volume amplification; in $n$ dimensions, $\abs{\det J}$ is an $n$-volume scaling. If $\det J = 0$ at a point, the map collapses dimensions locally (singularity — handle with care). If $\det J < 0$, $\Phi$ flips orientation; we take absolute value because area is unsigned in this context (signed forms appear in differential geometry).

Tensor-product quadrature (Simpson, Gauss, midpoint on a grid) in $d$ dimensions needs $n^d$ function evaluations for $n$ points per axis. Error decays as $n^{-p}$ where $p$ is the rule's order — so error vs cost: $\text{err} \sim N^{-p/d}$ for $N = n^d$ evaluations. For $d=6$, even Simpson ($p=4$) gives $N^{-2/3}$ — slow.

Monte Carlo: pick $N$ uniform random points in $R$, average $f$ values, multiply by $\text{vol}(R)$. By CLT, error decays as $1/\sqrt N$ regardless of $d$. So once $d \gtrsim 2p$, Monte Carlo wins. Bayesian inference routinely operates in $d = 100$–$10^6$ — quadrature is hopeless, Monte Carlo (especially MCMC/HMC) is the only option.

This is the famous "curse of dimensionality" for deterministic integration. Quasi-Monte Carlo (low-discrepancy sequences like Sobol/Halton) bridges the gap with $\sim N^{-1}(\log N)^d$ — better than MC for moderate $d \approx 10$–$50$ but still defeated eventually.

The Riemann integral in $\R^n$ has serious limitations: it can't handle wildly discontinuous integrands and behaves badly under limits ($\lim\int \neq \int\lim$ for many useful sequences). Lebesgue's reformulation: replace "partition the domain" with "partition the range". For a measurable $f\ge 0$, $\int_R f\dd\mu = \sup\{\int_R s\dd\mu : s \text{ simple, } 0\le s\le f\}$ where simple functions take finitely many values.

Fubini's theorem in measure-theoretic form (for product measures $\mu\times\nu$): $\int_{X\times Y} f\,d(\mu\times\nu) = \int_X\int_Y f(x,y)\,d\nu\,d\mu$ when $f$ is non-negative measurable (Tonelli) or when $\int\abs{f}<\infty$ (Fubini). This is the right framework for probability theory: an integral over $\R^n$ is the same kind of object as an expectation $\int_\Omega X\,dP$, and product measures correspond to independent random variables.

The change of variables formula generalizes: for $\Phi:U\to V$ a $C^1$-diffeomorphism between open sets in $\R^n$, $\int_V f\,d\lambda = \int_U (f\circ\Phi)\abs{\det D\Phi}\,d\lambda$ where $\lambda$ is Lebesgue measure. Same formula, just rigorously proven for measurable $f$ on Lebesgue-measurable sets.

Cylindrical $(r,\theta,z)$ with $\dd{V}=r\dd{r}\dd\theta\dd{z}$: choose when the region has rotational symmetry about an axis (typically the $z$-axis) and the bounds in $z$ are simple — e.g.\ paraboloid $z=r^2$ to plane $z=1$, or cylinder cores.

Spherical $(\rho,\phi,\theta)$ with $\dd{V}=\rho^2\sin\phi\dd\rho\dd\phi\dd\theta$ (here $\rho$ = distance from origin, $\phi$ = polar angle from $+z$ axis, $\theta$ = azimuth): choose when the region has radial symmetry from a point — balls, spherical shells, cones from the origin, gravitational/electric problems of point sources.

Quick test: does the integrand or region simplify when written in those coordinates? $x^2+y^2 \to r^2$ (cylindrical) or $\rho^2\sin^2\phi$ (spherical); $x^2+y^2+z^2 \to \rho^2$ (spherical, dramatically simpler). If yes, switch. Don't forget the Jacobian — that $r$ or $\rho^2\sin\phi$ is mandatory.