3
$\begingroup$

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?

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

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]

Figure 1

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

$\endgroup$
2
  • $\begingroup$ 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? $\endgroup$ – Davide Venturelli Feb 27 at 9:47
  • 1
    $\begingroup$ @DavideVenturelli Yes it is common problem. Some function visualization algorithms can generate errors when compute intermediate values with functions like y2[t]. $\endgroup$ – Alex Trounev Feb 27 at 11:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.