Solution (source code)

= Solution

Let $B_r(x_0)$ have closure contained in the domain of a <harmonic function> $u$. Put $\omega_n=|B_1(0)|$, so $|B_r|=\omega_nr^n$ and $|\partial B_r|=n\omega_nr^{n-1}$. The sphere and ball <mean value properties for harmonic functions> are
$$
\boxed{u(x_0)=\frac1{|\partial B_r|}\int_{\partial B_r(x_0)}u\,dS=\frac1{|B_r|}\int_{B_r(x_0)}u\,dx.}
$$
To prove the sphere formula, define its average using a fixed unit sphere:
$$
M(r)=\frac1{n\omega_n}\int_{S^{n-1}}u(x_0+r\theta)\,dS_\theta.
$$
Differentiating under the integral and applying the <divergence theorem> gives
$$
M'(r)=\frac1{n\omega_nr^{n-1}}\int_{\partial B_r(x_0)}\partial_\nu u\,dS=\frac1{n\omega_nr^{n-1}}\int_{B_r(x_0)}\Delta u\,dx=0.
$$
Continuity gives $M(r)\to u(x_0)$ as $r\downarrow0$, proving the spherical average identity. Integrating that identity over radii with the polar-coordinate weight yields
$$
\int_{B_r(x_0)}u=\int_0^r n\omega_nt^{n-1}M(t)\,dt=\omega_nr^n u(x_0),
$$
which proves the ball identity. For $n=1$, the sphere average is the average of the two endpoints, and the same formulas hold with the counting surface measure.