Integration of three dimensional function gives wrong answer!

I have the following integration

$$\int_{-\infty}^{\infty}d^{3}\mathbf{p}\nabla\cdot\frac{\mathbf{p}}{(p_{x}^{2}+p_{y}^{2}+p_{z}^{2}+m^{2})^{3/2}}$$ $$=\int_{-\infty}^{\infty}d^{3}\mathbf{p}\left(\partial_{x}\frac{p_{x}}{(p_{x}^{2}+p_{y}^{2}+p_{z}^{2}+m^{2})^{3/2}}\right.\\\partial_{y}\frac{p_{x}}{(p_{x}^{2}+p_{y}^{2}+p_{z}^{2}+m^{2})^{3/2}}\\\left.\partial_{z}\frac{p_{x}}{(p_{x}^{2}+p_{y}^{2}+p_{z}^{2}+m^{2})^{3/2}}\right)$$

When I put the first part into Mathematica, I found:

 Assuming[m^2 > 0,
Integrate[D[px (px^2 + py^2 + pz^2 + m^2)^(-3/2), px],
{px, -∞, +∞}, {py, -∞, +∞}, {pz, -∞, +∞}]]


result: $$4\pi$$.

However, if I put all the three terms:

Assuming[m^2 > 0,
Integrate[
D[px (px^2 + py^2 + pz^2 + m^2)^(-3/2), px] +
D[py (px^2 + py^2 + pz^2 + m^2)^(-3/2), py] +
D[pz (px^2 + py^2 + pz^2 + m^2)^(-3/2), pz],
{px, -∞, +∞}, {py, -∞, +∞}, {pz, -∞, +∞}]]


the result is also $$4\pi$$. This is very strange. Should it be $$12\pi$$?

Update

In other words, why should the integral depend on the order of the integration variables?

$$\int dp_{x}dp_{y}dp_{z}\left[\frac{1}{(p_{x}^{2}+p_{y}^{2}+p_{z}^{2}+m^{2})^{3/2}}-\frac{3p_{x}^{2}}{(p_{x}^{2}+p_{y}^{2}+p_{z}^{2}+m^{2})^{5/2}}\right]=4\pi$$ $$\neq\int dp_{y}dp_{z}dp_{x}\left[\frac{1}{(p_{x}^{2}+p_{y}^{2}+p_{z}^{2}+m^{2})^{3/2}}-\frac{3p_{x}^{2}}{(p_{x}^{2}+p_{y}^{2}+p_{z}^{2}+m^{2})^{5/2}}\right]=0$$

#

By the way, if I use Nintegrate the order does not affect the result, very strange!!

m=1.;
NIntegrate[
D[px (px^2 + py^2 + pz^2 + m^2)^(-3/2),
px], {px, -\[Infinity], +\[Infinity]}, {py, -\[Infinity], +\
\[Infinity]}, {pz, -\[Infinity], +\[Infinity]}]


result: 12.5565.

NIntegrate[
D[px (px^2 + py^2 + pz^2 + m^2)^(-3/2),
px], {py, -\[Infinity], +\[Infinity]}, {pz, -\[Infinity], +\
\[Infinity]}, {px, -\[Infinity], +\[Infinity]}]


result: 12.5565.

Very strange!!! And this is very dangerous for the numerical calculation.

• Did you set a value for m in NIntegrate? I get the same result setting m=1. The numerical value is 4 Pi. Commented Nov 16, 2019 at 13:34
• @mikado Yes, m=1. And the order of integration variables now does not affect the result. Commented Nov 16, 2019 at 13:35
• I think Fubini–Tonelli does not apply. It’s similar to the last counterexample here en.m.wikipedia.org/wiki/Fubini%27s_theorem Commented Nov 16, 2019 at 16:50
• "Very strange!!! And this is very dangerous for the numerical calculation." -- 1. NIntegrate does not use iterative integration. 2. It is assumed that the Fubini-Tonelli theorem applies to most of the integrals given to NIntegrate. 3. For the integrals in the question you are most likely going to get different results using Cartesian rules. Commented Nov 17, 2019 at 20:00

The integral over the subregion does not converge:

Integrate[
(m^2 - 2 x^2 + y^2 + z^2)/(m^2 + x^2 + y^2 + z^2)^(5/2),
{y, -Infinity, Infinity},
{z, -Infinity, Infinity},
{x, -Sqrt[1 + y^2 + z^2], Sqrt[1 + y^2 + z^2]},
Assumptions -> m > 0 && {x, y, z} \[Element] Reals]
(*  Infinity  *)


The triple integral does not equal the iterated integral, something that Integrate[] misses.

The surface m^2 - 2 x^2 + y^2 + z^2 == 0 divides space into a region over which the integral diverges to positive infinity and one over which the integral diverges to negative infinity. One could try to choose a principal value. One has to be aware that one can obtain any result. The surface m^2 - 2 x^2 + y^2 + z^2 == 0 was a convenient (and somewhat obvious) choice for analyzing the divergence of the integral. It is not necessarily for it to be used to define a principal value. A common choice is as follows. It has the appealing attraction of corresponding somewhat with the symmetry of the integral. Since over a ball $$B$$ centered at the origin we have by symmetry $$\textstyle \int_B \frac{x^2}{\left(m^2+x^2+y^2+z^2\right)^{5/2}} \; dV = \int_B \frac{y^2}{\left(m^2+x^2+y^2+z^2\right)^{5/2}} \; dV = \int_B \frac{z^2}{\left(m^2+x^2+y^2+z^2\right)^{5/2}} \; dV \,,$$ therefore we get some cancellation and $$\int_B \frac{m^2-2 x^2+y^2+z^2}{\left(m^2+x^2+y^2+z^2\right)^{5/2}} \; dV = \int_B \frac{m^2}{\left(m^2+x^2+y^2+z^2\right)^{5/2}} \; dV = \frac{4 \pi R^3}{3 \left(m^2+R^2\right)^{3/2}}$$ which converges to $$4\pi/3$$ as the radius $$R$$ goes to infinity.

But maybe its attraction is as a Siren leading sailors into a shipwreck.

• The integration limits you specify for z will be imaginary for some values of x and y. I don't know if this matters. Commented Nov 17, 2019 at 18:00

The integral is indeed equal to $$\iiint (f_x+f_y+f_z)=4\pi$$ where $$f_i=\partial_i(p_i/(p^2+m^2)^{3/2})$$. This is easy to prove using spherical symmetry and e.g. the Gauss theorem (the integral is basically the residue at infinity, and so independent of $$m$$).

The integral is perfectly convergent; indeed, it is easy to see that $$(f_x+f_y+f_z)\sim 1/r^5$$:

Div[{px, py, pz}/(px^2 + py^2 + pz^2 + m^2)^(3/2), {px, py, pz}] /. {px -> r Cos[θ] Sin[ϕ], py -> r Cos[θ] Cos[ϕ], pz -> r Sin[θ]} // FullSimplify
Series[%, {r, ∞, 4}]
(* O[1/r]^5 *)


The problem is that the partial integrals $$\iiint f_i$$ do not exist individually. Indeed, they are $$f_i\sim 1/r^3$$:

D[px/(px^2 + py^2 + pz^2 + m^2)^(3/2), px] /. {px -> r Cos[θ] Sin[ϕ], py -> r Cos[θ] Cos[ϕ], pz -> r Sin[θ]} // FullSimplify
Series[%, {r, ∞, 2}]
(* O[1/r]^3 *)


(This, together with $$\mathrm d\boldsymbol p=4\pi p^2\mathrm dp$$ means that the integrand is $$\sim 1/r$$, which is not integrable).

Unfortunately, Mathematica was not able to identify the divergence of the integral: the result it yields is just meaningless.

• So Fubinis theorem does not apply here? Interesting. Commented Nov 17, 2019 at 14:34

This appears to be a case where you genuinely cannot change the order of integration. I think it's a Mathematics problem not a Mathematica problem.

Define relevant assumptions

\$Assumptions = {px^2 > 0, py^2 > 0, pz^2 > 0};


Evaluate and simplify the integrand

expr =
D[px (px^2 + py^2 + pz^2 + m^2)^(-3/2), px] /. m -> 1 // FullSimplify
(* (1 - 2 px^2 + py^2 + pz^2)/(1 + px^2 + py^2 + pz^2)^(5/2) *)


The integral w.r.t. px is zero. (This can be verified easily by specifying numerical values for py and pz)

Integrate[expr, {px, -∞, ∞}]
(* 0 *)


Integrating w.r.t. py and pz

Integrate[expr, {py, -∞, ∞}, {pz, -∞, ∞}]
(* (2 π)/(1 + px^2)^(3/2) *)


Since the integrand is always positive, this is unsurprisingly non-zero

Integrate[%, {px, -∞, ∞}]
(* 4 π *)


Not a full answer, but we can see Mathematica returns a different answer, depending on the order in which the integration is performed.

Assuming[m^2 > 0,
Integrate[
D[px (px^2 + py^2 + pz^2 + m^2)^(-3/2),
px], {px, -∞, +∞}, {py, -∞, +∞}, {pz, -∞, +∞}]]
(* 4 π *)

Assuming[m^2 > 0,
Integrate[
D[px (px^2 + py^2 + pz^2 + m^2)^(-3/2),
px], {pz, -∞, +∞}, {px, -∞, +∞}, {py, -∞, +∞}]]
(* 0 *)

• why is the order of this integration matters? This is a very smooth function since m**2>0. I don't see any reasons why I can't change the order. Commented Nov 16, 2019 at 13:20