Theorem: let V be a compact and simply-connected region in \mathbb{R}^3 whose boundary, denoted \partial V, is a piecewise smooth, closed, and oriented surface. The orientation is taken to be outward, meaning the normal vector at any point on the surface points away from the interior of V.
Consider a vector field \mathbf{F}: U \to \mathbb{R}^3, where U is an open subset of \mathbb{R}^3 that contains the region V. If \mathbf{F} is continuously differentiable in U (i.e., its components have continuous first-order partial derivatives), the divergence theorem states that the volume integral of the divergence of \mathbf{F} over V is equal to the surface integral of the flux of \mathbf{F} across the boundary \partial V:
\iiint_V (\nabla \cdot \mathbf{F}) \, \mathrm{d}V = \oiint_{\partial V} \mathbf{F} \cdot \mathrm{d}^2\boldsymbol{\Sigma}
In this expression, \nabla \cdot \mathbf{F} represents the divergence of the vector field \mathbf{F}, and \mathrm{d}^2\boldsymbol{\Sigma} = \mathbf{n} \, \mathrm{d}S is the oriented surface element, with \mathbf{n} being the outward-pointing unit normal to the surface \partial V.
Proof: we proceed by establishing the theorem for a specific class of “simple” volumes and then extending the result to more general regions through decomposition. The strategy involves decomposing the vector field \mathbf{F} into its components and proving the theorem for each component separately.
Let the vector field be \mathbf{F}(x,y,z) = F_x(x,y,z)\mathbf{i} + F_y(x,y,z)\mathbf{j} + F_z(x,y,z)\mathbf{k}. The Divergence Theorem can be written as the sum of three independent equalities:
\begin{aligned} \iiint_V \frac{\partial F_x}{\partial x} \, \mathrm{d}V &= \oiint_{\partial V} (F_x\mathbf{i}) \cdot \mathbf{n} \, \mathrm{d}S \\ \iiint_V \frac{\partial F_y}{\partial y} \, \mathrm{d}V &= \oiint_{\partial V} (F_y\mathbf{j}) \cdot \mathbf{n} \, \mathrm{d}S \\ \iiint_V \frac{\partial F_z}{\partial z} \, \mathrm{d}V &= \oiint_{\partial V} (F_z\mathbf{k}) \cdot \mathbf{n} \, \mathrm{d}S \end{aligned}
We will provide a detailed proof for the third equality. The proofs for the other two are entirely analogous.
Let us consider a “z-simple” region V, which is a volume defined by:
V = \{ (x,y,z) \in \mathbb{R}^3 \mid (x,y) \in D, \ g_1(x,y) \le z \le g_2(x,y) \}
Here, D is the projection of V onto the xy-plane, and z=g_1(x,y) and z=g_2(x,y) are the functions describing the lower and upper bounding surfaces of V, respectively.
First, we evaluate the volume integral using iterated integration and the fundamental theorem of calculus:
\begin{aligned} \iiint_V \frac{\partial F_z}{\partial z} \, \mathrm{d}V &= \iint_D \left( \int_{g_1(x,y)}^{g_2(x,y)} \frac{\partial F_z}{\partial z} \, \mathrm{d}z \right) \mathrm{d}A \\ &= \iint_D F_z(x,y,z) \bigg|_{z=g_1(x,y)}^{z=g_2(x,y)} \mathrm{d}A \\ &= \iint_D \left[ F_z(x,y, g_2(x,y)) - F_z(x,y, g_1(x,y)) \right] \mathrm{d}A \end{aligned}
Next, we compute the flux integral \oiint_{\partial V} (F_z\mathbf{k}) \cdot \mathbf{n} \, \mathrm{d}S.
The boundary surface \partial V consists of three parts: the top surface S_2 (where z=g_2(x,y)), the bottom surface S_1 (where z=g_1(x,y)), and the side surface S_s.
For the top surface S_2, parametrized by \mathbf{r}(x,y) = (x, y, g_2(x,y)), the outward (upward) normal vector is given by \mathbf{n}_2 \mathrm{d}S = \left(-\frac{\partial g_2}{\partial x}\mathbf{i} - \frac{\partial g_2}{\partial y}\mathbf{j} + \mathbf{k}\right) \mathrm{d}A. The flux of F_z\mathbf{k} through S_2 is:
\iint_{S_2} (F_z\mathbf{k}) \cdot \mathbf{n}_2 \, \mathrm{d}S = \iint_D F_z(x,y, g_2(x,y)) \, \mathbf{k} \cdot \left(-\frac{\partial g_2}{\partial x}\mathbf{i} - \frac{\partial g_2}{\partial y}\mathbf{j} + \mathbf{k}\right) \mathrm{d}A = \iint_D F_z(x,y, g_2(x,y)) \, \mathrm{d}A
For the bottom surface S_1, parametrized by \mathbf{r}(x,y) = (x, y, g_1(x,y)), the standard parametrization yields an upward normal. However, the outward normal for the bottom surface must point downward. We then reverse its direction:
\mathbf{n}_1 \mathrm{d}S = - \left(-\frac{\partial g_1}{\partial x}\mathbf{i} - \frac{\partial g_1}{\partial y}\mathbf{j} + \mathbf{k}\right) \mathrm{d}A
The flux through S_1 is then:
\iint_{S_1} (F_z\mathbf{k}) \cdot \mathbf{n}_1 \, \mathrm{d}S = \iint_D F_z(x,y, g_1(x,y)) \, \mathbf{k} \cdot \left(\frac{\partial g_1}{\partial x}\mathbf{i} + \frac{\partial g_1}{\partial y}\mathbf{j} - \mathbf{k}\right) \mathrm{d}A = -\iint_D F_z(x,y, g_1(x,y)) \, \mathrm{d}A
For the side surface S_s, the outward normal vector \mathbf{n}_s is everywhere horizontal.
Consequently, its dot product with the vertical vector \mathbf{k} is zero, i.e., \mathbf{n}_s \cdot \mathbf{k} = 0. The flux through the side surface is therefore zero.
Summing the flux contributions from all parts of the boundary surface gives:
\oiint_{\partial V} (F_z\mathbf{k}) \cdot \mathbf{n} \, \mathrm{d}S = \iint_D F_z(x,y, g_2(x,y)) \, \mathrm{d}A - \iint_D F_z(x,y, g_1(x,y)) \, \mathrm{d}A
Comparing this result with the result from the volume integral, we see they are identical. This proves the third equality for z-simple regions. By symmetric arguments for x-simple and y-simple regions, the other two equalities are also established.
A general region can be subdivided into a finite number of simple regions. The sum of the volume integrals over the subregions is the volume integral over the entire region. For the surface integrals, the flux across any shared interior boundary between two subregions cancels out, as the outward normal vectors are in opposite directions. This leaves only the flux through the exterior boundary of the original region, and therefore, the theorem holds for any such general region.