Any possibility to speed up this code?

Question

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}];

Answer 1

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