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


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}$?

  • $\begingroup$ Integrate only supports integrating with respect to scalar variables, i.e. you need to integrate by the components of the vector. See the doc page on Integrate on how, and pay attention to spelling of names and capitalisation. (Integrate, E^x, and don't use N as a variable because it's a built-in symbol) $\endgroup$ – Szabolcs Feb 18 '13 at 15:59
  • $\begingroup$ You don't need Mathematica to compute this. Switch to the generalized spherical coordinates and integrate over r and the angles. This is exmplained in many texts on Quantum Field Theory, look up dimensional regularization. $\endgroup$ – Leonid Shifrin Feb 18 '13 at 16:04
  • 2
    $\begingroup$ Actually, in this case, it is even much simpler, your integral factorizes into a product of simple integrals. $\endgroup$ – Leonid Shifrin Feb 18 '13 at 16:11
  • $\begingroup$ I have voted to reopen because I believe there is a legitimate interpretation of the question, as suggested by my edit of it (and I recall seeing similar issues raised on this site but do not recall seeing a general answer). $\endgroup$ – whuber Feb 19 '13 at 20:26
  • $\begingroup$ @LeonidShifrin I'm fairly sure dimensional reguralization isn't necessary for multidimensional gaussian integrals :) It factorizes, also, you can go to generalized spherical coordinates. However this seems to really be asking "how do I set up multidimensional integrals". In any case I think the real answer is "you're better off with some sort of monte carlo scheme for large N and arbitrary integrands" $\endgroup$ – acl Feb 20 '13 at 2:20

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:

    (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:

 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:

  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\}$


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.

| improve this answer | |
  • $\begingroup$ For this problem, we can use the Expection and distribution to calculate it. But it is still not fast for the large dimension and complex gaussian form. Because we usuall have the anlytical result for the Gaussian integrals, can we find some package for dealing with this kind of problem? $\endgroup$ – Orders Mar 19 '13 at 8:21
  • $\begingroup$ @Orders I added some more speed comparisons. I would say that the package you're looking for is probably just the set of statistical capabilities built into Mathematica. Maybe there's other packages thatI'm not aware of. $\endgroup$ – Jens Mar 20 '13 at 0:14

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.

| improve this answer | |

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.

| improve this answer | |

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}]

enter image description here it seems this reproduces it quite well:

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


So the answer seems to be for arbitrary $N$

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

Of course this is not a proof ;-)

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.