14
$\begingroup$

FINAL EDIT: The "curious behavior" is, as explained by user21, due to the fact that ImplicitRegion decides the RegionDimension according to the form of the input. If you are interested in a workaround for the issue, see also the post of M. Stern. Mathematica support agrees that the questions surrounding this behavior have been answered by the contributors to this post. Thanks to all!

Can anyone explain the following curious behavior? I first integrate over a disk, and get the expected answer:

NIntegrate[1, {x, y} ∈  Disk[]]

3.14159

Now I construct a numerical region and recalculate:

Needs["NDSolve`FEM`"]
inDisk[x_?NumericQ, y_?NumericQ] := x^2 + y^2 <= 1
nrDisk = 
  ToNumericalRegion[ImplicitRegion[inDisk[x, y] == True, {{x, -1, 1}, {y, -1, 1}}]];
NIntegrate[1, {x, y} ∈  nrDisk]

6.28054

The answer is approximately twice what it should be. Any ideas?

Bonus points: If I define the disk another way, why does NIntegrate just spit it back?

r2Disk[x_?NumericQ, y_?NumericQ] := x^2 + y^2
nrDisk = 
  ToNumericalRegion[ImplicitRegion[r2Disk[x, y] <= 1, {{x, -1, 1}, {y, -1, 1}}]];
NIntegrate[1, {x, y} ∈  nrDisk]

NIntegrate[1, {x, y} ∈ nrDisk]

Improve this question
$\endgroup$
8
  • 1
    $\begingroup$ The problem seems to be due to the fact that ImplicitRegion[] is unable to handle pure Boolean-valued functions like your inDisk[]. $\endgroup$
    – J. M.'s slightly less busy ♦
    Nov 3, 2017 at 14:30
  • 1
    $\begingroup$ I get the expected result ( 10.1 ) $\endgroup$
    – george2079
    Nov 3, 2017 at 14:36
  • $\begingroup$ @J.M. Thanks, this and your posted answer seem like key insights. Yet ConstantRegionQ[ImplicitRegion[x^2 + y^2 <= 1, {{x, -1, 1}, {y, -1, 1}}]] returns True. What is the technical difference between the expression x^2 + y^2 <= 1 and a pure Boolean-valued function? $\endgroup$
    – Will.Mo
    Nov 3, 2017 at 14:37
  • $\begingroup$ @george2079 I am using Mathematica 11.1.1.0. $\endgroup$
    – Will.Mo
    Nov 3, 2017 at 14:38
  • 1
    $\begingroup$ Nothing, supposedly. I'd report this to support were I you. (Also, it's m_goldberg who posted an answer, not me. :)) $\endgroup$
    – J. M.'s slightly less busy ♦
    Nov 3, 2017 at 14:42

3 Answers 3

Reset to default
7
$\begingroup$

The circumference of a unit disk is 2π.

When you hide the internal of the function and the == is the only thing ImplicitRegion does see it tries to construct a region with RegionDimension 1:

Needs["NDSolve`FEM`"]
inDisk[x_?NumericQ, y_?NumericQ] := x^2 + y^2 <= 1
nrDisk = ToNumericalRegion[
   ImplicitRegion[inDisk[x, y] == True, {{x, -1, 1}, {y, -1, 1}}]];
2 Pi - NIntegrate[1, {x, y} \[Element] nrDisk]

0.00264963

What will work better in this case is to use:

mesh = ToElementMesh[
   ImplicitRegion[
    If[inDisk[x, y], 1, -1] > 0, {{x, -1, 1}, {y, -1, 1}}]];
Pi - NIntegrate[1, {x, y} \[Element] mesh]
7.120859435438831`*^-7

Note that since there is no symbolic region associated with the mesh it's not really necessary to put the mesh in a NumericalRegion.

You can do this:

nr = ToNumericalRegion[
   ImplicitRegion[
    If[inDisk[x, y], 1, -1] > 0, {{x, -1, 1}, {y, -1, 1}}]];
ToElementMesh[nr]
Pi - NIntegrate[1, {x, y} \[Element] nr]
3.9561815157185265`*^-7

I'd need to look at why one needs to help NIntegrate with the numerical region in this case by calling ToElementMesh prior to calling NIntegrate.

There is a suggestion that ImplicitRegion should be able handle boolean region specifications better but that's out of my control so I can not say when and if this will come.

Improve this answer
$\endgroup$
3
  • $\begingroup$ +1 for pointing out that it's the circumference. When you try different radii r it becomes obvious that the result is proportional to r instead of r^2. $\endgroup$
    – M. Stern
    Nov 14, 2017 at 1:21
  • $\begingroup$ Thanks @user21, I think this basically resolves the question. I found a similar workaround, but defined a separate function r2Disk that returns a smooth numerical value, which can be compared with 0. (I assumed that ImplicitRegion would appreciate a smooth function, but apparently it is just as happy with your step function.) By the way, is there any way you know of to "inform" ImplicitRegion that it should try to construct a 2D region (other than giving it an inequality as input)? $\endgroup$
    – Will.Mo
    Nov 14, 2017 at 8:56
  • 1
    $\begingroup$ @Will.Mo, no, there is no way to "inform" ImplicitRegion about that. $\endgroup$
    – user21
    Nov 14, 2017 at 22:55
14
$\begingroup$

According to the documentation of ToNumericalRegion:

The region r should be a constant region for which ConstantRegionQ gives True.

However, 

Needs["NDSolve`FEM`"]
inDisk[x_?NumericQ, y_?NumericQ] := x^2 + y^2 <= 1
r = ToNumericalRegion[ImplicitRegion[inDisk[x, y] == True, {{x, -1, 1}, {y, -1, 1}}]];
ConstantRegionQ @ r

False

so your region definition doesn't meet the requirements of ToNumericalRegion.

The following works.

r = ImplicitRegion[x^2 + y^2 <= 1, {{x, -1, 1}, {y, -1, 1}}];
ConstantRegionQ @ r

True

Area @ r

π

NIntegrate[1, {x, y} ∈ r]

3.14159

Improve this answer
$\endgroup$
10
  • $\begingroup$ one could add that also DiscretizeRegion gives a workaround: nrDisk = DiscretizeRegion[ ImplicitRegion[r2Disk[x, y] <= 1, {{x, -1, 1}, {y, -1, 1}}], MaxCellMeasure -> 10, AccuracyGoal -> 7]; $\endgroup$
    – M. Stern
    Nov 3, 2017 at 14:49
  • $\begingroup$ Thanks very much! However, the explicit boolean expression within ImplicitRegion is, unfortunately, what I'm trying to avoid. I'm interested in using a more complicated boolean function that cannot be entered this way. $\endgroup$
    – Will.Mo
    Nov 3, 2017 at 14:49
  • $\begingroup$ @M.Stern You cannot discretize a region that does not satisfy ConatantRegionQ. $\endgroup$
    – m_goldberg
    Nov 3, 2017 at 14:50
  • 1
    $\begingroup$ @m_goldberg Thanks also for your insight about ToNumericalRegion. The mystery seems to remain as to why ImplicitRegion is unable to produce a well-formed region (one that passes the ConatantRegionQ test) from a general Boolean function. I will forward this issue to support. $\endgroup$
    – Will.Mo
    Nov 4, 2017 at 8:20
  • 1
    $\begingroup$ @Will.Mo. You are aware, are you not, that ImplicitRegion[inDisk[x, y] == True is not a boolean function? $\endgroup$
    – m_goldberg
    Nov 5, 2017 at 3:54
7
$\begingroup$

I noticed that ContourPlot does an awesome job discretizing the region. I don't know how this scales and how complicated your real problem is, but for the disk it works quite well:

ContourPlot excludes complex numbers from the region by default, so a convenient way to plot the region is to set the function value to the imaginary unit if the point is outside the region. Then we can get the mesh by using DiscretizeGraphics.

inDisk[x_?NumericQ, y_?NumericQ] := x^2 + y^2 <= 1;
mr = DiscretizeGraphics@
   ContourPlot[
    Piecewise[{{0, inDisk[x, y]}}, I], {x, -1, 1}, {y, -1, 1}, 
    Mesh -> All, PlotPoints -> 15]

enter image description here

If we want we can simplify the mesh:

simplemr = 
  TriangulateMesh[mr, MaxCellMeasure -> Infinity, 
   PerformanceGoal -> "Speed"]

enter image description here

And then integrate over this region:

NIntegrate[1, {x, y} ∈ simplemr]

3.14027

Not perfect, but close.

Edit: A cleaner solution might be to use RegionFunction instead of the imaginary unit, for example:

inDisk[x_?NumericQ, y_?NumericQ, f_] := x^2 + y^2 <= 1;
mr = BoundaryDiscretizeGraphics@
   ContourPlot[0, {x, -1, 1}, {y, -1, 1}, PlotPoints -> 5, 
    MaxRecursion -> 6, RegionFunction -> inDisk];
NIntegrate[1, {x, y} ∈ mr]

3.14118

Note that I added an argument f to inDisk, because that's the form RegionFunction expects. Also BoundaryDiscretizeGraphics should suffice here. Larger values of MaxRecursion can make the mesh region more accurate.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Slick tricks! It seems yet another one way is RegionPlot, i.e. mr = DiscretizeGraphics[ RegionPlot[inDisk[x, y] == True, {x, -1, 1}, {y, -1, 1}]]. I think the key idea here is that Mathematica's plotting routines are more robust/efficient than ImplicitRegion, which is currently a bit delicate. $\endgroup$
    – Will.Mo
    Nov 4, 2017 at 7:52
  • $\begingroup$ Yes, here's another example of delicate behavior: Area[ImplicitRegion[(x^2 + y^2 <= 1) == True, {{x, -1, 1}, {y, -1, 1}}]] gives 0, without complaining. Interestingly while you would expect RegionPlot and ContourPlot to work similar internally, it's easier to approach the correct value of $\pi$ with ContourPlot... $\endgroup$
    – M. Stern
    Nov 4, 2017 at 13:15
  • $\begingroup$ Yet another solution a lá @M.Stern can be found if one can replace the condition with a smooth function. Define r2Disk[x_?NumericQ, y_?NumericQ] := x^2 + y^2. Then nrDisk = ToNumericalRegion[ ImplicitRegion[r2Disk[x, y] <= 1, {{x, -1, 1}, {y, -1, 1}}]]. Then meshDisk = ToElementMesh[nrDisk];. Then check the accuracy: Pi - NIntegrate[1, {x, y} ∈ meshDisk]. This yields $\sim 10^{-7}$ error, i.e. the same as a "perfect" Disk[]. $\endgroup$
    – Will.Mo
    Nov 5, 2017 at 6:36
  • $\begingroup$ That is, the same error as with a mesh generated from a Disk[]. NIntegrate directly over a Disk[] is better, $\sim 10^{-9}$. $\endgroup$
    – Will.Mo
    Nov 5, 2017 at 6:46

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.