# error with NIntegrate and RegionPlot

Bug introduced in 10.2 or earlier and persisting through 11.0.1 or later

Here's a simplified example of what I'm trying to do:

RegionPlot[
NIntegrate[PDF[NormalDistribution[0, 1], a], {a, 0, y}] >= 0.2,
{x, -1, 1}, {y, 0.1, 0.7}]


(yes, it doesn't depend on x)

This returns an inscrutable error:

Throw::nocatch: "Uncaught \!$$Throw[\(-HolonomicDifferentialRootReduceDumpy[NIntegrateLevinRuleDumpx]$$
+ \*SuperscriptBox[\"HolonomicDifferentialRootReduceDumpy\", \"\[Prime]\",
MultilineFunction->None][NIntegrateLevinRuleDumpx],
NIntegrateLevinRuleDumpFastLookupHolonomicDifferentialEquation]\) returned to top level."


I'm using Mathematica 10.2. Any ideas what's wrong?

EDIT: I'm not looking for algebraic simplifications or substituting Integrate for NIntegrate. The above code is just an example to reproduce the error. In the code I actually want to run, NIntegrate is the only option.

• Have you tried to figure out where the error come from? To izolate the problem try divide and conquer approach first... – mmal Sep 10 '15 at 20:02
• I can reproduce that error with 10.1, even defining a pure numeric function to pass to RegionPlot. The NIntegrate by itself works fine. B-G? – george2079 Sep 10 '15 at 20:34

## 4 Answers

This is a bug in RegionPlot. For a possible workaround, try the following undocumented option

RegionPlot[NIntegrate[PDF[NormalDistribution[0, 1], a], {a, 0, y}] >= 0.2,
{x, -1, 1}, {y, 0.1, 0.7}, "NumericalFunction" -> False] • I cannot find documentation of the option "NumericalFunction" for RegionPlot nor Graphics – Bob Hanlon Sep 11 '15 at 16:22
• It seems that 11.3 introduces additional error NIntegrate::nlim: a = y is not a valid limit of integration. and also the bug is still not solved, any idea that why the new error shows up? – luyuwuli Feb 23 '19 at 16:21

EDIT : Changed for your edited question

\$Version

(* "10.2.0 for Mac OS X x86 (64-bit) (July 7, 2015)" *)


Define a helper function that is defined only for numeric arguments

f[y_?NumericQ] :=
NIntegrate[
PDF[NormalDistribution[0, 1], a],
{a, 0, y}];

rgn = ImplicitRegion[
f[y] >= 0.2 && -1 <= x <= 1 && 0.1 <= y <= 0.7,
{x, y}];


However, this is very sloo...oow

RegionPlot[rgn, PlotRange -> {{-1, 1}, {0.1, 0.7}}] //
AbsoluteTiming // Column ContourPlot is much, much faster

f[0.62] >= 0.2

(* True *)

ContourPlot[f[y],
{x, -1, 1}, {y, 0.1, 0.7},
Contours -> {0.2},
Epilog -> Text["f[y] \[GreaterEqual] 0.2", {0, 0.62}]] //
AbsoluteTiming // Column The largest that your integral can be is for y = 0.7

dist = NormalDistribution[0, 1];

Integrate[PDF[dist, a], {a, 0, 0.7}]

(* 0.258036 *)


This is equivalent to

CDF[dist, 0.7] - CDF[dist, 0]

(* 0.258036 *)


Even if you were to integrate from -Infinity, the largest that the integral could be is

Integrate[PDF[dist, a], {a, -Infinity, 0.7}]

(* 0.758036 *)


or equivalently,

CDF[dist, 0.7] - CDF[dist, -Infinity]

(* 0.758036 *)


or more simply

CDF[dist, 0.7]

(* 0.758036 *)


Consequently, since you are looking for the region for which the integral is greater than or equal to 0.95, your region is empty. Note the use of Integrate rather than NIntegrate

RegionPlot[
Integrate[PDF[NormalDistribution[0, 1], a], {a, 0, y}] >= 0.95, {x, -1,
1}, {y, 0.1, 0.7}] If you reverse the inequality then

RegionPlot[
Integrate[PDF[dist, a], {a, 0, y}] < 0.95, {x, -1, 1}, {y, 0.1, 0.7}] Or the same result with

RegionPlot[CDF[dist, y] - CDF[dist, 0] < 0.95, {x, -1, 1}, {y, 0.1, 0.7}]

• be that as it may, we shouldn't expect an "inscrutable" error message. – george2079 Sep 11 '15 at 2:30
• Yes, it has nothing to do with NormalDistribution, CDF or PDF. Try RegionPlot[NIntegrate[Sin[a], {a, 0, y}] >= 0.95, {x, -1, 1}, {y, 0.1, 0.7}] vs RegionPlot[NIntegrate[a, {a, 0, y}] >= 0.95, {x, -1, 1}, {y, 0.1, 0.7}]. – mmal Sep 11 '15 at 7:04
• Thanks for this, though substituting Integrate for NIntegrate isn't an option for me. See my edited post. – mlubin Sep 11 '15 at 14:45

Try this small variation over your original request

   RegionPlot[
Integrate[PDF[NormalDistribution[0, 1], a], {a, 0, x}] > 0.3 //
Evaluate, {x, -10, 10}, {y, -1, 1}] • doesn't work on v9 :( – Dr. belisarius Sep 10 '15 at 20:53
• Release .. had to look that up, superseded in version 2 (1991!). (Evaluate works as well ) – george2079 Sep 11 '15 at 2:39
• @george2079 oops… I am old and I forget :-) Then against does the job! – chris Sep 11 '15 at 5:34

While the behavior has already been confirmed a bug -- an uncaught Throw from an internal function must always be one, right? -- here are a couple more workarounds.

How I analyze the problem of finding a workaround: From the context NIntegrateLevinRuleDump in the error message, one might infer that NIntegrate is trying to determine whether to use (or even using) "LevinRule" for the integration method. Since it's for oscillatory integrands, it doesn't seem that appropriate. Probably it is occurring during the "SymbolicProcessing" phase of NIntegrate. And if not, we could pick an integration rule manually.

RegionPlot[
NIntegrate[PDF[NormalDistribution[0, 1], a], {a, 0, y},
Method -> {"GlobalAdaptive", "SymbolicProcessing" -> 0}] >=
0.2, {x, -1, 1}, {y, 0.1, 0.7}]; // AbsoluteTiming RegionPlot[
NIntegrate[PDF[NormalDistribution[0, 1], a], {a, 0, y},
Method -> {"GlobalAdaptive", Method -> "GaussKronrodRule"}] >=
0.2, {x, -1, 1}, {y, 0.1, 0.7}]; // AbsoluteTiming Turning off symbolic processing is fastest, but you lose automatic method selection for certain types of integral such as when the integrand is oscillatory. By comparison, ilian's mystery option "NumericalFunction" -> False takes about as long as picking the "GaussKronrodRule" explicitly. For some reason, this method does not produce a spurious error message (which comes from RegionPlot evaluating the argument symbolically).

RegionPlot[
NIntegrate[PDF[NormalDistribution[0, 1], a], {a, 0, y}] >= 0.2,
{x, -1, 1}, {y, 0.1, 0.7},
"NumericalFunction" -> False]; // AbsoluteTiming All produce plots that look the same: • This is a good approach to analyze the problem (I upvoted). I would like to point out that singularity handling for one dimensional integrals is not switched off by "SymbolicProcessing"->0. That is done with the option setting "SingularityHandler" -> None in the integration strategy methods. – Anton Antonov Oct 3 '15 at 14:35
• @AntonAntonov Thanks for the correction. I hope the update is more accurate. (And thanks for the upvote.) – Michael E2 Oct 3 '15 at 23:45