# 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. – mikado Nov 16 at 13:34
• @mikado Yes, m=1. And the order of integration variables now does not affect the result. – ZHANG Juenjie Nov 16 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 – Michael E2 Nov 16 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. – Anton Antonov Nov 17 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. – mikado Nov 17 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. – lalala Nov 17 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. – ZHANG Juenjie Nov 16 at 13:20