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?
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
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
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.