As in the title I have to evaluate this triple integral:


I've been trying to solve this since a week ago.

The first thing I've done was understand the meaning of the integral. I think this integral represents the mass (as an example) of a unitary cube which contains materials of different density.

The values of the materials density are, point for point, the inverse of the square of spheres centered in the origin plus one.

The max value of the density is $1$ in the origin of the cube and the min value is $\frac{1}{16}$ on the opposite vertex.

I suppose that the value of the integral is $\frac{\pi^2}{32}$

I've tried to use simple substitutions without any results, so I tried to change the coordinates with spherical and cylindrical systems. The spherical coordinates give me an incredibly long sum of integrals and I doubt that they're all integrable as elementary functions.

The cylindrical gives me the following result


which I'm not able to solve.

My instinct tells me that there is a trick in some steps where I can observe that a difficult integral actually is exactly half of another one simpler but I can't figure out where.

I'll appreciate any kind of suggestions.

  • $\begingroup$ @Squirtle: actually I believe it is simpler to avoid spherical coordinates and just invoke an integral transform and the error function. $\endgroup$ Oct 14, 2018 at 23:05
  • $\begingroup$ I do not understand the Erf answer below. However, the indefinite integral is relatively easy to compute using spherical coordinates. If one can slink back to rectangular coordinates for the definite integral evaluation, that should work. I tried it, but the algebra is nasty. But in principle it should work. $\endgroup$
    – Ben W
    Oct 14, 2018 at 23:23
  • $\begingroup$ I've edited your question; normally you should state clearly that you have posted an attempt/solution if you have done so in an answer, otherwise it's easy to miss it. I've voted to reopen as well. $\endgroup$
    – user21820
    Sep 11, 2019 at 7:47
  • 1
    $\begingroup$ @Jack perhaps OP used wolfram differently to confirm the answer, not sure how posting the screenshot with an aproximation of the answer helps here. $\endgroup$
    – Zacky
    Sep 11, 2019 at 15:32
  • 1
    $\begingroup$ @Zacky: Indeed I was not able to reproduce the result OP claims in Mathematica, neither the exact result in Wolfram Alpha. The approximation could be one way to provide context (although seems to be redundant in this case) and "confirm" numerically the claimed result. A screen shot from OP would be more helpful. $\endgroup$
    – user9464
    Sep 11, 2019 at 15:46

4 Answers 4


Here's an intuitive solution involving surfaces of hypercubes.

Let $$v := \iiint_{[0,1]^3}\frac{dxdydz}{(1+x^2+y^2+z^2)^2}$$

be the value you're looking for.

Suppose a point-shaped light source shines uniformly from the origin of 4D space with intensity of $1$ watt per cubic radian. Let $p$ be a parametrization of a cubic piece $P$ of hyperplane, specifically $p(x,y,z) = (1,x,y,z)$ for $x, y, z \in [0,1]$.

It is clear that no part of this piece of hyperplane is in shadow. If we compute the irradiance $I(x,y,z)$ of arriving light in watts per cubic unit, we find

$$I(x,y,z) = \frac{1}{(1+x^2+y^2+z^2)^2}.$$

If we integrate this irradiance over the surface of $P$, we find that the total power arriving at $P$ is $v$ watts per cubic unit.

Now consider the surface of a hypercube with corners at $(\pm 1, \pm 1, \pm 1, \pm 1)$. We know that $P$ is an octant of one of the eight cells of this surface, and every other octant receives the same irradiance, so we know that the total power arriving at the hypercube is $64v$. However, the hypercube captures every bit of light emitted by our light source, so we know that $64v$ equals $2\pi^2$, the total wattage of light emitted by our light source.

Therefore $v = \frac{2\pi^2}{64} = \frac{\pi^2}{32}$.

  • 3
    $\begingroup$ It's very nice to see how integrals can solve "real life" problems. $\endgroup$
    – Zacky
    Jul 25, 2019 at 13:11
  • $\begingroup$ Technically it's the other way around, I found a "real life" situation and used it to solve an integral. $\endgroup$
    – Magma
    Jul 25, 2019 at 13:13
  • $\begingroup$ Ingenious! So basically, the result is just a piece of surface area of a hypersphere! $\endgroup$
    – orion
    Jul 25, 2019 at 14:59
  • $\begingroup$ Yeah, you could probably also formulate it like that. If you take the parametrization $(x, y, z) \mapsto \frac 1{\sqrt{1+x^2+y^2+z^2}}(1, x, y, z)$ of a bit of $3$-sphere, then $\frac {dx\,dy\,dz}{(1+x^2+y^2+z^2)^2}$ would be the surface element of the parametrization. $\endgroup$
    – Magma
    Jul 25, 2019 at 18:01

Let us denote the integral by $I$, $$I=\int_{0}^{1} \int_{0}^{1} \int_{0}^{1} \frac{dx dy dz}{(1+x^2+y^2+z^2)^2}=\int_{0}^{1} \int_{0}^{1} \int_{0}^{1} \int_{0}^{\infty} t \exp[-(1+x^2+y^2+z^2)t] dx dy dz dt$$ Due to the symmetry of $I$ in $x,y,z$; we can write $$I=\int_{0}^{\infty} dt ~t e^{-t} \left ( \int_{0}^{1} e^{-tx^2} dx \right)^3=\frac{\pi \sqrt{\pi}}{8} \int_{0}^{\infty} t e^{-t} \mbox{erf}~^3(\sqrt{t}) dt= \frac{\pi^2}{8} \int_{0}^{1} v^3 dv=\frac{\pi^2}{32}. $$ Here $\mbox{erf}(x)=\frac{2}{\sqrt{\pi}} \int_{0} ^{x} e^{-t^2} dt$ and we have used $\mbox{erf}(\sqrt{t})=v$.

I also post the screen shot of from Mathematica,

  • 8
    $\begingroup$ Why did you ask a question, then answer it yourself less than 10 minutes later in the answer section and not in your original question? $\endgroup$ Jul 25, 2019 at 12:36
  • 2
    $\begingroup$ @Mattos The question is still open for other solutions. $\endgroup$
    – Z Ahmed
    Jul 25, 2019 at 12:39
  • 1
    $\begingroup$ @Mattos We don't know how long has he being trying the problem before asking. $\endgroup$
    – ajotatxe
    Jul 25, 2019 at 12:47
  • 5
    $\begingroup$ Answering ownself question should be encouraged. However it would be nice to give some time to others to solve it before posting the answer. $\endgroup$
    – Zacky
    Jul 25, 2019 at 13:02
  • 2
    $\begingroup$ @Zacky Seeing multiple solutions is also a fun. $\endgroup$
    – Z Ahmed
    Jul 25, 2019 at 13:05

Since for any $a>0$ we have $\frac{1}{a^2}=\int_{0}^{+\infty} w e^{-aw}\,dw$, by Fubini's theorem the original integral can be written as

$$ \int_{0}^{+\infty} w e^{-w}\left(\frac{\sqrt{\pi}\,\text{Erf}(\sqrt{w})}{2\sqrt{w}}\right)^3\,dw\stackrel{w\mapsto w^2}{=}\frac{\pi\sqrt{\pi}}{4}\int_{0}^{+\infty}e^{-w^2}\text{Erf}^3(w)\,dw $$ and the RHS is clearly $$ \frac{\pi\sqrt{\pi}}{4}\left[\frac{\sqrt{\pi}}{8}\,\text{Erf}^4(w)\right]_{0}^{+\infty} =\color{red}{\frac{\pi^2}{32}}.$$

Interesting side-note: the same approach in dimension $2$ gives a relation between $\iint_{(0,1)^2}\frac{dx\,dy}{(1+x^2+y^2)^2}=\frac{1}{\sqrt{2}}\arctan\frac{1}{\sqrt{2}}$ and $\int_{0}^{+\infty}\left(1-\text{Erf}^3(w)\right)\,dw$.
In particular it proves that $$ \int_{0}^{+\infty}\left(1-\text{Erf}^3(x)\right)\,dx = \frac{6\sqrt{2}}{\pi\sqrt{\pi}}\,\arctan\frac{1}{\sqrt{2}} $$ which is not recognized by Mathematica (or at least by my version).

  • $\begingroup$ I'm interested in this question, but I do not understand your answer. How on earth does Fubini's Theorem give us the error function and all that jazz? $\endgroup$
    – Ben W
    Oct 14, 2018 at 23:25
  • $\begingroup$ Nevermind, I see it now. Cool! $\endgroup$
    – Ben W
    Oct 14, 2018 at 23:26
  • 3
    $\begingroup$ Just write it down: $$ \iiint_{(0,1)^3}\frac{d\mu}{(1+x^2+y^2+z^2)^2}=\int_{0}^{+\infty}\iiint_{(0,1)^3} w e^{-w(1+x^2+y^2+z^2)}\,dx\,dy\,dz\,dw = \int_{0}^{+\infty}we^{-w}\left(\int_{0}^{1}e^{-wu^2}\,du\right)^3\,dw.$$ There is just an unusual "dimensional tour" $3\mapsto 4 \mapsto 1$. $\endgroup$ Oct 14, 2018 at 23:27
  • 1
    $\begingroup$ You blowed up my mind, thank so much! That's a so simple solution but still so incredibly beautiful. You rock! $\endgroup$
    – Anonymous
    Oct 15, 2018 at 7:43

Another solution:

Let us use $x=r \cos \phi,~ y= r \sin \phi$ and write $$I=\int_{0}^{1}\int_{0}^{1} \int_{0}^{1} \frac{dx dy dz}{(1+x^2+y^2+z^2)^2}= 2\int_{0}^{1} \int_{0}^{\pi/4} \int_{0}^{\mbox{sec}\phi} \frac{r dr d\phi dz}{(1+r^2+z^2)^2}~~~(1).$$ $$ \Rightarrow I =\int_{0}^{1} \int_{0}^{\pi/4} \left(\frac{1}{1+z^2}-\frac{1}{(1+z^2+\mbox{sec}^2\phi)}\right) d\phi dx=\int_{0}^{1} \int_{0}^{1} \frac{1}{(1+z^2)(2+z^2+t^2)} dz~ dt,~~~~(2)$$ $t=\mbox{sec}\phi$ taken here. Next, interchange $t$ and $z$ to get $$I=\int_{0}^{1} \int_{0}^{1} \frac{dz dt}{(1+t^2)(2+z^2+t^2)}~~~~(3)$$ Adding (2) and (3), we get $$I=\frac{1}{2} \int_{0}^{1} \int_{0}^{1} \frac{1}{(2+z^2+t^2)} \left( \frac{1}{1+z^2}+\frac{1}{1+t^2} \right)dx~dt=\frac{1}{2} \int_{0}^{1} \int_{0}^{1} \frac{ dz~ dt}{(1+z^2)(1+t^2)}= \frac{1}{2} \left ( \int_{0}^{1}\frac {dz}{1+z^2} \right)^2=\frac{\pi^2}{32}.$$


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .