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, 2013 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$ Feb 18, 2013 at 16:04
  • 2
    $\begingroup$ Actually, in this case, it is even much simpler, your integral factorizes into a product of simple integrals. $\endgroup$ Feb 18, 2013 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, 2013 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, 2013 at 2:20

4 Answers 4


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.


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.

  • $\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, 2013 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, 2013 at 0:14

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.


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 ;-)


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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