# 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
Commented Sep 10, 2015 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? Commented Sep 10, 2015 at 20:34

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 Commented Sep 11, 2015 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? Commented Feb 23, 2019 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. Commented Sep 11, 2015 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
Commented Sep 11, 2015 at 7:04
• Thanks for this, though substituting Integrate for NIntegrate isn't an option for me. See my edited post. Commented Sep 11, 2015 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 :( Commented Sep 10, 2015 at 20:53
• Release .. had to look that up, superseded in version 2 (1991!). (Evaluate works as well ) Commented Sep 11, 2015 at 2:39
• @george2079 oops… I am old and I forget :-) Then against does the job! Commented Sep 11, 2015 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. Commented Oct 3, 2015 at 14:35
• @AntonAntonov Thanks for the correction. I hope the update is more accurate. (And thanks for the upvote.) Commented Oct 3, 2015 at 23:45