1
$\begingroup$

I have a double variable integral, so I can define a function which is the integration of the other variable, and finally, I integrate the other variable. here is the code:

y[x_?NumericQ] := NIntegrate[v^2 Exp[-x v^2], {v, 0, \[Infinity]}]
NIntegrate[y[x], {x, 20, 1000000}]

It takes a long time to run; also it gives two errors: NIntegrate::inumr and General::stop.

So, I thought doing the integration together at one time is faster, and indeed it is faster. Here is the code:

NIntegrate[v^2 Exp[-x v^2], {v, 0, \[Infinity]}, {x, 20, 1000000}]

It does not have errors and it is ten times faster. However, the final result is a different value. So, I have three questions:

1- Which one should I trust?

2- Why do they have different results?

3- How can I make the NIntegrate faster? The original problem is kind of complicated and since I am going to plot the contour plot of two variables, it takes an extremely long time to give me the plot. For each point, the second code takes about 0.2 seconds, which is for a ContourPlot is extremely high. I try to decrease the number of the points with PlotPoints, but still, it takes a long time and also I lose the quality of the graph. I tried to decrease the PrecisionGoal to 4 inside the NIntegate, it helps but not too much.

Thank you so much for your kind helps.

Improve this question
$\endgroup$
6
  • 2
    $\begingroup$ Show us your original problem. This one can be done easily with {int = Integrate[ v^2 Exp[-x v^2], {v, 0, \[Infinity]}, {x, 20, 1000000}], int // N} to yield {((-1 + 100 Sqrt[5]) Sqrt[\[Pi]])/2000, 0.19728} . For the original problem, try to do at least integration over one variable analytically. $\endgroup$
    – Akku14
    Commented May 7, 2022 at 18:10
  • 2
    $\begingroup$ The problem is undersampling at the low end of the x range. (The v domain being infinite is transformed.) This performs better NIntegrate[v^2 Exp[-x v^2] Dt[x, u] /. x -> Exp[u], {v, 0, \[Infinity]}, {u, Log[20], Log[1000000]}]. If your actual problem is different, then the problem may lie in that difference. Usually, unhappy code is unhappy in its own way. $\endgroup$
    – Michael E2
    Commented May 7, 2022 at 18:41
  • 1
    $\begingroup$ A strategic division of the integration region: NIntegrate[v^2 Exp[-x v^2], {v, 0, \[Infinity]}, {x, 20, 1000000}, Exclusions -> {x v^2 == 30, x v^2 == 40}] -- again focused the particular example. $\endgroup$
    – Michael E2
    Commented May 8, 2022 at 22:37
  • $\begingroup$ Things to do to test a result (stability is taken to indicate reliability): Increase MinPrecision -> 3 (for integrals of dimension $d$, the time increases by a factor of $n^d$). Increase WorkingPrecision to 16, 24, 32. Increase PrecisionGoal 6, 8, 10 (but as the goal reaches the working precision, you'll get errors because the difference is how much round-off/ill-conditioning the integration can tolerate). $\endgroup$
    – Michael E2
    Commented May 8, 2022 at 22:44
  • $\begingroup$ @MichaelE2: Can you kindly explain how do you come to the trick Exclusions -> {x v^2 == 30, x v^2 == 40}? TIA. $\endgroup$
    – user64494
    Commented May 9, 2022 at 5:16

1 Answer 1

Reset to default
1
$\begingroup$

Indeed, in version 13 on Windows 10

y[x_?NumericQ] := NIntegrate[v^2 Exp[-x v^2], {v, 0, \[Infinity]}]
NIntegrate[y[x], {x, 20, 1000000}]

0.19728

, but

NIntegrate[v^2 Exp[-x v^2], {v, 0, Infinity}, {x, 20, 1000000}]

0.0638567

Both results are produced without any warnings and error communications. The latter result is not correct (and this is a bug) in view of

Integrate[v^2 Exp[-x v^2], {v, 0, Infinity}, {x, 20, 1000000}]

((-1 + 100 Sqrt[5]) Sqrt[\[Pi]])/2000

N[%]

0.19728

There is a workaround:

NIntegrate[v^2 Exp[-x v^2], {v, 0, Infinity}, {x, 20, 1000000}, 
Method -> {"GaussKronrodRule", "GaussPoints" -> 6}]

0.19728

Improve this answer
$\endgroup$
7
  • 1
    $\begingroup$ NIntegrate[v^2 Exp[-x v^2], {v, 0, Infinity}, {x, 20, 1000000}, Method -> {"GaussKronrodRule", "GaussPoints" -> 4}] results in 0.190341. Therefore, any symbolic integration is not used. $\endgroup$
    – user64494
    Commented May 8, 2022 at 4:50
  • 1
    $\begingroup$ On a fresh kernel NIntegrate[v^2 Exp[-x v^2], {v, 0, Infinity}, {x, 20, 1000000}, Method -> {"GaussKronrodRule", "GaussPoints" -> 6}] // AbsoluteTiming produces {1.10956, 0.19728}. $\endgroup$
    – user64494
    Commented May 8, 2022 at 7:41
  • 1
    $\begingroup$ I would like to use this opportunity to report a bug in Integrate[v^2 Exp[-x *(v^2 + v)], {x, 20, 1000000}, {v, 0, Infinity}] which produces 1/4 E^5 Sqrt[\[Pi]/5] - (E^250000 Sqrt[\[Pi]])/2000 + 1/6 (3 \[Pi] (Erfi[500] - Erfi[Sqrt[5]]) + 15 HypergeometricPFQ[{1, 1}, {3/2, 2}, 5] - 750000 HypergeometricPFQ[{1, 1}, {3/2, 2}, 250000] + 5 HypergeometricPFQ[{1, 1}, {2, 5/2}, 5] - 250000 HypergeometricPFQ[{1, 1}, {2, 5/2}, 250000] + Log[8] - 15 Log[10]). Its numerical value equals 7.62467*10^108560 though it should be less than 0.19728. $\endgroup$
    – user64494
    Commented May 8, 2022 at 9:11
  • $\begingroup$ Thank you so much for your smart trick. It works, however, in the original problem I have to increase GausssPoints which makes it even slower than the first integral. So, I posted the original problem here (mathematica.stackexchange.com/questions/268049/…), I am so appreciative if you have any comments on it. $\endgroup$
    – Mehrdad
    Commented May 10, 2022 at 1:16
  • $\begingroup$ "[...] result is not correct (and this is a bug) in view of [...]" -- I am not sure is it a bug or "a feature." The multi-dimensional rule, "MultidimensionalRule", is somewhat sparse of points and can produce convergence quicker. So, providing more sampling points -- via "CartisianRule", or MinRecursion, etc. -- might produce more faithful results. And that is a fairly well-known phenomena. (Also documented.) $\endgroup$
    – Anton Antonov
    Commented Oct 7, 2022 at 16:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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