NDEigensystem: dependence on how the region is defined

Question

Excited to see version 11 advertised as solving partial differential eigenproblems - then disappointed when I tried to use NDEigensystem beyond simple cases. I wonder if anyone can explain why different definitions of the same region can lead to different behaviour.

My example is to solve a 2D Laplacian with a 4-lobed Dirichlet boundary defined by the plane-polar curve r=1+.3*Cos[4*phi]. The shape is unimportant, it's just an example. I defined this region in 3 different ways. The first two capture all points from the origin to r for all angles phi. The boundary r(x,y) along the radial line through x,y is expressed in cartesian coords either via phi using the ArcTan function (region1) or directly from x and y after some trig manipulation on a piece of paper (region2). The third way forms region3 rather artificially from region2 by taking 1/4 of it, making rotated copies and using RegionUnion:

boundary1[x_, y_] := 1. + 0.3*Cos[4*ArcTan[x, y]]; 
region1 = ImplicitRegion[x^2 + y^2 <= boundary1[x, y]^2, {x, y}]; 

boundary2[x_, y_] := 1. + 0.3*(1 - 8*x^2*(y^2/(x^2 + y^2)^2)); 
region2 = ImplicitRegion[x^2 + y^2 <= boundary2[x, y]^2, {x, y}]; 

regionSectorA = 
 RegionIntersection[Disk[{0, 0}, 2., {-Pi/4, Pi/4}], region2];
regionSectorB = 
  TransformedRegion[regionSectorA, RotationTransform[Pi/2]]; 
regionSectorC = 
  TransformedRegion[regionSectorA, RotationTransform[Pi]]; 
regionSectorD = 
  TransformedRegion[regionSectorA, RotationTransform[-Pi/2]]; 
region3 = 
  RegionUnion[regionSectorA, regionSectorB, regionSectorC, 
   regionSectorD];

They all seem to give the same, valid, region when plotted:

{RegionPlot[region1, PlotRange -> RegionBounds[region1]], RegionPlot[region2, PlotRange -> RegionBounds[region2]], RegionPlot[region3, PlotRange -> RegionBounds[region3]]}

but when I then try to find the first solution of the eigensystem using:

sol1 = NDEigensystem[{-Laplacian[u[x, y], {x, y}], DirichletCondition[u[x, y] == 0, True]}, u[x, y], Element[{x, y}, region1], 1]
sol2 = NDEigensystem[{-Laplacian[u[x, y], {x, y}], DirichletCondition[u[x, y] == 0, True]}, u[x, y], Element[{x, y}, region2], 1]
sol3 = NDEigensystem[{-Laplacian[u[x, y], {x, y}], DirichletCondition[u[x, y] == 0, True]}, u[x, y], Element[{x, y}, region3], 1]

the outcomes are quite different. region1 fails completely, with the unhelpful error message "Message text not found"! region2 succeeds without error. region3 complains that the bounds to the region are infinite (even though it's just plotted a finite region for me, above, as confirmed by RegionBounds) then delivers a different solution to the one from region2. Plotting the eigenfunctions for region2 and region3 using:

Quiet[Plot3D[-sol2[[2,1]], Element[{x, y}, region2], PlotRange -> RegionBounds[region3], PlotLabel -> Sqrt[sol2[[1,1]]]]]
Quiet[Plot3D[sol3[[2,1]], Element[{x, y}, region3], PlotRange -> RegionBounds[region3], PlotLabel -> Sqrt[sol3[[1,1]]]]]

shows that the region3 computation has arbitrarily truncated to the unit square before solving.

Can anyone please shed some light on what's going on? More than anything I'm looking for an understanding of what region definitions do and don't work, and why. (FYI the more-messy problem I really want to solve defines r(phi) for a Pi/n sector that then repeats periodically with phi. I can't think of a way to define this region without either ArcTan or RegionUnion.) I'm running version 11.0 on a 64 bit Win 8.1 PC.

Answer 1

Sorry to hear that you are disappointed. I'll get the messages issue fixed. You can directly generate the mesh and use that with NDEigensystem like so:

Needs["NDSolve`FEM`"]
sol1 = NDEigensystem[{-Laplacian[u[x, y], {x, y}], 
    DirichletCondition[u[x, y] == 0, True]}, u[x, y], 
   Element[{x, y}, ToElementMesh[region1]], 1][[1]]
{7.939178404359806`}

sol2 = NDEigensystem[{-Laplacian[u[x, y], {x, y}], 
    DirichletCondition[u[x, y] == 0, True]}, u[x, y], 
   Element[{x, y}, region2], 1][[1]]
{7.9390407461104076`}

sol3 = NDEigensystem[{-Laplacian[u[x, y], {x, y}], 
    DirichletCondition[u[x, y] == 0, True]}, u[x, y], 
   Element[{x, y}, ToElementMesh[region3, RegionBounds[region3]]], 1][[1]]
{7.940050668297508`}

The discrepancy in the result of the first eigenvalue is expected if you need a more accurate result you'd need a finer mesh, as detailed in the NDEigensystem ref page.

Here is what happens in the first case. The code needs to validate the PDE coefficients that are given and for that it needs a test coordinate. It is finding this test coordinate which takes long in this case and the code then gives up. You avoid this by specifying the mesh directly, then a test coordinate is just any coordinate from the mesh, which is easy (efficient) to extract.

In the third example the fast computation of the region bounds fails. Calling RegionBounds explicitly uses a different slower but more accurate algorithm.

One thing I am not quite sure about is what your definition of the region should be at {0,0} as

RegionMember[region1, {0, 0}]
ArcTan::indet
ArcTan::indet
0.` <= Indeterminate