Difficult numerical integration

Question

I have a complicated function which I'm trying to integrate numerically like this:

y2[t_, OptionsPattern[]] :=
  Module[
   {α = OptionValue[α], T = OptionValue[T], ν = OptionValue[ν], k = OptionValue[k], d = OptionValue[d], r = OptionValue[r], ϵ = OptionValue[ϵ], x0 = OptionValue[x0], n = OptionValue[n]},
   NIntegrate[ν/(2 π) E^(-ν k (t - s)) E^(-d q^α (r + q^2) (s - s′)) q Exp[-q^2 ϵ^2] Exp[-(T/k) (1 - E^(-ν k (s - s′)) - 1/2 (E^(-ν k s) - E^(-ν k s′))^2) q^2] Sin[q x0 (E^(-ν k s′) - E^(-ν k s))] (d q^α + ν T E^(-ν k (s - s′)) q^2/(r + q^2)), {s, 0, t}, {s′, 0, s}, {q, 0, ∞}, Method -> "QuasiMonteCarlo", MaxPoints -> n, AccuracyGoal -> ∞]];

Options[y2] = {α -> 0, T -> 1, ν -> 1, k -> 1, ϵ -> 0.001, x0 -> 1, d -> 1, r -> 1, n -> 10^5};

I've been almost forced to use "QuasiMonteCarlo" as integration method, because otherwise the integration takes ages. I'm interested in the result for values of the parameters around the ones I chose as default. I do eventually get a result, but with a lot of errors of the kind NIntegrate::maxp, that is

NIntegrate: The integral failed to converge after 100000 integrand evaluations.

The integrand function seems to behave reasonably well when I plot it. Could you help me understand what I'm doing wrong?

Answer 1

We can compute this integral numerically without messages by using Method -> "AdaptiveQuasiMonteCarlo" as follows (here s1=s′)

y2[t_, OptionsPattern[]] := 
  Module[{\[Alpha] = OptionValue[\[Alpha]], 
    T = OptionValue[T], \[Nu] = OptionValue[\[Nu]], 
    k = OptionValue[k], d = OptionValue[d], 
    r = OptionValue[r], \[Epsilon] = OptionValue[\[Epsilon]], 
    x0 = OptionValue[x0], n = OptionValue[n]}, 
   NIntegrate[\[Nu]/(2 \[Pi]) E^(-\[Nu] k (t - 
          s)) E^(-d q^\[Alpha] (r + q^2) (s - 
          s1)) q Exp[-q^2 \[Epsilon]^2] Exp[-(T/k) (1 - 
         E^(-\[Nu] k (s - s1)) - 
         1/2 (E^(-\[Nu] k s) - E^(-\[Nu] k s1))^2) q^2] Sin[
      q x0 (E^(-\[Nu] k s1) - 
         E^(-\[Nu] k s))] (d q^\[Alpha] + \[Nu] T E^(-\[Nu] k (s - 
             s1)) q^2/(r + q^2)), {s, 0, t}, {s1, 0, s}, {q, 
     0, \[Infinity]}, Method -> "AdaptiveQuasiMonteCarlo"]];

Options[y2] = {\[Alpha] -> 0, T -> 1, \[Nu] -> 1, 
   k -> 1, \[Epsilon] -> 10^-3, x0 -> 1, d -> 1, r -> 1};

To plot y2[t] we use

lst = Table[{t, y2[t]}, {t, 0, 2, .02}];

ListPlot[lst]

Note, that for $\alpha =0$ we can reduce dimension since integral over {q, 0, \[Infinity]} can be calculated exactly.