# Difficult numerical integration

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?

• What is the range for y2[t] plot? – Alex Trounev Feb 25 at 22:24
• Good point, usually t in [0,10] is enough – Davide Venturelli Feb 26 at 17:12

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,

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.

• Thanks. I realized that I still get errors if I plot y2[t] using Plot[], while I don't if I produce a list first and then use ListPlot, as you suggested. Is it common? – Davide Venturelli Feb 27 at 9:47
• @DavideVenturelli Yes it is common problem. Some function visualization algorithms can generate errors when compute intermediate values with functions like y2[t]. – Alex Trounev Feb 27 at 11:01