Understanding of the Mean Value Theorem in PDE I learned the following theorem in Folland's Introduction to Partial Differential Equations(p.69 Chapter 2):

Suppose $u$ is harmonic on an open set $\Omega\subset{\mathbb R}^n$. If $x\in\Omega$ and $r>0$  is small enough so that $\overline{B_r(x)}\subset\Omega$, then
  $$u(x)=\frac{1}{r^{n-1}\omega_n}\int_{S_r(x)}u(y)d\sigma(y)=\frac{1}{\omega_n}\int_{S_1(0)}u(x+ry)d\sigma(y)，$$
  where $$\omega_n=\frac{2\pi^{n/2}}{\Gamma(n/2)}.$$

I found that I could not immediately reconstruct a proof for the theorem. A key point is that one needs to use the Green's identity, which is a basic property of harmonic functions. But I don't see any "clue" that how people actually come up with this theorem and such proof. (Maybe this is the common problem, at least for me, for most of the textbooks.) A curious search in Google returns nothing satisfactory to me. Since this is a basic property of harmonic functions, I am wondering that if one needs to know this history of harmonic functions in order to know this theorem well. 
Here is my question:


*

*Can any one here come up with a motivation of this theorem in PDE? 


My second question may be more vague:


*

*How can I approach the proof of this theorem more "naturally" instead of just remembering bunch of facts? (In the language of Polya, any heuristics here?)

 A: I can  tell you my favourite proof of the mean value property, which I find more intuitive than the one via Green's identity. To make the proof rigorous, you need to know about integration over the manifold $O(n)$ of all orthogonal matrices, but you can just depict what I do as averaging over all orthogonal matrices. Let's take $x=0$ for simplicity. Heuristically, you'd expect that
$$
  \frac{1}{r^{n-1}\sigma_{n-1}}\int_{S_r(0)}u(y)\,d\sigma(y) = {\rlap{\;\bar{}}\int_{O(n)}} u(Az)\,d\sigma(A) =: f(z)
$$
for $z\in\mathbb R^n$ with $|z|=r$. On the left hand side you take the average of $u(y)$ over all $y$ with $|y|=r$, on the right hand side it's the average of $u(Az)$ over all the orthogonal matrices $A$, which should be and is indeed the same.
Now you can differentiate the right hand side under the integral sign: take the Laplacian with respect to $z$. Then, due to $\Delta u=0$, you see that $f$ is harmonic. Moreover, $f$ is radially symmetric, i.e., $f(z)$ depends on $|z|$ only. Finally, use the fact that a radially symmetric harmonic function defined on all of $\mathbb R^n$ is constant. This yields $f(z) = f(0) = u(0)$, and we're done.
A: Here's a slightly-less-than-rigorous heuristic of an infinitesimal version of the mean value theorem, which provides a sort of motivation for the macroscopic version. (You can probably make this precise and intuitive using non-standard analysis, or you can just make it precise using $\epsilon-\delta$. Here I only give the intuition.) 
Suppose we don't know anything about Harmonicity of a function, and we want to think about mean values. In particular, give a function $u$ and a point $x$, we want to think about how large $u(x)$ is compared with the average of its infinitesimal neighbours. So you consider
$$ \int_{B_{\delta}(x)} u(y) - u(x) dy $$
where $B_\delta$ is the ball of radius $\delta$. 
Of course, if $u$ is continuous, then as $\delta\to 0$ the above expression vanishes. So we need to renormalise by dividing by an appropriate factor of $\delta$. But forget that for the time being. Now, we can assume, by translation, that $x = 0$. And we assume $u$ is sufficiently smooth that we can Taylor expand $u(y)$
$$ u(y) = u(0) + \sum_{i = 1}^d y_i\partial_iu + \frac12 \sum_{i,j = 1}^d y_iy_j\partial^2_{ij} u + \ldots $$
Now, the $u(0) - u(0)$ term cancels out. The first order terms vanishes, because $\partial_iu(0)$ is just some constant, and you are integrating $y_i$, which is an odd function, over a symmetric domain. You also see that by the same token, the integral of $y_iy_j\partial^2_{ij}u(0)$ in $B_\delta$ is zero, if $i\neq j$. So you are left with that the lowest order term 
$$ \int_{B_\delta}u(y) - u(0) dy \sim \int_{B_\delta} \sum_{i = 1}^d y_i^2 \partial_i^2 u(0) dy  = \triangle u(0) \cdot \int_{B_\delta} y_1^2 dy$$
where, by spherical symmetry, the integral over the ball of $y_i^2$ is some fixed constant independent of $i$. So you have that the Laplacian of a function measures the infinitesimal deviation of a function from its mean. 
Once you have the infinitesimal version, the macroscopic version should be something that suggests itself as possibly being true. 
A: The following proof is (almost verbatim) from the freely available Harmonic Function Theory [pdf] by Axler, Bourdon, and Ramey, which is also a nice book I'd recommend for studying harmonic functions.
Let $n>2$. Define $\Omega = \{ \epsilon<\|x\|<1, x\in\mathbb{R}^n \}$ and $v(x)=\|x\|^{2-n}$, and let $u(x)$ be a harmonic function. Denote the unit sphere by $S$. Use Green's second identity on $u$ and $v$ to get
$$ 0 = (2-n) \int_S u\, ds - (2-n)\epsilon^{1-n} \int_{\epsilon S} u\, ds - \int_S \frac{\partial u}{\partial n} ds -\epsilon^{2-n} \int_{\epsilon S} \frac{\partial u}{\partial n} ds.$$
We can use the fact that $\oint \partial u / \partial n ds = 0 $ (again by Green's) to take out the latter two integrals. Then we have, after normalization,
$$ \frac{1}{\omega_n} \int_S u\, ds = \frac{1}{\epsilon^{n-1} \omega_n} \int_{\epsilon S} u\, ds = u(x) \text{ as } \epsilon\to 0. $$
This reasoning generalizes from the unit sphere to any sphere by appropriately scaling or translating $u$. For the planar case with $n=2$, define $v(x) = \ln \| x\|$ and use the same reasoning.
