# NIntegrate with vector variable

Here's a Gaussian integral

NIntegrate[E^-(x^2 + y^2) E^(Sqrt[3] x y), {x,y} ∈ FullRegion[2]]


which evaluates to 2π // N. This is fine for a 2D integral, but to extend this to higher dimensions, it would be more convenient to write it as

NIntegrate[E^-r.r E^(Sqrt[3] Most[r].Rest[r]), r ∈ FullRegion[2]]


This evaluates to the correct answer, but gives errors saying that Most and Rest have been applied to the atomic expression r. This is presumably happening because NIntegrate tries to evaluate its first argument before substituting in values for r.

What's the most sensible way to avoid these error messages? (I would rather not just silence them, in case it leads to unanticipated problems with a more complicated example.)

(The same is true with Integrate instead, but I suspect there will be different solutions for the two cases, and I'm more interested in NIntegrate.)

• It would be nice to learn what is behind the structure of the second factor of your integrand? Where does the integral come from? – Dr. Wolfgang Hintze Aug 17 '16 at 22:14
• @Dr.WolfgangHintze With appropriately chosen coefficients, this is the partition function of a set of harmonic oscillators coupled by springs. (The second term comes from expanding (x[i]-x[i+1])^2.) I just chose Sqrt[3] to make the integral give a nice result, but it was obviously a bad choice, as your answer shows. – Stephen Powell Aug 19 '16 at 9:59

First of all let us introduce the dimensionalty of the problem explicitly.

In the example case of three dimensions:

d = 3;


the integration variable is defined as the vector

r = Array[vr, d]

(* Out[2]= {vr[1], vr[2], vr[3]} *)


and the integral becomes

f3[a_] := Integrate[E^(-r.r) E^(a Most[r].Rest[r]), r \[Element] FullRegion[d]]


Here we have replaced Sqrt[3] by some real number "a", because, as we shall see, the integral for d = 3 diverges for a^2>=2.

Examples

f3[0]

(* Out[20]= \[Pi]^(3/2) *)

f3[1]

(* Out[21]= Sqrt[2] \[Pi]^(3/2) *)

f3[Sqrt[2]]

(* Out[22]= \[Infinity] *)


The more the the integral diverges for a = Sqrt[3]

Writing the integral more explicitly we find the explicit formula for f3[a]:

Integrate[E^(-r.r) E^(a Most[r].Rest[r]), {vr[
1], -\[Infinity], \[Infinity]}, {vr[2], -\[Infinity], \[Infinity]}, {vr[
3], -\[Infinity], \[Infinity]}]

(* Out[26]= ConditionalExpression[(Sqrt[2] \[Pi]^(3/2))/(
Sqrt[4 - a^2] Sqrt[(-2 + a^2)/(-4 + a^2)]), Re[(-2 + a^2)/(-4 + a^2)] > 0] *)


Hence

f3x[a_] = %[[1]]

(* Out[27]= (Sqrt[2] \[Pi]^(3/2))/(Sqrt[4 - a^2] Sqrt[(-2 + a^2)/(-4 + a^2)]) *)


To address the programming issue of restricting integrand evaluation to vector r, you need to define a function with pattern-restricted argument for your integrand.

Block[{i},
i[r_?VectorQ] := E^-r.r E^(Sqrt[3] Most[r].Rest[r]);
NIntegrate[i[r], r ∈ FullRegion[2]]
]
(*  6.28318  *)


Beware that this can come at a cost, although not in the OP's test case. The use of the pattern-restricted function i[r] essentially makes the integrand a numeric black box. NIntegrate cannot symbolically analyze the integrand and choose an appropriate integration method automatically. The responsibility for this falls on the user.

The use of ?VectorQ is similar to the use of ?NumericQ discussed here.

• @Miachel E2 Oh yes, this restriction to vectors was what I was looking for (and did "by hand"). – Dr. Wolfgang Hintze Aug 17 '16 at 22:15

Here's a possible solution (with Sqrt[3] replaced by 1):

With[{r = Array[r, 2]}, NIntegrate[E^-r.r E^(Most[r].Rest[r]), r ∈ FullRegion[Length[r]]]]