I have an integral that needs to be evaluated using NIntegrate:

∫100dx1∫100dx2(1+(x2−x1)2)−1/6e−∫x10dx′∫x10dx″f(x′−x″)−∫x20dx′∫x20dx″f(x′−x″)×∫x10dx′∫x10dx″f(x1−x′)f(x2−x″),\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}}.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}];

=================

3

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 at 17:33

@george2079 I hate the non-italic differential dd, 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 at 17:57

1

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 at 18:00

Sorry? There should be no f[0,..]. It is just f[x]. I have fixed that.

– titanium

Aug 11 at 18:03

1

@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 at 18:19

=================

2 Answers

2

=================

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]) *)

You are missing a minus sign in the exponential

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

f[x2 – z2], {x1, 0, 10}, {x2, 0, 10}, {z1, 0, x1}, {z2, 0,

x1}] // AbsoluteTiming

(* {14.2359, 0.355381} *)