Theorem: Gauss-Ostrogradsky Theorem

Let $U$ be a subset of $\R^3$ which is compact and has a piecewise smooth boundary $\partial U$.
Let $\mathbf V: \R^3 \to \R^3$ be a smooth vector field defined on a neighborhood of $U$.
Then:
:$\ds \iiint \limits_U \paren {\nabla \cdot \mathbf V} \rd v = \iint \limits_{\partial U} \mathbf V \cdot \mathbf n \rd S$
where $\mathbf n$ is the unit normal to $\partial U$, directed outwards.
That is, the sum of all expansions of $\mathbf V$ within $U$ equals the flux of $\mathbf V$ across the boundary of $U$.

Proof:
It suffices to prove the theorem for rectangular prisms.
The Riemann-sum nature of the triple integral then guarantees the theorem for arbitrary regions.
{{explain|Link to a result demonstrating the above, or explain it better here}}
Let:
:$R = \set {\tuple {x, y, z}: a_1 \le x \le a_2, b_1 \le y \le b_2, c_1 \le z \le c_2}$ 
and let $S = \partial R$, oriented outward. 
Then:
:$S = A_1 \cup A_2 \cup A_3 \cup A_4 \cup A_5 \cup A_6$
where:
:$A_1, A_2$ are those sides perpendicular to the $x$-axis
:$A_3, A_4$ are those sides perpendicular to the $y$ axis
and
:$A_5, A_6$ are those sides perpendicular to the $z$-axis
and in all cases the lower subscript indicates a side closer to the origin.
{{explain|A diagram here is required}}
Let:
:$\mathbf V = M \mathbf i + N \mathbf j + P \mathbf k$
where $M, N, P: \R^3 \to \R$.
Then:
{{begin-eqn}}
{{eqn | l = \iiint_R \nabla \cdot \mathbf V \rd v
      | r = \iiint_R \paren {\frac {\partial M} {\partial x} + \frac {\partial N} {\partial y} + \frac {\partial P} {\partial z} } \rd x \rd y \rd z
}}
{{eqn | r = \iiint_R \frac {\partial M} {\partial x} \rd x \rd y \rd z + \iiint_M \frac {\partial N} {\partial y} \rd x \rd y \rd z + \iiint_M \frac {\partial P} {\partial z} \rd x \rd y \rd z
}}
{{eqn | r = \int_{c_1}^{c_2} \int_{b_1}^{b_2} \paren {\map M {a_2, y, z} - \map M {a_1, y, z} } \rd y \rd z 
}}
{{eqn | o =
      | ro= +
      | r = \int_{c_1}^{c_2} \int_{a_1}^{a_2} \paren {\map N {x, b_2, z} - \map N {x, b_1, z} } \rd x \rd z 
}}
{{eqn | o =
      | ro= +
      | r = \int_{b_1}^{b_2} \int_{a_1}^{a_2} \paren {\map P {x, y, c_2} - \map P {x, y, c_1} } \rd x \rd y
}}
{{eqn | r = \iint_{A_2} M \rd y \rd z - \iint_{A_1} M \rd y \rd z
}}
{{eqn | o =
      | ro= +
      | r = \iint_{A_4} N \rd x \rd z - \iint_{A_3} N \rd x \rd z 
}}
{{eqn | o =
      | ro= +
      | r = \iint_{A_6} P \rd x \rd y - \iint_{A_5} P \rd x \rd y
}}
{{end-eqn}}
We turn now to examine $\mathbf n$:  
{{begin-eqn}}
{{eqn | l = \mathbf n
      | r = \tuple {-1, 0, 0}
      | c = on $A_1$
}}
{{eqn | l = \mathbf n
      | r = \tuple {1, 0, 0}
      | c = on $A_2$
}}
{{eqn | l = \mathbf n
      | r = \tuple {0, -1, 0}
      | c = on $A_3$
}}
{{eqn | l = \mathbf n
      | r = \tuple {0, 1, 0}
      | c = on $A_4$
}}
{{eqn | l = \mathbf n
      | r = \tuple {0, 0, -1}
      | c = on $A_5$
}}
{{eqn | l = \mathbf n
      | r = \tuple {0, 0, 1}
      | c = on $A_6$
}}
{{end-eqn}}
Hence:
{{begin-eqn}}
{{eqn | l = \mathbf V \cdot \mathbf n
      | r = -M
      | c = on $A_1$
}}
{{eqn | l = \mathbf V \cdot \mathbf n
      | r = M
      | c = on $A_2$
}}
{{eqn | l = \mathbf V \cdot \mathbf n
      | r = -N
      | c = on $A_3$
}}
{{eqn | l = \mathbf V \cdot \mathbf n
      | r = N
      | c = on $A_4$
}}
{{eqn | l = \mathbf V \cdot \mathbf n
      | r = -P
      | c = on $A_5$
}}
{{eqn | l = \mathbf V \cdot \mathbf n
      | r = P
      | c = on $A_6$
}}
{{end-eqn}}
We also have:
{{begin-eqn}}
{{eqn | l = \d S
      | r = \d y \rd z
      | c = on $A_1$ and $A_2$
}}
{{eqn | l = \d S
      | r = \d x \rd z
      | c = on $A_3$ and $A_4$
}}
{{eqn | l = \d S
      | r = \d x \rd y
      | c = on $A_5$ and $A_6$
}}
{{end-eqn}}
where $\d S$ is the area element.
This is true because each side is perfectly flat, and constant with respect to one coordinate.
Hence:
{{begin-eqn}}
{{eqn | l = \iint_{A_2} \mathbf V \cdot \mathbf n \rd S
      | r = \iint_{A_2} M \rd y \rd z
}}
{{end-eqn}}
and in general:
{{begin-eqn}}
{{eqn | l = \sum_{i \mathop = 1}^6 \iint_{A_i} \mathbf V \cdot \mathbf n \rd S
      | r = \iint_{A_2} M \rd y \rd z - \iint_{A_1} M \rd y \rd z + \iint_{A_4} N \rd x \rd z - \iint_{A_3} N \rd x \rd z + \iint_{A_6} P \rd x \rd y - \iint_{A_5} P \rd x \rd y
}}
{{end-eqn}}
Hence the result:
:$\ds \iiint_R \nabla \cdot \mathbf V \rd v = \iint_{\partial R} \mathbf V \cdot \mathbf n \rd S$
{{qed}}