# How do I integrate over this region?

I want to compute this integral: $$\int_{\Omega}e^{i\sum_{n=1}^{N}\frac{(p_n+p_{n-1})(q_n-q_{n-1})}{2}}dq_1\cdots dq_{N-1}dp_1\cdots dp_{N-1}\,,$$ where $$\Omega=\{|q_n-q_{n-1}|<\epsilon,\, |p_n-p_{n-1}|<\epsilon:1\le n\le N\}$$ and where $$p_0,q_0,p_N,q_N$$ are fixed. How would I implement this integral in Mathematica? I need $$N$$ to be arbitrary and so I have an arbitrary number of variables I need to define.

For the case $$N=2$$ and for arbitrarily chosen valued of $$p_0,q_0,p_2,q_2$$ this can be implemented with the code:

region = ImplicitRegion[Abs[b-0] < 1 && Abs[0.1-b] < 1 && Abs[y-0.5]<1 && Abs[1-y]<1, {y, b}]

f[y_,b_]=Exp[I*((1+y)*(0.1-b)/2+(y+0.5)*(b)/2)]

Integrate[f[y,b],y,b,{y,b}\[Element]region]


Here $$y,b$$ represents $$p_1,q_1\,.$$

• Posting the code for $N=2$ would be helpful.
– JimB
Sep 22, 2023 at 22:38
• Thanks for your advice. I just included the code for $N=2\,.$
– JLA
Sep 23, 2023 at 5:35
• Can you kindly precise the description of $\Omega$? Also don't switch to other variables than $p_i,q_j$. Sep 23, 2023 at 7:56
• It looks like it can be done analytically, or, at least written as a product of 2D integrals. Just change to the new difference variables. The Jacobian matrix is a block diagonal with just 2 nonzero diagonals (Toeplizt matrix). Sep 23, 2023 at 22:48
• @yarchik I'd like to do it analytically but I think it may be a bit more complicated than that, for example for $N=2$ the domain (in the $q$ variables) is $|q_1-q_0|<\epsilon, |q_2-q_1|<\epsilon\,.$ I can let $q'=q_1-q_0$ and the new domain (in the $q$ variables) is $|q'|<\epsilon, |q_2-q'-q_0|<\epsilon.$ And then the $p$ variables appear as $p_n+p_{n-1}$ in the integrand, not $p_n-p_{n-1}\,.$ Though maybe you see something that I don't (which I'd be interested in, I'll continue thinking about it).
– JLA
Sep 23, 2023 at 23:45

Change Integrate - command to Integrate[f[y, b], {y, b} \[Element] region] and use symbolic numbers ( 1/2instead of .5,...)

region = ImplicitRegion[Abs[b - 0] < 1 && Abs[1/10 - b] < 1 && Abs[y- 1/2] < 1 &&Abs[1 - y] < 1, {y, b}]

f[y_, b_] = Exp[I*((1 + y)*(1/10 - b)/2 + (y + 1/2)*(b)/2)]

Integrate[f[y, b], {y, b} \[Element] region]
(*80 E^(-(I/5)) (-1 + E^((3 I)/40) - 2 I E^((41 I)/80) Sin[3/80])*)
N [%] (*2.81469 + 0.211498 I*)


In the general case this can be done as follows (Your notations are somewhat changed because N is reserved in WL.).

n = 5; p[0] = 1/2; q[0] = 0; q[n] = 1/10; p[n] = 1; \[Epsilon] = 1;
integrand = Exp[I*Sum[(p[j] + p[j - 1])*(q[j] - q[j - 1])/2, {j, 1, n}]];
reg1 = ImplicitRegion[ Reduce[Table[RealAbs[p[j] - p[j - 1]] < \[Epsilon], {j, 1, n}],
Table[p[j], {j, 1, n - 1}], Reals] // FullSimplify, Evaluate[Table[p[j], {j, 1, n - 1}]]];
reg2 = ImplicitRegion[ Reduce[Table[RealAbs[q[j] - q[j - 1]] < \[Epsilon], {j, 1, n}],
Table[q[j], {j, 1, n - 1}], Reals] // FullSimplify,
Evaluate[Table[q[j], {j, 1, n - 1}]]];
NIntegrate[integrand, Table[p[j], {j, 1, n - 1}] \[Element] reg1,
Table[q[j], {j, 1, n - 1}] \[Element] reg2, PrecisionGoal -> 2, AccuracyGoal -> 2]


69.6059 + 6.12144 I

Replacing n=5; by n=2;, one obtains 2.81469 + 0.211498 I. Some comments: the set of the integration is the Cartesian product of two sets, but the integral does not split in the product of two integrals because of (p[j] + p[j - 1])*(q[j] - q[j - 1]) in the exponent.

• This is a nice answer, thank you. I think it's the unfortunate case that mathematica can't compute the integral for large $N$ values (like $N>10$), or analytically.
– JLA
Sep 23, 2023 at 23:55
• @JLA: For $n=10$ we deal with 18 dimensions and NIntegrate produces "DiscretizeRegion::cdim: The region given at position 1 in DiscretizeRegion[RegionProduct[ImplicitRegion[(-(1/2)<p[1]<=0&&((Inequality[<<5>>]&&Or[<<2>>])||(Inequality[<<5>>]&&Or[<<2>>])))||(<<1>>),{p[1],<<8>>}],<<1>>]] is in dimension 18. DiscretizeRegion only supports dimensions 1 through 3." and "NIntegrate::femrdim: FiniteElement method can only be applied to regions of embedding dimension 1, 2, or 3.". Try Method->"MonteCarlo" option in this case. Good luck! Sep 24, 2023 at 5:36

The integral goes over the product of parallels on both sides of diagonals

$$q_{n}-q_{n-1<\epsilon}$$

The exponential is a product, the volume 2n-form, too.

Formulate the integral in the difference variable over all of the 2n-dimensional symmetric hypercube using the difference variables

$$Q_n = q_{n}-q_{n-1}, \ P_n = p_n- p_{n-1}$$

      With[{\[CurlyEpsilon] = 1/4},
RegionPlot[((a < b && a + \[CurlyEpsilon] > b) || (a > b &&
b + \[CurlyEpsilon] > a)), {a, -1, 1}, {b, -1, 1}]]


and replace the domain of integration by multiplication of the characteristic functions

  IndicatorFunction[x_Real, {a_Real, b_Real}] :=
HeavisidePi[(x - (a + b)/2)/(b - a) ]

\[Chi]= IndicatorFunction


$$\text{dint}=\prod _{k=1}^3 dq_k \prod _{l=1}^3 dp_l \left(\prod _{n=1}^3 e^{\frac{1}{2} i \left(p_{n-1}+p_n\right) \left(q_n-q_{n-1}\right)}\right) \left(\prod _{m=2}^3 \chi \left(q_m-q_{m-1},\{-\varepsilon ,\varepsilon \}\right)\right) \left(\prod _{l=2}^3 \chi \left(p_l-p_{l-1},\{-\varepsilon ,\varepsilon \}\right)\right)$$

\begin{align}& dp_1 dp_2 dp_3 dq_1 dq_2 dq_3 \\& * \exp \left(\frac{1}{2} i \left(p_0+p_1\right) \left(q_1-q_0\right)+\frac{1}{2} i \left(p_1+p_2\right) \left(q_2-q_1\right)+\frac{1}{2} i \left(p_2+p_3\right) \left(q_3-q_2\right)\right) \\& * \chi \left(p_2-p_1,\{-\varepsilon ,\varepsilon \}\right) \chi \left(p_3-p_2,\{-\varepsilon ,\varepsilon \}\right) \chi \left(q_2-q_1,\{-\varepsilon ,\varepsilon \}\right) \chi \left(q_3-q_2,\{-\varepsilon ,\varepsilon \}\right)\end{align}