# How to improve the speed of this nested NIntegrate?

I have an integral that needs to be evaluated using NIntegrate: $$\int_0^{10}\mathrm{d}x_1\int_0^{10}\mathrm{d}x_2\left(1+(x_2-x_1)^2\right)^{-1/6}\mathrm{e}^{-\int_0^{x_1}\mathrm{d}x'\int_0^{x_1}\mathrm{d}x''f(x'-x'')-\int_0^{x_2}\mathrm{d}x'\int_0^{x_2}\mathrm{d}x''f(x'-x'')}\times\int_0^{x_1}\mathrm{d}x'\int_0^{x_1}\mathrm{d}x''f(x_1-x')f(x_2-x''),$$ where $$f(x)=\frac{1}{(1+x^2)^{1/3}}.$$

If I follow the instruction given in this link, the computational time it takes is too much. Is there any way to improve the computational time?

Here is the code:

f[x_]:=1/(1+x^2)^(1/3);

i1[x1_?NumericQ, x2_?NumericQ] :=  i1[x1, x2] =
NIntegrate[f[z1 - z2], {z1, 0, x1}, {z2, 0, x1}];

i2[x1_?NumericQ, x2_?NumericQ] :=
i2[x1, x2] =
NIntegrate[f[z1 - z2], {z1, 0, x2}, {z2, 0, x2}];

NIntegrate[((1+(x2-x1)^2)^(-1/6)Exp[(i1[x1, x2] + i2[x1, x2])]f[x1-z1]*
f[x2-z2], {x1, 0, 10}, {x2, 0, 10},{z1,0,x1},{z2,0,x1}];

• Am I the only one who hates the notation with the differential operator in front of the integrand? This looks like four integrals multiplied together. Please post the code that takes too much time. – george2079 Aug 11 '16 at 17:33
• @george2079 I hate the non-italic differential $d$, which seems to me to be a recent invention of LaTeX-addled minds (compare with plain TeX, TeXBook, p. 169, as well as the few centuries of printing pre-2000). The "wrong" order is standard in physics, I believe. – Michael E2 Aug 11 '16 at 17:57
• A few issues, what is f[0,..] ? You have unbalanced parenthesis in last expression. If you are trying to use memoization the second := should just be = – george2079 Aug 11 '16 at 18:00
• Sorry? There should be no f[0,..]. It is just f[x]. I have fixed that. – titanium Aug 11 '16 at 18:03
• @MichaelE2 But according to tug.org/TUGboat/Articles/tb18-1/tb54becc.pdf, you have to set it in Roman. But I agree: I hate that rule. – QuantumDot Aug 11 '16 at 18:19

there is one more analytic integral we can pull out:

(using @BobHanlon's i[x]):

g[x1_, x2_] =
Assuming[{x1 > 0, x2 > 0},
Integrate[f[x1 - z1]*f[x2 - z2], {z1, 0, x1}, {z2, 0, x1}]];
NIntegrate[(1 + (x2 - x1)^2)^(-1/6) Exp[-(i[x1] + i[x2])]
g[x1, x2] , {x1, 0, 10}, {x2, 0, 10}] // AbsoluteTiming


{1.30975, 0.355381}

BTW all the other x1/x2 symmetry in the formulation makes me wonder if one of those x1 integration limits should be an x2 (of course I have no idea what the actual calculation is)

f[x_] = (1 + x^2)^(-1/3);


Use Integrate rather than NIntegrate for i1 and i2.

Clear[i, i1, i2]

i[x_] = Assuming[{x > 0},
Integrate[f[z1 - z2], {z1, 0, x}, {z2, 0, x}]]

(*  1/2 (3 - 3 (1 + x^2)^(2/3) + 4 x^2 Hypergeometric2F1[1/3, 1/2, 3/2, -x^2])  *)


NIntegrate[(1 + (x2 - x1)^2)^(-1/6) Exp[-(i[x1] + i[x2])] f[x1 - z1]*