What is the fastest way of doing one integral in Mathematica?

Question

I want to compute the following numerical integral in Mathematica

$\int_0^L dy \int_0^y d \bar{y} f(y,\bar{y}),$

where $f(y,\bar{y})$ is a very complicated function.

I show my code below where the important line is the last one where I define the integral that I want to make. My problem is that when I try to compute for example

t1[4, L, co, QS]

it takes much more time than I should expect. So my question is: What is the fastest numerical integral method for dealing with this kind of integrals?

gaussexp[pol_, arg_] := Coefficient[pol, arg, 1]^2/(-4 Coefficient[pol, arg, 2]) + Coefficient[pol, arg, 0]
gaussfac[pol1_, pol2_, arg_] := -Pi / Coefficient[pol2, arg, 2] (Coefficient[pol1, arg, 0] + Coefficient[pol1, arg, 1] Coefficient[pol2, arg, 1]/(-2 Coefficient[pol2, arg, 2]) + Coefficient[pol1, arg, 2] (1 + Coefficient[pol2, arg, 1]^2/(-2 Coefficient[pol2, arg, 2]))/(-2 Coefficient[pol2, arg, 2]))
G1[y_, yb_, xp_, yp_] := Exp[-I k (y - yb) - Qs^2/4 (y - yb)^2 (xp - yp)/lp]
G2[x_, y_, yb_] := (-I Qs^2/(2 Pi a1)) Exp[I Qs^2/(2 a1) (x - y)^2 + I Qs^2/8 a1 1/180 (15 (y - yb)^2 + 15 (y - yb) (x - y) + 4 (x - y)^2) - Qs^2/4 ((y - yb)^2 + (y - yb) (x - y) + 1/3 (x - y)^2)]

der = Simplify[D[G1[u, yb, lp, xb], yb] D[G2[u, y, yb], y]] Exp[-Qs^2/4 (y - yb)^2 a3];
exp2 = Exponent[der, E];
fac2 = der/Exp[exp2];
solexp = gaussexp[exp2, u];
solfac = gaussfac[fac2, exp2, u];
expt1a = gaussexp[solexp + I q1 y - I q2 yb, y];
fact1a = gaussfac[solfac, solexp + I q1 y - I q2 yb, y];

η = 2;
L = 30/Cosh[η];
QS = 1;
μt = 2/10;
co = Sqrt[2] Exp[-η] QS L;

exp = Simplify[
   expt1a /. {yb -> 0, lp -> l, a1 -> e (xb - x)/l, a3 -> x/l, k -> ke, q1 -> qe}];
fac = Simplify[
   fact1a /. {yb -> 0, lp -> l, a1 -> e (xb - x)/l, a3 -> x/l, k -> kf, q1 -> qf}];
exp2 = exp /. {qe^2 -> q^2, qe ke -> q k Cos[t], ke^2 -> k^2};
fac2 = fac /. {qf^2 -> q^2, kf qf -> q k Cos[t], kf^2 -> k^2};

t1[k1_, l1_, e1_, Qs1_] := 
 2 Re[NIntegrate[
    k1^2 e1^2/(Qs1^4 l1^2) fac2 Exp[exp2] /. {k -> k1, l -> l1, e -> e1, Qs -> Qs1}, 
       {q, 0, Infinity}, {t, 0, 2 Pi}, {x, 0, l1}, {xb, 0, x}]
   ]

Answer 1

High dimensional integrals are often hard to evaluate accurately. A dimension greater than one is often hard, and greater than two, harder. A quick method with poor convergence is a Monte Carlo method:

integrand = k1^2 e1^2/(Qs1^4 l1^2) fac2 Exp[exp2];

ClearAll[t1];
t1[k1_, l1_, e1_, Qs1_] := With[{mp = mp},
   2 Re[NIntegrate[
      k1^2 e1^2/(Qs1^4 l1^2) fac2 Exp[exp2] /.
        {k -> k1, l -> l1, e -> e1, Qs -> Qs1} // Evaluate,
      {q, 0, Infinity}, {t, 0, 2 Pi}, {x, 0, l1}, {xb, 0, x},
      Method -> {"QuasiMonteCarlo", "MaxPoints" -> mp}, 
      PrecisionGoal -> 3]]
   ];

PrintTemporary@Dynamic@{Clock[Infinity]};
mp = 10^5;
t1[4, L, co, QS] // AbsoluteTiming
mp = 10^6;
t1[4, L, co, QS] // AbsoluteTiming

NIntegrate::maxp: The integral failed to converge after 100000 integrand evaluations. NIntegrate obtained 40.1598 +503.395 I and 22.13462812471233` for the integral and error estimates.

(*  {1.09656, 80.3196}  *)

NIntegrate::maxp: The integral failed to converge after 1000000 integrand evaluations. NIntegrate obtained 42.1303 +500.87 I and 8.240494084903718` for the integral and error estimates.

(*  {10.7654, 84.2605}  *)

Generally, the time is proportional to the number of points mp and the error is inversely proportion to the square root Sqrt[mp].

Update:

The default rule applied here is the Levin rule. The Multidimensional rule turns out to be more effective. This is easier than the original improvement below, though roughly equivalent in terms of time and accuracy.

t1[k1_, l1_, e1_, Qs1_] := 2 Re[NIntegrate[
    k1^2 e1^2/(Qs1^4 l1^2) fac2 Exp[exp2] /. {k -> k1, l -> l1, 
       e -> e1, Qs -> Qs1} // Evaluate,
    {q, 0, Infinity}, {t, 0, 2 Pi}, {x, 0, l1}, {xb, 0, x},
    Method -> "MultidimensionalRule", PrecisionGoal -> 3]]

PrintTemporary@Dynamic@{Clock[Infinity]};
t1[4, L, co, QS] // AbsoluteTiming
(*  {16.3445, 81.5715}  *)

Original workaround:

Greater accuracy requires analysis; that is, it takes work to figure if anything is possible. The integrand is periodic in t, which means it is amenable to the "Trapezoidal" method. The three-dimensional subintegral is still not easy.

ClearAll[t1];
t1[k0_?NumericQ, l0_, e0_, Qs0_] := Module[{t2},
   t2[k1_, l1_, e1_, Qs1_, t1_Real] :=
    NIntegrate[
     k1^2 e1^2/(Qs1^4 l1^2) fac2 Exp[exp2] /.
        {k -> k1, 
         l -> l1, e -> e1, Qs -> Qs1, t -> tt} // Re // Evaluate,
     {q, 0, Infinity}, {x, 0, l1}, {xb, 0, x},
     Method -> "MultidimensionalRule", PrecisionGoal -> 4, 
     AccuracyGoal -> 3];

   2 Re[NIntegrate[
      t2[k0, l0, e0, Qs0, t],
      {t, 0, 2 Pi}, Method -> {"Trapezoidal"}, PrecisionGoal -> 3]]
   ];

PrintTemporary@Dynamic@{Clock[Infinity]};
t1[4, L, co, QS] // AbsoluteTiming

(*  {32.1898, 81.5715}  *)

The integrand is also symmetric in t, so one can save half the evaluations and half the time:

ClearAll[t1];
t1[k0_?NumericQ, l0_, e0_, Qs0_] := Module[{t2, t3},
   t2[k1_, l1_, e1_, Qs1_, t1_Real] := t3[k1, l1, e1, Qs1, t1];
   t3[k1_, l1_, e1_, Qs1_, tt_] :=
    t3[k1, l1, e1, Qs1, t0_ /; t0 == tt] =           (* save tt value *)
     t3[k1, l1, e1, Qs1, t0_ /; t0 == 2 Pi - tt] =   (* save 2 Pi - tt value *)
      NIntegrate[
       k1^2 e1^2/(Qs1^4 l1^2) fac2 Exp[exp2] /.
          {k -> k1, 
           l -> l1, e -> e1, Qs -> Qs1, t -> tt} // Re // Evaluate,
       {q, 0, Infinity}, {x, 0, l1}, {xb, 0, x},
       Method -> "MultidimensionalRule", PrecisionGoal -> 4, 
       AccuracyGoal -> 3];

   times = {}; err = 0;
   2 Re[NIntegrate[
      t2[k0, l0, e0, Qs0, t],
      {t, 0, 2 Pi}, Method -> {"Trapezoidal"}, PrecisionGoal -> 3]]
   ];

PrintTemporary@Dynamic@{Clock[Infinity]};
t1[4, L, co, QS] // AbsoluteTiming

(*  {16.4071, 81.5715}  *)

Answer 2

On my laptop the evaluation of t1[4, L, co QS] is almost instantaneous:

In[24]:= t1[4, L, co QS] // AbsoluteTiming Out[24]= {6.*10^-6, t1[4, 30 Sech[2], (30 Sqrt[2] Sech[2])/E^2]}

Also, you should try harder to show a minimal example of your problem.