# (Numerical) calculation of nested integrals

The problem

I'd like to (numerically) compute nested integrals of the type $$I(a;k) \equiv \left( \prod_{i=1}^k \int_0^\infty \mathrm{d}x_i \right) \exp \left( - \frac{1}{4a} \sum_{i=1}^k \left( x_i - x_{i+1} \right)^2 - \frac{a}{2} \sum_{i=1}^k (x_i + x_{i+1}) \right)$$ where $$a > 0$$ is a parameter and $$x_{k+1} \equiv x_1$$.

Analytic result for $$k=1,2$$

For $$k = 1$$ and $$k = 2$$ the integrals can be done analytically with results $$I(a;1) = 1/a, ~~~ I(a;2) = \sqrt{\frac{\pi}{2a}} ~ e^{a^3/2} ~ \text{Erfc}\left( \frac{a^{3/2}}{\sqrt{2}} \right),$$ but as far as I'm aware no closed-form expression is known for $$k > 2$$. One suggestion discussed in the comments was to write the argument of the exponential as a quadratic form $$\boldsymbol{x}^\top A \boldsymbol{x} + \boldsymbol{J}^\top \boldsymbol{x}$$ with $$A/\boldsymbol{J}$$ some $$a$$-dependent symmetric matrix/vector and attempt to use the same logic as in the derivation of the standard multiple-Gaussian integral (that is, orthogonally diagonalizing $$A = O^\top D O$$ and changing variables $$\boldsymbol{y} = O \boldsymbol{x}$$). This is complicated however due to the range of integration being $$[0,\infty)^k$$ and not $$(-\infty,\infty)^k$$ -- these are not Gaussian integrals.

My code so far and why it's wrong

So far I have written the following NIntegrate command:

int[a_, k_] :=
NIntegrate[
Exp[-1/(4 a) (Sum[(x[i] - x[i + 1])^2, {i, 1, k - 1}] + (x[k] -
x[1])^2) -
a/2 (Sum[(x[i] + x[i + 1]), {i, 1, k - 1}] + (x[k] + x[1]))],
Sequence @@ Table[{x[i], 0, \[Infinity]}, {i, k}] // Evaluate,
AccuracyGoal -> 3];


For values $$a$$ of order unity the integral should be of order unity, hence the unambitious AccuracyGoal. An accuracy of $$10^{-3}$$ for $$a = O(1)$$ would suffice for my purposes though it would be nice if it could be systematically improved. It seems like the integrals become more difficult to do when $$a$$ is small -- let's say we'd like $$10^{-1} < a < 1$$ but again ways to systematically enlarge the range would be nice.

To see that the above code gives the wrong answer we can compare this numerics with the exact result for the $$k=2$$ case:$k=2$." />

For $$a \lesssim 0.3$$ the numerics clearly start failing beyond the error bar allowed by the AccuracyGoal setting, though Mathematica produces no error message. For $$k=3,4$$ an error message (a NIntegrate::slwcon) is produced when $$a$$ drops below some threshold value, but because of what we saw above it is doubtful that the code computes the correct value for the integral within the error bar even when $$a$$ is slightly larger and no error message is produced. I've tried increasing the WorkingPrecision but that gives an NIntegrate::precw error and doesn't solve the $$k=2$$ discrepancy.

Check on the correct solution

The numerical solution for $$k > 2$$ should be checked by computing also $$\mathrm{d} I(a;k)/\mathrm{d}a = \left( \prod_{i=1}^k \int_0^\infty \mathrm{d}x_i \right) \left[ \frac{1}{4a^2} \sum_{i=1}^k \left( x_i - x_{i+1} \right)^2 - \frac{1}{2} \sum_{i=1}^k (x_i + x_{i+1}) \right] \exp \left( - \frac{1}{4a} \sum_{i=1}^k \left( x_i - x_{i+1} \right)^2 - \frac{a}{2} \sum_{i=1}^k (x_i + x_{i+1}) \right),$$ computing its (numerical) integral from some reference value $$a_0$$ to a variable upper limit $$b$$ and showing that it agrees with $$I(b;k)-I(a_0;k)$$ (for some range of $$b$$ and some moderate values of $$k$$).

• I guess there is an analytic expression for integer $k \geq 3$, because it is just a Gaussian integral. Dec 10, 2020 at 2:53
• I think you're right -- let me think about this. I'll post the solution via that route if I find it. Anyhow it would be nice to have a numerical algorithm, because in slightly different scenarios that might still work while there might not be an analytic solution. We could check the numerical algorithm against the exact answer too. Dec 10, 2020 at 14:13
• Just some warn message. Dec 10, 2020 at 14:48
• Sorry @ΑλέξανδροςΖεγγ we were wrong. Like you suggested it is possible to write the argument of the exponential as a quadratic form, schematically x^T A x + J^T x with A a symmetric matrix and J,x vectors. We can orthogonally diagonalize A like A = P^T D P with P^T P = P P^T = 1. There seems to always be a single zero eigenvalue, but that's not necessarily a problem because the eigenvectors of P still form an orthonormal basis and we can change variables like y = Px. This doesn't change the measure factor because det(P) = 1. But I think it does change the integration domain in some weird way. Dec 11, 2020 at 19:13
• The integrals are not Gaussian because the integration domain is over [0,\infty)^k not (-\infty,\infty)^k. Dec 11, 2020 at 19:13

Here we use

k=4;
ListCorrelate[{1, -1}, Array[x, k], 1]


to make x[k+1] just equal to x[1]

{x[1] - x[2], x[2] - x[3], x[3] - x[4], -x[1] + x[4]}

And we use  Method -> "LocalAdaptive"

iint[a_, k_] :=
NIntegrate[
Exp[-1/(4 a)*(ListCorrelate[{1, -1}, Array[x, k], 1]^2 // Total) -
a/2 (ListCorrelate[{1, 1}, Array[x, k], 1] // Total)],
Sequence @@ Table[{x[i], 0, Infinity}, {i, k}] // Evaluate,
Table[iint[1/5, k], {k, 1, 5}]


{5., 2.61318, 2.00599, 1.79427, 1.73224}

iiint[a_, k_] :=
Integrate[
Exp[-1/(4 a)*(ListCorrelate[{1, -1}, Table[Subscript[x, i], {i, k}],
1]^2 // Total) -
a/2 (ListCorrelate[{1, 1}, Table[Subscript[x, i], {i, k}], 1] //
Total)],
Sequence @@ Table[{Subscript[x, i], 0, Infinity}, {i, k}] //
Evaluate]
ParallelTable[iiint[1/5, i] // N, {i, 1, 2}]


{5., 2.61319}

Appendix

We can also test the antiderivative.

integrate[a_, k_] :=
Integrate[
Exp[-1/(4 a)*(ListCorrelate[{1, -1}, Table[Subscript[x, i], {i, k}],
1]^2 // Total) -
a/2 (ListCorrelate[{1, 1}, Table[Subscript[x, i], {i, k}], 1] //
Total)], Sequence @@ Table[Subscript[x, i], {i, k}]]
integrate[a, 2]

• Sorry @cvgmt, the integrand you get with ListCorrelate E^(-(1/2) a (x[1] + 2 x[2] + 2 x[3] + x[4]) - ((x[1] - x[2])^2 + (x[2] - x[3])^2 + (x[3] - x[4])^2 + (-x[1] + x[4])^2)/( 4 a))  is not the integrand the OP provides E^(-(1/2) a (2 x[1] + 2 x[2] + 2 x[3] + 2 x[4]) - ((x[1] - x[2])^2 + (x[2] - x[3])^2 + (x[3] - x[4])^2 + (-x[1] + x[4])^2)/( 4 a))  (I didn't check which is right according to LATEx formula) Dec 10, 2020 at 5:25
• @Akku14 Thanks! Yes, ListCorrelate[{1, 1}, Array[x, k], 1] instead of  ListCorrelate[{1, 1}, Array[x, k]] Dec 10, 2020 at 6:17
• ParallelTable[iint[1/5, k], {k, 1, 10}] get {5., 2.44203, 1.83674, 1.72395, 1.69251, 1.70752, 1.77206, 1.86193, 2.05917, 2.4763}` Dec 10, 2020 at 14:53