Integral of a polynomial over a three-dimensional ball 
Let $f$ be a polynomial of total degree at most three in
  $(x,y,z)\in\mathbb{R}^3$. Prove that
  $$\int\limits_{x^2+y^2+z^2\leq1}f(x,y,z)\,dx\,dy\,dz = \frac{4\pi
 f((0,0,0))}{3} + \frac{2\pi(\Delta f)((0,0,0))}{15}$$ Here $\Delta =
 \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} +
 \frac{\partial^2}{\partial z^2}$ is the Laplacian operator on
  $\mathbb{R}^3$.

I even don't know how to approach this problem. Since the LHS of the equality is given by a volume integral, I thought about applying the divergence theorem, but wasn't able to do it. Any help?
 A: This probably isn't the slickest way of presenting the answer, but here goes. We're going to repeatedly use the change of variables theorem in what follows. Let $B$ denote the unit ball in $\Bbb{R}^3$. Then, by the change of variables formula, if $\alpha,\beta,\gamma$ are non-negative integers, then
\begin{equation}
\int_B x^{\alpha}y^{\beta}z^{\gamma} \, dV = - \int_B x^{\alpha}y^{\beta}z^{\gamma} \, dV = 0,
\end{equation}
provided that atleast one of $\alpha,\beta,\gamma$ is odd. For example if $\alpha$ is odd, use the change of variables $(x,y,z) \mapsto (-x,y,z)$. If $\beta$ is odd, use $(x,y,z) \mapsto (x,-y,z)$, and use $(x,y,z) \mapsto (x,y,-z)$ if $\gamma$ is odd. For example, $\int_B x \, dV = \int_B yz^2 \, dV = 0$, etc. This works because the region of integration is still the unit ball $B$, the absolute value of the determinant is $1$, but the actual integrand picks up a minus sign. 
This observation immediately implies that we can ignore "cubic" terms, the "cross-quadratic terms" (like $xy$) and linear terms in $f(x,y,z)$, because these terms all integrate to $0$. So, if we write
\begin{equation}
f(x,y,z) = \text{cubic terms} + a_1x^2 + a_2y^2 + a_3 z^2 + \text{cross quadratic terms} + \text{linear terms} +  c,
\end{equation}
then to quickly integrate such an expression, notice that (by symmetry in change of variables) that
\begin{equation}
\int_B x^2 \, dV = \int_B y^2 \, dV = \int_B z^2 \, dV = \dfrac{1}{3} \int_B (x^2 + y^2 + z^2) \, dV
\end{equation}
The last expression can be easily integrated using spherical coordinates, and the answer is $\dfrac{4\pi}{15}.$
Lastly integrating the constant term means we simply multiply by the volume of the unit ball. Hence, putting this all together, we find that
\begin{align}
\int_B f \, dV = (a_1 + a_2 + a_3) \cdot \dfrac{4 \pi}{15} + c \cdot \dfrac{4\pi}{3} \tag{$*$}
\end{align}
What this shows is that if $f$ is a polynomial of degree at most $3$, then integrating over the unit ball only depends on the constant term, and the coefficients of the "pure quadratic terms" (like $x^2,y^2,z^2$). It is easy to verify that
\begin{align}
a_1+a_2+a_3 = \dfrac{\Delta f(0,0,0)}{2} \quad \text{and} \quad c= f(0,0,0). \tag{$**$}
\end{align}
So, substituting $(**)$ into $(*)$, we find that
\begin{equation}
\int_B f\, dV = \dfrac{4 \pi}{3} \cdot f(0,0,0) + \dfrac{2 \pi}{15} \cdot \Delta f (0,0,0).
\end{equation}

Additional Remarks:
I found this question pretty interesting, so I tried to generalize this result to $n$-dimensions, and here's what I came up with so far. I'll prove that 

If $f: \Bbb{R}^n \to \Bbb{R}$ is a polynomial of degree at most $3$, $B = \{\xi \in \Bbb{R}^n: \lVert \xi\rVert^2 \leq 1\}$ is the closed unit ball, and $dV$ denotes the $n$-dimensional volume element, then
  \begin{align}
\int_B f \, dV &= f(0) \cdot \text{vol}(B) + \dfrac{\text{trace}(D^2f_0)}{2}\cdot \lambda   \\
&= f(0) \cdot \text{vol}(B) + \dfrac{\Delta f(0)}{2}\cdot \lambda,
\end{align}
  where $D^2f_0$ is the second differential of $f$ at $0$ (a symmetric bilinear form), and $\Delta = \sum_{i=1}^n \frac{\partial^2}{\partial x_i^2}$ is the Laplacian, and $\lambda$ is a constant, computed by
  \begin{equation}
\lambda := \dfrac{1}{n} \int_B \lVert \xi \rVert^2 \, dV
\end{equation}

The proof of this is very similar to the one I gave above in the special case. First, note that since $f$ is a polynomial by assumption, it equals its own Taylor polynomial:
\begin{equation}
f(\xi) = f(0) + Df_0(\xi) + \dfrac{1}{2}D^2f_0(\xi)^2 + \dfrac{1}{6}D^3f_0(\xi)^3
\end{equation}
where $D^kf_0$ is a symmetric $k$-linear map from $\Bbb{R}^n \times \cdots \times \Bbb{R}^n$ into $\Bbb{R}$ and $(\xi)^k$ denotes the element $(\xi,\dots, \xi) \in \Bbb{R}^n \times \cdots \times \Bbb{R}^n$ (k-products).
Hence to compute $\int_B f \, dV$, we have to sum up $4$-terms. By an almost identical argument I gave above, the cubic and linear terms all vanish after integration:
\begin{equation}
\int_B Df_0(\xi) \, dV = \int_V D^3f_0(\xi)^3 \, dV = 0.
\end{equation}
So, we have that
\begin{align}
\int_B f \, dV &= \int_B \left [f(0) + \dfrac{1}{2} D^2f_0(\xi)^2  \right] \, dV \\
&= f(0)\cdot \text{vol}(B) + \dfrac{1}{2} \int_B D^2f_0(\xi)^2 \, dV
\end{align}
In the second term, $D^2f_0(\xi)^2$ is a sum of terms of the form $\xi_i \xi_j \cdot (\partial_i \partial_j f)(0)$. But now notice that (again change of variables) if $i \neq j$ then
\begin{equation}
\int_B \xi_i \xi_j \, dV = - \int_B \xi_i \xi_j \, dV = 0
\end{equation}
Hence, the only contribution to the integral comes from terms where $i=j$ "the diagonal terms". More precisely,
\begin{align}
\int_B D^2f_0(\xi,\xi) \, dV &= \sum_{i=1}^n (\partial_i^2 f)(0) \cdot \left( \int_B (\xi_i)^2 \, dV \right) \tag{$\ddot{\smile}$}
\end{align}
But now notice that (again symmetry in change of variables) that
\begin{align}
\int_B (\xi_1)^2 \, dV  = \dots = \int_B (\xi_n)^2 \, dV = \dfrac{1}{n} \int_B \lVert \xi \rVert^2 \, dV =: \lambda \tag{$\ddot{\smile} \ddot{\smile}$}
\end{align}
Substituting $\ddot{\smile} \ddot{\smile}$ into $\ddot{\smile}$ yields the result
\begin{align}
\int_B D^2f_0(\xi)^2 \, dV &= \lambda \cdot \sum_{i=1}^n (\partial_i^2f)(0) \\
&= \lambda \cdot \Delta f(0) = \lambda \cdot \text{trace}(D^2f_0) 
\end{align}
This proves that

\begin{equation}
\int_B f \, dV = f(0) \cdot \text{vol}(B) + \dfrac{\text{trace}(D^2f_0)}{2}\cdot \lambda 
\end{equation}

In the case $n=3$, everything was nice because we could easily compute $\text{vol}(B)$ and $\lambda$ explicitly using spherical coordinates. In higher dimensions, this will necessarily be more complicated, and I think it will involve the use of gamma functions and stuff. Also, if we allow for higher order polynomials, I'm pretty sure we can get a formula involving $f(0),D^2f_0, D^4f_0,D^6f_0 \dots$, although things will probably get more messy. 
