3N-dimensional integral

Question

I'm trying to compute a multidimensional integral with a variable number of dimensions.

The integral is as follows: $$ \int d^{3N}\!p~e^{-\frac{\beta}{2m}\vec p^2}. $$

I have tried this

Integrate[e^(-a*{p1,p2,p3}^2),{{p1,p2,p3}^N,-Infinity,Infinity}]

but it's not working at all. I guess I have to determine the vector before, or?

More generally, to what extent can Mathematica compute such integrals for arbitrary $N$ and arbitrary integrands in place of $e^{-\frac{\beta}{2m}\vec p^2}$?

Answer 1

For integrands where Fubini's Theorem is applicable (and that would be the great majority), we may iterate the integral. The following solution--intended to work with fixed limits of integration in all dimensions--does this with two definitions: one to perform the integration over the last variable and another to get everything started.

integrateND[f_, limits_] /; Length[limits] > 1 := Module[{x, n = Length[limits]},
   Integrate[integrateND[f[##, x] &, Most[limits]], {x, limits[[-1, 1]], limits[[-1, 2]]}]];
integrateND[f_, limits_]  := Module[{x}, Integrate[f[x], {x, limits[[1, 1]], limits[[1, 2]]}]];

The arguments of integrateND are a function of an arbitrary number $n$ variables and a list of $n$ limits (lower and upper, each given as a list). For example, the integration in the question with $N=2$ (which is a six-dimensional integral) can now be performed as

With[{n = 2}, 
 Assuming[\[Beta] / (2 m) > 0, 
  integrateND[Exp[-\[Beta] / (2 m) {##}.{##}] &, ConstantArray[{-Infinity, Infinity}, 3 n]]]]

$\frac{8 \pi ^3 m^3}{\beta ^3}$

High-dimensional multiple integrals should be invoked judiciously: they can take a very long time to compute or not be computable at all.

Answer 2

Although the original question is a little too narrow, I'm going to interpret it as asking about general multidimensional integrals in which Gaussians appear together with arbitrary factors in the integrand:

$$\iiint f(\,\vec{r}\,)\,\exp(\,-\frac{1}{2}\vec{r}^T\Sigma^{-1}\vec{r}\,)\,dx\,dy\,dz$$

Here $\Sigma^{-1}$ could be a symmetric positive definite matrix, or in the simplest case just $\sigma^2$ times the identity matrix. The dimension can be anything, but here I'll just specialize to three.

Then you may be interested in this observation which speeds up such integrals tremendously in comparison with straightforward multiple integration:

Edit: An integrand that doesn't factor

Start by defining the Gaussian in three dimensions:

f = PDF[MultinormalDistribution[{0, 0, 
    0}, \[Sigma] IdentityMatrix[3]], {x, y, z}]

$\frac{e^{\frac{1}{2} \left(-\frac{x^2}{\sigma }-\frac{y^2}{\sigma }-\frac{z^2}{\sigma }\right)}}{2 \sqrt{2} \pi ^{3/2} \sqrt{\sigma ^3}}$

Try the following as a triple integral:

Timing[Integrate[
  Integrate[Integrate[
    (x^2 + y^3)/(1 + z^2)
      f, {x, -Infinity, Infinity}], {y, -Infinity, \
Infinity}], {z, -Infinity, Infinity}]
 ]

$\left\{1.48519,\text{ConditionalExpression}\left[\frac{\sqrt{\frac{ \pi }{2}} e^{\frac{1}{2 \sigma }} \sqrt{\sigma ^3} \text{erfc}\left(\frac{1}{\sqrt{ 2} \sqrt{\sigma }}\right)}{\sigma },\Re(\sigma )>0\right]\right\}$

Now do the same using Expectation, keeping in mind the different normalization of the result:

Timing[Expectation[(x^2 + y^3)/(
  1 + z^2), {x, y, z} \[Distributed] 
   MultinormalDistribution[{0, 0, 0}, \[Sigma] IdentityMatrix[3]]]
 ]

$\left\{0.775868,\sqrt{\frac{\pi }{2}} e^{\frac{1}{2 \sigma }} \sqrt{\sigma } \text{erfc}\left(\frac{1}{\sqrt{ 2} \sqrt{\sigma }}\right)\right\}$

So the result is obtained about twice as fast using Expectation.

Fourier transforms as a special case

My initial example was the Fourier transform of a Gaussian,

$$\iiint \exp(- i\, \vec{k}\cdot \vec{r}\,)\,\exp(\,-\frac{1}{2}\vec{r}^T\Sigma^{-1}\vec{r}\,)\,dx\,dy\,dz$$

where the speed advantage of Expectation compared to multiple Integrate with infinite integration boundaries is more significant than above:

Expectation[
 Exp[-I {kx, ky, kz}.{x, y, z}], {x, y, z} \[Distributed] 
  MultinormalDistribution[{0, 0, 0}, σ IdentityMatrix[3]]]

$e^{-\frac{1}{2} \sigma \left(\text{kx}^2+\text{ky}^2+\text{kz}^2\right)}$

This works for symbolic integrands (note I didn't specify a value for $\sigma$ or the components of $\vec{k}$), and you don't have to specify the list of integration boundaries explicitly either. The timing for the above is about 7 times faster than doing the calculation using three nested Integrates, where I also have to add assumptions for the parameters (specifically $\sigma>0$, $\vec{k}$ real).

Although my initial claim is correct that this is faster than the straightforward triple integral, one can do a lot better in this special case by just using FourierTransform directly:

Timing[FourierTransform[f, {x, y, z}, {kx, ky, kz}]]

$\left\{0.007332,\frac{e^{-\frac{1}{2} \sigma \left(\text{kx}^2+\text{ky}^2+\text{kz}^2\right)} }{2 \sqrt{2} \pi ^{3/2} \left(\frac{1}{\sigma }\right)^{3/2} \sqrt{\sigma ^3}}\right\}$

Interestingly, though, we can even beat this really fast result using the multinormal distribution again:

Timing[CharacteristicFunction[
  MultinormalDistribution[{0, 0, 0}, \[Sigma] IdentityMatrix[3]], {kx,
    ky, kz}]
 ]

$\left\{0.000185,e^{\frac{1}{2} \left(\text{kx}^2 (-\sigma )-\text{ky}^2 \sigma -\text{kz}^2 \sigma \right)}\right\}$

Conclusion

For Gaussian integrals over all space (or momentum space, as in the question), the approach using MultinormalDistribution is complementary to whuber's solution: general Gaussian integrals can be evaluated by using Expectation and similar tools for probability distributions, such as CharacteristicFunction.

Answer 3

For your particular example, you could use FullRegion[n] as the integration domain:

Integrate[Exp[-(β/(2m)) p.p], p ∈ FullRegion[6], Assumptions->β>0 && m>0]

(8 m^3 π^3)/β^3

which is the same as whuber's answer.

Answer 4

Stealing from @carl's answer, we can do it for a few dimensions and 'guess' the trend.

tt = Table[
  Integrate[Exp[-(\[Beta]/(2 m)) p.p], p \[Element] FullRegion[i], 
   Assumptions -> \[Beta] > 0 && m > 0], {i, 2, 10}]

it seems this reproduces it quite well:

tt2 = Table[(\[Beta]/(2 m \[Pi]))^(-i/2) , {i, 2, 10}];
tt/tt2 // PowerExpand;

{1,1,1,1,1,1,1,1,1}

So the answer seems to be for arbitrary $N$

((\[Beta]/(2 m \[Pi]))^(-N/2))

Of course this is not a proof ;-)