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

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


## 2 Answers

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

• One might try MultidimensionalRule on the 4D integral. I thought it was the default but realized NIntegrate probably was not using it just when I got called away. I’ll try it out later when I get back to a computer. But feel free to steal a march on me. :) – Michael E2 Mar 23 at 18:45

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

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

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

• Um, that output suggests that the integral has not been calculated, it has just returned the input unevaluated. – KraZug Mar 22 at 20:31
• There's a missing comma. I've fixed it in the OP. – Michael E2 Mar 22 at 20:34