Just as the gradient is "the direction of steepest ascent", and the divergence is "amount of stuff created at a point", is there a nice interpretation of the Laplacian Operator (a.k.a. divergence of gradient)?
Assume the function $f$ is $C^2$ in a neighbourhood of ${\bf 0}\in{\mathbb R}^n$. Using the Taylor expansion of $f$ at ${\bf 0}$ one can prove the following:
$$\Delta f({\bf 0}) =\lim_{r\to 0+}\ {2n\over r^2}{1\over \omega(S_r)} \int\nolimits_{S_r}\bigl(f({\bf x})-f({\bf 0})\bigr)\ {\rm d}\omega({\bf x})\ .\qquad (*)$$
This formula says that $\Delta f({\bf 0})$ is "essentially" (i.e., up to the scaling factor ${2n \over r^2}$) equal to the average difference $f({\bf x})-f({\bf 0})$ over small spheres around ${\bf 0}$.

Using this interpretation one gets, e.g., an intuitive understanding of the heat equation
$${\partial u\over\partial t}=a^2\ \Delta u\ ,$$
namely: If averaged over small spheres around a point ${\bf p}$ it is hotter than at ${\bf p}$ itself, then in the next second the temperature at ${\bf p}$ will rise.

Given the interest in the above formula $(*)$, here are some hints for the proof: By Taylor's theorem one has
$$f({\bf x})-f({\bf 0})= \sum_{i=1}^n f_{.i} x_i +{1\over2}\sum_{i,k} f_{.ik} x_i x_k + o(|{\bf x}|^2)\qquad({\bf x}\to{\bf 0}) .$$
Here the $f_{.i}$ and the $f_{.ik}$ are the partial derivatives of $f$ evaluated at ${\bf 0}$, whence constants. Now we integrate this over $S_r$ with respect to the surface measure ${\rm d}\omega$  and  obtain
$$\int\nolimits_{S_r}\bigl(f({\bf x})-f({\bf 0})\bigr)\ {\rm d}\omega({\bf x})={1\over2}\sum_{i}f_{.ii}\int\nolimits_{S_r}x_i^2\ {\rm d}\omega({\bf x}) +o\bigl(r^{2+(n-1)}\bigr) \qquad(r\to0+)\ ,$$
because all other terms are odd in at least one variable. The integrals $\int\nolimits_{S_r}x_i^2\ {\rm d}\omega({\bf x})$ are all equal; therefore we have 
$$\int\nolimits_{S_r}x_i^2\ {\rm d}\omega({\bf x})={1\over n}\int\nolimits_{S_r}\sum_k x_k^2\ {\rm d}\omega({\bf x})={r^2\over n}\omega(S_r)\qquad(1\leq i\leq n)\ .$$
Putting it all together we obtain
$$\int\nolimits_{S_r}\bigl(f({\bf x})-f({\bf 0})\bigr)\ {\rm d}\omega({\bf x})={r^2\over 2n}\omega(S_r)\Delta f({\bf 0}) +o(r^{n+1})\qquad(r\to 0+)\ ,$$
and solving for $\Delta f({\bf 0})$ we get the stated formula.

For a proof using Gauss' theorem see here:

Nice way of thinking about the Laplace operator... but what's the proof?
