The multiple integral
Over a box R = [a_1, b_1] × ... × [a_n, b_n], the multiple integral is built exactly as in one variable. Cut R into sub-boxes by a partition, form a Riemann sum of (function value) × (sub-box volume), and ask whether these sums converge as the mesh shrinks. If they do, f is Riemann integrable and the limit is the integral. A bounded f on a box is integrable exactly when its set of discontinuities has measure zero — the Lebesgue criterion.
Fubini: integrate one variable at a time
Computing a multidimensional limit directly is hopeless. Fubini's theorem rescues us: a double (or n-fold) integral equals an iterated integral, computed one variable at a time, and the order of integration may be swapped. Each inner integral is an ordinary one-variable integral you already know how to do.
Fubini on a box (f continuous, hence integrable):
integral_R f dV = integral_a^b ( integral_c^d f(x, y) dy ) dx
= integral_c^d ( integral_a^b f(x, y) dx ) dy.
Example: f(x, y) = x y over R = [0, 1] x [0, 2].
inner: integral_0^2 x y dy = x * [ y^2/2 ]_0^2 = x * 2 = 2x.
outer: integral_0^1 2x dx = [ x^2 ]_0^1 = 1.
Swap the order to double-check:
inner: integral_0^1 x y dx = y * [ x^2/2 ]_0^1 = y/2.
outer: integral_0^2 (y/2) dy = [ y^2/4 ]_0^2 = 1. Same answer.Change of variables: the Jacobian measures volume
The one-variable substitution rule grows a Jacobian factor in n dimensions. Under a C^1 invertible map T, a tiny box maps to a tiny parallelepiped whose volume is scaled by |det DT|. So change of variables inserts that absolute Jacobian determinant as the local volume-distortion factor.
Change of variables. If T : U -> V is a C^1 bijection with C^1 inverse,
integral_V f(y) dy = integral_U f( T(u) ) * | det DT(u) | du.
Polar coordinates T(r, theta) = ( r cos theta, r sin theta ):
DT = [ cos theta -r sin theta ; sin theta r cos theta ]
det DT = r cos^2 theta + r sin^2 theta = r, so | det DT | = r.
Compute the Gaussian integral via this. Let I = integral_{R^2} e^{-(x^2+y^2)} dA.
I = integral_0^{2pi} integral_0^{infinity} e^{-r^2} * r dr dtheta
= 2pi * [ -1/2 e^{-r^2} ]_0^{infinity}
= 2pi * (1/2) = pi.
Therefore integral_{-inf}^{inf} e^{-x^2} dx = sqrt(I) = sqrt(pi).Look back at the whole ladder: the total derivative gave us a linear approximation, its determinant is exactly the local volume scale, and that is precisely what appears in the integral. Differentiation and integration in several variables meet in the Jacobian.