# Any possibility to speed up this code?

I am running a Mathematica code to evaluate an integral. The probelm is this code is running too slow. Is there any possibility to speed this up? Here is the code:

xi[R_, z_] := 9*(0.00001/(R^2 + z^2))^(1/6);
i1[a_] := N[2.37783 a^(2 - 1/3)];
Dz[R_, z_] := (R^2 + z^2)^(1/3);
i3[R_?NumericQ, z1_?NumericQ, z2_?NumericQ] :=
i3[R, z1, z2] =
NIntegrate[xi[R, x1 - x2], {x1, 0, z1}, {x2, 0, z2}];

i4[R_, z1_?NumericQ, z2_?NumericQ] :=
i4[R, z1, z2] =
1.98153 (z1^(2/3))*NIntegrate[xi[R, z2 - x8], {x8, 0, z1}];

i5[R_, z1_?NumericQ, z2_?NumericQ] :=
i5[R, z1, z2] =
1.98153 (z2^(2/3))*NIntegrate[xi[R, z1 - x7], {x7, 0, z2}];

i6[R_, z1_?NumericQ, z2_?NumericQ] :=
i6[R, z1, z2] =
1.98153 (z1^(2/3)) *1.98153 (z2^(2/3)) +
NIntegrate[xi[R, z1 - x7]*xi[R, z2 - x8], {x7,
0, z1}, {x8, 0, z2}];

NIntegrate[Exp[-4*(z1 + z2)]* Exp[(i1[z1] + i1[z2])] (1/Sqrt[Dz[0, z1 - z2]]
Exp[i3[0, z1, z2]]*(16 + xi[0, z1 - z2] + (i4[0, z1, z2] + i5[0, z1, z2] +
i6[0, z1, z2])) -
1/Sqrt[Dz[0.01, z1 - z2]]Exp[(i3[0.01, z1, z2])]*(16 + xi[0.01,
z1 - z2] + (i4[0.01, z1, z2] + i5[0.01, z1, z2] +
i6[0.01, z1, z2]))), {z1, 0, 0.25}, {z2, 0, 0.25},
Method -> {"GlobalAdaptive", "SymbolicProcessing" -> 0}];
• i3 can be solved symbolically with Integrate, and it may be possible to Integrate i4, i5, and i6 as well. Nov 24, 2016 at 2:43

To make my comment above more concrete, the integrals defining i3, i4, i5, and i6 can be evaluated symbolically, saving a significant amount of time. Define

f = 9*(10^-5/(R^2 + z^2))^(1/6);

Then, the i3 integral is

Integrate[f /. z -> x1 - x2, {x1, 0, z1}, {x2, 0, z2},
Assumptions -> z1 > 0 && z2 > 0 && R > 0]

(* (1/(5 10^(5/6) R^(1/3))) 9 (3 R^(1/3) (R^(5/3) - (R^2 + z1^2)^(5/6) +
(R^2 + (z1 - z2)^2)^(5/6) - (R^2 + z2^2)^(5/6)) +
5 z1^2 Hypergeometric2F1[1/6, 1/2, 3/2, -(z1^2/R^2)] -
5 (z1 - z2)^2 Hypergeometric2F1[1/6, 1/2, 3/2, -((z1 - z2)^2/R^2)] +
5 z2^2 Hypergeometric2F1[1/6, 1/2, 3/2, -(z2^2/R^2)]) *)

The i4 integral is

1.98153 (z1^(2/3))*Integrate[f /. z -> z2 - x8, {x8, 0, z1},
Assumptions -> z1 > 0 && z2 > 0 && R > 0]
(* (2.61764 z1^(2/3) ((z1 - z2) Hypergeometric2F1[1/6, 1/2, 3/2, -((z1 - z2)^2/R^2)] +
z2 Hypergeometric2F1[1/6, 1/2, 3/2, -(z2^2/R^2)]))/R^(1/3) *)

The i5 integral is

1.98153 (z2^(2/3))*Integrate[f /. z -> z1 - x7, {x7, 0, z2},
Assumptions -> z1 > 0 && z2 > 0 && R > 0]
(* (2.61764 z2^(2/3) (z1 Hypergeometric2F1[1/6, 1/2, 3/2, -(z1^2/R^2)] + (-z1 + z2)
Hypergeometric2F1[1/6, 1/2, 3/2, -((z1 - z2)^2/R^2)]))/R^(1/3) *)

Finally, the i6 integral is

1.98153 (z1^(2/3))*1.98153 (z2^(2/3)) +
Integrate[(f /. z -> z1 - x7) (f /. z -> z2 - x8), {x7, 0, z1}, {x8, 0, z2},
Assumptions -> z1 > 0 && z2 > 0 && R > 0]
(* 3.92646 z1^(2/3) z2^(2/3) + (81 z1 z2 Hypergeometric2F1[1/6, 1/2, 3/2, -(z1^2/R^2)]
Hypergeometric2F1[1/6, 1/2, 3/2, -(z2^2/R^2)])/(10 10^(2/3) R^(2/3)) *)
• Sorry for disturbing you @bbgodfrey, but this is my only way to contact you. Could you please see my post. Thank you in advance. Jul 24, 2019 at 10:52